This tutorial can be downloaded link.

Intro Tutorial 10: Computing Excitation Energies Using Time-Dependent Density Functional Theory (TDDFT)

This tutorial shows how to compute neutral excitation energies of the negatively charged nitrogen-vacancy defect (NV\(^-\)) in diamond using linear-response time-dependent density functional theory (TDDFT) within the Tamm-Dancoff approximation. In TDDFT, the excitation energies are computed as the eigenvalues of the Liouville superoperator, with the corresponding eigenvectors representing the transition density. Tutorial 11 shows how to compute absorption spectra using TDDFT. Tutorial 12 shows how to compute analytical nuclear forces and perform excited state geometry optimization using TDDFT. More details about the TDDFT implementation in WEST can be found in Y. Jin, V. Yu, M. Govoni, A. Xu, and G. Galli, J. Chem. Theory Comput. 19, 8689–8705 (2023).

In this tutorial, the TDDFT calculations are conducted using both the semi-local PBE functional and the dielectric-dependent hybrid (DDH) functional.

Step 1: TDDFT calculation with the PBE functional

A TDDFT calculation with a local or semi-local functional contains three steps:

  • 1.1 Performing a mean-field calculation, using density functional theory (DFT), as the starting point of TDDFT

  • 1.2 Computing the neutral excitation energies of triplet excited states using spin-conserving TDDFT

  • 1.3 Computing the neutral excitation energies of singlet excited states using spin-flip TDDFT

Step 1.1: Mean-field starting point

We first perform the DFT calculation using the Quantum ESPRESSO code.

Download the following files to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/pbe_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 pbe_pw.in file, the input for the pw.x code:

[1]:
%%bash
cat pbe_pw.in
&CONTROL
calculation = 'scf'
prefix = 'pbe'
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 2 cores.

[ ]:
%%bash
mpirun -n 2 pw.x -i pbe_pw.in > pbe_pw.out

Step 1.2: Spin-conserving TDDFT calculation

Now we run the wbse.x executable to compute the excitation energies of triplet states of NV\(^-\) in diamond using spin-conserving TDDFT.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/pbe_sc_wbse.in

Let us inspect the pbe_sc_wbse.in file:

[2]:
%%bash
cat pbe_sc_wbse.in
input_west:
  qe_prefix: pbe
  west_prefix: pbe_sc
  outdir: ./

wbse_init_control:
  wbse_init_calculation: S
  solver: TDDFT

wbse_control:
  wbse_calculation: D
  n_liouville_eigen: 2
  trev_liouville: 0.00000001
  trev_liouville_rel: 0.000001
  l_preconditioning: True

The solver: TDDFT keyword instructs the code to perform a TDDFT calculation. 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.

We run wbse.x on 2 cores:

[ ]:
%%bash
mpirun -n 2 wbse.x -i pbe_sc_wbse.in > pbe_sc_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:

[ ]:
%%bash
mkdir -p pbe_sc.wbse.save
wget -N -q https://west-code.org/doc/training/nv_tddft/pbe_sc_summary.json -O pbe_sc.wbse.save/summary.json

The calculated excitation energys (in Rydberg) can be found in the file pbe_sc.wbse.save/summary.json.

[3]:
%%bash
cat pbe_sc.wbse.save/summary.json
{
  "plep": {
    "n_plep_eigen": 2,
    "eigenval": [
      0.14406126316112416,
      0.1441942044061179
    ]
  }
}

These two states are the \(m_S = +1\) sublevels of the \(^3E\) excited states of NV\(^-\) in diamond.

Step 1.3: Spin-flip TDDFT calculation

Now we run the wbse.x executable to compute the excitation energies of the singlet states of NV\(^-\) in diamond using spin-flip TDDFT.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/pbe_sf_wbse.in

Let us inspect the pbe_sf_wbse.in file:

[4]:
%%bash
cat pbe_sf_wbse.in
input_west:
  qe_prefix: pbe
  west_prefix: pbe_sf
  outdir: ./

wbse_init_control:
  wbse_init_calculation: S
  solver: TDDFT

wbse_control:
  wbse_calculation: D
  n_liouville_eigen: 4
  trev_liouville: 0.00000001
  trev_liouville_rel: 0.000001
  l_preconditioning: True
  l_spin_flip: True
  l_spin_flip_kernel: True
  l_spin_flip_alda0: False

The l_spin_flip: True and l_spin_flip_kernel: True keywords instruct the code to perform a spin-flip TDDFT calculation with the spin-flip kernel.

We run wbse.x on 2 cores:

[ ]:
%%bash
mpirun -n 2 wbse.x -i pbe_sf_wbse.in > pbe_sf_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:

[ ]:
%%bash
mkdir -p pbe_sf.wbse.save
wget -N -q https://west-code.org/doc/training/nv_tddft/pbe_sf_summary.json -O pbe_sf.wbse.save/summary.json

The calculated excitation energys (in Rydberg) can be found in the file pbe_sf.wbse.save/summary.json.

[5]:
%%bash
cat pbe_sf.wbse.save/summary.json
{
  "plep": {
    "n_plep_eigen": 4,
    "eigenval": [
      0.006752572854356471,
      0.043495765189432344,
      0.043678785251997566,
      0.10011245174283542
    ]
  }
}

The four states are the \(m_S = 0\) sublevel of the \(^3A_2\) state, the doubly degenerate \(^1E\) states, and the \(^1A_1\) state of NV\(^-\) in diamond.

Step 2: TDDFT calculation with the DDH functional

A TDDFT calculation with a hybrid functional contains five steps:

  • 2.1 Performing a self-consistent DFT calculation as the starting point of TDDFT

  • 2.2 Performing a non-self-consistent DFT including more empty bands to construct the adaptively compressed exchange (ACE) operator, a low-rank approximation to the exact exchange operator

  • 2.3 Computing the Coulomb integrals for pairs of occupied Kohn-Sham orbitals

  • 2.4 Computing the neutral excitation energies of triplet excited states using spin-conserving TDDFT

  • 2.5 Computing the neutral excitation energies of singlet excited states using spin-flip TDDFT

Step 2.1: Mean-field starting point

We first perform the DFT calculation using the Quantum ESPRESSO code.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_pw.in

We can now inspect the ddh_pw.in file, the input for the pw.x code:

[6]:
%%bash
cat ddh_pw.in
&CONTROL
calculation = 'scf'
prefix = 'ddh'
pseudo_dir = './'
/
&SYSTEM
ibrav = 0
ecutwfc = 50
tot_charge = -1.0
nspin = 2
tot_magnetization = 2.0
nbnd = 150
nat = 63
ntyp = 2
input_dft = 'pbe0'
ecutfock = 100.0
exx_fraction = 0.18
/
&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

Compared to the file pbe_pw.in, the additional keywords input_dft = 'pbe0' and exx_fraction = 0.18 instruct the DFT calculation to use the PBE0 functional, incorporating an exact exchange fraction of 0.18. This approach is termed the dielectric-dependent hybrid (DDH) functional, as the fraction is derived from the inverse of the macroscopic dielectric constant of diamond.

We run pw.x on 2 cores, and make a copy of the save directory generated by pw.x.

[ ]:
%%bash
mpirun -n 2 pw.x -i ddh_pw.in > ddh_pw.out
cp -r ddh.save ddh_copy.save

Step 2.2: Non-self-consistent-field calculation for the ACE operator

We perform the non-self-consistent-field (nscf) calculation for the ACE operator that is used in the TDDFT calculation in the case of hybrid functionals.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_nscf.in

We can now inspect the ddh_nscf.in file:

[7]:
%%bash
cat ddh_nscf.in
&CONTROL
calculation = 'nscf'
prefix = 'ddh'
pseudo_dir = './'
/
&SYSTEM
ibrav = 0
ecutwfc = 50
tot_charge = -1.0
nspin = 2
tot_magnetization = 2.0
nbnd = 512
nat = 63
ntyp = 2
input_dft = 'pbe0'
ecutfock = 100.0
exx_fraction = 0.18
/
&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 content of ddh_nscf.in is almost identical to that of ddh_pw.in, except that in ddh_nscf.in we request a non-self-consistent-field (nscf) calculation with a larger value of nbnd = 512, which is the number of Kohn-Sham orbitals used to construct the ACE operator.

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 again on 2 cores, and copy the file that stores the ACE operator to the directory saved in the previous self-consistent run.

[ ]:
%%bash
mpirun -n 2 pw.x -i ddh_nscf.in > ddh_nscf.out
cp ddh.save/ace* ddh_copy.save/
rm -rf ddh.save ddh.wfc*
mv ddh_copy.save ddh.save

Step 2.3: Coulomb integrals

We compute the Coulomb integrals for pairs of occupied Kohn-Sham orbitals using the wbse_init.x executable.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_wbse_init.in

Let us inspect the ddh_wbse_init.in file:

[8]:
%%bash
cat ddh_wbse_init.in
input_west:
  qe_prefix: ddh
  west_prefix: ddh_west
  outdir: ./

wbse_init_control:
  wbse_init_calculation: S
  solver: TDDFT
  localization: W
  overlap_thr: 0.001

The localization: W keyword instructs the code to use Wannier functions, which are a localized representation of the Kohn-Sham orbitals. If the overlap between two Wannier functions is below overlap_thr: 0.001, the evaluation of the corresponding Coulomb integral is skipped, thus reducing the computational cost.

We run wbse_init.x on 2 cores:

[ ]:
%%bash
mpirun -n 2 wbse_init.x -i ddh_wbse_init.in > ddh_wbse_init.out

Step 2.4: Spin-conserving TDDFT calculation

Now we run the wbse.x executable to compute the excitation energies of triplet states of NV\(^-\) in diamond using spin-conserving TDDFT.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_sc_wbse.in

Let us inspect the ddh_sc_wbse.in file:

[9]:
%%bash
cat ddh_sc_wbse.in
input_west:
  qe_prefix: ddh
  west_prefix: ddh_sc
  outdir: ./

wbse_init_control:
  wbse_init_calculation: S
  solver: TDDFT
  localization: W
  overlap_thr: 0.001

wbse_control:
  wbse_calculation: D
  n_liouville_eigen: 2
  trev_liouville: 0.00000001
  trev_liouville_rel: 0.000001
  l_preconditioning: True
  n_exx_lowrank: 512

The n_exx_lowrank: 512 keyword means that the ACE operator is constructed with 512 Kohn-Sham orbitals.

We run wbse.x on 2 cores:

[ ]:
%%bash
ln -s ddh_west.wbse_init.save ddh_sc.wbse_init.save
mpirun -n 2 wbse.x -i ddh_sc_wbse.in > ddh_sc_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:

[ ]:
%%bash
mkdir -p ddh_sc.wbse.save
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_sc_summary.json -O ddh_sc.wbse.save/summary.json

The calculated excitation energys (in Rydberg) can be found in the file ddh_sc.wbse.save/summary.json.

[10]:
%%bash
cat ddh_sc.wbse.save/summary.json
{
  "plep": {
    "n_plep_eigen": 2,
    "eigenval": [
      0.16518277920006938,
      0.16529209364139377
    ]
  }
}

These two states are the \(m_S = +1\) sublevel of the \(^3E\) excited states of NV\(^-\) in diamond.

Step 2.5: Spin-flip TDDFT calculation

Now we run the wbse.x executable to compute the excitation energies of singlet states of NV\(^-\) in diamond using spin-flip TDDFT.

Download the following file to your working directory:

[ ]:
%%bash
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_sf_wbse.in

Let us inspect the ddh_sf_wbse.in file:

[11]:
%%bash
cat ddh_sf_wbse.in
input_west:
  qe_prefix: ddh
  west_prefix: ddh_sf
  outdir: ./

wbse_init_control:
  wbse_init_calculation: S
  solver: TDDFT
  localization: W
  overlap_thr: 0.001

wbse_control:
  wbse_calculation: D
  n_liouville_eigen: 4
  trev_liouville: 0.00000001
  trev_liouville_rel: 0.000001
  l_preconditioning: True
  l_spin_flip: True
  l_spin_flip_kernel: True
  l_spin_flip_alda0: False
  n_exx_lowrank: 512

We run wbse.x on 2 cores:

[ ]:
%%bash
ln -s ddh_west.wbse_init.save ddh_sf.wbse_init.save
mpirun -n 2 wbse.x -i ddh_sf_wbse.in > ddh_sf_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:

[ ]:
%%bash
mkdir -p ddh_sf.wbse.save
wget -N -q https://west-code.org/doc/training/nv_tddft/ddh_sf_summary.json -O ddh_sf.wbse.save/summary.json

The calculated excitation energys (in Rydberg) can be found in the file ddh_sf.wbse.save/summary.json.

[12]:
%%bash
cat ddh_sf.wbse.save/summary.json
{
  "plep": {
    "n_plep_eigen": 4,
    "eigenval": [
      0.010238564822618666,
      0.05862318898692176,
      0.058865528147033344,
      0.1472998272595149
    ]
  }
}

The four states are the \(m_S = 0\) sublevel of the \(^3A_2\) state, the doubly degenerate \(^1E\) states, and the \(^1A_1\) state of NV\(^-\) in diamond.

Step 3: Summary of results

We extract the excitation energies from the JSON files.

[13]:
import json
import numpy as np

ry2ev = 13.605662285137

with open("pbe_sc.wbse.save/summary.json", "r") as f:
    data = json.load(f)
    pbe_sc_vee = np.array(data["plep"]["eigenval"])
    pbe_sc_vee *= ry2ev

with open("pbe_sf.wbse.save/summary.json", "r") as f:
    data = json.load(f)
    pbe_sf_vee = np.array(data["plep"]["eigenval"])
    pbe_sf_vee -= pbe_sf_vee[0]
    pbe_sf_vee *= ry2ev

with open("ddh_sc.wbse.save/summary.json", "r") as f:
    data = json.load(f)
    ddh_sc_vee = np.array(data["plep"]["eigenval"])
    ddh_sc_vee *= ry2ev

with open("ddh_sf.wbse.save/summary.json", "r") as f:
    data = json.load(f)
    ddh_sf_vee = np.array(data["plep"]["eigenval"])
    ddh_sf_vee -= ddh_sf_vee[0]
    ddh_sf_vee *= ry2ev

Now we compare the computed excitation energies with the experimental zero-phonon line (ZPL):

[14]:
import pandas as pd

data = {
    "State": ["$^3A_2$", "$^1E$", "$^1A_1$", "$^3E$"],
    "TDDFT@PBE (eV)": [0.0, pbe_sf_vee[1], pbe_sf_vee[3], pbe_sc_vee[0]],
    "TDDFT@DDH (eV)": [0.0, ddh_sf_vee[1], ddh_sf_vee[3], ddh_sc_vee[0]],
    "Expt. ZPL (eV)": [0.0, "0.34 - 0.43", "1.51 - 1.60", 1.945]
}

df = pd.DataFrame(data)
pd.options.display.float_format = "{:.3f}".format
df
[14]:
State TDDFT@PBE (eV) TDDFT@DDH (eV) Expt. ZPL (eV)
0 $^3A_2$ 0.000 0.000 0.000
1 $^1E$ 0.500 0.658 0.34 - 0.43
2 $^1A_1$ 1.270 1.865 1.51 - 1.60
3 $^3E$ 1.960 2.247 1.945

The TDDFT@DDH calculations employed two numerical approximations:

  • the ACE approximation to the exact exchange operator;

  • the nearsightedness principle to restrict the Coulomb integrals to overlapping Wannier functions.

These approximations introduce minimal errors in the calculated excitation energies, as demonstrated below.

[16]:
data = {
    "State": ["$^3A_2$", "$^1E$", "$^1A_1$", "$^3E$"],
    "TDDFT@DDH, w/ approximations (eV)": [0.0, ddh_sf_vee[1], ddh_sf_vee[3], ddh_sc_vee[0]],
    "TDDFT@DDH, w/o approximations (eV)": [0.0, 0.659, 1.866, 2.250],
}

df = pd.DataFrame(data)
pd.options.display.float_format = "{:.3f}".format
df
[16]:
State TDDFT@DDH, w/ approximations (eV) TDDFT@DDH, w/o approximations (eV)
0 $^3A_2$ 0.000 0.000
1 $^1E$ 0.658 0.659
2 $^1A_1$ 1.865 1.866
3 $^3E$ 2.247 2.250