4.1 Stochastic Calculation of Band Gap Renormalization.
In this example we will learn to do a stochastic Monte Carlo calculation for the band gap renormalization. For this type of calculation we need a phonon restart/checkpoint xml file using a Cartesian phonon calulcation or Normal Mode phonon calculation.
This calulation is very similar to a normal-mode phonon calculation. First we generate the displacements and then using qbox we compute the total electronic energies and band energies for all displaced coordinates. Last step is very simple. We compute band gaps for each displaced coordinate average over all configurations. This would be our electron-phonon renormalized band gap at that particular temperature.
We need to import pyepfd_io module to read the input file and we need coord_util module to create the displaced coordinate using the Monte Carlo / Stochastic method. Lets import them first.
[1]:
from pyepfd.coord_util import *
from pyepfd.pyepfd_io import *
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.0
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2023-05-09 20:51:39.451446
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 sd_inp (stochastic displacement input) 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]:
sd_inp = read_pyepfd_info(file_path='../1_cartesian_phonon/fdphonon.xml')
Time spent on read_pyepfd_info class: 0.0004673004150390625 s.
Now we will define the input dynamical matrix for stochastic displacement. If refined dynamical matrix (i.e. after performing acoustic sum rule) exists, it would take that one, otherwise it would take dynamical matrix object.
[3]:
if sd_inp.ref_dynmatrix is not None:
inp_dynmat = sd_inp.ref_dynmatrix
else:
inp_dynmat = sd_inp.dynmatrix
Now we will initiate an xyz class in writing mode. Atoms are obtained from atom objects in sd_inp class
[4]:
sd_xyz = xyz(file_path='mcap.xyz',io='w',atoms=sd_inp.atoms)
Time spent on xyz class: 0.00021219253540039062 s.
Like displacing the atoms during a (E)NMFD phonon calculation, we instantiate an object named sdmoves (stochastic displacement moves) of ionic_mover class. Here we would use mode = ‘SD’, i.e. stochastic displacement.
Moreover, we need to provide an algorithm (algo keyword) for stochastic displacement. 6 algorithms are available: (1)oneshot (‘os’),
(2)oneshot with antithetic pair (‘osap’),
(3)one shot with random sample of signs (osr), this is equivalent to thermal line sampling, see https://journals.aps.org/prb/abstract/10.1103/PhysRevB.93.014302
Same like ‘osr’ but with antithetic pairs.
(5)Monte Carlo Sampling of Normal Modes (‘mc’)
Same as ‘mc’ but with antithetic pairs of each chosen random sample.
Usually algorithm with antithetic pairs are recommended. Below example we are choosing MCAP.
Here ngrid is the number of independent samples to be drawn for algorithms 3-6 described above. If an algorithm has antithetic pairs in it then actual numbers of configuration would be 2 times ngrid. In the below example chosen value of ngrid is 250 and so it would generate 500 configurations.
We also need to chose a temperature. In the below example it is 100 K
Therefore, the definition of sdmoves class would read as:
[5]:
sdmoves = ionic_mover( atoms = sd_inp.atoms,
opt_coord = sd_inp.coord,
mode = 'SD',
algo ='mcap',
ngrid = 250,
dynmat = inp_dynmat,
asr = sd_inp.asr,
temperature = 100)
Time spent on ionic_mover class: 0.02138376235961914 s.
Now all the displaced coordinates are stored within sdmoves object. Next step is to write them into an xyz file using the sd_xyz class we initiated before.
[6]:
for j in range(sdmoves.disp_coord.shape[1]):
sd_xyz.write(cell = sd_inp.cell,coord = sdmoves.disp_coord[:,j])
del sd_xyz
[7]:
%%bash
ls mcap.xyz
mcap.xyz
Note that the last step was a Monte Carlo simulation that involves pseudo random numbers. Therefore each time you run the sdmoves class it would create a different set of displaced coordinates. Only one OS and OSAP algorithms would produce the same set of coordinates because it does not involves any pseudo random numbers. Let us restart the kernel for completion of the file writing.
The next step is to prepare the qbox input files for computation of electronic properties, e.g., band energies for these configurations. This step is very similar to the previous examples 1.2 and 2.2. We will use xyz2qbox tool. In this case we do not need to save the wave functions (unlike example 2.2).
[8]:
%run ../../../utils/qbox_utils/xyz2qbox.py mcap.xyz
File mcap.xyz has 3 number of atoms in each frame and a total 500 number of frames.
No of species: 2
Present species are:
O_species
C_species
The file mcap.xyz is an i-PI output with pos_unit: angstrom and cell_unit: angstrom
Enter start frame number: (Default = 1)
Enter last frame number: (Default = 500)
Enter step frame: (Default = 1)
Enter qbox command for 1st iteration: (Default = ' randomize_wf, run -atomic_density 0 100 10')
Hint: more than one commands should be separated by a comma.
Enter qbox run command: (Default = ' run 0 60 10')
Hint: more than one commands should be separated by a comma.
Enter qbox plot command: (Default = ' None')
Hints:
(1) Example: plot options filename (without .cube extension)
(2) More than one commands should be separated by a comma.
Enter qbox spectrum command: (Default = ' None')
Hints:
(1) Example: spectrum options filename (without .dat extension)
Do you want to save sample for each configurations (y/n)? (Default = n)
Enter filename of qbox input: (Default = qbox.i) mcap-1.i
Enter xc: (Default = PBE)B3LYP
Enter wf_dyn: (Default = JD)
Enter ecut (Ry): (Default = 50.0)
Enter scf_tol (Ry): (Default = 1e-08) 1e-9
Enter nempty: (Default = 100) 8
Enter nspin (1/2): (Default = 1)
Enter delta_spin: (Default = None)
Enter net_charge: (Default = 0)
Enter pseudopotential: (Default = ONCV_PBE-1.0)
Warning!
Files:
O_ONCV_PBE-1.0.xml
C_ONCV_PBE-1.0.xml
must be present in the following pseudo-potential directory.
../pseudos/
This created a qbox input file with name mcap-1.i. Just like example 1.2, here also multiple input file can be (must for a large system) created. For that purpose the bash wrapper many_xyz2qbox.sh located at util/qbox_util folder can be used instead of running xyz2qbox.py multiple times. The next step is to submit those jobs and wait until all calculations finished properly.
Here we already computed the HOMO-LUMOs for a set of MC trajectory (saved_example.xyz) generated for the same temperature and same algorithm. The related qbox input file is saved_example-1.i and the output file is saved_example-1.r. Taking this output file as an example, I will show how to compute renormalized bad gap.
First we will extract the Kon-Sham eigenvalues for HOMO (band 7,8) and LUMO.
[9]:
%%bash
source ../../../env.sh
rm -rf Eigenvalues #Removing the previous existing Eigenvalues folder
extract_eigval.sh -prefix saved_example -seq 1 1 -orb 7 9
[10]:
%%bash
head -n 25 orbital_kp_0_0_0_is_0.dat
#Orbital-7 #Orbital-8 #Orbital-9
-10.487170364117 -10.477806462848 -1.023844849688
-10.491930996223 -10.476740223367 -1.056185101151
-10.485337121956 -10.475223504189 -1.081698326149
-10.465087173793 -10.443357728981 -1.267443801201
-10.472053419490 -10.465358595320 -0.970187929921
-10.512543172977 -10.512527824290 -0.898384711847
-10.527799679100 -10.517207132901 -0.923443428935
-10.478123219864 -10.462103393953 -1.148608321486
-10.339938214676 -10.339440662568 -0.954741724093
-10.423617068754 -10.401157571192 -1.244167491779
-10.532306461367 -10.529566959494 -0.894040197086
-10.409028400597 -10.399754091388 -1.076863441171
-10.477285216996 -10.465041191181 -1.139326536596
-10.448040183680 -10.431506465550 -1.111545941111
-10.478714643717 -10.467120526323 -1.114443749835
-10.491228021287 -10.489006381255 -0.958756414210
-10.487713978646 -10.474126658401 -0.999063507444
-10.486938627166 -10.481050436941 -1.025233757119
-10.421068833951 -10.403397269651 -1.193176789321
-10.375467136826 -10.371551646788 -1.016743984958
-10.403133578003 -10.399024563790 -0.977891719837
-10.499967797745 -10.497579430661 -0.966816043298
-10.454634021664 -10.454160120113 -1.003325292033
-10.528476645895 -10.524321592394 -0.906084170125
Generated orbital_kp_0_0_0.dat file lists energies for bands 7,8 (HOMO) and 9 (LUMO) for all 500 sampled configurations. We will compute the average HOMO energy for each configurations by taking the average energies of band 7 and 8. Then we compute the band gap for each configurations by subtracting the average HOMO energy from LUMO energy (band 9).
This can be easily done using an “awk” command.
[11]:
%%bash
awk '{print $3-($1+$2)/2}' orbital_kp_0_0_0_is_0.dat | sed '1d' > gap.dat
head -n 25 gap.dat
9.45864
9.42815
9.39858
9.18678
9.49852
9.61415
9.59906
9.3215
9.38495
9.16822
9.6369
9.32753
9.33184
9.32823
9.35847
9.53136
9.48186
9.45876
9.21906
9.35677
9.42319
9.53196
9.45107
9.62031
9.61705
To obtain the renormalized gap, the next step is to average over all 500 configurations.
[12]:
from pyepfd.epfd import mc_convergence
import matplotlib.pyplot as plt
co2_100K = mc_convergence(file_path='gap.dat',algo='mcap', usecols=(0))
mc = co2_100K.avg[:,0] #Our gap was in column 0 of the file
emc = 2*co2_100K.sd_mean[:,0] # Error is 2 times standard deviation of mean
plt.plot( [i+1 for i in range(len(mc))], mc, label='MC')
plt.fill_between([i+1 for i in range(len(mc))], mc + emc, mc - emc, alpha = 0.3)
plt.xlabel("# of independent samples")
plt.ylabel("Renormalized Gap (eV)")
plt.show()
/home/arpan/.local/lib/python3.8/site-packages/pyepfd/epfd.py:290: RuntimeWarning: invalid value encountered in double_scalars
sd_mean = np.sqrt( np.var(data_chunk, axis = axis)/float(data_chunk.shape[0] - 1) )
[ ]: