Skip to content

Commit

Permalink
improve convergence of GW bandstructure calculation & improve manual
Browse files Browse the repository at this point in the history
  • Loading branch information
@JWilhelm
JWilhelm authored Jul 30, 2024
1 parent feaa26d commit 1b6f62e
Show file tree
Hide file tree
Showing 13 changed files with 662 additions and 574 deletions.
74 changes: 35 additions & 39 deletions docs/methods/properties/bandstructure_gw.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
energies of molecules. In this tutorial, we describe the band gap problem of DFT to motivate the
usage of *GW*, we briefly discuss the theoretical framework of *GW* and the *GW* implementation in
CP2K. We also provide input and output files of *GW* calculations performed with CP2K. The *GW*
theory of this tutorial is taken from [](#Golze2019) and many useful details on *GW* can be found in
this review.
theory of this tutorial is taken from \[[](#Golze2019)\] and many useful details on *GW* can be
found in this review.

## 1. The band gap problem of DFT

Expand All @@ -25,12 +25,12 @@ electronic band structure of a solid. This approximation comes with limitations:

- When using one of the common GGA exchange-correlation (xc) functionals, the band gap in the KS-DFT
band structure $\varepsilon_{n\mathbf{k}}^\text{DFT}$ is much too small compared to experimental
band gaps (Fig. 26 in [](#Golze2019)). Even with the exact xc functional, the band gap in the
band gaps (Fig. 26 in \[[](#Golze2019))\]. Even with the exact xc functional, the band gap in the
KS-DFT band structure $\varepsilon_{n\mathbf{k}}^\text{DFT}$ will be too small due to the
derivative discontinuity.

- The KS-DFT band structure $\varepsilon_{n\mathbf{k}}^\text{DFT}$ is insensitive to screening by
the environment. As an example, the KS eigenvalue gap
- The GGA band structure $\varepsilon_{n\mathbf{k}}^\text{DFT}$ is insensitive to screening by the
environment. As an example, the GGA eigenvalue gap
$\varepsilon_{n=\text{LUMO}}^\text{DFT}-\varepsilon_{n=\text{HOMO}}^\text{DFT}$ of a molecule is
almost identical for the molecule in gas phase and the molecule being on a surface. In experiment
however, the surface induces an image-charge effect which can reduce the HOMO-LUMO gap of the
Expand Down Expand Up @@ -80,9 +80,10 @@ $\varepsilon_{n\mathbf{k}}^\text{DFT}$ and we add the xc contribution from
$\braket{\psi_{n\mathbf{k}}|
\Sigma^{G_0W_0}(\varepsilon_{n\mathbf{k}}^{G_0W_0}) |\psi_{n\mathbf{k}}}$.

The DFT xc functional can influence the *G*<sub>0</sub>*W*<sub>0</sub> quasiparticle energies. For
molecules, it is recommended to start the *G*<sub>0</sub>*W*<sub>0</sub> calculation from the PBE0
functional and for solids, from the PBE functional.
The DFT xc start functional of *G*<sub>0</sub>*W*<sub>0</sub> can influence the
*G*<sub>0</sub>*W*<sub>0</sub> quasiparticle energies. For molecules, it is recommended to start the
*G*<sub>0</sub>*W*<sub>0</sub> calculation from the PBE0 functional and for solids, from the PBE
functional.

CP2K also allows to perform eigenvalue-selfconsistency in $G$ (ev*GW*<sub>0</sub>) and
eigenvalue-selfconsistency in $G$ and in $W$ (ev*GW*). For solids, it has been shown that band gaps
Expand All @@ -91,18 +92,19 @@ from ev*GW*$_0$@PBE can be in better agreement with experimental band gaps than

CP2K contains three different *GW* implementations:

- *GW* for molecules [](#Wilhelm2016) (Sec. 3)
- *GW* for molecules \[[](#Wilhelm2016)\] (Sec. 3)
- *GW* for computing the band structure of a solid with small unit cell with *k*-point sampling in
DFT (publication in preparation TODO: insert reference, Sec. 4)
- *GW* for computing the band structure of a large cell in a Γ-only approach [](#Graml2024) (Sec. 5)
- *GW* for computing the band structure of a large cell in a Γ-only approach \[[](#Graml2024)\]
(Sec. 5)

In the following, we will discuss details and usage of these *GW* implementations.

## 3. *GW* for molecules

For starting a *G*<sub>0</sub>*W*<sub>0</sub>, ev*GW*<sub>0</sub> or ev*GW* calculation, one needs
to set the [RUN_TYPE](#CP2K_INPUT.GLOBAL.RUN_TYPE) to `ENERGY` and one needs to put the following
section:
For starting a *G*<sub>0</sub>*W*<sub>0</sub>, ev*GW*<sub>0</sub> or ev*GW* calculation for a
molecule, one needs to set the [RUN_TYPE](#CP2K_INPUT.GLOBAL.RUN_TYPE) to `ENERGY` and one needs to
put the following section:

```
&XC
Expand All @@ -124,19 +126,21 @@ section:
In the upper *GW* section, the following keywords have been used:

- [QUADRATURE_POINTS](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.QUADRATURE_POINTS): Number
of imaginary-frequency points used for computing the self-energy. 100 points are usually enough
for converging quasiparticle energies within 10 meV.
of imaginary-frequency points used for computing the self-energy (Eq. (21) in
\[[](#Wilhelm2016))\]. 100 points are usually enough for converging quasiparticle energies within
10 meV.

- [SELF_CONSISTENCY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.SELF_CONSISTENCY):
Determines which *GW* self-consistency (*G*<sub>0</sub>*W*<sub>0</sub>, ev*GW*<sub>0</sub> or
ev*GW*) is used to calculate the *GW* quasiparticle energies.

The numerical precision of the *GW* implementation is 10 meV compared to reference calculations, for
example the *GW*100 test set [](#vanSetten2015). To help new users get familiar with the *GW*
example the *GW*100 test set \[[](#vanSetten2015)\]. To help new users get familiar with the *GW*
implementation in CP2K, we recommend starting by reproducing the HOMO and LUMO
*G*<sub>0</sub>*W*<sub>0</sub>@PBE quasiparticle energies of the H<sub>2</sub>O molecule from the
GW100 test set. The reference values are $\varepsilon_\text{HOMO}^{G_0W_0\text{@PBE}}$ = -11.97 eV
and $\varepsilon_\text{LUMO}^{G_0W_0\text{@PBE}}$ = 2.37 eV; input and output files are available
and $\varepsilon_\text{LUMO}^{G_0W_0\text{@PBE}}$ = 2.37 eV; CP2K input and output files for the
*G*<sub>0</sub>*W*<sub>0</sub>@PBE calculation of the H<sub>2</sub>O molecule are available
[here](https://github.com/cp2k/cp2k-examples/tree/master/bandstructure_gw/1_H2O_GW100).

The following settings from DFT will also have an influence on *GW* quasiparticle energies:
Expand All @@ -145,19 +149,11 @@ The following settings from DFT will also have an influence on *GW* quasiparticl
*G*<sub>0</sub>*W*<sub>0</sub>, ev*GW*<sub>0</sub> or ev*GW* calculation. For molecules, we
recommend either ev*GW*<sub>0</sub>@PBE or *G*<sub>0</sub>*W*<sub>0</sub>@PBE0. For further
guidance on selecting an appropriate DFT starting functional and self-consistency scheme for your
system, you may consult [](#Golze2019) or the following video:

```{youtube} 1vUuethWhbs
---
url_parameters: ?start=5563
align: center
privacy_mode:
---
```
system, you may consult \[[](#Golze2019)\].

- [BASIS_SET](#CP2K_INPUT.FORCE_EVAL.SUBSYS.KIND.BASIS_SET): The basis set is of Gaussian type and
can have strong influence on the quasiparticle energies. For computing quasiparticle energies, a
basis set extrapolation is necessary, see Fig. 2a in [](#Wilhelm2016) and we recommend
basis set extrapolation is necessary, see Fig. 2a in \[[](#Wilhelm2016)\] and we recommend
all-electron GAPW calculations with correlation-consistent basis sets from the
<a href="https://www.basissetexchange.org/" target="_blank">EMSL database</a>. For computing the
HOMO-LUMO gap from *GW*, we recommend augmented basis sets, for example
Expand Down Expand Up @@ -247,11 +243,11 @@ details in forthcoming publication (TODO Link):
radius of truncated Coulomb metric in Å. A larger cutoff leads to faster the RI basis set
convergence, but also computational cost increases. A cutoff of 7 Å is an accurate choice.
- [&SOC](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.BANDSTRUCTURE.SOC): Activates spin-orbit coupling (SOC
from Hartwigsen-Goedecker-Hutter pseudopotentials \[[Hartwigsen1998](#Hartwigsen1998)\]. SOC also
from Hartwigsen-Goedecker-Hutter pseudopotentials \[[Hartwigsen1998](#Hartwigsen1998)\]). SOC also
needs `POTENTIAL_FILE_NAME GTH_SOC_POTENTIALS`.
- [&BANDSTRUCTURE_PATH](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.BANDSTRUCTURE.BANDSTRUCTURE_PATH): Specify
the *k*-path in the Brillouin zone for computing the band structure. Relative *k*-coordinates are
needed which you can retrieve for your crystal structure from [](#Setyawan2010).
needed which you can retrieve for your crystal structure from \[[](#Setyawan2010)\].

We recommend TZVP-MOLOPT basis sets together with GTH/HGH pseudopotentials, see basis set
convergence study in (TODOref). At present, 2d-periodic boundary conditions are supported, 1d- and
Expand All @@ -271,21 +267,21 @@ hours on 20 large-memory nodes).

## 5. *GW* for large cells in Γ-only approach

For a large unit cell, a Γ-only *GW* algorithm is available in CP2K [](#Graml2024). The requirement
on the cell is that elements of the density matrix decay by several orders of magnitude when the two
basis functions of the matrix element have a distance similar to the cell size. As rule of thumb,
for a 2d material, a 9x9 unit cell is large enough for the Γ-only algorithm, see Fig. 1 in
[](#Graml2024).
For a large unit cell, a Γ-only *GW* algorithm is available in CP2K \[[](#Graml2024)\]. The
requirement on the cell is that elements of the density matrix decay by several orders of magnitude
when the two basis functions of the matrix element have a distance similar to the cell size. As rule
of thumb, for a 2d material, a 9x9 unit cell is large enough for the Γ-only algorithm, see Fig. 1 in
\[[](#Graml2024)\].

The input file for a Γ-only *GW* calculation is identical as for *GW* for small cells with *k*-point
sampling except that the `&KPOINTS` section in DFT needs to be avoided. An exemplary input and
sampling except that the `&KPOINTS` section in DFT needs to be removed. An exemplary input and
output is available
\[<a href="https://github.com/cp2k/cp2k-examples/tree/master/bandstructure_gw/5_9x9_supercell_2d_MoS2" target="_blank">here</a>\]
\[<a href="https://github.com/cp2k/cp2k-examples/tree/master/bandstructure_gw/5_9x9_supercell_2d_MoS2" target="_blank">here</a>\].
Running the input file requires access to a large computer (the calculation took 2.5 hours on 32
nodes on Noctua2 cluster in Paderborn). The computational parameters from this input file reach
numerical convergence of the band gap within ~ 50 meV (TZVP basis set, 10 time and frequency
points). Detailed convergence tests are available in the SI, Table S1 of [](#Graml2024) We recommend
the numerical parameters from the input file for large-scale GW calculations. The code prints
restart files with ending .matrix that can be used to restart a crashed calculation.
points). Detailed convergence tests are available in the SI, Table S1 of \[[](#Graml2024)\] We
recommend the numerical parameters from the input file for large-scale GW calculations. The code
prints restart files with ending .matrix that can be used to restart a crashed calculation.

In case anything does not work, please feel free to contact jan.wilhelm (at) ur.de.
108 changes: 20 additions & 88 deletions docs/methods/properties/optical/bethe-salpeter.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,12 @@ The following ingredients are necessary for computing $\Omega^{(n)}$:
- Occupied Kohn-Sham (KS) orbitals $\varphi_i(\mathbf{r})$ and empty KS orbitals
$\varphi_a(\mathbf{r})$ from a DFT calculation, where $i=1,\ldots,N_\mathrm{occ}$ and
$a=N_\mathrm{occ}+1,\ldots,N_\mathrm{occ}+N_\mathrm{empty}$,
- $GW$ eigenvalues $\varepsilon_i^{GW}$ and $\varepsilon_a^{GW}$ of corresponding KS orbitals.
- *GW* eigenvalues $\varepsilon_i^{GW}$ and $\varepsilon_a^{GW}$ of corresponding KS orbitals.

It is possible to use $G_0W_0$, ev$GW_0$ or ev$GW$ eigenvalues, see details in [GW] and in Ref.
\[[](#Golze2019)\], i.e. we perform BSE@$G_0W_0$/ev$GW_0$/ev$GW$@DFT. Thus, also input parameters
for a DFT and $GW$ calculation influence the BSE calcualation (see full input in Sec.
It is possible to use *G*<sub>0</sub>*W*<sub>0</sub>, ev*GW*<sub>0</sub> or ev*GW* eigenvalues, see
details in [GW] and in Ref. \[[](#Golze2019)\], i.e. we perform
BSE@*G*<sub>0</sub>*W*<sub>0</sub>/ev*GW*<sub>0</sub>/ev*GW*@DFT. Thus, also input parameters for a
DFT and *GW* calculation influence the BSE calcualation (see full input in Sec.
[3.1](#header-input-file)). We obtain optical properties from BSE solving the following generalized
eigenvalue problem that involves the block matrix $ABBA$:

Expand Down Expand Up @@ -86,8 +87,9 @@ For starting a BSE calculation one needs to set the [RUN_TYPE](#CP2K_INPUT.GLOBA
In the upper GW/BSE section, the following keywords have been used:

- [SELF_CONSISTENCY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.SELF_CONSISTENCY):
Determines which GW self-consistency ($G_0W_0$, ev$GW_0$ or ev$G_0W_0$) is used to calculate the
single-particle GW energies $\varepsilon_p^{GW}$ needed in the BSE calculation.
Determines which *GW* self-consistency (*G*<sub>0</sub>*W*<sub>0</sub>, ev*GW*<sub>0</sub> or
ev*GW*) is used to calculate the single-particle *GW* energies $\varepsilon_p^{GW}$ needed in the
BSE calculation.

- [TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA): Three options available:

Expand Down Expand Up @@ -123,8 +125,8 @@ The following settings from DFT will also have an influence on the BSE excitatio
- [XC_FUNCTIONAL](#CP2K_INPUT.FORCE_EVAL.DFT.XC.XC_FUNCTIONAL): Choose between one of the available
xc-functionals. The starting point can have a profound influence on the excitation energies
\[[Bruneval2015](#Bruneval2015), [](#Gui2018), [](#Jacquemin2017)\]. We either recommend the PBE
functional as DFT starting point when using BSE@ev$GW_0$@PBE or the PBE0 functional, when using
BSE@$G_0W_0$@PBE0 (see also
functional as DFT starting point when using BSE@ev*GW*<sub>0</sub>@PBE or the PBE0 functional,
when using BSE@*G*<sub>0</sub>*W*<sub>0</sub>@PBE0 (see also
[SELF_CONSISTENCY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.SELF_CONSISTENCY)).

- [BASIS_SET](#CP2K_INPUT.FORCE_EVAL.SUBSYS.KIND.BASIS_SET): Specify the basis set, which affects
Expand All @@ -143,7 +145,7 @@ We have benchmarked the numerical precision of our BSE implementation and we fou
agreement within only 10 meV compared to the BSE implementation in FHI aims \[[](#Liu2020)\].

The current BSE implementation in CP2K works for molecules. The inclusion of periodic boundary
conditions in a $\Gamma$-only approach and with full $k$-point sampling is work in progress.
conditions in a Γ-only approach and with full *k*-point sampling is work in progress.

(header-example)=

Expand All @@ -153,7 +155,7 @@ conditions in a $\Gamma$-only approach and with full $k$-point sampling is work

### 3.1 Input file

In this section, we provide a minimal example of a BSE calculation on $\mathrm{H}_2$. For the
In this section, we provide a minimal example of a BSE calculation on H<sub>2</sub>. For the
calculation you need the input file BSE_H2.inp and the aug-cc-DZVP basis
([Download](https://github.com/cp2k/cp2k-examples/tree/master/bethe-salpeter/H2)).

Expand All @@ -163,83 +165,13 @@ Please copy both files into your working directory and run CP2K by
mpirun -n 1 cp2k.psmp BSE_H2.inp
```

which requires 12 GB RAM and takes roughly 2 minutes on 1 core. You can find the output file also in
the [Download](https://github.com/cp2k/cp2k-examples/tree/master/bethe-salpeter/H2).

```none
&GLOBAL
PROJECT H2
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS-aug ! Custom Basis set file (aug-cc-pVDZ and aug-cc-pVDZ-RIFIT
POTENTIAL_FILE_NAME POTENTIAL ! from the Basis Set Exchange library - www.basissetexchange.org/)
&QS
METHOD GAPW ! All electron calculation
EPS_DEFAULT 1.0E-16
EPS_PGF_ORB 1.0E-16
&END QS
&POISSON
PERIODIC NONE
PSOLVER MT
&END
&SCF
SCF_GUESS RESTART
EPS_SCF 1e-7
&END SCF
&XC
&XC_FUNCTIONAL PBE0 ! Choice of functional has a profound influence on the results
&END XC_FUNCTIONAL
&WF_CORRELATION
&RI_RPA ! In the RI_RPA and the GW section, additional numerical parameters, e.g.
&GW
SELF_CONSISTENCY G0W0 ! can be changed to EV_GW0 or EV_GW
&BSE
TDA TDA+ABBA ! Diagonalizing ABBA and A
SPIN_CONFIG SINGLET ! or TRIPLET
NUM_PRINT_EXC 20 ! Number of printed excitations
ENERGY_CUTOFF_OCC -1 ! Set to positive numbers (eV) to
ENERGY_CUTOFF_EMPTY -1 ! truncate matrices A_ia,jb and B_ia,jb
&END BSE
&END GW
&END RI_RPA
&END WF_CORRELATION
&END XC
&END DFT
&SUBSYS
&CELL
ABC 20 20 20
PERIODIC NONE
&END CELL
&COORD
H 0.0000 0.0000 0.0000 ! H2 molecule geometry from GW100 Paper
H 0.0000 0.0000 0.74144
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY
&KIND H
BASIS_SET ORB aug-cc-pVDZ ! For production runs, the basis set should be checked for convergence.
BASIS_SET RI_AUX aug-cc-pVDZ-RIFIT ! In general, pVDZ should be a solid choice.
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET ORB aug-cc-pVDZ
BASIS_SET RI_AUX aug-cc-pVDZ-RIFIT
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
```

The basis sets `aug-cc-pVDZ` and `aug-cc-pVDZ-RIFIT` in `BASIS-aug` can be obtained from the Basis
Set Exchange Library:
which requires 12 GB RAM and takes roughly 2 minutes on 1 core. You can find the input and output
file [here](https://github.com/cp2k/cp2k-examples/tree/master/bethe-salpeter/H2). We use the basis
sets `aug-cc-pVDZ` and `aug-cc-pVDZ-RIFIT` from the file `BASIS-aug`. These basis sets can be
obtained from the Basis Set Exchange Library:
<a href="https://www.basissetexchange.org/basis/aug-cc-pvdz/format/cp2k/?version=1&elements=1" target="_blank">aug-cc-pVDZ</a>,
<a href="https://www.basissetexchange.org/basis/aug-cc-pvdz-rifit/format/cp2k/?version=1&elements=1" target="_blank">aug-cc-pVDZ-RIFIT</a>.
The geometry for $\mathrm{H}_2$ was taken from [](#vanSetten2015).
The geometry for H<sub>2</sub> was taken from \[[](#vanSetten2015)\].

### 3.2 Output

Expand Down Expand Up @@ -277,14 +209,14 @@ printed:
BSE| 4 1 => 5 -full- 0.7076
```

In the case of the $\mathrm{H}_2$, the lowest excitation *n* = 1 is mainly built up by a transition
from the HOMO (i=1) to the LUMO (a=2), what is apparent from
In the case of H<sub>2</sub>, the lowest excitation *n* = 1 is mainly built up by a transition from
the HOMO (i=1) to the LUMO (a=2), what is apparent from
$X_{i=\text{HOMO},a=\text{LUMO}}^{(n=1)}=0.663$, and also contains a considerable contribution from
the 1=>4 (HOMO=>LUMO+2) transition ($X_{i=\text{HOMO},a=\text{LUMO+2}}^{(n=1)}=0.2549$ ).

### 3.3 Large scale calculations

Going to larger systems is a challenge for a $GW$+BSE-calculation, since the memory consumption
Going to larger systems is a challenge for a *GW*+BSE-calculation, since the memory consumption
increases with $N_\mathrm{occ}^2 N_\mathrm{empty}^2$. The used $N_\mathrm{occ}$, $N_\mathrm{empty}$
and the required memory of a calculation are printed in the output file to estimate the memory
consumption. We recommend to use several nodes to provide the required memory. In the following, we
Expand Down
Loading

0 comments on commit 1b6f62e

Please sign in to comment.