Fix the subroutine that identifies if a phonon mode is Raman active.

[mesonepigreco]

Sometimes, the subroutine computing the Raman active modes gives a false negative if no Raman tensor is provided.

dyn.GetRamanActive()

The reason is that a random Raman tensor is generated, and then, the Raman activity is measured for all modes.
It is possible that some modes are, by chance, almost orthogonal to the Raman tensor, thus resulting in a false Raman inactive mode.

This could be either fixed by lowering the threshold (dangerous, needs to be tested), or by generating a Raman tensor explicitly to have a nonzero projection on every normal mode before symmetrization (probably better as it would be deterministic).

[Shinleon]

Hello! I was just surfing through good first issues and was wondering if I could help with this? I was reading through this description of Raman matrices along with your code in Phonon.py but was wondering what it meant for a tensor to have a nonzero projection on every normal mode?

[mesonepigreco]

Hi,
So, a Raman tensor in cartesian space is a 4-dimensional array 3x3x3xN, that quantifies how much the polarizability (which is a 3x3 matrix ) of material changes if we move an atom (a 3xN vector indicating the position of the atoms, where N is the number of atoms, 3 per cartesian component).
The Raman tensor projected on the phonon modes is equation 2 of this work ($A^\nu$) Lazzeri, Mauri, Physical Review Letters, 90, 036401, 2005. In particular, this is the relation:

Equation 1 allows us to compute the Raman intensity of the mode $\nu$, determined by the Raman tensor $A^\nu$ (exactly the 3x3 matrix of your reference).

The target of the mentioned subroutine is to determine whether some Raman modes are active based on symmetries analysis. This is done by generating a random Cartesian Raman tensor ( $\frac{\partial^3 \mathcal E}{\partial E_a\partial E_b \partial u_{kl}}$)

self.raman_tensor = np.zeros((3,3, 3*self.structure.N_atoms), dtype = np.double)
self.raman_tensor[:,:,:] = np.random.uniform( size = self.raman_tensor.shape)

enforcing symmetries on it,

# Get the symmetries
qe_sym = symmetries.QE_Symmetry(self.structure)
if use_spglib:
    qe_sym.SetupFromSPGLIB()
else:
    qe_sym.SetupQPoint()

# Symmetrize the effective charges
qe_sym.ApplySymmetryToRamanTensor(self.raman_tensor)

then calculating the trace of $A^\nu$ for each mode, and checking whether it is zero.

res = np.zeros(self.structure.N_atoms * 3, dtype = np.double)
for i, j in itertools.product(range(3) , range(3)):
    pol_in = np.zeros(3)
    pol_in[i] = 1
    pol_out = np.zeros(3)
    pol_out[j] = 1

    res += self.GetRamanResponce(pol_in, pol_out)

is_raman_active = res > 1e-5

However, as mentioned in the issue, it could happen with a random generation of the original Raman Tensor that the evaluated $A^\nu$ is very small by chance, thus resulting in a false identification that a mode should be Raman invisible by symmetries.

To fix this, we need to replace the generation of the random Raman tensor $\frac{\partial^3 \mathcal E}{\partial E_a\partial E_b \partial u_{kl}}$ with generating a nonzero $A^\nu$ for every $\nu$ mode (the identity 3x3 matrix, for example), then invert equation 2 to get the full Raman tensor $\frac{\partial^3 \mathcal E}{\partial E_a\partial E_b \partial u_{kl}}$. (then the method should proceed as before:

• apply symmetries on $\frac{\partial^3 \mathcal E}{\partial E_a\partial E_b \partial u_{kl}}$
• compute the intensity of every mode
• check whether the intensity is above a threshold

I marked this as a good first issue as it is very well isolated, as one only needs to work on a specific subroutine. The solution requires changing a few lines of code to that subroutine (maybe less than 10). However, it requires a bit of understanding of the problem.