1.2.1 Phonon Calculation: Displacing the Atoms along Normal Modes.
Once a Cartesian Phonon Calculation is performed and Normal modes are obtained, phonons can be recomputed by displacing the atoms along those normal modes. A phonon calculated along normal mode displacements have less numerical noise, particularly for low frequency modes and these frequencies give numerically more accurate free energies, see for example the following two papers: (1) https://doi.org/10.1021/ct500291x & (2) https://doi.org/10.1021/jacs.6b08646.
To calculate electron phonon coupling energies (EPCE), which are defined as the second derivative of the band energies (Kohn-Sham eigenvalues) with respect to a particular phonon mode, see Eq. S2 of the following paper: https://doi.org/10.1103/PhysRevMaterials.5.L070801. The sum of EPCEs over all phonon modes is the zero-point renormalization (ZPR) of that band. To compute EPCEs or ZPR using a frozen phonon method, it is necessary to recompute phonons by displacing the atoms along normal modes. In this exercise we will learn that.
The first step is to prepare the displaced coordinates. For that purpose again we import the following coordinate utilities tools, electron-phonon classes and pyepfd input/output tools from pyepfd.
[1]:
from pyepfd.coord_util import *
from pyepfd.elph_classes import *
from pyepfd.pyepfd_io import *
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.0
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2023-05-09 20:08:35.845891
Next, we will read all information from the checkpoint/restart XML file created in the previous step, i.e., Cartesian phonon calculation. We would instantiate an object named enmfd_inp of class read_pyepfd_info that belongs to coord_util object. The class read_pyepfd_info has several objects such as acoustic sum rule (asr), refined dynamical matrix (ref_dynmatrix), optimized coordinates(coord), atoms and cell.
[2]:
enmfd_inp = read_pyepfd_info(\
file_path='../1_cartesian_phonon/fdphonon.xml')
Time spent on read_pyepfd_info class: 0.0005269050598144531 s.
Now we would save these informations into different variables.
[3]:
asr = enmfd_inp.asr
inp_dynmat = enmfd_inp.ref_dynmatrix
opt_coord = enmfd_inp.coord
atoms = enmfd_inp.atoms
cell = enmfd_inp.cell
Like Cartesian phonon calculation, we again instantiate an object named fdmoves of ionic_mover class. Here we would use mode = ‘ENMFD’, i.e. energy-scaled normal mode finite displacement. Here we also need to provide the dynamical matrix (key: dynmat) and target energy displacement (deltae). Our target energy displacement is 0.001 Ha.
[4]:
fdmoves = ionic_mover( atoms = atoms,\
opt_coord = opt_coord,\
mode = 'ENMFD', \
dynmat = inp_dynmat,\
deltax = 0.005,\
deltae = 0.001 )
#Mode = 1 Disp-step(au) = 81.7641 Disp-step(Freq-scaled) = nan
#Config Disp(au) Disp(Freq-scaled)
2 81.7641 nan
3 -81.7641 nan
#Mode = 2 Disp-step(au) = 84.5714 Disp-step(Freq-scaled) = nan
#Config Disp(au) Disp(Freq-scaled)
4 84.5714 nan
5 -84.5714 nan
#Mode = 3 Disp-step(au) = 85.3889 Disp-step(Freq-scaled) = nan
#Config Disp(au) Disp(Freq-scaled)
6 85.3889 nan
7 -85.3889 nan
#Mode = 4 Disp-step(au) = 81.7642 Disp-step(Freq-scaled) = 0.0162
#Config Disp(au) Disp(Freq-scaled)
8 81.7642 0.0162
9 -81.7642 -0.0162
#Mode = 5 Disp-step(au) = 82.5022 Disp-step(Freq-scaled) = 0.0239
#Config Disp(au) Disp(Freq-scaled)
10 82.5022 0.0239
11 -82.5022 -0.0239
#Mode = 6 Disp-step(au) = 15.3008 Disp-step(Freq-scaled) = 0.8272
#Config Disp(au) Disp(Freq-scaled)
12 15.3008 0.8272
13 -15.3008 -0.8272
#Mode = 7 Disp-step(au) = 15.3008 Disp-step(Freq-scaled) = 0.8272
#Config Disp(au) Disp(Freq-scaled)
14 15.3008 0.8272
15 -15.3008 -0.8272
#Mode = 8 Disp-step(au) = 7.2918 Disp-step(Freq-scaled) = 0.5710
#Config Disp(au) Disp(Freq-scaled)
16 7.2918 0.5710
17 -7.2918 -0.5710
#Mode = 9 Disp-step(au) = 4.1770 Disp-step(Freq-scaled) = 0.4322
#Config Disp(au) Disp(Freq-scaled)
18 4.1770 0.4322
19 -4.1770 -0.4322
Time spent on ionic_mover class: 0.0017631053924560547 s.
/home/arpan/.local/lib/python3.8/site-packages/pyepfd/coord_util.py:676: RuntimeWarning: invalid value encountered in sqrt
freq_scaling = np.sqrt(nmfd.omega[imode])
First 5 modes are rigid body translation and rotation mode and in frequency scaled coordinate displacement are imaginary. This is okay.
Now all the displaced coordinates are stored within fdmoves object. Next step is to write them into an xyz file. To do that we instantiate an object named out_xyz of class xyz with input/output mode write(w).
[5]:
out_xyz = xyz(file_path = 'enmfd_phonon.xyz',\
io = 'w', atoms = atoms)
Time spent on xyz class: 0.0002453327178955078 s.
Next we write the enmfd_phonon.xyz file step by step applying a loop.
[6]:
ndisp = fdmoves.disp_coord.shape[1]
for i in range(ndisp):
out_xyz.write(cell = cell,coord = fdmoves.disp_coord[:,i])
del out_xyz
This generated the coordinates by displacing the atoms along normal modes and wrote it into enmfd_phonon.xyz file. Let us have a look at this file.
[7]:
%%bash
cat enmfd_phonon.xyz
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 0 positions{angstrom} cell{angstrom}
O 0 0 1.1571328
O 0 0 -1.1571313
C 0 0 -1.5875317e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 1 positions{angstrom} cell{angstrom}
O -2.1644539e-08 -1.3546308e-09 1.3098932
O -1.0825663e-08 -8.5003505e-09 -1.0043709
C -1.623509e-08 -4.927498e-09 0.15275873
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 2 positions{angstrom} cell{angstrom}
O 2.1644539e-08 1.3546308e-09 1.0043725
O 1.0825663e-08 8.5003505e-09 -1.3098916
C 1.623509e-08 4.927498e-09 -0.1527619
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 3 positions{angstrom} cell{angstrom}
O 1.5349446e-07 0.23739929 1.1571328
O -6.6299837e-07 -0.09108796 -1.1571313
C -2.5475279e-07 0.073155326 -1.5879651e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 4 positions{angstrom} cell{angstrom}
O -1.5349446e-07 -0.23739929 1.1571328
O 6.6299837e-07 0.09108796 -1.1571313
C 2.5475279e-07 -0.073155326 -1.5870984e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 5 positions{angstrom} cell{angstrom}
O 0.1862744 -5.6314584e-07 1.1571328
O -0.18790509 6.5617577e-08 -1.1571312
C -0.00081572774 -2.4876349e-07 -1.5828014e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 6 positions{angstrom} cell{angstrom}
O -0.1862744 5.6314584e-07 1.1571328
O 0.18790509 -6.5617577e-08 -1.1571313
C 0.00081572774 2.4876349e-07 -1.5922621e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 7 positions{angstrom} cell{angstrom}
O -0.15367421 -3.37338e-07 1.1571328
O -0.15184276 2.6639613e-07 -1.1571313
C -0.15275848 -3.5470314e-08 -1.6037902e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 8 positions{angstrom} cell{angstrom}
O 0.15367421 3.37338e-07 1.1571329
O 0.15184276 -2.6639613e-07 -1.1571312
C 0.15275848 3.5470314e-08 -1.5712733e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 9 positions{angstrom} cell{angstrom}
O -1.9778686e-07 -0.052928795 1.1571328
O 2.1686863e-08 -0.22031741 -1.1571313
C -8.8049775e-08 -0.13662328 -1.593362e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 10 positions{angstrom} cell{angstrom}
O 1.9778686e-07 0.052928795 1.1571328
O -2.1686863e-08 0.22031741 -1.1571312
C 8.8049775e-08 0.13662328 -1.5817015e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 11 positions{angstrom} cell{angstrom}
O 0.016641655 0.0054577579 1.1571328
O 0.016641724 0.0054577804 -1.1571312
C -0.044336641 -0.01454054 -1.5961057e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 12 positions{angstrom} cell{angstrom}
O -0.016641655 -0.0054577579 1.1571329
O -0.016641724 -0.0054577804 -1.1571313
C 0.044336641 0.01454054 -1.5789578e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 13 positions{angstrom} cell{angstrom}
O 0.0054577478 -0.016641625 1.1571328
O 0.0054577703 -0.016641693 -1.1571312
C -0.014540513 0.044336559 -1.5914266e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 14 positions{angstrom} cell{angstrom}
O -0.0054577478 0.016641625 1.1571328
O -0.0054577703 0.016641693 -1.1571313
C 0.014540513 -0.044336559 -1.5836369e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 15 positions{angstrom} cell{angstrom}
O 5.1422293e-09 1.285583e-09 1.1731096
O 5.1422505e-09 1.2855883e-09 -1.1731078
C -1.369991e-08 -3.425046e-09 -1.9321051e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 16 positions{angstrom} cell{angstrom}
O -5.1422293e-09 -1.285583e-09 1.141156
O -5.1422505e-09 -1.2855883e-09 -1.1411547
C 1.369991e-08 3.425046e-09 -1.2429584e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 17 positions{angstrom} cell{angstrom}
O -9.5921236e-10 1.054293e-10 1.1619138
O -9.5921631e-10 1.0542974e-10 -1.15235
C 2.5555303e-09 -2.8088439e-10 -0.012739473
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 18 positions{angstrom} cell{angstrom}
O 9.5921236e-10 -1.054293e-10 1.1523518
O 9.5921631e-10 -1.0542974e-10 -1.1619125
C -2.5555303e-09 2.8088439e-10 0.012736298
The next step is to convert them into qbox input. This can be done using xyz2qbox available in util/qbox_utils folder. This is shown in example 1.2.2
[ ]: