Source code for jarvis.analysis.phonon.ir

"""
Modules for analyzing infrared intensities.

Please find more details in https://doi.org/10.1038/s41524-020-0337-2 .
"""
from jarvis.core.spectrum import Spectrum
import numpy as np
import os


[docs]def normalize_vecs(phonon_eigenvectors, masses): """ Return the eigenvectors after division of each component by sqrt(mass). Adapted from https://github.com/JMSkelton/Phonopy-Spectroscopy/ TODO: include LO-TO splitting. """ nmodes = len(phonon_eigenvectors) nmasses = len(masses) natoms = nmasses sqrt_masses = np.sqrt(masses) eigendisplacements = np.zeros((nmodes, natoms, 3), dtype="float64") for i in range(0, nmodes): eigenvector = phonon_eigenvectors[i] for j in range(0, natoms): eigendisplacements[i, j, :] = np.divide( eigenvector[j], sqrt_masses[j] ) return eigendisplacements
[docs]def ir_intensity( phonon_eigenvectors=[], phonon_eigenvalues=[], masses=[], born_charges=[], factor=33.35641, nac=True, epsilon=[], enforce_positive_freqs=True, smoothen=True, ): """Calculate IR intensity using DFPT.""" # TODO:add non-analytical correction eigendisplacements = normalize_vecs(phonon_eigenvectors, masses) becDim1, becDim2, becDim3 = np.shape(born_charges) freq = [] ir_ints = [] for i, v in zip(eigendisplacements, phonon_eigenvalues): eigendisplacement = i eigDim1, eigDim2 = np.shape(eigendisplacement) irIntensity = 0.0 for a in range(0, 3): sumTemp1 = 0.0 for j in range(0, eigDim1): sumTemp2 = 0.0 for b in range(0, 3): sumTemp2 += born_charges[j][a][b] * eigendisplacement[j][b] sumTemp1 += sumTemp2 irIntensity += sumTemp1 ** 2 if enforce_positive_freqs: if v > 0: freq.append(v * factor) # Thz to cm-1 ir_ints.append(irIntensity) else: freq.append(v * factor) # Thz to cm-1 ir_ints.append(irIntensity) if smoothen: freq, ir_ints = Spectrum(x=freq, y=ir_ints).smoothen_spiky_spectrum() return freq, ir_ints
[docs]def ir_intensity_phonopy( run_dir=".", vasprun="vasprun.xml", BornFileName="BORN", PoscarName="POSCAR", ForceConstantsName="FORCE_CONSTANTS", supercell=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], nac=True, symprec=1e-5, primitive=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], degeneracy_tolerance=1e-5, vector=[0, 0, 0], smoothen=False, ): """Calculate IR intensity using DFPT and phonopy.""" from phonopy import Phonopy from phonopy.interface.vasp import read_vasp from phonopy.file_IO import ( parse_BORN, parse_FORCE_CONSTANTS, ) import shutil from phonopy.units import VaspToCm # from phonopy.phonon.degeneracy import ( # degenerate_sets as get_degenerate_sets, # ) # adapted from https://github.com/JaGeo/IR # TODO: Make directory indepndent cwd = os.getcwd() print("Directory:", cwd) os.chdir(run_dir) if not os.path.exists(vasprun): shutil.copy2(vasprun, "vasprun.xml") cmd = str("phonopy --fc vasprun.xml") os.system(cmd) born_file = os.path.join(os.getcwd(), BornFileName) cmd = str("phonopy-vasp-born > ") + str(born_file) os.system(cmd) from jarvis.io.vasp.outputs import Vasprun v = Vasprun(vasprun) strt = v.all_structures[0] strt.write_poscar(PoscarName) unitcell = read_vasp(PoscarName) phonon = Phonopy( unitcell, supercell_matrix=supercell, primitive_matrix=primitive, factor=VaspToCm, symprec=symprec, ) natoms = phonon.get_primitive().get_number_of_atoms() force_constants = parse_FORCE_CONSTANTS(filename=ForceConstantsName) phonon.set_force_constants(force_constants) masses = phonon.get_primitive().get_masses() phonon.set_masses(masses) BORN_file = parse_BORN(phonon.get_primitive(), filename=BornFileName) BORN_CHARGES = BORN_file["born"] # print ('born_charges2',BORN_CHARGES) if nac: phonon.set_nac_params(BORN_file) frequencies, eigvecs = phonon.get_frequencies_with_eigenvectors(vector) # frequencies=VaspToTHz*frequencies/VaspToCm # x, y = ir_intensity( # phonon_eigenvectors=np.real(eigvecs), # phonon_eigenvalues=frequencies, # masses=masses, #np.ones(len(masses)), # born_charges=born_charges, # smoothen=smoothen, # ) NumberOfBands = len(frequencies) EigFormat = {} for alpha in range(NumberOfBands): laufer = 0 for beta in range(natoms): for xyz in range(0, 3): EigFormat[beta, alpha, xyz] = eigvecs[laufer][alpha] laufer = laufer + 1 Intensity = {} intensities = [] for freq in range(len(frequencies)): Intensity[freq] = 0 tmp = 0 for alpha in range(3): asum = 0 for n in range(natoms): for beta in range(3): asum = asum + BORN_CHARGES[n, alpha, beta] * np.real( EigFormat[n, freq, beta] ) / np.sqrt(masses[n]) tmp += asum Intensity[freq] = Intensity[freq] + np.power(np.absolute(asum), 2) intensities.append(Intensity[freq]) os.chdir(cwd) return frequencies, intensities
""" if __name__ == "__main__": from jarvis.io.vasp.outputs import Vasprun, Outcar out = Outcar("../../tests/testfiles/io/vasp/OUTCAR.JVASP-39") vrun = Vasprun("../../tests/testfiles/io/vasp/vasprun.xml.JVASP-39") phonon_eigenvectors = vrun.dfpt_data["phonon_eigenvectors"] vrun_eigs = vrun.dfpt_data["phonon_eigenvalues"] phonon_eigenvalues = out.phonon_eigenvalues masses = vrun.dfpt_data["masses"] born_charges = vrun.dfpt_data["born_charges"] x, y = ir_intensity( phonon_eigenvectors=phonon_eigenvectors, phonon_eigenvalues=phonon_eigenvalues, masses=masses, born_charges=born_charges, smoothen=False, ) for i, j in zip(x, y): if j > 0.1: print("first", i, j) print() x, y = ir_intensity_phonopy() for i, j in zip(x, y): if j > 0.1: print("later", i, j) # for i, j in zip(phonon_eigenvalues, vrun_eigs): # print(i, j) """