This tutorial can be downloaded link.

Advanced Tutorial 3: Embedded Bethe-Salpeter Equation (eBSE) Calculations¶

In this tutorial, we show how to perform embedded BSE (eBSE) calculations with the WEST code. Within eBSE, we formulate an effective Hamiltonian for a selected number of orbitals in the form

\[H^{\mathrm{eBSE}}_{vc, v'c'} = \left( \epsilon^{GW}_{c} - \epsilon^{GW}_{v} \right)\delta_{vv'}\delta_{cc'} - W_{cv, c'v'} + V_{cv, c'v'},\]

where \(v\) and \(v'\) include all occupied states in the active space, \(c\) and \(c'\) all unoccupied states. \(\epsilon^{GW}\) are the quasiparticle energies in the \(\mathrm{G_0W_0}\) approximation.

Matrix elements of the direct interaction \(W_{cv, c'v'}\) are defined as

\[w_{cv,c'v'} = \int d^3\mathbf{r} d^3\mathbf{r}' \phi_{v}(\mathbf{r})\phi_{v'}(\mathbf{r})W_0(\mathbf{r},\mathbf{r}',\omega=0) \phi_{c}(\mathbf{r}')\phi_{c'}(\mathbf{r}'),\]

where \(W_0(\mathbf{r},\mathbf{r}',\omega)\) is the screened Coulomb potential in the random-phase approximation (RPA). These matrix elements describe the Coulomb attraction between the excited electrons and the valence holes.

Matrix elements of the repulsive exchange interaction \(V_{cv, c'v'}\) are given by

\[V_{cv,c'v'} = \int d^3\mathbf{r} d^3\mathbf{r}' \phi_{v}(\mathbf{r})\phi_{c}(\mathbf{r})W^R_0(\mathbf{r},\mathbf{r}',\omega=0) \phi_{v'}(\mathbf{r}')\phi_{c'}(\mathbf{r}'),\]

where \(W^R_0(\mathbf{r},\mathbf{r}',\omega)\) is the Coulomb potential in the constrained random-phase approximation (cRPA).

The eigenvalues and -states of the eBSE Hamiltonian yield the neutral excitations of the system, i.e.

\[H^{\mathrm{eBSE}}X^\lambda = E^\lambda X^\lambda,\]

where \(E^\lambda\) are the neutral excitation energies, and \(X^\lambda\) are the corresponding excited many-body states.

eBSE calculations are performed in two separate steps: 1. In an initial WEST calculation, the quasiparticle energies and matrix elements \(W\) and \(V\) are calculated and stored in a JSON file. 2. The eBSE Hamiltonian is constructed from the parameters in the JSON file and diagonalized. This step is performed by the WESTpy Python package.

Step 1: Parameters of the eBSE Hamiltonian¶

Step 1.1: Mean-field starting point¶

As a starting point for the many-body perturbation theory (MBPT) calculation, we perform a density functional theory (DFT) calculation using Quantum ESPRESSO.

Note that embedded BSE calculations require spin-polarized DFT calculations as a starting point. This is distinct from Quantum Defect Embedding Theory (QDET), which generally starts from spin-unpolarized DFT calculations.

Download the following files in your working directory:

[1]:

%%bash
wget -N -q https://west-code.org/doc/training/nv_diamond_63_spinpol/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

Let us inspect the pw.in file:

[2]:

%%bash
cat pw.in

&CONTROL
calculation = 'scf'
pseudo_dir = './'
/
&SYSTEM
input_dft = 'PBE'
ibrav = 0
ecutwfc = 50
tot_charge = -1
nspin = 2
tot_magnetization = 2
nbnd = 176
nat = 63
ntyp = 2
/
&ELECTRONS
diago_full_acc = .true.
/
K_POINTS gamma
CELL_PARAMETERS angstrom
7.136012  0.000000  0.000000
0.000000  7.136012  0.000000
0.000000  0.000000  7.136012
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.99996000  0.99996000  0.99996000
C    0.12495000  0.12495000  0.12495000
C    0.99905000  0.25039000  0.25039000
C    0.12350000  0.37499000  0.37499000
C    0.25039000  0.99905000  0.25039000
C    0.37499000  0.12350000  0.37499000
C    0.25039000  0.25039000  0.99905000
C    0.37499000  0.37499000  0.12350000
C    0.00146000  0.00146000  0.50100000
C    0.12510000  0.12510000  0.62503000
C    0.00102000  0.24944000  0.74960000
C    0.12614000  0.37542000  0.87402000
C    0.24944000  0.00102000  0.74960000
C    0.37542000  0.12614000  0.87402000
C    0.24839000  0.24839000  0.49966000
C    0.37509000  0.37509000  0.61906000
C    0.00146000  0.50100000  0.00146000
C    0.12510000  0.62503000  0.12510000
C    0.00102000  0.74960000  0.24944000
C    0.12614000  0.87402000  0.37542000
C    0.24839000  0.49966000  0.24839000
C    0.37509000  0.61906000  0.37509000
C    0.24944000  0.74960000  0.00102000
C    0.37542000  0.87402000  0.12614000
C    0.99883000  0.50076000  0.50076000
C    0.12502000  0.62512000  0.62512000
C    0.99961000  0.74983000  0.74983000
C    0.12491000  0.87493000  0.87493000
C    0.25216000  0.50142000  0.74767000
C    0.37740000  0.62659000  0.87314000
C    0.25216000  0.74767000  0.50142000
C    0.37740000  0.87314000  0.62659000
C    0.50100000  0.00146000  0.00146000
C    0.62503000  0.12510000  0.12510000
C    0.49966000  0.24839000  0.24839000
C    0.61906000  0.37509000  0.37509000
C    0.74960000  0.00102000  0.24944000
C    0.87402000  0.12614000  0.37542000
C    0.74960000  0.24944000  0.00102000
C    0.87402000  0.37542000  0.12614000
C    0.50076000  0.99883000  0.50076000
C    0.62512000  0.12502000  0.62512000
C    0.50142000  0.25216000  0.74767000
C    0.62659000  0.37740000  0.87314000
C    0.74983000  0.99961000  0.74983000
C    0.87493000  0.12491000  0.87493000
C    0.74767000  0.25216000  0.50142000
C    0.87314000  0.37740000  0.62659000
C    0.50076000  0.50076000  0.99883000
C    0.62512000  0.62512000  0.12502000
C    0.50142000  0.74767000  0.25216000
C    0.62659000  0.87314000  0.37740000
C    0.74767000  0.50142000  0.25216000
C    0.87314000  0.62659000  0.37740000
C    0.74983000  0.74983000  0.99961000
C    0.87493000  0.87493000  0.12491000
N    0.48731000  0.48731000  0.48731000
C    0.49191000  0.76093000  0.76093000
C    0.62368000  0.87476000  0.87476000
C    0.76093000  0.49191000  0.76093000
C    0.87476000  0.62368000  0.87476000
C    0.76093000  0.76093000  0.49191000
C    0.87476000  0.87476000  0.62368000

We can now run pw.x on 2 cores.

[ ]:

%%bash
mpirun -n 2 pw.x -i pw.in > pw.out

Step 1.2: Calculation of the dielectric screening¶

As for GW and QDET calculations within WEST, we first have to determine the dielectric screening before we can proceed to calculate the excitations of the spin defect. In WEST, the dielectric screening is obtained from the projective dielectric eigendecomposition (PDEP) technique. The calculation with wstat.x requires an input file wstat.in.

[3]:

%%bash
wget -N -q https://west-code.org/doc/training/nv_diamond_63_spinpol/wstat.in

Once again, we can have a look at the input file:

[4]:

%%bash
cat wstat.in

wstat_control:
  wstat_calculation: S
  n_pdep_eigen: 512
  trev_pdep: 0.00001

As we can see, there are no input parameters in wstat.in specific to QDET. We can now execute wstat.x with the following command.

[ ]:

%%bash
mpirun -n 2 wstat.x -i wstat.in > wstat.out

Step 1.3: eBSE calculation¶

The calculation of matrix elements for eBSE within WEST is identical to that of matrix elements for QDET, with the exception of a single keyword: By adding the keyword l_qdet_verbose: true, we enable the write-out of both matrix elements of the RPA-screened and cRPA-screened matrix elements (typically only matrix elements of the cRPA-screened potential are written to JSON file).

You can download the input file with the following command:

[5]:

%%bash
wget -N -q https://west-code.org/doc/training/nv_diamond_63_spinpol/wfreq.in

We see that the wfreq.in file looks exactly like the input file used for the QDET calculation in Tutorial 5.

[6]:

%%bash
cat wfreq.in

wstat_control:
  wstat_calculation: S
  n_pdep_eigen: 512
  trev_pdep: 0.00001

wfreq_control:
  wfreq_calculation: XWGQH
  macropol_calculation: C
  l_qdet_verbose: true
  l_enable_off_diagonal: true
  n_pdep_eigen_to_use: 512
  qp_bands: [87, 122, 123, 126, 127, 128]
  n_refreq: 300
  ecut_refreq: 2.0

The wfreq.x calculation is performed with the following command

[ ]:

%%bash
mpirun -n 2 wfreq.x -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:

[7]:

%%bash
mkdir west.wfreq.save
wget -N -q https://west-code.org/doc/training/nv_diamond_63_spinpol/wfreq.json -O west.wfreq.save/wfreq.json

Step 2: Diagonalization of the eBSE Hamiltonian¶

[8]:

from westpy.qdet import eBSEResult

# construct object for effective Hamiltonian
ebse = eBSEResult(filename='west.wfreq.save/wfreq.json', spin_flip_=True)

# diagonalize Hamiltonian
solution = ebse.solve()


 _    _ _____ _____ _____
| |  | |  ___/  ___|_   _|
| |  | | |__ \ `--.  | |_ __  _   _
| |/\| |  __| `--. \ | | '_ \| | | |
\  /\  / |___/\__/ / | | |_) | |_| |
 \/  \/\____/\____/  \_/ .__/ \__, |
                       | |     __/ |
                       |_|    |___/

WEST version     :  5.5.0
Today            :  2024-02-02 14:37:36.123063
===============================================================
Solving eBSE Hamiltonian...
===============================================================

			diag[1RDM - 1RDM(GS)]
	E [eV]	char	87	122	123	126	127	128
0	0.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000
1	0.422	2.867	-0.000	-0.004	-0.003	-0.014	0.000	0.022
2	0.422	1.334	-0.000	-0.004	-0.003	-0.014	0.023	-0.001
3	1.054	1.000	-0.000	-0.000	-0.000	-0.000	-0.001	0.001
4	2.242	2.210	-0.000	-0.004	-0.007	-0.472	0.482	0.000
5	2.243	2.262	-0.000	-0.004	-0.007	-0.472	0.000	0.482
6	4.026	2.232	-0.001	-0.013	-0.473	-0.010	0.497	0.000
7	4.026	2.240	-0.001	-0.013	-0.473	-0.010	0.000	0.497
8	5.179	2.233	-0.001	-0.477	-0.017	-0.003	0.498	0.000
9	5.179	2.239	-0.001	-0.477	-0.017	-0.003	0.000	0.498
10	9.272	2.236	-0.498	-0.001	-0.000	-0.000	0.500	0.000
11	9.273	2.236	-0.498	-0.001	-0.000	-0.000	-0.000	0.500

-----------------------------------------------------

Calling the function solve() of the object eBSEResult writes the excitation energies (in eV), spin multiplicity and relative occupation (compared to the occupation of the ground state) to screen in a manner similar to the output of a QDET calculation. See also Tutorial 5.

Let us consider the first excitation (# 1) as an example: It has an energy of 0.427 eV above the ground state and has a multiplicity of 1.081. That means that it is a singlet, which should have a multiplicity of 1. The difference is called the spin contamination, describing that the spin state is note pure.

For a more detailed analysis, one can access the full output stored in the solution dictionary:

[9]:

print([key for key in solution.keys()])

['hamiltonian', 'evs_au', 'evs', 'evcs', 'rdm1s', 'mults', 'excitations']

Comparison between eBSE and QDET¶

We can easily compare the results from eBSE and QDET as we can perform both calculations from the same output:

[10]:

from westpy.qdet import QDETResult

# construct object for effective Hamiltonian
qdet = QDETResult(filename='west.wfreq.save/wfreq.json')

# diagonalize Hamiltonian
qdet_solution = qdet.solve(nelec=(5,5))

-----------------------------------------------------
===============================================================
Building effective Hamiltonian...
nspin: 2
occupations: [[1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 0. 0.]]
===============================================================

			diag[1RDM - 1RDM(GS)]
	E [eV]	char	87	122	123	126	127	128
0	0.000	3-	0.000	0.000	0.000	0.000	0.000	0.000
1	0.419	1-	-0.000	-0.004	-0.008	-0.035	0.000	0.048
2	0.419	1-	-0.000	-0.004	-0.008	-0.035	0.048	-0.000
3	1.205	1-	-0.001	-0.009	-0.010	-0.035	0.027	0.028
4	1.843	3-	-0.001	-0.006	-0.061	-0.431	0.497	0.002
5	1.843	3-	-0.001	-0.006	-0.061	-0.431	0.002	0.497
6	2.800	1-	-0.000	-0.022	-0.017	-0.427	0.465	0.000
7	2.801	1-	-0.000	-0.022	-0.017	-0.427	0.000	0.466
8	4.490	1-	-0.003	-0.030	-0.082	-0.846	0.480	0.480
9	4.957	3-	-0.007	-0.341	-0.118	-0.033	0.499	0.000

-----------------------------------------------------

As we can see, both approaches yield similar energies and qualitatively the same level diagram (triplet ground state followed by 3 singlet excitations), but the spin multiplicity is different. Generally, we observe more pronounced spin contamination in eBSE than in QDET.