1.2.3 Post-processing Qbox Outputs to compute dynamical matrix.
Here we would learn how to post-process the qbox output(s) enmfd-1.r after all the job(s) regarding the displaced coordinates have finished.
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:16:12.115427
Like example 2.1, we again read all information from the checkpoint/restart XML file created in the 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.0005414485931396484 s.
Now we would save these informations into different variable, similar to example 2.1.
[3]:
asr = enmfd_inp.asr
inp_dynmat = enmfd_inp.ref_dynmatrix
opt_coord = enmfd_inp.coord
atoms = enmfd_inp.atoms
cell = enmfd_inp.cell
cell_v = abc2h(cell)
Now, we would read the qbox output(s) and store them within a list object. Let us first initialize the list object named qbouts. Then we would read the qbox outputs using the qbox class available within coord_utils.
[4]:
qbouts = []
for i in range(1):
qbouts.append(qbox(file_path = 'enmfd-'+str(i+1)+'.r',io = 'r'))
Reording atoms, forces, and coordinates.
Time spent on qbox class: 0.017906665802001953 s.
Note: during qbox calculation sometimes it rearranges the atom order. Therefore PyEPFD again reverting back to input order.
First, we would read the qbox outputs and store them within a list object. Let us first initialize the list object named qbouts.
[5]:
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.(Note as here we have only one qbox output file the next step is unnecessary, however usually we will have more than one qbox output to harness the parallel submission of of multiple jobs)
[6]:
for i in range(1,len(qbouts)):
coords = np.vstack((coords,qbouts[i].coords))
forces = np.vstack((forces,qbouts[i].forces))
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.
[7]:
phonon = phonon_calculator(forces = forces,\
mass = mass,\
dynmat = inp_dynmat,\
mode = 'ENMFD',\
deltax = 0.005,\
deltae = 0.001,\
ngrid = 1)
Please remember to use the same mode, deltax and deltae that was used to prepare the input files, i.e., on tutorial enmfd_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.
[8]:
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.
[9]:
restart = write_pyepfd_info(inp_dynmat = inp_dynmat,\
dynmat = dm.dynmatrix,\
ref_dynmat = dm.refdynmatrix,\
mass = mass,\
atoms = atoms,\
opt_coord = opt_coord,\
cell = cell,\
file_name='enmfdphonon.xml',\
mode='enmfd',\
deltax=0.005,\
deltae= 0.001,\
asr = asr)
# Deleting the restart object to finishing the file writing
del restart
Time spent at write_info 0.0017 s
[10]:
%%bash
cat enmfdphonon.xml
<pyepfd>
<phonon mode='enmfd'>
<ngrid> 1 </ngrid>
<deltax> 0.005 </deltax>
<deltae> 0.001 </deltae>
<asr> lin </asr>
<inp_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 ]
</inp_dynmat>
<dynmat shape='(9, 9)'>
[ 1.52149567e-06, -2.91302688e-12, 4.40786142e-14, 1.00370699e-06, -1.38621426e-12,
-4.18653074e-12, -2.98797305e-06, -2.22351831e-09, 3.40829809e-12, -2.91302688e-12,
1.45170088e-06, 3.17891915e-12, -9.32830194e-13, 1.08651713e-06, 8.91803023e-13,
2.23137564e-09, -2.96824894e-06, 5.97604638e-12, 4.40786142e-14, 3.17891915e-12,
3.46788383e-05, -2.04652743e-13, 1.37362634e-11, -2.89580040e-06, 2.83697844e-13,
7.15144431e-12, -3.66635112e-05, 1.00370699e-06, -9.32830194e-13, -2.04652743e-13,
1.51822645e-06, 5.27884977e-13, -4.67908801e-12, -2.99049244e-06, -2.22795549e-09,
2.84124381e-12, -1.38621426e-12, 1.08651713e-06, 1.37362634e-11, 5.27884977e-13,
1.33558110e-06, 8.89315199e-12, 2.22812614e-09, -3.02818958e-06, 2.72793760e-11,
-4.18653074e-12, 8.91803023e-13, -2.89580040e-06, -4.67908801e-12, 8.89315199e-12,
3.46791765e-05, 1.04296605e-11, 8.69163359e-12, -3.66638958e-05, -2.98797305e-06,
2.23137564e-09, 2.83697844e-13, -2.99049244e-06, 2.22812614e-09, 1.04296605e-11,
7.14882744e-06, -2.95473689e-12, -6.77145994e-12, -2.22351831e-09, -2.96824894e-06,
7.15144431e-12, -2.22795549e-09, -3.02818958e-06, 8.69163359e-12, -2.95473689e-12,
7.14956898e-06, 1.59883079e-11, 3.40829809e-12, 5.97604638e-12, -3.66635112e-05,
2.84124381e-12, 2.72793760e-11, -3.66638958e-05, -6.77145994e-12, 1.59883079e-11,
8.46401151e-05 ]
</dynmat>
<ref_dynmat shape='(9, 9)'>
[ 1.31717371e-06, -7.25798794e-13, 6.22844237e-13, 1.31717914e-06, -7.25801781e-13,
-3.76082417e-12, -3.04047638e-06, 1.67538576e-12, 3.62174484e-12, -7.25798794e-13,
1.31717570e-06, 7.02167317e-13, -7.25801781e-13, 1.31718113e-06, -7.55927237e-13,
1.67538576e-12, -3.04048097e-06, 6.20477885e-14, 6.22844237e-13, 7.02167317e-13,
3.46719008e-05, 6.22846800e-13, 7.02170207e-13, -2.90273971e-06, -1.43773230e-12,
-1.62083643e-12, -3.66668359e-05, 1.31717914e-06, -7.25801781e-13, 6.22846800e-13,
1.31718456e-06, -7.25804769e-13, -3.76083964e-12, -3.04048889e-06, 1.67539266e-12,
3.62175975e-12, -7.25801781e-13, 1.31718113e-06, 7.02170207e-13, -7.25804769e-13,
1.31718655e-06, -7.55930348e-13, 1.67539266e-12, -3.04049348e-06, 6.20480439e-14,
-3.76082417e-12, -7.55927237e-13, -2.90273971e-06, -3.76083964e-12, -7.55930348e-13,
3.46722354e-05, 8.68123688e-12, 1.74493226e-12, -3.66672221e-05, -3.04047638e-06,
1.67538576e-12, -1.43773230e-12, -3.04048889e-06, 1.67539266e-12, 8.68123688e-12,
7.01843385e-06, -3.86734929e-12, -8.36019539e-12, 1.67538576e-12, -3.04048097e-06,
-1.62083643e-12, 1.67539266e-12, -3.04049348e-06, 1.74493226e-12, -3.86734929e-12,
7.01844445e-06, -1.43226996e-13, 3.62174484e-12, 6.20477885e-14, -3.66668359e-05,
3.62175975e-12, 6.20480439e-14, -3.66672221e-05, -8.36019539e-12, -1.43226996e-13,
8.46395618e-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 is similar to the restart file we created at the end of example 1.1.3, i.e., the one we used as input for normal mode phonon calculations except here is one additional dynmacila matric element named inp_dynamt that is the initial dynamical matrix used to define the normal modes along which displacements are performed.
We can visualize the normal modes as we did in example 1.1.3. The below example shows how to do it using VMD file as an example.
[11]:
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='nmd', #Options are: 'axsf','nmd','molden'
file_path='refdynmat')
# Deleting the nmodes object to finish printing the information into the file
del nmodes
This will create an nmd file named refdynmat.nmd cintaining the normal mode vectors.
[12]:
%%bash
cat refdynmat.nmd
title PyEPFD-nmodes
names O O C
coordinates 0.0 0.0 1.157132840007336 0.0 0.0 -1.1571312524755888 0.0 0.0 -1.587531747e-06
mode 1 -44347.6 -0.000245977 -0.00198492 7.05801e-05 -0.00495603 0.000405116 7.05801e-05 -0.00260101 -0.000789901 7.05801e-05
mode 2 -214174 -0.00109769 -0.0020006 0.000150957 0.000918096 -0.00449894 0.000150957 -8.97939e-05 -0.00324977 0.000150957
mode 3 329802 -0.00223329 -0.00398869 -6.44854e-05 0.0020469 0.00299984 -6.44854e-05 -9.31928e-05 -0.000494414 -6.44854e-05
mode 4 101100 -0.00483253 0.00239793 -4.45383e-05 6.49351e-05 -0.00038747 -4.45383e-05 -0.00238379 0.00100523 -4.45383e-05
mode 5 71131.4 -4.99698e-05 8.27297e-05 0.00352578 9.81602e-05 0.000234485 0.00352578 2.40954e-05 0.000158607 0.00352578
mode 6 0.00146652 0.00205646 0.00067057 -1.0886e-09 0.00205647 0.000670573 1.53405e-09 -0.00547883 -0.00178653 -5.93389e-10
mode 7 0.00146652 0.00067057 -0.00205646 7.61687e-11 0.000670573 -0.00205647 6.06901e-11 -0.00178653 0.00547883 -1.82309e-10
mode 8 0.000743305 6.50042e-10 2.16218e-10 0.00414051 6.50045e-10 2.16219e-10 -0.00414049 -1.73184e-09 -5.76047e-10 -2.34076e-08
mode 9 0.000422301 -2.3297e-10 -3.99146e-12 0.00216302 -2.32971e-10 -3.99148e-12 0.00216305 6.20679e-10 1.0634e-11 -0.00576275
This file can be visualized using VMD. If we choose fmt=‘axsf’ / ‘molden’ that would create .axsf / .molden file respectively. They can be visualized with XCrysden and Molden, respectively.
[ ]: