1.1.3 Post-processing Qbox Outputs to compute dynamical matrix.
Here we would learn how to post-process the qbox outputs after all the jobs regarding the displaced coordinates have finished. Here the qbox output files are fd-1.r and fd-2.r.
For that purpose we need to import the following objects 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:03:02.680114
First, we would read the qbox outputs and store them within a list object. Let us first initialize the list object named qbouts.
[2]:
qbouts = []
Now we would read the qbox outputs using the qbox class available within coord_utils.
[3]:
for i in range(2):
qbouts.append(qbox(file_path = 'fd-'+str(i+1)+'.r',io = 'r'))
Time spent on qbox class: 0.012288093566894531 s.
Time spent on qbox class: 0.010625839233398438 s.
We are running the loop twice as we have only two outputs. Change it according to your needs. qbox class have objects that stores coordinates (coords), atoms, mass, forces etc. Now we would store them into separate variables from the list of qbox objects. For that purpose first we initialize these objects with the first element of the qbouts.
[4]:
coords = qbouts[0].coords
atoms = qbouts[0].atoms
mass = qbouts[0].mass
forces = qbouts[0].forces
Next we will append other elements upto the length of the qbouts.
[5]:
for i in range(1,len(qbouts)):
coords = np.vstack((coords,qbouts[i].coords))
forces = np.vstack((forces,qbouts[i].forces))
The optimized geometry is the first set of coordinates. So we define it accordingly,
[6]:
opt_coord = coords[0]
The cell parameters are fixed. Therefore we can take the first cell_parameters and conver them into vectors using abc2h function available within pyepfd.coord_utils.
[7]:
cell_v = abc2h(qbouts[0].cell)
Now we compute the force constant matrix and diagonaliz it to obtained normal-modes and its frequencies. This is done using the phonon_calculator class available within elph_classes.
[8]:
phonon = phonon_calculator(forces = forces,\
mass = mass,\
mode = 'FD',\
deltax = 0.005)
Please remember to use the same mode and deltax that was used to prepare the input files, i.e., on tutorial cart_phonon_1.
Next, we would take the computed dynamical matrix and instantiate a dynamica matrix (dm) class available in elph_classes. This class has the methods to refine the dynamical matrix using the suitable acoustic sum rule (asr). Here we have a linear molecule so we will choose asr = ‘lin’. For a non-linear polyatomic molecule and crystal, one should use asr = ‘lin’ & asr = ‘crystal’, respectively.
[9]:
dm = dm(dynmat = phonon.dynmat, mass = mass)
dm.apply_asr(opt_coord = opt_coord, asr = 'lin')
Now it is time to write all information into a checkpoint/restart file that can be used for future reference. This file would be the starting point for calculating a normal-mode finite difference method, stochastic method, etc.
[10]:
restart_file = write_pyepfd_info(inp_dynmat = None,\
dynmat = dm.dynmatrix,\
ref_dynmat = dm.refdynmatrix,\
mass = mass,atoms = atoms,\
opt_coord = opt_coord,\
cell = qbouts[0].cell,\
file_name='fdphonon.xml',\
mode='fd',\
deltax=0.005,\
asr = 'lin')
# deleting the restart file object to finish writing
del restart_file
Time spent at write_info 0.0011 s
Let us have a look at the restart file fdphonon.xml.
[11]:
%%bash
cat fdphonon.xml
<pyepfd>
<phonon mode='fd'>
<ngrid> 1 </ngrid>
<deltax> 0.005 </deltax>
<deltae> 0.001 </deltae>
<asr> lin </asr>
<dynmat shape='(9, 9)'>
[ 2.12445535e-07, 3.42875298e-11, 1.71437649e-11, 1.33933949e-06, 0.00000000e+00,
-1.71437649e-11, -3.07800297e-06, 1.97867238e-11, 3.95734476e-11, 3.42875298e-11,
2.12445535e-07, 1.71437649e-11, 3.42875298e-11, 1.33932235e-06, -1.71437649e-11,
1.97867238e-11, -3.07800297e-06, 1.97867238e-11, 1.71437649e-11, 1.71437649e-11,
3.45437948e-05, 0.00000000e+00, 0.00000000e+00, -3.07092832e-06, 0.00000000e+00,
0.00000000e+00, -3.67259006e-05, 1.33933949e-06, 3.42875298e-11, 0.00000000e+00,
2.12274097e-07, 0.00000000e+00, -3.42875298e-11, -3.07790403e-06, 1.97867238e-11,
-3.95734476e-11, 0.00000000e+00, 1.33932235e-06, 0.00000000e+00, 0.00000000e+00,
2.12342672e-07, 1.71437649e-11, 0.00000000e+00, -3.07788425e-06, -1.97867238e-11,
-1.71437649e-11, -1.71437649e-11, -3.07092832e-06, -3.42875298e-11, 1.71437649e-11,
3.45451663e-05, 0.00000000e+00, 0.00000000e+00, -3.67272857e-05, -3.07800297e-06,
1.97867238e-11, 0.00000000e+00, -3.07790403e-06, 0.00000000e+00, 0.00000000e+00,
5.83324319e-06, 0.00000000e+00, -2.28371329e-11, 1.97867238e-11, -3.07800297e-06,
0.00000000e+00, 1.97867238e-11, -3.07788425e-06, 0.00000000e+00, 0.00000000e+00,
5.83328886e-06, 0.00000000e+00, 3.95734476e-11, 1.97867238e-11, -3.67259006e-05,
-3.95734476e-11, -1.97867238e-11, -3.67272857e-05, -2.28371329e-11, 0.00000000e+00,
8.22001133e-05 ]
</dynmat>
<ref_dynmat shape='(9, 9)'>
[ 1.16570877e-06, -1.27450455e-12, 1.77449774e-12, 1.16571357e-06, -1.27450979e-12,
-7.58294525e-12, -2.69084475e-06, 2.94198171e-12, 6.70390356e-12, -1.27450455e-12,
1.16571224e-06, 1.48889657e-12, -1.27450979e-12, 1.16571704e-06, -8.50449734e-13,
2.94198171e-12, -2.69085276e-06, -7.36872634e-13, 1.77449774e-12, 1.48889657e-12,
3.44489401e-05, 1.77450504e-12, 1.48890270e-12, -3.16584535e-06, -4.09613279e-12,
-3.43687002e-12, -3.61058355e-05, 1.16571357e-06, -1.27450979e-12, 1.77450504e-12,
1.16571837e-06, -1.27451504e-12, -7.58297646e-12, -2.69085583e-06, 2.94199381e-12,
6.70393115e-12, -1.27450979e-12, 1.16571704e-06, 1.48890270e-12, -1.27451504e-12,
1.16572184e-06, -8.50453234e-13, 2.94199381e-12, -2.69086383e-06, -7.36875666e-13,
-7.58294525e-12, -8.50449734e-13, -3.16584535e-06, -7.58297646e-12, -8.50453234e-13,
3.44501870e-05, 1.75039675e-11, 1.96312172e-12, -3.61072746e-05, -2.69084475e-06,
2.94198171e-12, -4.09613279e-12, -2.69085583e-06, 2.94199381e-12, 1.75039675e-11,
6.21136742e-06, -6.79107529e-12, -1.54748460e-11, 2.94198171e-12, -2.69085276e-06,
-3.43687002e-12, 2.94199381e-12, -2.69086383e-06, 1.96312172e-12, -6.79107529e-12,
6.21138590e-06, 1.70094788e-12, 6.70390356e-12, -7.36872634e-13, -3.61058355e-05,
6.70393115e-12, -7.36875666e-13, -3.61072746e-05, -1.54748460e-11, 1.70094788e-12,
8.33458035e-05 ]
</ref_dynmat>
</phonon>
<cell shape='(6)'>
[ 2.00000000e+01, 2.00000000e+01, 2.00000000e+01, 9.00000000e+01, 9.00000000e+01,
9.00000000e+01 ]
</cell>
<mass shape='(9)'>
[ 2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04, 2.91651223e+04,
2.91651223e+04, 2.18941669e+04, 2.18941669e+04, 2.18941669e+04 ]
</mass>
<opt_coord shape='(9)'>
[ 0.00000000e+00, 0.00000000e+00, 2.18666400e+00, 0.00000000e+00, 0.00000000e+00,
-2.18666100e+00, 0.00000000e+00, 0.00000000e+00, -3.00000000e-06 ]
</opt_coord>
<atoms shape='(3)'>
[ O, O, C ]
</atoms>
</pyepfd>
We see this is an xml file.It can be used to gather all information for further pyEPFD calculations, for example a normal-mode Hessian calculation or stochastic electron-phonon calculation.
Sometime we would like to visualize the normal modes. For that purpose common visualization codes can be used and pyEPFD has methods to write normal-modes into an ‘axsf’, ‘nmd’ or ‘molden’ files that can be visualized with XCRYSDEN, VMD or MOLDEN programs, respectively. The below example shows how to do it using XCRYSDEN file as an example.
[12]:
nmodes = write_nmode(atoms = atoms,
cell_v = cell_v,
opt_coord = opt_coord,
mode_v = dm.refV, # refined normal mode vectors after ASR
mode_freq = dm.refomega, #refined frequencies
fmt='axsf', #Options are: 'axsf','nmd','molden'
file_path='refdynmat')
# Deleting the nmodes object to finish printing the information into the file
del nmodes
[13]:
%%bash
cat refdynmat.axsf
ANIMSTEPS 9
CRYSTAL
PRIMVEC
10.583544980 0.000000000 0.000000000
0.000000000 10.583544980 0.000000000
0.000000000 0.000000000 10.583544980
PRIMCOORD 1
3 1
O 0 0 1.15713 0.00543872 -0.000154545 -7.63158e-07
O 0 0 -1.15713 -0.00078962 9.78154e-06 -7.63158e-07
C 0 0 -1.58753e-06 0.00232455 -7.23814e-05 -7.63158e-07
PRIMCOORD 2
3 1
O 0 0 1.15713 -1.5501e-05 0.00180638 -0.000434643
O 0 0 -1.15713 0.00502252 -0.000506981 -0.000434643
C 0 0 -1.58753e-06 0.00250351 0.000649697 -0.000434643
PRIMCOORD 3
3 1
O 0 0 1.15713 0.000169877 0.00511544 1.96058e-05
O 0 0 -1.15713 -0.00179793 -0.000373095 1.96058e-05
C 0 0 -1.58753e-06 -0.000814027 0.00237117 1.96058e-05
PRIMCOORD 4
3 1
O 0 0 1.15713 -4.82234e-07 0.000348649 0.000149111
O 0 0 -1.15713 -0.00032832 -0.0054022 0.000149111
C 0 0 -1.58753e-06 -0.000164401 -0.00252678 0.000149111
PRIMCOORD 5
3 1
O 0 0 1.15713 1.6699e-06 -0.000180755 -0.0035005
O 0 0 -1.15713 -0.000647508 -0.00016926 -0.0035005
C 0 0 -1.58753e-06 -0.00032292 -0.000175007 -0.0035005
PRIMCOORD 6
3 1
O 0 0 1.15713 0.00205537 0.000673933 -2.22525e-09
O 0 0 -1.15713 0.00205537 0.000673936 3.0202e-09
C 0 0 -1.58753e-06 -0.0054759 -0.00179549 -1.05895e-09
PRIMCOORD 7
3 1
O 0 0 1.15713 0.000673933 -0.00205537 -8.26929e-12
O 0 0 -1.15713 0.000673936 -0.00205537 3.69328e-10
C 0 0 -1.58753e-06 -0.00179549 0.0054759 -4.80965e-10
PRIMCOORD 8
3 1
O 0 0 1.15713 1.33266e-09 3.33171e-10 0.00414054
O 0 0 -1.15713 1.33266e-09 3.33172e-10 -0.00414047
C 0 0 -1.58753e-06 -3.55046e-09 -8.87632e-10 -8.92961e-08
PRIMCOORD 9
3 1
O 0 0 1.15713 -4.33957e-10 4.76973e-11 0.00216297
O 0 0 -1.15713 -4.33959e-10 4.76975e-11 0.0021631
C 0 0 -1.58753e-06 1.15615e-09 -1.27075e-10 -0.00576275
The refdynmat.axsf containes the normal modes after refining them applying acoustic sum rule. This file can be visualized with XCRYSDEN. If we choose fmt=‘nmd’ / ‘molden’ that would create .nmd / .molden file respectively. They can be visualized with VMD and Molden, respectively.
[ ]: