1.4.1 EPCE and ZPR using FP: Preparing the input files
After a normal-mode phonon calculation we need to post-process the data to compute electron-phonon coupling energies (EPCE), zero-point renormalization (ZPR) as well as electron-phonon renormalization at finite temperatures using a frozen phonon (harmonic) approach.
The first step is to compute the overlap integral matrix for degenerate bands. For the CO2 example, HOMO is 2-fold degenerate, therefore, we need to compute the overlap integral matrix between the positive and negative displaced coordinates along each normal mode.
For this purpose first we would import the mode_overlap class.
[1]:
from pyepfd.overlap import mode_overlap
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.0
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2023-05-09 20:36:41.809705
A CO2 molecule has 9 normal modes. We would compute the overlap matrix elements for all 9 modes. Although, LUMO is singly degenerate, just to show how we can include several degenerate space in a single calculation, we would include also LUMO (orbital index 9). This may take a few minutes on a laptop.
[2]:
overlap = mode_overlap(nmode=9,\
orbital_space=[[7,8],[9]],\
directory='../2_normal_mode_phonon/',
cube_prefix='wf')
# Deleting the overlap object to complete file writing
del overlap
[3]:
%%bash
ls -lrth *.overlap
-rw-rw-r-- 1 arpan arpan 899 May 9 20:38 orbital-0.overlap
-rw-rw-r-- 1 arpan arpan 401 May 9 20:38 orbital-1.overlap
We see that two overlap files are created. Let us look at the first file: orbital-0.overlap.
[4]:
%%bash
cat orbital-0.overlap
#Orbitals indices: [7, 8]
#normal mode index = 1
0.58352 -0.646224
-0.646224 -0.58352
#normal mode index = 2
-0.571502 -1.50436e-06
-2.04339e-06 -0.846159
#normal mode index = 3
-0.824038 -9.34836e-07
7.68441e-07 0.49162
#normal mode index = 4
-0.642178 -4.97481e-08
4.85663e-08 0.871273
#normal mode index = 5
-0.659458 3.49759e-06
2.89872e-06 0.872965
#normal mode index = 6
0.998597 0.000107763
0.000107612 0.981333
#normal mode index = 7
0.998597 -0.000107678
-0.000107666 0.981333
#normal mode index = 8
-0.978867 -0.196916
0.196916 -0.978867
#normal mode index = 9
0.185788 0.973269
0.973269 -0.185788
We see that there are 9 (2 x 2) overlap matrices. Each matrix element is the overlap integral between the negative and positive displaced structure along a normal mode. The diagonal elements are the overlap integrals between the 7th and 8th orbital indices. Off-diagonal elements are overlap integrals between 7th and 8th orbital indices.
This file is necessary to calculate the second derivative of band energies when degeneracies are present.
As expected orbital-1.overlap has only one orbital (9th) in the orbital space there are only one value for each mode. We do not need this file to compute the second derivative of LUMO, as it is singly degenerate.
[5]:
%%bash
cat orbital-1.overlap
#Orbitals indices: [9]
#normal mode index = 1
-0.979298
#normal mode index = 2
0.982993
#normal mode index = 3
0.981145
#normal mode index = 4
-0.986205
#normal mode index = 5
0.985829
#normal mode index = 6
-0.921256
#normal mode index = 7
0.921263
#normal mode index = 8
0.998094
#normal mode index = 9
0.999743
For large system with many normal modes the calculation of overlap matrices may take long. Therefore, this step is parallelized over number of modes using MPI4Py. In that case it is easier to run it on a cluster. For that purpose first create a file calc_overlap.py, that should read as:
[6]:
%%bash
cat calc_overlap.py
#!/usr/bin/env python3
from pyepfd.overlap import mode_overlap
overlap = mode_overlap(nmode=9,\
orbital_space=[[7,8],[9]],\
directory='../2_normal_mode_phonon/',
cube_prefix='wf')
Now we can make the file as an executable using the following command
[7]:
%%bash
chmod 755 calc_overlap.py
Now this can be run using mpirun as shown below, here I am using only 2 mpi process, but you can use n-number of MPI processes where n is less than equalto 0.5*total number of modes.
[1]:
%%bash
mpirun -np 2 python3 calc_overlap.py
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.0
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2023-05-09 20:39:14.453256
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.0
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2023-05-09 20:39:14.454862
[2]:
%%bash
ls -lrth *.overlap
-rw-rw-r-- 1 arpan arpan 899 May 9 20:40 orbital-0.overlap
-rw-rw-r-- 1 arpan arpan 401 May 9 20:40 orbital-1.overlap
Besides overlap matrices, we need total electronic energies for displaced coordinates and band(orbital) energies. The electronic energies can be extracted from the qbox outputs using grep tool as below.
[5]:
%%bash
grep '<etotal>' ../2_normal_mode_phonon/enmfd-?.r | awk '{print $2}' > etotal.dat
This would write the electronic energies sequentially in a file named etotal.dat. There are total 19 configurations, so there will be 19 numbers.
[6]:
%%bash
cat etotal.dat
-37.52801335
-37.52805984
-37.52805985
-37.52754654
-37.52754654
-37.52715354
-37.52715354
-37.52808833
-37.52808833
-37.52805442
-37.52805442
-37.52690255
-37.52690255
-37.52690255
-37.52690255
-37.52702210
-37.52700642
-37.52700281
-37.52700337
Next step is to extract the band (orbital) energies. We will do it in two steps using the script extract_eigval.sh that can be found at util/qbox_tools folder of the pyepfd distribution. First step is to put add the path of the util/qbox_tools in the PATH variable. Then we would extract the homo energies (in this case orbital indices 7 8) after extracting copy them into a file named homo.eigval.dat. Then we would extract the lumo energies and copy them into a file lumo.eigval.dat.
[8]:
%%bash
source ../../../env.sh
cd ../2_normal_mode_phonon/
rm -rf Eigenvalues #Removes old directory
extract_eigval.sh -prefix enmfd -seq 1 1 -orb 7 8
Note, here we have only one qbox output file enmfd-1.r. Therefore, our prefix would be enmfd while the range would be 1 to 1. So we used “-seq 1 1”. If we had 12 output files then we had to use “-seq 1 12”.
Our HOMO are orbitals 7 and 8, so the orbital range is 7 to 8. So we used “-orb 7 8” option.
This will create a file named orbital_kp_0_0_0.dat Let us look into that file.
[9]:
%%bash
cd ../2_normal_mode_phonon/
cat orbital_kp_0_0_0_is_0.dat
#Orbital-7 #Orbital-8
-10.509136120898 -10.509136120808
-10.509182319751 -10.509182319484
-10.509182131230 -10.509182130746
-10.501317018104 -10.501194144059
-10.501317355303 -10.501194482823
-10.498513779430 -10.497985889827
-10.498513685829 -10.497985796856
-10.510094365923 -10.508908852667
-10.510094389441 -10.508908878769
-10.508022276177 -10.507268484282
-10.508021558974 -10.507267763318
-10.501391720266 -10.494668036040
-10.501391575875 -10.494667892401
-10.501391574934 -10.494667891409
-10.501391589185 -10.494667905889
-10.495537912275 -10.495537911992
-10.522408878497 -10.522408878494
-10.494621364128 -10.494621364091
-10.494629033591 -10.494629033587
We see that the energies of orbitals 7 and 8 (in eV) are written. Now we will move this file to current directory and rename it homo.eigval.dat.
[10]:
%%bash
mv ../2_normal_mode_phonon/orbital_kp_0_0_0_is_0.dat homo.eigval.dat
Similarly we will extract the LUMO orbital energies and copy it to the current directory.
[11]:
%%bash
source ../../../env.sh
cd ../2_normal_mode_phonon/
extract_eigval.sh -prefix enmfd -seq 1 1 -orb 9 9
Directory Eigenvalues exists
File eigvals.enmfd-1 exists
Skipping extracting eigenvalues from enmfd-1.r
Directory Eigenvalues/Orbital_energy_time_evol_kp_0_0_0_is_0 exists
Just like before we move the orbital_kp_0_0_0.dat into current directory as lumo.eigval.dat
[12]:
%%bash
mv ../2_normal_mode_phonon/orbital_kp_0_0_0_is_0.dat lumo.eigval.dat
We also need the restart/checkpoint file enmfdphonon.xml that we prepared in example 2.3. Lets copy it to the current directory.
[14]:
%%bash
cp ../2_normal_mode_phonon/enmfdphonon.xml .
[15]:
%%bash
ls
1_epce_zpr.ipynb
2_epce_zpr.ipynb
calc_overlap.py
enmfdphonon.xml
epfd_out
etotal.dat
homo.eigval.dat
lumo.eigval.dat
orbital-0.overlap
orbital-1.overlap
We need the following files for EPCE and ZPR calculation. (1)etotal.dat –> Total electronic energies for all displaced configurations. (2)enmfdphonon.xml –> All information regarding the phonon calculation. (3)homo.eigval.dat –> HOMO/VBM (for CO2 2 fold degenerate, so band indices 7 & 8) energies for all displaced configurations. (4)lumo.eigval.dat –> LUMO/CBM energies for all displaced configurations. (5)orbital-0.overlap –> Overlap matrices for degenerate HOMO level. If LUMO were degenerate we would also need orbital-1.overlap.
Next step would be to post-process these files that would be the subject of example 4.2.
[ ]: