Eigenvector normalization problem in Raman workflow #762
Description
Hi,
I think the eigenvector normalization in Raman workflow has a bug that In practice can make small errors in Raman tensors and intensities. For example, MoS2 is a well-known system and there are a lot of reports about that in literature. It has only two Raman-active modes, but when calculated with atomate, an extra mode appeared in the Raman spectrum.
The issue raised from WriteNormalmodeDisplacedPoscar class in atomate/vasp/firetasks/write_inputs.py and RamanTensorToDb class in atomate/vasp/firetasks/parse_outputs.py
Eigenvector is size [natoms,3]. atomate default approach calculates the norm over row (x,y,z) for each atom, and gets the vector of size [natoms,1], and then normalizes each row. It means that each atom is displaced by the step size (default = 0.005 Å). That makes no sense if the eigenvector is such that, e.g., S should dominate but Mo also has a small non-zero contribution (due to noise, for example). Also, the default approach missed dividing eigenvectors by sqrt(mass).
The correct way to normalize is to divide each row of eigenvector by sqrt(mass) and then normalize with the norm of the whole eigenvector (a single number).
`class WriteNormalmodeDisplacedPoscar(FiretaskBase):
required_params = ["mode", "displacement"]
def run_task(self, fw_spec):
mode = self["mode"]
disp = self["displacement"]
structure = Structure.from_file("POSCAR")
masses = np.array([site.specie.data['Atomic mass'] for site in structure])
nm_eigenvecs = np.array(fw_spec["normalmodes"]["eigenvecs"])
nm_eigenvec_sqm = nm_eigenvecs[mode, :, :] / np.sqrt(masses[:, np.newaxis])
nm_norms = np.linalg.norm(nm_eigenvec_sqm)
nm_displacement = (nm_eigenvec_sqm * disp / nm_norms)
for i, vec in enumerate(nm_displacement):
structure.translate_sites(i, vec, frac_coords=False)
structure.to(fmt="poscar", filename="POSCAR")`
`
class RamanTensorToDb(FiretaskBase):
optional_params = ["db_file"]
def run_task(self, fw_spec):
nm_eigenvecs = np.array(fw_spec["normalmodes"]["eigenvecs"])
nm_eigenvals = np.array(fw_spec["normalmodes"]["eigenvals"])
structure = fw_spec["normalmodes"]["structure"]
masses = np.array([site.specie.data["Atomic mass"] for site in structure])
# To get the actual eigenvals, the values read from vasprun.xml must be multiplied by -1.
# frequency_i = sqrt(-e_i)
# To convert the frequency to THZ: multiply sqrt(-e_i) by 15.633
# To convert the frequency to cm^-1: multiply sqrt(-e_i) by 82.995
nm_frequencies = np.sqrt(np.abs(nm_eigenvals)) * 82.995 # cm^-1
d = {
"structure": structure.as_dict(),
"formula_pretty": structure.composition.reduced_formula,
"normalmodes": {
"eigenvals": fw_spec["normalmodes"]["eigenvals"],
"eigenvecs": fw_spec["normalmodes"]["eigenvecs"],
},
"frequencies": nm_frequencies.tolist(),
}
# store the displacement & epsilon for each mode in a dictionary
mode_disps = fw_spec["raman_epsilon"].keys()
modes_eps_dict = defaultdict(list)
for md in mode_disps:
modes_eps_dict[fw_spec["raman_epsilon"][md]["mode"]].append(
[
fw_spec["raman_epsilon"][md]["displacement"],
fw_spec["raman_epsilon"][md]["epsilon"],
]
)
# raman tensor = finite difference derivative of epsilon wrt displacement.
raman_tensor_dict = {}
**scale = structure.volume / 4.0 / np.pi**
for k, v in modes_eps_dict.items():
**nm_eigenvec_sqm = nm_eigenvecs[k, :, :] / np.sqrt(masses[:, np.newaxis])**
**nm_norms = np.linalg.norm(nm_eigenvec_sqm)**
raman_tensor = (np.array(v[0][1]) - np.array(v[1][1])) / (v[0][0] - v[1][0])
# frequency in cm^-1
omega = nm_frequencies[k]
if nm_eigenvals[k] > 0:
logger.warning(f"Mode: {k} is UNSTABLE. Freq(cm^-1) = {-omega}")
**raman_tensor = scale * raman_tensor * nm_norms**
raman_tensor_dict[str(k)] = raman_tensor.tolist()
d["raman_tensor"] = raman_tensor_dict
d["state"] = "successful"
# store the results
db_file = env_chk(self.get("db_file"), fw_spec)
if not db_file:
with open("raman.json", "w") as f:
f.write(json.dumps(d, default=DATETIME_HANDLER))
else:
db = VaspCalcDb.from_db_file(db_file, admin=True)
db.collection = db.db["raman"]
db.collection.insert_one(d)
logger.info("Raman tensor calculation complete.")
return FWAction()
`