This tutorial can be downloaded link.
Intro Tutorial 13: Computing Excitation Energies by Solving the Bethe-Salpeter Equation (BSE) Using Density Matrix Perturbation Theory (DMPT)¶
This tutorial shows how to compute neutral excitation energies of the negatively charged nitrogen-vacancy defect (NV\(^-\)) in diamond by solving the Bethe-Salpeter Equation (BSE) using the density matrix perturbation theory (DMPT). In BSE, the excitation energies are computed as the eigenvalues of the Liouville superoperator, with the corresponding eigenvectors representing the transition density.
WEST implements two distinct methods to compute the dielectric screening for the BSE solver:
the projective dielectric eigenpotentials (PDEP) technique, which builds a low rank representation of the static dielectric matrix that enters the screened Coulomb interaction, the random phase approximation (RPA) is used;
the finite field (FF) method, which allows for efficient calculations within and beyond the RPA.
Both methods circumvent the explicit computation of empty electronic states and the storage and inversion of large matrices. In addition, the computational cost can be reduced by using a localized representation of the ground state Kohn-Sham wavefunctions.
This tutorial focuses on the PDEP-based method as described in Rocca et al., J. Chem. Phys. 133, 164109 (2010) and Rocca et al., Phys. Rev. B 85, 045116 (2012). The FF-based method is covered in Tutorial 8. Tutorial 7 shows how to compute absorption spectra by solving the BSE.
Step 1: Mean-field starting point¶
We first perform the mean-field electronic structure calculation within DFT using the Quantum ESPRESSO code.
Download the following files to your working directory:
[1]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_bse/pw.in
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/C_ONCV_PBE-1.2.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/N_ONCV_PBE-1.2.upf
We can now inspect the pw.in file, the input for the pw.x code:
[3]:
%%bash
cat pw.in
&CONTROL
calculation = 'scf'
pseudo_dir = './'
/
&SYSTEM
ibrav = 0
ecutwfc = 50
tot_charge = -1.0
nspin = 2
tot_magnetization = 2.0
nbnd = 150
nat = 63
ntyp = 2
/
&ELECTRONS
diago_full_acc = .true.
/
K_POINTS gamma
CELL_PARAMETERS angstrom
7.136000 0.000000 0.000000
0.000000 7.136000 0.000000
0.000000 0.000000 7.136000
ATOMIC_SPECIES
C 12.01099968 C_ONCV_PBE-1.2.upf
N 14.00699997 N_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS crystal
C 0.1269815510 0.3733371195 0.3775718026
C 0.1253081282 0.3763002936 0.8746894487
C 0.1260060117 0.8738380396 0.3754459361
C 0.1250639790 0.8750661683 0.8749376784
C 0.6248201544 0.3811585813 0.3751813598
C 0.6224309085 0.3733360188 0.8730196975
C 0.6250110086 0.8765296624 0.3749876809
C 0.6245505450 0.8738390861 0.8739957865
C 0.0000259749 0.0000259633 -0.0000256868
C 0.0014358428 -0.0014358100 0.5010495568
C 0.0014328181 0.4989469740 0.0014351388
C 0.0012018772 0.4992148548 0.5007857943
C 0.4989470397 -0.0014328307 0.0014351823
C 0.4992149137 0.0012018289 0.5007857375
C 0.4992116247 0.4992116312 -0.0012034708
C 0.1250667730 0.1250667914 0.1249294365
C 0.1253032977 0.1253032524 0.6237019642
C 0.1260062437 0.6245495725 0.1261632799
C 0.1269771899 0.6224340478 0.6266613449
C 0.6245496185 0.1260062840 0.1261633007
C 0.6224340177 0.1269771772 0.6266613352
C 0.6250092594 0.6250092159 0.1234699036
C 0.6248234427 0.6248234639 0.6188442442
C 0.0009910041 0.2504149699 0.2494871226
C 0.0003716741 0.2501475881 0.7498513998
C 0.0009645933 0.7495989686 0.2503986681
C 0.0009941125 0.7505129665 0.7495871392
C 0.4985394814 0.2524198182 0.2522377911
C 0.5078677486 0.2391591291 0.7608354182
C 0.5003280630 0.7516382226 0.2483618688
C 0.4985401685 0.7477659930 0.7475823593
C 0.3748654285 0.3748654245 0.1250552701
C 0.3749783337 0.8749127142 0.1250879590
C 0.3748650726 0.8749450711 0.6251340225
C 0.8749127198 0.3749782927 0.1250879634
C 0.8749450684 0.3748650501 0.6251340257
C 0.8750493463 0.8750493517 0.1249501816
C 0.8749120516 0.8749120381 0.6250201881
C 0.2504149894 -0.0009909944 0.2494870713
C 0.2501475446 0.0003715462 0.7498514576
C 0.2524198093 0.4985394738 0.2522377661
C 0.2391592245 0.5078678650 0.7608354474
C 0.7495989932 0.0009646155 0.2503986775
C 0.7505130382 -0.0009940706 0.7495871127
C 0.7516382084 0.5003280124 0.2483619045
C 0.7477659502 0.4985401253 0.7475823293
C 0.3733371729 0.1269814887 0.3775717513
C 0.3763002438 0.1253079989 0.8746895222
C 0.3811585422 0.6248202431 0.3751812920
C 0.3733360636 0.6224310126 0.8730197732
C 0.8738379990 0.1260060464 0.3754459661
C 0.8750662653 0.1250640377 0.8749376416
C 0.8765296848 0.6250109436 0.3749877093
C 0.8738390803 0.6245504485 0.8739957771
C 0.2501515174 0.2501515049 -0.0003784384
C 0.2391603565 0.2391602431 0.4921393755
C 0.2504126363 0.7505106105 0.0009953037
C 0.2524180055 0.7477644305 0.5014586422
C 0.7505106094 0.2504126558 0.0009952558
C 0.7477643646 0.2524180264 0.5014586712
C 0.7495988945 0.7495988379 -0.0009639807
C 0.7516394023 0.7516393858 0.4996700487
N 0.5130505433 0.5130506248 0.4869539927
The calculation is spin polarized (i.e., nspin = 2 and tot_magnetization = 2.0), representing the \(m_S = +1\) sublevel of the \(^3A_2\) many-body state of NV\(^-\) in diamond.
We run pw.x on 32 cores.
[ ]:
%%bash
mpirun -n 32 pw.x -nk 2 -i pw.in > pw.out
Step 2: Calculation of dielectric screening¶
As detailed in Tutorial 1, the static dielectric screening is computed using the projective dielectric eigendecomposition (PDEP) technique.
Download the following file to your working directory:
[5]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_bse/wstat.in
Let us inspect the wstat.in file:
[8]:
%%bash
cat wstat.in
wstat_control:
wstat_calculation: S
n_pdep_eigen: 512
trev_pdep: 0.00001
We run wstat.x on 512 cores:
[ ]:
%%bash
mpirun -n 512 wstat.x -ni 16 -nk 2 -i wstat.in > wstat.out
Step 3: Calculation of quasiparticle corrections¶
The GW electronic structure is computed treating the frequency integration of the correlation part of the self energy with the contour deformation technique and by computing the dielectric screening at multipole frequencies with Lanczos iterations.
Download the following file to your working directory:
[10]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_bse/wfreq.in
Let us inspect the wfreq.in file:
[13]:
%%bash
cat wfreq.in
wstat_control:
wstat_calculation: S
n_pdep_eigen: 512
trev_pdep: 0.00001
wfreq_control:
wfreq_calculation: XWGQO
macropol_calculation: C
n_pdep_eigen_to_use: 512
qp_bandrange: [121, 140]
The wfreq_calculation: XWGQO keyword instructs the code to compute quasiparticle corrections (XWGQ) for bands specified in the qp_bandrange: [121, 140] keyword, and to compute optical properties (XWO) including the dielectric function.
We run wfreq.x on 512 cores:
[ ]:
%%bash
mpirun -n 512 wfreq.x -ni 4 -nk 2 -nb 4 -i wfreq.in > wfreq.out
If the reader does NOT have the computational resources to run the calculation, the WEST output file needed for the next step can be directly downloaded as:
[15]:
%%bash
mkdir -p west.wfreq.save
wget -N -q https://west-code.org/doc/training/nv_bse/optics.json -O west.wfreq.save/optics.json
The quasiparticle corrections and optical properties are stored in the files west.wfreq.save/wfreq.json and west.wfreq.save/optics.json, respectively. From the latter, we can extract the macroscopic dielectric constant.
[17]:
import json
with open('west.wfreq.save/optics.json') as json_file :
data = json.load(json_file)
eps_infty = data['eps1'][0]
print(eps_infty)
7.750227417643485
Step 4: BSE calculation¶
Step 4.1: BSE initialization¶
We perform an initialization step to compute the screened exchange integrals using the wbse_init.x executable.
Download the following file to your working directory:
[19]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_bse/wbse_init.in
Let us inspect the wbse_init.in file:
[22]:
%%bash
cat wbse_init.in
input_west:
outdir: ./
wbse_init_control:
wbse_init_calculation: S
bse_method: PDEP
n_pdep_eigen_to_use: 512
localization: W
overlap_thr: 0.001
The bse_method: PDEP keyword selects the PDEP method. A total number of n_pdep_eigen_to_use: 512 PDEPs is used to represent the static dielectric matrix. The localization: W keyword instructs the code to use Wannier functions, which are a localized representation of the Kohn-Sham wavefunctions. If the overlap between two Wannier functions is below overlap_thr: 0.001, the evaluation of the corresponding screened exchange integral is skipped, thus reducing the computational cost.
We run wbse_init.x on 512 cores:
[ ]:
%%bash
mpirun -n 512 wbse_init.x -ni 16 -nk 2 -i wbse_init.in > wbse_init.out
Step 4.2: BSE excitation energies¶
Now we run the wbse.x executable to compute the lowest excitation energies of NV\(^-\).
Download the following file to your working directory:
[24]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_bse/wbse.in
Let us inspect the wbse.in file:
[27]:
%%bash
cat wbse.in
input_west:
outdir: ./
wbse_init_control:
wbse_init_calculation: S
bse_method: PDEP
n_pdep_eigen_to_use: 512
localization: W
overlap_thr: 0.001
wbse_control:
wbse_calculation: D
qp_correction: west.wfreq.save/wfreq.json
wbse_epsinfty: 7.75
n_liouville_eigen: 2
The wbse_calculation: D keyword instructs the code to compute the neutral excitation energies using the Davidson algorithm. The n_liouville_eigen: 2 keyword specifies that two lowest excitation energies are computed. The qp_correction keyword indicates the name of the JSON file that contains the quasiparticle correction computed by wfreq.x. The wbse_epsinfty keyword specifies the macroscopic relative dielectric constant of the material. It can be the experimental value if
known, or in this case, we use the value computed by wfreq.x.
We run wbse.x on 512 cores:
[ ]:
%%bash
mpirun -n 512 wbse.x -nk 2 -nb 16 -i wbse.in > wbse.out
If the reader does NOT have the computational resources to run the calculation, the WEST output file needed for the next step can be directly downloaded as:
[29]:
%%bash
mkdir -p west.wbse.save
wget -N -q https://west-code.org/doc/training/nv_bse/summary.json -O west.wbse.save/summary.json
The calculated excitation energys (in Rydberg) can be found in the file west.wbse.save/summary.json.
[32]:
%%bash
cat west.wbse.save/summary.json
{
"plep": {
"n_plep_eigen": 2,
"eigenval": [
0.16843124342116073,
0.1685732039135792
]
}
}
Step 5: Analysis of excited states¶
Finally, we run the westpp.x executable to analyze the composition of the excited states.
Download the following file to your working directory:
[34]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_bse/westpp.in
Let us inspect the westpp.in file:
[37]:
%%bash
cat westpp.in
input_west:
outdir: ./
westpp_control:
westpp_calculation: C
westpp_n_liouville_to_use: 2
westpp_range: [1, 2]
We run westpp.x on 32 cores:
[ ]:
%%bash
mpirun -n 32 westpp.x -i westpp.in > westpp.out
The westpp_calculation: C keyword instructs the code to perform a decomposition of the excited states into transitions from occupied to empty Kohn-Sham wavefunctions. Here we consider two excited states as specified by the westpp_range: [1, 2] keyword.
The output file westpp.out would include the following:
*-------------* THE PRINCIPLE PROJECTED COMPONENTS *-------------*
# Exciton : 1 | Excitation energy : 0.168431
# Transition_from | Transition_to | Coeffcient
2 126 | 2 127 | 0.905626
2 126 | 2 128 | 0.406721
# Exciton : 2 | Excitation energy : 0.168573
# Transition_from | Transition_to | Coeffcient
2 126 | 2 127 | 0.406721
2 126 | 2 128 | -0.905573
The first excited state has an excitation energy of 0.168431 Ry. This excitation has a major contribution by a transition from the 126th band (a\(_1\) defect orbital) in the second spin channel to the 127th band (one of two degenerate e defect orbitals) in the same spin channel, and a minor contribution by a transition from the 126th band (a\(_1\) defect orbital) in the second spin channel to the 128th band (the other e defect orbital) in the same spin channel. The second excited
state has an excitation energy of 0.168573 Ry. These two excited states are the \(m_S = +1\) sublevel of the \(^3E\) excited states of NV\(^-\) in diamond. By using westpp.x and westpp_calculation: X, we can generate cube files of the excited states, which can then be visualized with any software capable of plotting cube files.
We computed spin-conserving excitations in this tutorial. We can compute and analyze spin-flip excitations following the same steps. An example is given in Tutorial 10 for time-dependent density functional theory (TDDFT), which can be adapted to BSE.