comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

BSE: new features and update manual (#3781)

Browse source
This commit is contained in:
Max Graml 2024-11-20 16:34:36 +01:00 • committed by GitHub
parent 52af7ed493
commit 4dc10b3491
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
23 changed files with 3600 additions and 876 deletions
Download patch file Download diff file Expand all files Collapse all files

513
docs/methods/properties/optical/bethe-salpeter.md
View file

@ -10,60 +10,277 @@ input file of a BSE calculation and the corresponding output. For reviews on BSE
## 1. Theory and implementation of BSE
A central goal of a BSE calculation is to compute electronic excitation energies
$\Omega^{(n)}, n=1,2,\ldots$.
### 1.1 Electronic excitation energies and the excitation wave function
The following ingredients are necessary for computing $\Omega^{(n)}$:
A central goal of a BSE calculation is to compute excited state properties such as the electronic
excitation energies $\Omega^{(n)} (n=1,2,\ldots)$.
The following ingredients are necessary for such a computation:
- 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.
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
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 calculation (see full input in Sec. [3.1](#header-input-file)). We
obtain the properties of the excited states from the BSE by solving the following generalized
eigenvalue problem that involves the block matrix $ABBA$:
$$\left( \begin{array}{cc}A & B\\B & A\end{array} \right)\left( \begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array} \right) = \Omega^{(n)}\left(\begin{array}{cc}1&0\\0&-1\end{array}\right)\left(\begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array}\right) \quad .$$
$$
\left( \begin{array}{cc}A & B\\B & A\end{array} \right)\left( \begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array} \right) = \Omega^{(n)}\left(\begin{array}{cc}1&0\\0&-1\end{array}\right)\left(\begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array}\right) \quad .
$$
We abbreviate $A$ and $B$ as matrices with index $A_{ia,jb}$, i.e. they have
$N_\mathrm{occ}N_\mathrm{empty}$ rows and $N_\mathrm{occ}N_\mathrm{empty}$ columns. The entries of
$A$ and
$$ \begin{align}
$$
\begin{align}
A_{ia,jb} &= (\varepsilon_a^{GW}-\varepsilon_i^{GW})\delta_{ij}\delta_{ab} + \alpha^\mathrm{S/T}
v_{ia,jb} - W_{ij,ab}(\omega=0) \quad ,\\
B_{ia,jb} &= \alpha^\mathrm{(S/T)} v_{ia,bj} - W_{ib,aj}(\omega=0) \quad .
\end{align}$$
\end{align}
$$
where $\delta_{ij}$ is the Kronecker delta. The user sets $\alpha^S=2$ for computing singlet
excitations and $\alpha^T=0$ for computing triplet excitations. $v_{pq,rs}$ is the bare Coulomb
interaction and $W_{pq,rs}(\omega=0)$ the statically ($\omega=0$) screened Coulomb interaction ,
where $p,q,r,s \in [ 1, N_\mathrm{occ}+N_\mathrm{empty}]$ are KS orbital indices.
$(\mathbf{X}^{(n)},\mathbf{Y}^{(n)})$ with elements $X_{ia}^{(n)}$ and $Y_{ia}^{(n)}$ are the
eigenvectors of the excitation $n$ which relate to the wavefunction of the electronic excitation
\[[](#Blase2020)\],
interaction and $W_{pq,rs}(\omega=0)$ the statically ($\omega=0$) screened Coulomb interaction,
where $p,q,r,s \in [ 1, N_\mathrm{occ}+N_\mathrm{empty}]$ are KS orbital indices. Note here, that
the screened Coulomb interaction is always computed from DFT quantitites, i.e. $W_{0}(\omega=0)$
enters the BSE. $(\mathbf{X}^{(n)},\mathbf{Y}^{(n)})$ with elements $X_{ia}^{(n)}$ and
$Y_{ia}^{(n)}$ are the eigenvectors of the excitation $n$, which relate to the wave function of the
electronic excitation \[[](#Blase2020)\],
$$ \begin{align}
$$
\begin{align}
\Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) = \sum_{ia} X_{ia}^{(n)} \varphi_i(\mathbf{r}_h) \varphi_a(\mathbf{r}_e) + Y_{ia}^{(n)} \varphi_i(\mathbf{r}_e) \varphi_a(\mathbf{r}_h) \quad ,
\end{align}$$
\end{align}
$$
i.e. $X_{ia}^{(n)}$ and $Y_{ia}^{(n)}$ describe the transition amplitude between occupied orbital
$\varphi_i$ and empty orbital $\varphi_a$ of the $n$-th excitation.
The Tamm-Dancoff approximation (TDA) is also implemented, which constrains $B=0$. In case $A$ is
positive definite, and excitation energies $\Omega^{(n)}>0$, we have $\mathbf{Y}=0$ and $\mathbf{X}$
can be computed from the Hermitian eigenvalue problem
The Tamm-Dancoff approximation (TDA) is also implemented, which constrains $B=0$ and
$\mathbf{Y}^{(n)}=0$. In TDA, $\Omega^{(n)}_\mathrm{TDA}$ and $\mathbf{X}^{(n)}_\mathrm{TDA}$ can be
computed from the Hermitian eigenvalue problem
$$ A \mathbf{X}^{(n)}_\mathbf{TDA} = \Omega^{(n)}_\mathbf{TDA} \mathbf{X}^{(n)}_\mathbf{TDA} \quad .$$
$$
A \mathbf{X}^{(n)}_\mathrm{TDA} = \Omega^{(n)}_\mathrm{TDA} \mathbf{X}^{(n)}_\mathrm{TDA} \quad .
$$
Diagonalizing $A$ in TDA, or the full block-matrix $ABBA$, takes in the order of
$(N_\mathrm{occ} N_\mathrm{empty})^3$ floating point operations. This translates to a computational
scaling of $O(N^6)$ in the system size $N$.
### 1.2 Optical properties
The BSE further allows the investigation of optical properties. For example, the optical absorption
spectrum can be computed as the imaginary part of the dynamical dipole polarizability tensor
$\alpha_{\mu,\mu'}(\omega) $ with $(\mu,\mu'\in\{x,y,z\})$:
$$
\begin{align}
\alpha_{\mu,\mu'}(\omega)
= \sum_n \frac{2 \Omega^{(n)} d^{(n)}_{\mu} d^{(n)}_{\mu'}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2}
\quad ,
\end{align}
$$
where we have introduced an artificial broadening $\eta$. The transition moments $d^{(n)}_{\mu}$ are
computed in the length gauge $(\mu\in\{x,y,z\})$ as
$$
\begin{align}
d^{(n)}_{\mu} = \sqrt{2} \sum_{i,a} \langle \varphi_i|\hat{\mu}| \varphi_a \rangle (X_{ia}^{(n)} + Y_{ia}^{(n)})
\quad .
\end{align}
$$
When the molecules are not aligned, e.g. for gas phase and liquids, the spatial average is
sufficient, i.e. the optical absorption spectrum can be computed as
$$
\begin{align}
\mathrm{Im}\left[\bar{\alpha}(\omega)\right] = \frac{1}{3} \sum_{\mu\in\{x,y,z\}} \mathrm{Im}\left[\alpha_{\mu,\mu}(\omega)\right]
\quad .
\end{align}
$$
We can rewrite the last equation as
$$
\begin{align}
\mathrm{Im}\left[\bar{\alpha}(\omega)\right]
= \mathrm{Im}\left[
\sum_n \frac{f^{(n)}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2}
\right]
\quad .
\end{align}
$$
where we introduced the oscillator strengths $f^{(n)}$, which are defined by
$$
\begin{align}
f^{(n)} = \frac{2}{3} \Omega^{(n)} \sum_{\mu\in\{x,y,z\}} | d^{(n)}_{\mu} |^2
\quad .
\end{align}
$$
### 1.3 Visualizing excited states using Natural Transition Orbitals (NTOs)
In order to represent the exciton wave function independent of a specific choice of the molecular
orbitals $\varphi_p(\mathbf{r})$, we can rewrite it as
$$
\begin{align}
\Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) =
\sum_I {\lambda_I^{(n)}} \phi_I^{(n)}(\mathbf{r}_e) \chi_I^{(n)}(\mathbf{r}_h)
\quad .
\end{align}
$$
in terms of the natural transitions orbitals (NTOs) $\phi_I^{(n)}(\mathbf{r}_e) $ and
$\chi_I^{(n)}(\mathbf{r}_h)$. Here, we introduce the idea of electron and holes: In the excitation
process, an electron leaves a hole behind in the (formerly) occupied orbitals
$\chi_I^{(n)}(\mathbf{r}_h)$ and is afterwards located in the empty orbitals
$\phi_I^{(n)}(\mathbf{r}_e) $. We call $\phi_I^{(n)}(\mathbf{r}_e) $ and
$\chi_I^{(n)}(\mathbf{r}_h)$ a NTO pair of electron and hole with a weight $\lambda_I^{(n)}$, where
$I=1,\cdots,N_b$ with $N_b=N_\mathrm{occ}+N_\mathrm{empty}$.
The weights $\lambda_I^{(n)}$ are sorted in decreasing order, i.e.:
$$\lambda_1^{(n)} \ge \lambda_2^{(n)} \ge \cdots \ge 0$$
For many cases, we observe $\lambda_1^{(n)} \gg \lambda_2^{(n)}$, which indicates a dominant
single-particle excitation character for this $n$ and allows a visual analysis of a small number of
NTO pairs.
Assuming $\lambda_1^{(n)} = 1$ and $\lambda_{I\neq 1}=0$, the excitation wave function simply is
given as a product
$$
\begin{align}
\Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) =
\phi_1^{(n)}(\mathbf{r}_e) \chi_1^{(n)}(\mathbf{r}_h)
\quad .
\end{align}
$$
In this case, the electron is excited from the occupied NTO $\chi_1^{(n)}(\mathbf{r}_h)$ to the
empty state $\phi_1^{(n)}(\mathbf{r}_e)$. This process leaves a hole at $\mathbf{r}_h$ and creates
an excited electron at $\mathbf{r}_e$.
Computationally, the NTO pairs $\phi_I^{(n)}(\mathbf{r}_e) $ and $\chi_I^{(n)}(\mathbf{r}_h)$ as
well as $\lambda_I^{(n)}$ are obtained by a singular value decomposition of the transition density
matrix
$$
T = \begin{pmatrix}
0 & {X}^{(n)}\\
\left({Y}^{(n)} \right)^T & 0
\end{pmatrix} \quad ,
$$
i.e.:
$$
\begin{align}
{T}^{(n)} &=
{U}^{(n)}
{\Lambda^{(n)}}
\left({V}^{(n)}\right)^T
\\
\phi_I^{(n)}(\mathbf{r}_e) &= \sum_{p=1}^{N_b} \varphi_p(\mathbf{r}_e) V_{p,I}^{(n)} \quad ,
\\
\chi_I^{(n)}(\mathbf{r}_h) &= \sum_{q=1}^{N_b} \varphi_q(\mathbf{r}_h) U_{q,I}^{(n)} \quad .
\end{align}
$$
### 1.4 Measures for the size of an excited state
In this subsection, we introduce several measures, which allow us to analyse the wave function of an
excited state $\Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h)$ and its spatial properties,
following Ref. \[[](#Mewes2018)\].
To that end, we define the exciton expectation value with respect to a generic operator $\hat{O}$ as
$$
\begin{align}
\langle \hat{O} \rangle _\text{exc}^{(n)}
=
\frac{
\langle \Psi_\text{excitation}^{(n)} | \hat{O} | \Psi_\text{excitation}^{(n)}\rangle
}{
\langle \Psi_\text{excitation}^{(n)} | \Psi_\text{excitation}^{(n)}\rangle
}
\quad ,
\end{align}
$$
where we drop the excitation index $n$ from now on for better readability.
For each excitation level $n$, we have then several quantities, which allow us to quantify the
spatial extend of the electron, the hole and their combined two-particle character as combined
electron-hole pair, i.e. as an exciton. By that, we can often determine the type of the excited
state, i.e. distinguish between, e.g., valence, Rydberg or charge-transfer states
\[[](#Mewes2018)\].
First, we define the distance between electron and hole as
$$
\begin{align}
d_{h \rightarrow e} = | \langle \mathbf{r}_h - \mathbf{r}_e \rangle_\mathrm{exc} | \quad ,
\end{align}
$$
which can be used to distinguish different classes of excitations: For example in a charge-transfer
state, electron and hole sit on different parts of the molecule and therefore have a non-vanishing
electron-hole distance $d_{h \rightarrow e}$.
Further, we can measure the size of electron and hole, respectively, as
$$
\begin{align}
\sigma_{e/h} = \sqrt{
\langle \mathbf{r}_{e/h}^2 \rangle_\mathrm{exc}
- \langle \mathbf{r}_{e/h} \rangle_\mathrm{exc} ^2
} \quad ,
\end{align}
$$
which allow us to distinguish between Rydberg states, where $\sigma_h \ll \sigma_e$, and valence
states, where $\sigma_h \approx \sigma_e$.
Closely related to these quantities, we can also define the exciton size
$$
\begin{align}
d_\mathrm{exc} = \sqrt{ \langle |\mathbf{r}_h - \mathbf{r}_e|^2 \rangle_\mathrm{exc} } \quad .
\end{align}
$$
which quantifies the spatial extent of the combined electron-hole pair. As one would expect, the
exciton size $d_\mathrm{exc}$ increases when $d_{h \rightarrow e}$, $\sigma_{e}$ or $\sigma_{h}$
increase.
Finally, we quantify the correlation of electron and hole by the electron-hole correlation
coefficient
$$
\begin{align}
R_{eh} = \frac{1}{\sigma_e \sigma_h} \left( \langle \mathbf{r}_h \cdot \mathbf{r}_e \rangle_\mathrm{exc}
- \langle \mathbf{r}_h \rangle_\mathrm{exc} \cdot \langle \mathbf{r}_e \rangle_\mathrm{exc} \right)
\quad ,
\end{align}
$$
which allows us to distinguish between correlated ($R_{eh}>0$) motion, i.e. bound excitons, and
anticorrelated ($R_{eh}<0$) motion, where electron and hole try to avoid each other.
(header-input)=
## 2. BSE input
@ -73,13 +290,22 @@ For starting a BSE calculation one needs to set the [RUN_TYPE](#CP2K_INPUT.GLOBA
```
&GW
SELF_CONSISTENCY G0W0 ! can be changed to EV_GW0 or EV_GW
SELF_CONSISTENCY evGW0 ! We strongly recommend to use evGW0 (and PBE) for BSE runs
&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
TDA TDA+ABBA ! Diagonalizing ABBA and A
SPIN_CONFIG SINGLET ! or TRIPLET
NUM_PRINT_EXC 15 ! 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
NUM_PRINT_EXC_DESCR -1 ! Number of printed exciton descriptors, -1 to default to NUM_PRINT_EXC
&BSE_SPECTRUM ! Activates computation and output of optical absorption spectrum
ETA_LIST 0.01 0.02 ! Multiple broadenings can be specified within one run
&END BSE_SPECTRUM
&NTO_ANALYSIS
STATE_LIST 1 2 ! List of states for which NTOs are computed
CUBE_FILES T ! Write cube files for NTOs
STRIDE 6 6 6 ! Coarse grid for smaller example cube files
&END NTO_ANALYSIS
&END BSE
&END GW
```
@ -89,7 +315,9 @@ 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*<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.
BSE calculation. We recommend to always use ev*GW*<sub>0</sub> for BSE runs. The screened coulomb
interaction is always computed from the DFT level, i.e. $W_0(\omega=0)$ is used for the
$ABBA$-matrix, including BSE@ev*GW*@DFT.
- [TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA): Three options available:
@ -120,29 +348,61 @@ In the upper GW/BSE section, the following keywords have been used:
empty states in the interval
$\varepsilon_a\in[\varepsilon_{a=\text{LUMO}}^{GW},\varepsilon_{a=\text{LUMO}}^{GW}+E_\text{cut}^\text{empty}]$.
The following settings from DFT will also have an influence on the BSE excitation energies:
- [NUM_PRINT_EXC_DESCR](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.NUM_PRINT_EXC_DESCR):
Number of excitations, for which the exciton descriptors are printed.
- [BSE_SPECTRUM](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM): Activates
computation and printing of the optical absorption spectrum. For each chosen option in
[TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA), a file is printed with
columns frequency $\omega$, $\mathrm{Im}\left[\bar{\alpha}(\omega)\right]$ and the imaginary part
of the elements of the dynamical dipole polarizability tensor
$\mathrm{Im}\left[{\alpha_{\mu,\mu'}}(\omega)\right]$. The frequency range, step size and the
broadening $\eta$ can be specified by the user (cf. keywords in
[BSE_SPECTRUM](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM)). Further,
multiple broadenings $\eta$ can be given for one cp2k run
([ETA_LIST](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM.ETA_LIST)),
resulting in one file per specified $\eta$ and
[TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA) (see Sec.
[3.2](#header-output)).
- [NTO_ANALYSIS](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.NTO_ANALYSIS): Activates
computation and printing of the natural transition orbitals. Within this section, the user can
specify the NTO pairs, which should be printed, e.g. by the excitation indices of the desired
levels $n$, their minimum oscillator strength $f^{(n)}$ or by the minimum weight of their NTO
coefficient $\lambda_I$. For each NTO pair $\phi_I^{(n)}(\mathbf{r}_e) $ and
$\chi_I^{(n)}(\mathbf{r}_h)$ of a certain excitation $n$, which fulfills these requirements, and
the chosen options in [TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA), the
NTO coefficients $\lambda_I$ are computed and the corresponding NTOs are printed to separate files
if
[CUBE_FILES](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.NTO_ANALYSIS.CUBE_FILES)
is activated.
The following settings from DFT will also have an influence on the BSE excitation energies and
optical properties:
- [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*<sub>0</sub>@PBE or the PBE0 functional,
when using BSE@*G*<sub>0</sub>*W*<sub>0</sub>@PBE0 (see also
\[[](#Knysh2024)\]. Motivated by the discussion in \[[](#Schambeck2024)\], we strongly recommend
to use BSE@ev*GW*<sub>0</sub>@PBE, i.e. the PBE functional as DFT starting point (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
$N_\mathrm{empty}$ and thus the size of the matrices $A_{ia,jb}$ and $B_{ia,jb}$. The
<a href="https://www.basissetexchange.org/basis/aug-cc-pvdz/format/cp2k/?version=1&elements=1" target="_blank">`aug-cc-pVDZ`</a>
basis set should be sufficient for most calculations, but needs to be checked regarding
convergence, e.g. using
<a href="https://www.basissetexchange.org/basis/aug-cc-pvtz/format/cp2k/?version=1&elements=1" target="_blank">`aug-cc-pVTZ`</a>.
<a href="https://www.basissetexchange.org/basis/aug-cc-pvtz/format/cp2k/?version=1&elements=1" target="_blank">`aug-cc-pVTZ`</a>
.
The memory consumption of the BSE algorithm is large, it is approximately
$100 \cdot N_\mathrm{occ}^2 N_\mathrm{empty}^2$ Bytes. You can see $N_\mathrm{occ}$,
$N_\mathrm{empty}$ and the estimated memory consumption from the BSE output. The BSE implementation
is well parallelized, i.e. you can use several nodes that can provide the memory.
We have benchmarked the numerical precision of our BSE implementation and we found excellent
agreement within only 10 meV compared to the BSE implementation in FHI aims \[[](#Liu2020)\].
We have benchmarked the numerical precision of our BSE implementation and compared its results to
the BSE implementation in FHI aims \[[](#Liu2020)\]. For our recommended settings, i.e.
BSE@ev*GW*<sub>0</sub>@PBE with the aug-cc-pVDZ basis set, we have found excellent agreement with
less than 5 meV mean absolute deviation over the first 10 excitation levels and the 28 molecules in
*Thiel's set*.
The current BSE implementation in CP2K works for molecules. The inclusion of periodic boundary
conditions in a Γ-only approach and with full *k*-point sampling is work in progress.
@ -156,7 +416,7 @@ conditions in a Γ-only approach and with full *k*-point sampling is work in pro
### 3.1 Input file
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
calculation you need the input file BSE_H2.inp and the aug-cc-pVDZ basis
([Download](https://github.com/cp2k/cp2k-examples/tree/master/bethe-salpeter/H2)).
Please copy both files into your working directory and run CP2K by
@ -165,82 +425,165 @@ 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 input and output
which requires 5 GB RAM and takes roughly 45 seconds 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>.
<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 H<sub>2</sub> was taken from \[[](#vanSetten2015)\].
(header-output)=
### 3.2 Output
The BSE calculation outputs the excitation energies $\Omega^{(n)}$ for *n* = 1, ...,
$N_\text{exc}^\text{print}$:
```none
BSE| Excitation energies from solving the BSE without the TDA:
BSE|
BSE| Excitation n Multiplet TDA/full BSE Excitation energy Ω (eV)
BSE| 1 Singlet State -full- 11.8621
BSE| 2 Singlet State -full- 12.8264
BSE| 3 Singlet State -full- 15.6415
BSE| 4 Singlet State -full- 15.6415
BSE| Excitation energies from solving the BSE without the TDA:
BSE|
BSE| Excitation n Spin Config TDA/ABBA Excitation energy Ω^n (eV)
BSE| 1 Singlet -ABBA- 12.0361
BSE| 2 Singlet -ABBA- 13.0327
BSE| 3 Singlet -ABBA- 15.8729
BSE| 4 Singlet -ABBA- 15.8729
```
Afterwards, the single-particle transitions, i.e. the eigenvector elements $X_{ia}^{(n)}$, with
$|X_{ia}^{(n)}|$ > [EPS_X](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.EPS_X), are
printed:
Afterwards, the single-particle transitions, i.e. the eigenvector elements $X_{ia}^{(n)}$ and
$Y_{ia}^{(n)}$, with $|X_{ia}^{(n)}|$ >
[EPS_X](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.EPS_X) or $|Y_{ia}^{(n)}|$ >
[EPS_X](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.EPS_X), are printed. The third
column is either
- $\Rightarrow$: for excitations, i.e. entries of X_ia^n, or
- $\Leftarrow$: for deexcitations, i.e. entries of Y_ia^n.
```none
BSE| Excitations are built up by the following single-particle transitions,
BSE| neglecting contributions where |X_ia^n| < 0.10 :
BSE| -- Quick reminder: HOMO i = 1 and LUMO a = 2 --
BSE|
BSE| Excitation n i => a TDA/full BSE |X_ia^n|
BSE|
BSE| 1 1 => 2 -full- 0.6630
BSE| 1 1 => 4 -full- 0.2549
BSE|
BSE| 2 1 => 3 -full- 0.7059
BSE|
BSE| 3 1 => 6 -full- 0.7076
BSE|
BSE| 4 1 => 5 -full- 0.7076
BSE| The following single-particle transitions have significant contributions,
BSE| i.e. |X_ia^n| > 0.100 or |Y_ia^n| > 0.100, respectively :
BSE| -- Quick reminder: HOMO i = 1 and LUMO a = 2 --
BSE|
BSE| Excitation n i =>/<= a TDA/ABBA |X_ia^n|/|Y_ia^n|
BSE|
BSE| 1 1 => 2 -ABBA- 0.6687
BSE| 1 1 => 4 -ABBA- 0.2431
BSE|
BSE| 2 1 => 3 -ABBA- 0.7061
BSE|
BSE| 3 1 => 5 -ABBA- 0.7077
BSE|
BSE| 4 1 => 6 -ABBA- 0.7077
```
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$ ).
$X_{i=\text{HOMO},a=\text{LUMO}}^{(n=1)}= 0.6687$, 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.2431$ ).
In the next section, the optical properties, i.e. the transition moments $d^{(n)}_{r}$ and the
oscillator strengths $f^{(n)}$, are printed:
```none
BSE| Optical properties from solving the BSE without the TDA:
BSE|
BSE| Excitation n TDA/ABBA d_x^n d_y^n d_z^n Osc. strength f^n
BSE| 1 -ABBA- 0.000 -0.000 -1.191 0.418
BSE| 2 -ABBA- 0.000 0.000 0.000 0.000
BSE| 3 -ABBA- -0.836 0.713 -0.000 0.469
BSE| 4 -ABBA- -0.713 -0.836 -0.000 0.469
```
We can see that $n=1,3,4$ are optically active since they have a nonzero oscillator strength
$f^{(n)}$.
In case [BSE_SPECTRUM](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM) is
invoked, the frequency-resolved optical absorption spectrum is printed in addition. For each
specified $\eta$ in
[ETA_LIST](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM.ETA_LIST), a
separate file, e.g. `BSE-ABBA-eta=0.010.spectrum`, is created, which starts with
```none
# Imaginary part of polarizability α_{µ µ'}(ω) = \sum_n [2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] from Bethe Salpeter equation for method -ABBA-
# Frequency (eV) Im{α_{avg}(ω)} Im{α_xx(ω)} Im{α_xy(ω)} Im{α_xz(ω)} Im{α_yx(ω)} Im{α_yy(ω)} Im{α_yz(ω)} Im{α_zx(ω)} Im{α_zy(ω)} Im{α_zz(ω)}
0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.10000000 0.00005283 0.00003298 -0.00000000 0.00000000 -0.00000000 0.00003298 0.00000000 0.00000000 0.00000000 0.00009251
```
The next table summarizes some of the exciton descriptors of each excitation level $n$:
```none
BSE| Exciton descriptors from solving the BSE without the TDA:
BSE|
BSE| n c_n d_eh [Å] σ_e [Å] σ_h [Å] d_exc [Å] R_eh
BSE| 1 1.026 0.0000 2.0542 0.8576 2.2521 -0.0176
BSE| 2 1.004 0.0000 2.9922 0.8578 3.1127 -0.0000
BSE| 3 1.005 0.0000 1.5859 0.8568 1.8134 -0.0077
BSE| 4 1.005 0.0000 1.5859 0.8568 1.8134 -0.0077
```
where we obtain the normalization factor
$c_n = \langle \Psi_\text{excitation}^{(n)} | \Psi_\text{excitation}^{(n)}\rangle $, the vectorial
distance between electron and hole $d_{e\rightarrow h}$, sizes of electron $\sigma_e$ and hole
$\sigma_h$, the exciton size $d_\mathrm{exc}$ and the electron-hole correlation coefficient
$R_{eh}$.
Further, the expectation values $\langle \mathbf{r}_e \rangle_\mathrm{exc}$ and
$\langle \mathbf{r}_h \rangle_\mathrm{exc}$ as well as the reference point of origin and the atomic
geometry are printed.
The last section summarizes the characteristics of the NTO analysis:
```none
BSE| Number of excitations, for which NTOs are computed 2
BSE| Threshold for oscillator strength f^n ---
BSE| Threshold for NTO weights (λ_I)^2 0.010
BSE|
BSE| NTOs from solving the BSE without the TDA:
BSE|
BSE| Excitation n TDA/ABBA Index of NTO I NTO weights (λ_I)^2
BSE|
BSE| 1 -ABBA- 1 1.01320
BSE| 1 -ABBA- 2 0.01320
BSE|
BSE| 2 -ABBA- 1 1.00210
```
In the input file, we have specified a list of three excitation levels, no constraint on the
oscillator strengths $f^{(n)}$ and used the default threshold on the weights, given by
$\lambda_I^2 \leq 0.010$, which is repeated at the beginning of the section. Therefore, we obtain
two NTO pairs for $n=1$ with significant weight and one NTO pair for $n=2$ and $n=8$, respectively.
For $n=1$, we can verify that the NTO weights satisfy
$1.01320 + 0.01320 \approx \sum_I {\left( \lambda_I^{(n)}\right)}^2 = c_n \approx 1.026$ when
comparing to the output of the exciton descriptor section.
### 3.3 Large scale calculations
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
provide a sample output of a BSE calculation on a nanographene with 206 atoms which has a memory
requirement of 2.2 TB RAM:
consumption. [Here](https://github.com/cp2k/cp2k-examples/tree/master/bethe-salpeter/Nanographene),
you can find a sample output of a BSE@evGW0@PBE calculation on a nanographene with 206 atoms, which
has a peak memory requirement of 2.5 TB RAM:
```none
BSE| Cutoff occupied orbitals [eV] 80.000
BSE| Cutoff empty orbitals [eV] 10.000
BSE| First occupied index 155
BSE| Last empty index (not MO index!) 517
BSE| Energy of first occupied index [eV] -24.774
BSE| Energy of last empty index [eV] 7.863
BSE| Energy difference of first occupied index to HOMO [eV] 19.487
BSE| Energy difference of last empty index to LUMO [eV] 9.639
BSE| Number of GW-corrected occupied MOs 400
BSE| Number of GW-corrected empty MOs 600
BSE|
BSE| Total peak memory estimate from BSE [GB] 2221.591
BSE| Peak memory estimate per MPI rank from BSE [GB] 4.628
BSE| Cutoff occupied orbitals [eV] 80.000
BSE| Cutoff empty orbitals [eV] 10.000
BSE| First occupied index 155
BSE| Last empty index (not MO index!) 324
BSE| Energy of first occupied index [eV] -25.080
BSE| Energy of last empty index [eV] 9.632
BSE| Energy difference of first occupied index to HOMO [eV] 20.289
BSE| Energy difference of last empty index to LUMO [eV] 10.549
BSE| Number of GW-corrected occupied MOs 334
BSE| Number of GW-corrected empty MOs 324
BSE|
BSE| Total peak memory estimate from BSE [GB] 8.725E+02
BSE| Peak memory estimate per MPI rank from BSE [GB] 5.453
```
To enable BSE calculations on large molecules, we recommend to use a large supercomputer and setting
the keywords
To enable BSE calculations on large molecules, we recommend to use large clusters with increased RAM
and explicitly setting the keywords
[ENERGY_CUTOFF_OCC](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.ENERGY_CUTOFF_OCC)
and
[ENERGY_CUTOFF_EMPTY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.ENERGY_CUTOFF_EMPTY),

2
src/CMakeLists.txt
View file

@ -72,6 +72,8 @@ list(
bse_main.F
bse_full_diag.F
bse_iterative.F
bse_print.F
bse_properties.F
bse_util.F
bsse.F
cell_methods.F

313
src/bse_full_diag.F
View file

@ -13,11 +13,18 @@
! **************************************************************************************************
MODULE bse_full_diag
USE bse_print, ONLY: print_excitation_energies,&
print_exciton_descriptors,&
print_optical_properties,&
print_output_header,&
print_transition_amplitudes
USE bse_properties, ONLY: calculate_NTOs,&
exciton_descr_type,&
get_exciton_descriptors,&
get_oscillator_strengths
USE bse_util, ONLY: comp_eigvec_coeff_BSE,&
compute_absorption_spectrum,&
fm_general_add_bse,&
get_oscillator_strengths,&
print_transition_probabilities
get_multipoles_mo
USE cp_blacs_env, ONLY: cp_blacs_env_create,&
cp_blacs_env_release,&
cp_blacs_env_type
@ -39,7 +46,6 @@ MODULE bse_full_diag
USE message_passing, ONLY: mp_para_env_type
USE mp2_types, ONLY: mp2_type
USE parallel_gemm_api, ONLY: parallel_gemm
USE physcon, ONLY: evolt
USE qs_environment_types, ONLY: qs_environment_type
#include "./base/base_uses.f90"
@ -100,12 +106,12 @@ CONTAINS
CALL timeset(routineN, handle)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A'
END IF
!Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
SELECT CASE (mp2_env%bse%bse_spin_config)
CASE (bse_singlet)
alpha = 2.0_dp
CASE (bse_triplet)
@ -138,7 +144,7 @@ CONTAINS
matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
matrix_c=fm_A)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb'
END IF
@ -147,7 +153,7 @@ CONTAINS
matrix_a=fm_mat_S_bar_ij_bse, matrix_b=fm_mat_S_ab_bse, beta=0.0_dp, &
matrix_c=fm_W)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
END IF
@ -227,12 +233,12 @@ CONTAINS
CALL timeset(routineN, handle)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
END IF
! Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
SELECT CASE (mp2_env%bse%bse_spin_config)
CASE (bse_singlet)
alpha = 2.0_dp
CASE (bse_triplet)
@ -247,7 +253,7 @@ CONTAINS
CALL cp_fm_create(fm_W, fm_struct_v, name="fm_W_ibaj")
CALL cp_fm_set_all(fm_W, 0.0_dp)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
END IF
! v_ia,jb = \sum_P B^P_ia B^P_jb
@ -335,7 +341,7 @@ CONTAINS
CALL cp_fm_create(fm_work_product, fm_A%matrix_struct)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
END IF
@ -385,7 +391,7 @@ CONTAINS
fm_C)
CALL cp_fm_release(fm_work_product)
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
END IF
@ -469,6 +475,7 @@ CONTAINS
matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
matrix_c=fm_mat_eigvec_transform_sum)
CALL cp_fm_release(fm_sqrt_A_minus_B)
! This normalizes the eigenvectors
CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_sum, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
@ -479,6 +486,8 @@ CONTAINS
matrix_c=fm_mat_eigvec_transform_diff)
CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
CALL cp_fm_release(fm_eigvec_Z)
! This normalizes the eigenvectors
CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_diff, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
! Now, we add the two equations to obtain X_n
@ -497,7 +506,7 @@ CONTAINS
CALL cp_fm_release(fm_mat_eigvec_transform_diff)
CALL cp_fm_release(fm_mat_eigvec_transform_sum)
CALL success_message(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
CALL postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
homo, virtual, homo_irred, unit_nr, &
.FALSE., fm_eigvec_Y)
@ -553,7 +562,7 @@ CONTAINS
CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
CPASSERT(diag_info == 0)
CALL success_message(Exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
CALL postprocess_bse(Exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
homo, virtual, homo_irred, unit_nr, .TRUE.)
CALL cp_fm_release(fm_eigvec)
@ -577,7 +586,7 @@ CONTAINS
!> \param flag_TDA ...
!> \param fm_eigvec_Y ...
! **************************************************************************************************
SUBROUTINE success_message(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
homo, virtual, homo_irred, unit_nr, &
flag_TDA, fm_eigvec_Y)
@ -590,18 +599,23 @@ CONTAINS
LOGICAL, OPTIONAL :: flag_TDA
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message'
CHARACTER(LEN=*), PARAMETER :: routineN = 'postprocess_bse'
CHARACTER(LEN=10) :: info_approximation, multiplet
INTEGER :: handle, i_exc
INTEGER :: handle, idir, n_moments_di, &
n_moments_quad
REAL(KIND=dp) :: alpha
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dipole_exc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
TYPE(exciton_descr_type), ALLOCATABLE, &
DIMENSION(:) :: exc_descr
CALL timeset(routineN, handle)
!Prepare variables for printing
IF (mp2_env%ri_g0w0%bse_spin_config == 0) THEN
IF (mp2_env%bse%bse_spin_config == 0) THEN
multiplet = "Singlet"
alpha = 2.0_dp
ELSE
@ -617,190 +631,97 @@ CONTAINS
info_approximation = "-ABBA-"
END IF
!Get information about eigenvector
n_moments_di = 3
n_moments_quad = 9
! Compute BSE dipoles and oscillator strengths - Keep in memory for later usage
! Need dipoles also for spatial expectation values, which are well-defined also for triplets
ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
ALLOCATE (ref_point_multipole(3))
! Obtain dipoles in MO basis
CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
qs_env, mo_coeff, ref_point_multipole, 1, &
homo, virtual, fm_eigvec_X%matrix_struct%context)
! Compute exciton descriptors from these multipoles
IF (mp2_env%bse%num_print_exc_descr > 0) THEN
! Obtain quadrupoles in MO basis
ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
qs_env, mo_coeff, ref_point_multipole, 2, &
homo, virtual, fm_eigvec_X%matrix_struct%context)
CALL get_exciton_descriptors(exc_descr, fm_eigvec_X, &
fm_dipole_ij_trunc, fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
homo, virtual, unit_nr, &
mp2_env, qs_env, &
fm_eigvec_Y)
END IF
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA) *'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
'the BSE within the TDA:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb X_jb^n ) = Ω^n X_ia^n'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Full ("ABBA") Bethe Salpeter equation (BSE) (i.e. without TDA) *'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
'the BSE without the TDA:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n| |1 0| |X^n|'
WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', ' sum_jb ( A_ia,jb X_jb^n + B_ia,jb Y_jb^n ) = Ω^n X_ia^n'
WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', '- sum_jb ( B_ia,jb X_jb^n + A_ia,jb Y_jb^n ) = Ω^n Y_ia^n'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab'
IF (.NOT. flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
END IF
IF (.NOT. flag_TDA) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', 'The BSE is solved for Ω^n and X_ia^n as a hermitian problem, e.g. Eq.(42)'
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
! WRITE (unit_nr, '(T2,A4,T7,A69)') 'BSE|', 'C_ia,jb = sum_kc,ld ((A-B)^0.5)_ia,kc (A+B)_kc,ld ((A-B)^0.5)_ld,jb'
! WRITE (unit_nr, '(T2,A4)') 'BSE|'
! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X+Y)_ia,n = sum_jb,m ((A-B)^0.5)ia,jb Z_jb,m (Ω_m)^-0.5'
! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X-Y)_ia,n = sum_jb,m ((A-B)^-0.5)ia,jb Z_jb,m (Ω_m)^0.5'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A7,T19,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
WRITE (unit_nr, '(T2,A4,T7,A7,T19,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
WRITE (unit_nr, '(T2,A4,T7,A3,T19,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
WRITE (unit_nr, '(T2,A4,T7,A6,T18,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
WRITE (unit_nr, '(T2,A4,T7,A6,T18,A27)') 'BSE|', 'W_...:', 'Screened direct interaction'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T11,A12,T26,A11,T44,A8,T55,A27)') 'BSE|', &
'Excitation n', "Spin Config", 'TDA/ABBA', 'Excitation energy Ω^n (eV)'
IF (mp2_env%bse%bse_spin_config == 0) THEN
CALL get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
trans_mom_bse, oscill_str, polarizability_residues, &
mp2_env, homo, virtual, unit_nr, &
fm_eigvec_Y)
END IF
!prints actual energies values
IF (unit_nr > 0) THEN
DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
WRITE (unit_nr, '(T2,A4,T7,I16,T30,A7,T46,A6,T59,F22.4)') &
'BSE|', i_exc, multiplet, info_approximation, Exc_ens(i_exc)*evolt
END DO
END IF
!prints single particle transitions
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A61)') &
'BSE|', "Single-particle transitions are built up by (de-)excitations,"
WRITE (unit_nr, '(T2,A4,T7,A18)') &
'BSE|', "which we denote by"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T20,A2,T30,A40)') &
'BSE|', "=>", "for excitations, i.e. entries of X_ia^n,"
WRITE (unit_nr, '(T2,A4,T20,A2,T30,A42)') &
'BSE|', "<=", "for deexcitations, i.e. entries of Y_ia^n."
WRITE (unit_nr, '(T2,A4)') &
'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A73)') &
'BSE|', "The following single-particle transitions have significant contributions,"
WRITE (unit_nr, '(T2,A4,T7,A16,F5.3,A15,F5.3,A16)') &
'BSE|', "i.e. |X_ia^n| < ", mp2_env%ri_g0w0%eps_x, " or |Y_ia^n| < ", &
mp2_env%ri_g0w0%eps_x, ", respectively :"
WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
homo_irred, ' and LUMO a =', homo_irred + 1, " --"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A12,T30,A1,T32,A5,T42,A1,T49,A8,T64,A17)') &
"BSE|", "Excitation n", "i", "=>/<=", "a", 'TDA/ABBA', "|X_ia^n|/|Y_ia^n|"
! Prints basic definitions used in BSE calculation
CALL print_output_header(homo, virtual, homo_irred, flag_TDA, &
multiplet, alpha, unit_nr)
! Prints excitation energies up to user-specified number
CALL print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
info_approximation, mp2_env, unit_nr)
! Print single particle transition amplitudes, i.e. components of eigenvectors X and Y
CALL print_transition_amplitudes(fm_eigvec_X, homo, virtual, homo_irred, &
info_approximation, mp2_env, unit_nr, fm_eigvec_Y)
! Prints optical properties, if state is a singlet
CALL print_optical_properties(Exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
homo, virtual, homo_irred, flag_TDA, &
info_approximation, mp2_env, unit_nr)
! Print exciton descriptors if keyword is invoked
IF (mp2_env%bse%num_print_exc_descr > 0) THEN
CALL print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, mp2_env, qs_env)
END IF
DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
! Compute and print excitation wavefunctions
IF (mp2_env%bse%do_nto_analysis) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A47)') &
'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
!Iterate through eigenvector and print values above threshold
CALL print_transition_probabilities(fm_eigvec_X, "=>", info_approximation, &
i_exc, virtual, homo, homo_irred, &
unit_nr, mp2_env)
IF (PRESENT(fm_eigvec_Y)) THEN
CALL print_transition_probabilities(fm_eigvec_Y, "<=", info_approximation, &
i_exc, virtual, homo, homo_irred, &
unit_nr, mp2_env)
END IF
END DO
! Compute BSE dipoles and oscillator strengths
CALL get_oscillator_strengths(fm_eigvec_X, Exc_ens, &
dipole_exc, oscill_str, &
mp2_env, qs_env, mo_coeff, &
homo, virtual, unit_nr, &
fm_eigvec_Y)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A63)') &
'BSE|', "Dipoles d_r^n and oscillator strength f^n of excitation level n"
WRITE (unit_nr, '(T2,A4,T7,A17)') &
'BSE|', "are obtained from"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T12,A51)') &
'BSE|', "d_r^n = sqrt(2) sum_ia < ψ_i | µ | ψ_a > X_ia^n"
ELSE
WRITE (unit_nr, '(T2,A4,T12,A63)') &
'BSE|', "d_r^n = sqrt(2) sum_ia < ψ_i | µ | ψ_a > ( X_ia^n + Y_ia^n )"
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T20,A44)') &
'BSE|', "f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A19)') &
'BSE|', "where we introduced"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
'BSE|', "ψ_i:", "occupied molecular orbitals,"
WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
'BSE|', "ψ_a:", "empty molecular orbitals and"
WRITE (unit_nr, '(T2,A4,T7,A3,T15,A16)') &
'BSE|', "µ:", "dipole operator."
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A28)') &
'BSE|', "prelim Ref.: Eqs. (23), (24)"
WRITE (unit_nr, '(T2,A4,T7,A71)') &
'BSE|', "in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T35,A8,T45,A8,T55,A8,T65,A16)') 'BSE|', &
'Excitation n', "TDA/ABBA", "Dipole X", "Dipole Y", "Dipole Z", 'Osc. strength'
DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T35,F8.3,T45,F8.3,T55,F8.3,T65,F16.3)') &
'BSE|', i_exc, info_approximation, dipole_exc(1, 1, i_exc), dipole_exc(2, 1, i_exc), &
dipole_exc(3, 1, i_exc), oscill_str(i_exc)
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
CALL calculate_NTOs(fm_eigvec_X, fm_eigvec_Y, &
mo_coeff, homo, virtual, &
info_approximation, &
oscill_str, &
qs_env, unit_nr, mp2_env)
END IF
! Compute and print the absorption spectrum to external file
IF (mp2_env%ri_g0w0%bse_print_spectrum) THEN
CALL compute_absorption_spectrum(oscill_str, Exc_ens, &
info_approximation, unit_nr, mp2_env)
DO idir = 1, n_moments_di
CALL cp_fm_release(fm_dipole_ai_trunc(idir))
CALL cp_fm_release(fm_dipole_ij_trunc(idir))
CALL cp_fm_release(fm_dipole_ab_trunc(idir))
END DO
IF (mp2_env%bse%num_print_exc_descr > 0) THEN
DO idir = 1, n_moments_quad
CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
END DO
DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
DEALLOCATE (exc_descr)
END IF
DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
DEALLOCATE (ref_point_multipole)
IF (mp2_env%bse%bse_spin_config == 0) THEN
DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
END IF
DEALLOCATE (oscill_str, dipole_exc)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
@ -808,6 +729,6 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE success_message
END SUBROUTINE postprocess_bse
END MODULE bse_full_diag

24
src/bse_iterative.F
View file

@ -89,15 +89,15 @@ CONTAINS
num_res_unconverged = -1
num_exact_en_unconverged = -1
bse_davidson_abort_cond = mp2_env%ri_g0w0%davidson_abort_cond
num_exc_en = mp2_env%ri_g0w0%num_exc_en
num_add_start_z_space = mp2_env%ri_g0w0%num_add_start_z_space
fac_max_z_space = mp2_env%ri_g0w0%fac_max_z_space
num_new_t = mp2_env%ri_g0w0%num_new_t
num_davidson_iter = mp2_env%ri_g0w0%num_davidson_iter
eps_res = mp2_env%ri_g0w0%eps_res
eps_exc_en = mp2_env%ri_g0w0%eps_exc_en
z_space_energy_cutoff = mp2_env%ri_g0w0%z_space_energy_cutoff
bse_davidson_abort_cond = mp2_env%bse%davidson_abort_cond
num_exc_en = mp2_env%bse%num_exc_en
num_add_start_z_space = mp2_env%bse%num_add_start_z_space
fac_max_z_space = mp2_env%bse%fac_max_z_space
num_new_t = mp2_env%bse%num_new_t
num_davidson_iter = mp2_env%bse%num_davidson_iter
eps_res = mp2_env%bse%eps_res
eps_exc_en = mp2_env%bse%eps_exc_en
z_space_energy_cutoff = mp2_env%bse%z_space_energy_cutoff
num_Z_vectors_init = num_exc_en + num_add_start_z_space
@ -244,7 +244,7 @@ CONTAINS
END IF
IF (davidson_converged == 1) THEN
CALL success_message(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, &
CALL postprocess_bse(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, &
bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, &
num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space, max_res_norm, &
max_en_diff, num_en_unconverged, bse_davidson_abort_cond_string, &
@ -469,7 +469,7 @@ CONTAINS
!> \param success_abort_string ...
!> \param z_space_energy_cutoff ...
! **************************************************************************************************
SUBROUTINE success_message(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, &
SUBROUTINE postprocess_bse(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, &
bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, &
num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space, max_res_norm, &
max_en_diff, num_en_unconverged, bse_davidson_abort_cond_string, &
@ -487,7 +487,7 @@ CONTAINS
CHARACTER(LEN=10) :: success_abort_string
REAL(KIND=dp) :: z_space_energy_cutoff
CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message'
CHARACTER(LEN=*), PARAMETER :: routineN = 'postprocess_bse'
CHARACTER(LEN=10) :: multiplet
INTEGER :: handle, i

17
src/bse_main.F
View file

@ -20,10 +20,11 @@ MODULE bse_main
diagonalize_C
USE bse_iterative, ONLY: do_subspace_iterations,&
fill_local_3c_arrays
USE bse_util, ONLY: deallocate_matrices_bse,&
USE bse_print, ONLY: print_BSE_start_flag
USE bse_util, ONLY: adapt_BSE_input_params,&
deallocate_matrices_bse,&
estimate_BSE_resources,&
mult_B_with_W,&
print_BSE_start_flag,&
truncate_BSE_matrices
USE cp_fm_types, ONLY: cp_fm_release,&
cp_fm_type
@ -116,7 +117,7 @@ CONTAINS
my_do_tda = .FALSE.
my_do_abba = .FALSE.
!Method: Iterative or full diagonalization
SELECT CASE (mp2_env%ri_g0w0%bse_diag_method)
SELECT CASE (mp2_env%bse%bse_diag_method)
CASE (bse_iterdiag)
my_do_iterat_diag = .TRUE.
!MG: Basics of the Davidson solver are implemented, but not rigorously checked.
@ -125,7 +126,7 @@ CONTAINS
my_do_fulldiag = .TRUE.
END SELECT
!Approximation: TDA and/or full ABBA matrix
SELECT CASE (mp2_env%ri_g0w0%flag_tda)
SELECT CASE (mp2_env%bse%flag_tda)
CASE (bse_tda)
my_do_tda = .TRUE.
CASE (bse_abba)
@ -140,7 +141,7 @@ CONTAINS
! Link BSE debug flag against debug print level
logger => cp_get_default_logger()
IF (logger%iter_info%print_level == debug_print_level) THEN
mp2_env%ri_g0w0%bse_debug_print = .TRUE.
mp2_env%bse%bse_debug_print = .TRUE.
END IF
CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync()
@ -168,6 +169,8 @@ CONTAINS
gd_array, color_sub, para_env)
END IF
CALL adapt_BSE_input_params(homo_red, virtual_red, unit_nr, mp2_env)
IF (my_do_fulldiag) THEN
! Quick estimate of memory consumption and runtime of diagonalizations
CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
@ -220,12 +223,12 @@ CONTAINS
CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
fm_mat_Q_static_bse_gemm)
fm_mat_Q_static_bse_gemm, mp2_env)
DEALLOCATE (Eigenval_reduced)
IF (my_do_iterat_diag) THEN
! Contains untested Block-Davidson algorithm
CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
B_iaQ_bse_local, homo(1), virtual(1), mp2_env%ri_g0w0%bse_spin_config, unit_nr, &
B_iaQ_bse_local, homo(1), virtual(1), mp2_env%bse%bse_spin_config, unit_nr, &
Eigenval(:, 1, 1), para_env, mp2_env)
! Deallocate local 3c-B-matrices
DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)

793
src/bse_print.F Normal file
View file

@ -0,0 +1,793 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Routines for printing information in context of the BSE calculation
!> \par History
!> 10.2024 created [Maximilian Graml]
! **************************************************************************************************
MODULE bse_print
USE atomic_kind_types, ONLY: get_atomic_kind
USE bse_properties, ONLY: compute_absorption_spectrum,&
exciton_descr_type
USE bse_util, ONLY: filter_eigvec_contrib
USE cp_fm_types, ONLY: cp_fm_get_info,&
cp_fm_type
USE kinds, ONLY: dp
USE mp2_types, ONLY: mp2_type
USE particle_types, ONLY: particle_type
USE physcon, ONLY: angstrom,&
evolt
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_print'
PUBLIC :: print_BSE_start_flag, fm_write_thresh, print_excitation_energies, &
print_output_header, print_transition_amplitudes, print_optical_properties, &
print_exciton_descriptors
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param bse_tda ...
!> \param bse_abba ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE print_BSE_start_flag(bse_tda, bse_abba, unit_nr)
LOGICAL, INTENT(IN) :: bse_tda, bse_abba
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_BSE_start_flag'
INTEGER :: handle
CALL timeset(routineN, handle)
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) ' '
WRITE (unit_nr, '(T2,A79)') '*******************************************************************************'
WRITE (unit_nr, '(T2,A79)') '** **'
WRITE (unit_nr, '(T2,A79)') '** Bethe Salpeter equation (BSE) for excitation energies **'
IF (bse_tda .AND. bse_abba) THEN
WRITE (unit_nr, '(T2,A79)') '** solved with and without Tamm-Dancoff approximation (TDA) **'
ELSE IF (bse_tda) THEN
WRITE (unit_nr, '(T2,A79)') '** solved with Tamm-Dancoff approximation (TDA) **'
ELSE
WRITE (unit_nr, '(T2,A79)') '** solved without Tamm-Dancoff approximation (TDA) **'
END IF
WRITE (unit_nr, '(T2,A79)') '** **'
WRITE (unit_nr, '(T2,A79)') '*******************************************************************************'
WRITE (unit_nr, *) ' '
END IF
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param homo ...
!> \param virtual ...
!> \param homo_irred ...
!> \param flag_TDA ...
!> \param multiplet ...
!> \param alpha ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE print_output_header(homo, virtual, homo_irred, flag_TDA, &
multiplet, alpha, unit_nr)
INTEGER, INTENT(IN) :: homo, virtual, homo_irred
LOGICAL, INTENT(IN) :: flag_TDA
CHARACTER(LEN=10), INTENT(IN) :: multiplet
REAL(KIND=dp), INTENT(IN) :: alpha
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_output_header'
INTEGER :: handle
CALL timeset(routineN, handle)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA) *'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
'the BSE within the TDA:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb X_jb^n ) = Ω^n X_ia^n'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Full ("ABBA") Bethe Salpeter equation (BSE) (i.e. without TDA) *'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
'the BSE without the TDA:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n| |1 0| |X^n|'
WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', ' sum_jb ( A_ia,jb X_jb^n + B_ia,jb Y_jb^n ) = Ω^n X_ia^n'
WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', '- sum_jb ( B_ia,jb X_jb^n + A_ia,jb Y_jb^n ) = Ω^n Y_ia^n'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab'
IF (.NOT. flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
END IF
IF (.NOT. flag_TDA) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', 'The BSE is solved for Ω^n and X_ia^n as a hermitian problem, e.g. Eq.(42)'
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A7,T19,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
WRITE (unit_nr, '(T2,A4,T7,A7,T19,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
WRITE (unit_nr, '(T2,A4,T7,A3,T19,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
WRITE (unit_nr, '(T2,A4,T7,A6,T18,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
WRITE (unit_nr, '(T2,A4,T7,A6,T18,A27)') 'BSE|', 'W_...:', 'Screened direct interaction'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param Exc_ens ...
!> \param homo ...
!> \param virtual ...
!> \param flag_TDA ...
!> \param multiplet ...
!> \param info_approximation ...
!> \param mp2_env ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
info_approximation, mp2_env, unit_nr)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
INTEGER, INTENT(IN) :: homo, virtual
LOGICAL, INTENT(IN) :: flag_TDA
CHARACTER(LEN=10), INTENT(IN) :: multiplet, info_approximation
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_excitation_energies'
INTEGER :: handle, i_exc
CALL timeset(routineN, handle)
IF (unit_nr > 0) THEN
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T11,A12,T26,A11,T44,A8,T55,A27)') 'BSE|', &
'Excitation n', "Spin Config", 'TDA/ABBA', 'Excitation energy Ω^n (eV)'
END IF
!prints actual energies values
IF (unit_nr > 0) THEN
DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
WRITE (unit_nr, '(T2,A4,T7,I16,T30,A7,T46,A6,T59,F22.4)') &
'BSE|', i_exc, multiplet, info_approximation, Exc_ens(i_exc)*evolt
END DO
END IF
CALL timestop(handle)
END SUBROUTINE print_excitation_energies
! **************************************************************************************************
!> \brief ...
!> \param fm_eigvec_X ...
!> \param homo ...
!> \param virtual ...
!> \param homo_irred ...
!> \param info_approximation ...
!> \param mp2_env ...
!> \param unit_nr ...
!> \param fm_eigvec_Y ...
! **************************************************************************************************
SUBROUTINE print_transition_amplitudes(fm_eigvec_X, homo, virtual, homo_irred, &
info_approximation, mp2_env, unit_nr, fm_eigvec_Y)
TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
INTEGER, INTENT(IN) :: homo, virtual, homo_irred
CHARACTER(LEN=10), INTENT(IN) :: info_approximation
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
INTEGER, INTENT(IN) :: unit_nr
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_transition_amplitudes'
INTEGER :: handle, i_exc
CALL timeset(routineN, handle)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A61)') &
'BSE|', "Single-particle transitions are built up by (de-)excitations,"
WRITE (unit_nr, '(T2,A4,T7,A18)') &
'BSE|', "which we denote by"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T20,A2,T30,A40)') &
'BSE|', "=>", "for excitations, i.e. entries of X_ia^n,"
WRITE (unit_nr, '(T2,A4,T20,A2,T30,A42)') &
'BSE|', "<=", "for deexcitations, i.e. entries of Y_ia^n."
WRITE (unit_nr, '(T2,A4)') &
'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A73)') &
'BSE|', "The following single-particle transitions have significant contributions,"
WRITE (unit_nr, '(T2,A4,T7,A16,F5.3,A15,F5.3,A16)') &
'BSE|', "i.e. |X_ia^n| > ", mp2_env%bse%eps_x, " or |Y_ia^n| > ", &
mp2_env%bse%eps_x, ", respectively :"
WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
homo_irred, ' and LUMO a =', homo_irred + 1, " --"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A12,T30,A1,T32,A5,T42,A1,T49,A8,T64,A17)') &
"BSE|", "Excitation n", "i", "=>/<=", "a", 'TDA/ABBA', "|X_ia^n|/|Y_ia^n|"
END IF
DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
!Iterate through eigenvector and print values above threshold
CALL print_transition_amplitudes_core(fm_eigvec_X, "=>", info_approximation, &
i_exc, virtual, homo, homo_irred, &
unit_nr, mp2_env)
IF (PRESENT(fm_eigvec_Y)) THEN
CALL print_transition_amplitudes_core(fm_eigvec_Y, "<=", info_approximation, &
i_exc, virtual, homo, homo_irred, &
unit_nr, mp2_env)
END IF
END DO
CALL timestop(handle)
END SUBROUTINE print_transition_amplitudes
! **************************************************************************************************
!> \brief ...
!> \param Exc_ens ...
!> \param oscill_str ...
!> \param trans_mom_bse ...
!> \param polarizability_residues ...
!> \param homo ...
!> \param virtual ...
!> \param homo_irred ...
!> \param flag_TDA ...
!> \param info_approximation ...
!> \param mp2_env ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE print_optical_properties(Exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
homo, virtual, homo_irred, flag_TDA, &
info_approximation, mp2_env, unit_nr)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens, oscill_str
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: trans_mom_bse, polarizability_residues
INTEGER, INTENT(IN) :: homo, virtual, homo_irred
LOGICAL, INTENT(IN) :: flag_TDA
CHARACTER(LEN=10), INTENT(IN) :: info_approximation
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_optical_properties'
INTEGER :: handle, i_exc
CALL timeset(routineN, handle)
! Discriminate between singlet and triplet, since triplet state can't couple to light
! and therefore calculations of dipoles etc are not necessary
IF (mp2_env%bse%bse_spin_config == 0) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A60)') &
'BSE|', "Transition moments d_r^n (with r∈(x,y,z), in atomic units)"
WRITE (unit_nr, '(T2,A4,T7,A67)') &
'BSE|', "and oscillator strength f^n of excitation level n are obtained from"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T10,A)') &
'BSE|', "d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > X_ia^n"
ELSE
WRITE (unit_nr, '(T2,A4,T10,A)') &
'BSE|', "d_r^n = sum_ia sqrt(2) < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n )"
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T14,A)') &
'BSE|', "f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A19)') &
'BSE|', "where we introduced"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
'BSE|', "ψ_i:", "occupied molecular orbitals,"
WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
'BSE|', "ψ_a:", "empty molecular orbitals and"
WRITE (unit_nr, '(T2,A4,T9,A2,T14,A18)') &
'BSE|', "r:", "position operator."
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A28)') &
'BSE|', "prelim Ref.: Eqs. (23), (24)"
WRITE (unit_nr, '(T2,A4,T7,A71)') &
'BSE|', "in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A55)') 'BSE|', &
'Optical properties from solving the BSE within the TDA:'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', &
'Optical properties from solving the BSE without the TDA:'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T38,A5,T48,A5,T58,A5,T64,A17)') 'BSE|', &
'Excitation n', "TDA/ABBA", "d_x^n", "d_y^n", "d_z^n", 'Osc. strength f^n'
DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T35,F8.3,T45,F8.3,T55,F8.3,T65,F16.3)') &
'BSE|', i_exc, info_approximation, trans_mom_bse(1, 1, i_exc), trans_mom_bse(2, 1, i_exc), &
trans_mom_bse(3, 1, i_exc), oscill_str(i_exc)
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Check for Thomas-Reiche-Kuhn sum rule'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T35,A)') 'BSE|', &
'N_e = Σ_n f^n'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A24,T65,I16)') 'BSE|', &
'Number of electrons N_e:', homo_irred*2
WRITE (unit_nr, '(T2,A4,T7,A,T66,F16.3)') 'BSE|', &
'Sum over oscillator strengths Σ_n f^n :', SUM(oscill_str)
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN
CALL cp_warn(__LOCATION__, &
"Accuracy of TRK sum rule might suffer from cutoffs.")
END IF
END IF
! Compute and print the absorption spectrum to external file
IF (mp2_env%bse%bse_print_spectrum) THEN
CALL compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_ens, &
info_approximation, unit_nr, mp2_env)
END IF
ELSE
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
CALL cp_warn(__LOCATION__, &
"Requested triplet excitation cannot couple to light. "// &
"Skipping calculation of transition moments, "// &
"oscillator strengths, and spectrum.")
END IF
END IF
CALL timestop(handle)
END SUBROUTINE print_optical_properties
! **************************************************************************************************
!> \brief ...
!> \param fm_eigvec ...
!> \param direction_excitation ...
!> \param info_approximation ...
!> \param i_exc ...
!> \param virtual ...
!> \param homo ...
!> \param homo_irred ...
!> \param unit_nr ...
!> \param mp2_env ...
! **************************************************************************************************
SUBROUTINE print_transition_amplitudes_core(fm_eigvec, direction_excitation, info_approximation, &
i_exc, virtual, homo, homo_irred, &
unit_nr, mp2_env)
TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
CHARACTER(LEN=2), INTENT(IN) :: direction_excitation
CHARACTER(LEN=10), INTENT(IN) :: info_approximation
INTEGER :: i_exc, virtual, homo, homo_irred
INTEGER, INTENT(IN) :: unit_nr
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_transition_amplitudes_core'
INTEGER :: handle, k, num_entries
INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
CALL timeset(routineN, handle)
CALL filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
i_exc, virtual, num_entries, mp2_env)
! direction_excitation can be either => (means excitation; from fm_eigvec_X)
! or <= (means deexcitation; from fm_eigvec_Y)
IF (unit_nr > 0) THEN
DO k = 1, num_entries
WRITE (unit_nr, '(T2,A4,T14,I5,T26,I5,T35,A2,T38,I5,T51,A6,T65,F16.4)') &
"BSE|", i_exc, homo_irred - homo + idx_homo(k), direction_excitation, &
homo_irred + idx_virt(k), info_approximation, ABS(eigvec_entries(k))
END DO
END IF
DEALLOCATE (idx_homo)
DEALLOCATE (idx_virt)
DEALLOCATE (eigvec_entries)
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param exc_descr ...
!> \param ref_point_multipole ...
!> \param unit_nr ...
!> \param mp2_env ...
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, mp2_env, qs_env)
TYPE(exciton_descr_type), ALLOCATABLE, &
DIMENSION(:) :: exc_descr
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
INTENT(IN) :: ref_point_multipole
INTEGER, INTENT(IN) :: unit_nr
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_exciton_descriptors'
CHARACTER(LEN=1), DIMENSION(3) :: array_direction_str
INTEGER :: handle, i_dir, i_exc
REAL(KIND=dp) :: d_eh_dir, d_exc_dir, sigma_e_dir, &
sigma_h_dir
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, particle_set=particle_set)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Exciton descriptors for excitation level n are given by'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'd_eh = | <r_h - r_e>_exc |'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'σ_e = sqrt( <r_e^2>_exc - <r_e>_exc^2 )'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'σ_h = sqrt( <r_h^2>_exc - <r_h>_exc^2 )'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'COV_eh = <r_e r_h>_exc - <r_e>_exc <r_h>_exc'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'd_exc = sqrt( | < |r_h - r_e|^2 >_exc )'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
' = sqrt( d_eh^2 + σ_e^2 + σ_h^2 - 2 * COV_eh )'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'R_eh = COV_eh / (σ_e * σ_h)'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'where the expectation values <.>_exc are taken with respect to the '
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'exciton wavefunction of excitation n:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (exc_descr(1)%flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T20,A)') 'BSE|', &
'𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e) ,'
ELSE
WRITE (unit_nr, '(T2,A4,T20,A)') 'BSE|', &
'𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e)'
WRITE (unit_nr, '(T2,A4,T40,A)') 'BSE|', &
'+ Y_ia^n ψ_a(r_h) ψ_i(r_e) ,'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'i.e.'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T20,A)') 'BSE|', &
'< O >_exc = < 𝚿_n | O | 𝚿_n > / < 𝚿_n | 𝚿_n > ,'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (exc_descr(1)%flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'where c_n = < 𝚿_n | 𝚿_n > = 1 within TDA.'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'where c_n = < 𝚿_n | 𝚿_n > deviates from 1 without TDA.'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Here, we introduced'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
'BSE|', "ψ_i:", "occupied molecular orbitals,"
WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
'BSE|', "ψ_a:", "empty molecular orbitals and"
WRITE (unit_nr, '(T2,A4,T9,A2,T14,A)') &
'BSE|', "r:", "position operator."
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'prelim Ref.: Eqs. (15)-(22)'
WRITE (unit_nr, '(T2,A4,T7,A,A)') 'BSE|', &
'JCTC 2018, 14, 710-725; ', &
'http://doi.org/10.1021/acs.jctc.7b01145'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (exc_descr(1)%flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Exciton descriptors from solving the BSE within the TDA:'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Exciton descriptors from solving the BSE without the TDA:'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T10,A1,6X,A3,1X,4X,A10,5X,A10,5X,A10,3X,A11,8X,A4)') 'BSE|', &
'n', 'c_n', 'd_eh [Å]', 'σ_e [Å]', 'σ_h [Å]', 'd_exc [Å]', 'R_eh'
DO i_exc = 1, mp2_env%bse%num_print_exc_descr
WRITE (unit_nr, '(T2,A4,T7,I4,4X,F5.3,1X,5(2X,F10.4))') &
'BSE|', i_exc, exc_descr(i_exc)%norm_XpY, &
exc_descr(i_exc)%diff_r_abs*angstrom, &
exc_descr(i_exc)%sigma_e*angstrom, exc_descr(i_exc)%sigma_h*angstrom, &
exc_descr(i_exc)%diff_r_sqr*angstrom, exc_descr(i_exc)%corr_e_h*angstrom
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
! For debug runs, print first d_exc separately to allow the regtests to read in
IF (mp2_env%bse%bse_debug_print) THEN
IF (exc_descr(1)%flag_TDA) THEN
WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
'Exciton descriptor d_exc with TDA for n=1 is', exc_descr(1)%diff_r_sqr*angstrom
ELSE
WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
'Exciton descriptor d_exc without TDA for n=1 is', exc_descr(1)%diff_r_sqr*angstrom
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
! Print exciton descriptor resolved per direction
IF (mp2_env%bse%print_directional_exc_descr) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'We can restrict the exciton descriptors to a specific direction,'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'e.g. the x-components are:'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'd_eh^x = | <x_h - x_e>_exc |'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'σ_e^x = sqrt( <x_e^2>_exc - <x_e>_exc^2 )'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'σ_h^x = sqrt( <x_h^2>_exc - <x_h>_exc^2 )'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'COV_eh^x = <x_e x_h>_exc - <x_e>_exc <x_h>_exc'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'd_exc^x = sqrt( | < |x_h - x_e|^2 >_exc )'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
' = sqrt( (d_eh^x)^2 + (σ_e^x)^2'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
' + (σ_h^x)^2 - 2 * (COV_eh^x) )'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (exc_descr(1)%flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Exciton descriptors per direction from solving the BSE within the TDA:'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'Exciton descriptors per direction from solving the BSE without the TDA:'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T12,A1,2X,A9,5X,A12,5X,A12,5X,A12,3X,A13)') 'BSE|', &
'n', 'r = x/y/z', 'd_eh^r [Å]', 'σ_e^r [Å]', 'σ_h^r [Å]', 'd_exc^r [Å]'
DO i_exc = 1, mp2_env%bse%num_print_exc_descr
DO i_dir = 1, 3
array_direction_str = (/"x", "y", "z"/)
d_eh_dir = ABS(exc_descr(i_exc)%r_h(i_dir) - exc_descr(i_exc)%r_e(i_dir))
sigma_e_dir = SQRT(exc_descr(i_exc)%r_e_sq(i_dir) - exc_descr(i_exc)%r_e(i_dir)**2)
sigma_h_dir = SQRT(exc_descr(i_exc)%r_h_sq(i_dir) - exc_descr(i_exc)%r_h(i_dir)**2)
d_exc_dir = SQRT(d_eh_dir**2 + sigma_e_dir**2 + sigma_h_dir**2 &
- 2*exc_descr(i_exc)%cov_e_h(i_dir))
WRITE (unit_nr, '(T2,A4,T9,I4,10X,A1,1X,4(4X,F10.4))') &
'BSE|', i_exc, array_direction_str(i_dir), &
d_eh_dir*angstrom, &
sigma_e_dir*angstrom, sigma_h_dir*angstrom, &
d_exc_dir*angstrom
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
! Print the reference atomic geometry for the exciton descriptors
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'With the center of charge as reference point r_0,'
WRITE (unit_nr, '(T2,A4,T15,A7,F10.4,A2,F10.4,A2,F10.4,A1)') 'BSE|', &
'r_0 = (', ref_point_multipole(1)*angstrom, ', ', ref_point_multipole(2)*angstrom, ', ', &
ref_point_multipole(3)*angstrom, ')'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'we further obtain r_e and r_h from solving the BSE within the TDA'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') 'BSE|', &
'Excitation n', 'x_e [Å]', 'y_e [Å]', 'z_e [Å]'
DO i_exc = 1, mp2_env%bse%num_print_exc_descr
WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
'BSE|', i_exc, &
exc_descr(i_exc)%r_e_shift(:)*angstrom
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') 'BSE|', &
'Excitation n', 'x_h [Å]', 'y_h [Å]', 'z_h [Å]'
DO i_exc = 1, mp2_env%bse%num_print_exc_descr
WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
'BSE|', i_exc, &
exc_descr(i_exc)%r_h_shift(:)*angstrom
END DO
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'The reference atomic geometry for these values is given by'
END IF
CALL write_qs_particle_coordinates_bse(particle_set, unit_nr)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
CALL timestop(handle)
END SUBROUTINE print_exciton_descriptors
! **************************************************************************************************
!> \brief Debug function to write elements of a full matrix to file, if they are larger than a given threshold
!> \param fm ...
!> \param thresh ...
!> \param header ...
!> \param unit_nr ...
!> \param abs_vals ...
! **************************************************************************************************
SUBROUTINE fm_write_thresh(fm, thresh, header, unit_nr, abs_vals)
TYPE(cp_fm_type), INTENT(IN) :: fm
REAL(KIND=dp), INTENT(IN) :: thresh
CHARACTER(LEN=*), INTENT(IN) :: header
INTEGER, INTENT(IN) :: unit_nr
LOGICAL, OPTIONAL :: abs_vals
CHARACTER(LEN=*), PARAMETER :: my_footer = " | ENDING WRITING OF MATRIX", &
routineN = 'fm_write_thresh'
INTEGER :: handle, i, j, ncol_local, nrow_local
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
LOGICAL :: my_abs_vals
CALL timeset(routineN, handle)
IF (PRESENT(abs_vals)) THEN
my_abs_vals = abs_vals
ELSE
my_abs_vals = .FALSE.
END IF
CALL cp_fm_get_info(matrix=fm, &
nrow_local=nrow_local, &
ncol_local=ncol_local, &
row_indices=row_indices, &
col_indices=col_indices)
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) header
END IF
IF (my_abs_vals) THEN
DO i = 1, nrow_local
DO j = 1, ncol_local
IF (ABS(fm%local_data(i, j)) > thresh) THEN
WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
ABS(fm%local_data(i, j))
END IF
END DO
END DO
ELSE
DO i = 1, nrow_local
DO j = 1, ncol_local
IF (ABS(fm%local_data(i, j)) > thresh) THEN
WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
fm%local_data(i, j)
END IF
END DO
END DO
END IF
CALL fm%matrix_struct%para_env%sync()
IF (unit_nr > 0) THEN
WRITE (unit_nr, *) my_footer
END IF
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief Write the atomic coordinates to the output unit.
!> \param particle_set ...
!> \note Adapted from particle_methods.F [MG]
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE write_qs_particle_coordinates_bse(particle_set, unit_nr)
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_particle_coordinates_bse'
CHARACTER(LEN=2) :: element_symbol
INTEGER :: handle, iatom, natom
CALL timeset(routineN, handle)
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T13,A7,16X,A7,15X,A7,15X,A7)') 'BSE|', &
'Element', 'x [Å]', 'y [Å]', 'z [Å]'
natom = SIZE(particle_set)
DO iatom = 1, natom
CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
element_symbol=element_symbol)
WRITE (unit_nr, '(T2,A4,T8,A12,1X,3(5X,F15.4))') &
'BSE|', element_symbol, particle_set(iatom)%r(1:3)*angstrom
END DO
END IF
CALL timestop(handle)
END SUBROUTINE write_qs_particle_coordinates_bse
END MODULE bse_print

981
src/bse_properties.F Normal file
View file

@ -0,0 +1,981 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Routines for computing excitonic properties, e.g. exciton diameter, from the BSE
!> \par History
!> 10.2024 created [Maximilian Graml]
! **************************************************************************************************
MODULE bse_properties
USE bse_util, ONLY: fm_general_add_bse,&
print_bse_nto_cubes,&
reshuffle_eigvec,&
trace_exciton_descr
USE cp_files, ONLY: close_file,&
open_file
USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
cp_fm_trace
USE cp_fm_diag, ONLY: cp_fm_svd
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_submatrix,&
cp_fm_release,&
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_to_fm_submat,&
cp_fm_to_fm_submat_general,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_unit_nr,&
cp_logger_type
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE mp2_types, ONLY: mp2_type
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_methods, ONLY: write_qs_particle_coordinates
USE particle_types, ONLY: particle_type
USE physcon, ONLY: evolt
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_mo_types, ONLY: allocate_mo_set,&
deallocate_mo_set,&
init_mo_set,&
mo_set_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_properties'
PUBLIC :: exciton_descr_type
PUBLIC :: get_exciton_descriptors, get_oscillator_strengths, compute_absorption_spectrum, &
calculate_NTOs
! TYPE definitions for exciton wavefunction descriptors
TYPE exciton_descr_type
REAL(KIND=dp), DIMENSION(3) :: r_e, r_h, &
r_e_sq, r_h_sq, r_e_h, &
r_e_shift, r_h_shift, &
cov_e_h
REAL(KIND=dp) :: sigma_e, sigma_h, &
cov_e_h_sum, corr_e_h, &
diff_r_abs, diff_r_sqr, &
norm_XpY
LOGICAL :: flag_TDA
END TYPE exciton_descr_type
CONTAINS
! **************************************************************************************************
!> \brief Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n )
!> and oscillator strengths f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2
!> Prelim Ref.: Eqs. (23), (24)
!> in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290
!> \param fm_eigvec_X ...
!> \param Exc_ens ...
!> \param fm_dipole_ai_trunc ...
!> \param trans_mom_bse BSE dipole vectors in real space per excitation level
!> \param oscill_str Oscillator strength per excitation level
!> \param polarizability_residues Residues of polarizability ("tensorial oscillator strength")
!> per excitation level
!> \param mp2_env ...
!> \param homo_red ...
!> \param virtual_red ...
!> \param unit_nr ...
!> \param fm_eigvec_Y ...
! **************************************************************************************************
SUBROUTINE get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
trans_mom_bse, oscill_str, polarizability_residues, &
mp2_env, homo_red, virtual_red, unit_nr, &
fm_eigvec_Y)
TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
INTENT(IN) :: Exc_ens
TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_ai_trunc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
INTENT(OUT) :: trans_mom_bse
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
INTENT(OUT) :: oscill_str
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
INTENT(OUT) :: polarizability_residues
TYPE(mp2_type), INTENT(IN) :: mp2_env
INTEGER, INTENT(IN) :: homo_red, virtual_red, unit_nr
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_oscillator_strengths'
INTEGER :: handle, idir, jdir, n
TYPE(cp_fm_struct_type), POINTER :: fm_struct_dipole_MO_trunc_reordered, &
fm_struct_trans_mom_bse
TYPE(cp_fm_type) :: fm_eigvec_XYsum
TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_MO_trunc_reordered, &
fm_dipole_per_dir, fm_trans_mom_bse
CALL timeset(routineN, handle)
CALL cp_fm_struct_create(fm_struct_dipole_MO_trunc_reordered, fm_eigvec_X%matrix_struct%para_env, &
fm_eigvec_X%matrix_struct%context, 1, homo_red*virtual_red)
CALL cp_fm_struct_create(fm_struct_trans_mom_bse, fm_eigvec_X%matrix_struct%para_env, &
fm_eigvec_X%matrix_struct%context, 1, homo_red*virtual_red)
! Include excitonic amplitudes in dipoles, i.e. obtain "BSE dipoles":
! \vec{D}_n = sqrt(2) * sum_{i,a} \vec{D}_ai (X_{ai}^{(n)} + Y_{ai}^{(n)})
! Reorder dipoles in order to execute the sum over i and a by parallel gemm
DO idir = 1, 3
CALL cp_fm_create(fm_dipole_MO_trunc_reordered(idir), matrix_struct=fm_struct_dipole_MO_trunc_reordered, &
name="dipoles_mo_reordered")
CALL cp_fm_set_all(fm_dipole_MO_trunc_reordered(idir), 0.0_dp)
CALL fm_general_add_bse(fm_dipole_MO_trunc_reordered(idir), fm_dipole_ai_trunc(idir), 1.0_dp, &
1, 1, &
1, virtual_red, &
unit_nr, (/2, 4, 3, 1/), mp2_env)
CALL cp_fm_release(fm_dipole_per_dir(idir))
END DO
DO idir = 1, 3
CALL cp_fm_create(fm_trans_mom_bse(idir), matrix_struct=fm_struct_trans_mom_bse, &
name="excitonic_dipoles")
CALL cp_fm_set_all(fm_trans_mom_bse(idir), 0.0_dp)
END DO
! If TDA is invoked, Y is not present as it is simply 0
CALL cp_fm_create(fm_eigvec_XYsum, matrix_struct=fm_eigvec_X%matrix_struct, name="excit_amplitude_sum")
CALL cp_fm_set_all(fm_eigvec_XYsum, 0.0_dp)
CALL cp_fm_to_fm(fm_eigvec_X, fm_eigvec_XYsum)
IF (PRESENT(fm_eigvec_Y)) THEN
CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_XYsum, 1.0_dp, fm_eigvec_Y)
END IF
DO idir = 1, 3
CALL parallel_gemm('N', 'N', 1, homo_red*virtual_red, homo_red*virtual_red, SQRT(2.0_dp), &
fm_dipole_MO_trunc_reordered(idir), fm_eigvec_XYsum, 0.0_dp, fm_trans_mom_bse(idir))
END DO
! Get oscillator strengths themselves
ALLOCATE (oscill_str(homo_red*virtual_red))
ALLOCATE (trans_mom_bse(3, 1, homo_red*virtual_red))
ALLOCATE (polarizability_residues(3, 3, homo_red*virtual_red))
trans_mom_bse(:, :, :) = 0.0_dp
! Sum over all directions
DO idir = 1, 3
CALL cp_fm_get_submatrix(fm_trans_mom_bse(idir), trans_mom_bse(idir, :, :))
END DO
DO n = 1, homo_red*virtual_red
DO idir = 1, 3
DO jdir = 1, 3
polarizability_residues(idir, jdir, n) = 2.0_dp*Exc_ens(n)*trans_mom_bse(idir, 1, n)*trans_mom_bse(jdir, 1, n)
END DO
END DO
oscill_str(n) = 2.0_dp/3.0_dp*Exc_ens(n)*SUM(ABS(trans_mom_bse(:, :, n))**2)
END DO
CALL cp_fm_struct_release(fm_struct_dipole_MO_trunc_reordered)
CALL cp_fm_struct_release(fm_struct_trans_mom_bse)
DO idir = 1, 3
CALL cp_fm_release(fm_dipole_MO_trunc_reordered(idir))
CALL cp_fm_release(fm_trans_mom_bse(idir))
CALL cp_fm_release(fm_dipole_ai_trunc(idir))
END DO
CALL cp_fm_release(fm_eigvec_XYsum)
CALL timestop(handle)
END SUBROUTINE get_oscillator_strengths
! **************************************************************************************************
!> \brief Computes and returns absorption spectrum for the frequency range and broadening
!> provided by the user.
!> Prelim Ref.: C. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications
!> (Oxford University Press, Oxford, 2012), Eq. 7.51
!> \param oscill_str ...
!> \param polarizability_residues ...
!> \param Exc_ens ...
!> \param info_approximation ...
!> \param unit_nr ...
!> \param mp2_env ...
! **************************************************************************************************
SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_ens, &
info_approximation, unit_nr, mp2_env)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
INTENT(IN) :: oscill_str
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
INTENT(IN) :: polarizability_residues
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
INTENT(IN) :: Exc_ens
CHARACTER(LEN=10) :: info_approximation
INTEGER, INTENT(IN) :: unit_nr
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_absorption_spectrum'
CHARACTER(LEN=10) :: eta_str, width_eta_format_str
CHARACTER(LEN=30) :: file_name_spectrum
INTEGER :: handle, i, idir, j, jdir, k, num_steps, &
unit_nr_file, width_eta
REAL(KIND=dp) :: eta, freq_end, freq_start, freq_step, &
omega
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: abs_spectrum
REAL(KIND=dp), DIMENSION(:), POINTER :: eta_list
TYPE(cp_logger_type), POINTER :: logger
CALL timeset(routineN, handle)
freq_step = mp2_env%bse%bse_spectrum_freq_step_size
freq_start = mp2_env%bse%bse_spectrum_freq_start
freq_end = mp2_env%bse%bse_spectrum_freq_end
eta_list => mp2_env%bse%bse_eta_spectrum_list
! Calculate number of steps to fit given frequency range
num_steps = NINT((freq_end - freq_start)/freq_step) + 1
DO k = 1, SIZE(eta_list)
eta = eta_list(k)
! Some magic to get a nice formatting of the eta value in filenames
width_eta = MAX(1, INT(LOG10(eta)) + 1) + 4
WRITE (width_eta_format_str, "(A2,I0,A3)") '(F', width_eta, '.3)'
WRITE (eta_str, width_eta_format_str) eta*evolt
! Filename itself
file_name_spectrum = 'BSE'//TRIM(ADJUSTL(info_approximation))//'eta='//TRIM(eta_str)//'.spectrum'
! First column is frequency in eV, second column is imaginary part of the trace of the polarizability
! The following 9 columns are the entries of the polarizability tensor
ALLOCATE (abs_spectrum(num_steps, 11))
abs_spectrum(:, :) = 0.0_dp
! Calculate the imaginary part of the mean dipole polarizability α_{avg}(ω)
! which is given by (cf. C. Ullrichs Book on TDDFT, Eq. 7.51)
! α_{avg}(ω) = \sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}
! and then the imaginary part is (in the limit η -> 0)
! Im[α_{avg}(ω)] = \sum_{n=1}^{N_exc} f_n * η / ((ω - Ω^n)² + η²)
! where f_n are the oscillator strengths and E_exc the excitation energies
! For the full polarizability tensor, we have
! α_{µ µ'}(ω) = \sum_n [2 Ω^n d^n_µ d^n_µ'] / [(ω+iη)^2- (Ω^n)^2]
DO i = 1, num_steps
omega = freq_start + (i - 1)*freq_step
abs_spectrum(i, 1) = omega
DO j = 1, SIZE(oscill_str)
abs_spectrum(i, 2) = abs_spectrum(i, 2) - oscill_str(j)* &
AIMAG(1/((omega + CMPLX(0.0, eta, kind=dp))**2 - Exc_ens(j)**2))
DO idir = 1, 3
DO jdir = 1, 3
abs_spectrum(i, 2 + (idir - 1)*3 + jdir) = abs_spectrum(i, 2 + (idir - 1)*3 + jdir) &
- polarizability_residues(idir, jdir, j)* &
AIMAG(1/((omega + CMPLX(0.0, eta, kind=dp))**2 - Exc_ens(j)**2))
END DO
END DO
END DO
END DO
! Print it to file
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
unit_nr_file = cp_logger_get_default_unit_nr()
ELSE
unit_nr_file = -1
END IF
IF (unit_nr_file > 0) THEN
CALL open_file(file_name_spectrum, unit_number=unit_nr_file, &
file_status="UNKNOWN", file_action="WRITE")
WRITE (unit_nr_file, '(A,A6)') "# Imaginary part of polarizability α_{µ µ'}(ω) = \sum_n "// &
"[2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] from Bethe Salpeter equation for method ", &
TRIM(ADJUSTL(info_approximation))
WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", "Im{α_{avg}(ω)}", "Im{α_xx(ω)}", &
"Im{α_xy(ω)}", "Im{α_xz(ω)}", "Im{α_yx(ω)}", "Im{α_yy(ω)}", "Im{α_yz(ω)}", "Im{α_zx(ω)}", &
"Im{α_zy(ω)}", "Im{α_zz(ω)}"
DO i = 1, num_steps
WRITE (unit_nr_file, '(11(F20.8,1X))') abs_spectrum(i, 1)*evolt, abs_spectrum(i, 2:11)
END DO
CALL close_file(unit_nr_file)
END IF
DEALLOCATE (abs_spectrum)
END DO
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A50,A)') &
'BSE|', "Printed optical absorption spectrum to local file ", file_name_spectrum
WRITE (unit_nr, '(T2,A4,T7,A50)') &
'BSE|', "using the Eq. 7.51 from C. Ullrichs Book on TDDFT:"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T10,A73)') &
'BSE|', "Im{α_{avg}(ω)} = Im{\sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}}"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') &
'BSE|', "or for the full polarizability tensor:"
WRITE (unit_nr, '(T2,A4,T10,A)') &
'BSE|', "α_{µ µ'}(ω) = \sum_n [2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²]"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') &
'BSE|', "with transition moments d_µ^n, oscillator strengths f_n,"
WRITE (unit_nr, '(T2,A4,T7,A)') &
'BSE|', "and excitation energies Ω^n."
WRITE (unit_nr, '(T2,A4)') 'BSE|'
END IF
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param fm_X ...
!> \param fm_Y ...
!> \param mo_coeff ...
!> \param homo ...
!> \param virtual ...
!> \param info_approximation ...
!> \param oscill_str ...
!> \param qs_env ...
!> \param unit_nr ...
!> \param mp2_env ...
! **************************************************************************************************
SUBROUTINE calculate_NTOs(fm_X, fm_Y, &
mo_coeff, homo, virtual, &
info_approximation, &
oscill_str, &
qs_env, unit_nr, mp2_env)
TYPE(cp_fm_type), INTENT(IN) :: fm_X
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_Y
TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
INTEGER, INTENT(IN) :: homo, virtual
CHARACTER(LEN=10) :: info_approximation
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, INTENT(IN) :: unit_nr
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_NTOs'
CHARACTER(LEN=20), DIMENSION(2) :: nto_name
INTEGER :: handle, homo_irred, i, i_nto, info_svd, &
j, n_exc, n_nto, nao_full, nao_trunc
INTEGER, DIMENSION(:), POINTER :: stride
LOGICAL :: append_cube, cube_file
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval_svd_squ
REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_svd
TYPE(cp_fm_struct_type), POINTER :: fm_struct_m, fm_struct_mo_coeff, &
fm_struct_nto_holes, &
fm_struct_nto_particles, &
fm_struct_nto_set
TYPE(cp_fm_type) :: fm_eigvl, fm_eigvr_t, fm_m, fm_mo_coeff, fm_nto_coeff_holes, &
fm_nto_coeff_particles, fm_nto_set, fm_X_ia, fm_Y_ai
TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
TYPE(section_vals_type), POINTER :: bse_section, input, nto_section
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, &
input=input)
bse_section => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%GW%BSE")
nao_full = qs_env%mos(1)%nao
nao_trunc = homo + virtual
! This is not influenced by the BSE cutoff
homo_irred = qs_env%mos(1)%homo
! M will have a block structure and is quadratic in homo+virtual, i.e.
! occ virt
! | 0 X_i,a | occ = homo
! M = | Y_a,i 0 | virt = virtual
!
! X and Y are here not the eigenvectors X_ia,n - instead we fix n and reshape the combined ia index
! Notice the index structure of the lower block, i.e. X is transposed
CALL cp_fm_struct_create(fm_struct_m, &
fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
nao_trunc, nao_trunc)
CALL cp_fm_struct_create(fm_struct_mo_coeff, &
fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
nao_full, nao_trunc)
CALL cp_fm_struct_create(fm_struct_nto_holes, &
fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
nao_full, nao_trunc)
CALL cp_fm_struct_create(fm_struct_nto_particles, &
fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
nao_full, nao_trunc)
CALL cp_fm_create(fm_mo_coeff, matrix_struct=fm_struct_mo_coeff, &
name="mo_coeff")
! Here, we take care of possible cutoffs
! Simply truncating the matrix causes problems with the print function
! Therefore, we keep the dimension, but set the coefficients of truncated indices to 0
CALL cp_fm_to_fm_submat_general(mo_coeff(1), fm_mo_coeff, &
nao_full, nao_trunc, &
1, homo_irred - homo + 1, &
1, 1, &
mo_coeff(1)%matrix_struct%context)
! Print some information about the NTOs
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'The Natural Transition Orbital (NTO) pairs φ_I(r_e) and χ_I(r_h) for a fixed'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'excitation index n are obtained by singular value decomposition of T'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
' = (0 X)'
IF (PRESENT(fm_Y)) THEN
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'T = (Y^T 0)'
ELSE
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'T = (0 0)'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'T = U Λ V^T'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'φ_I(r_e) = \sum_p V_pI ψ_p(r_e)'
WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
'χ_I(r_h) = \sum_p U_pI ψ_p(r_e)'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'where we have introduced'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
'BSE|', "ψ_p:", "occupied and virtual molecular orbitals,"
WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
'BSE|', "φ_I(r_e):", "NTO state for the electron,"
WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
'BSE|', "χ_I(r_h):", "NTO state for the hole,"
WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
'BSE|', "Λ:", "diagonal matrix of NTO weights λ_I,"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
"The NTOs are calculated with the following settings:"
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Number of excitations, for which NTOs are computed', &
mp2_env%bse%num_print_exc_ntos
IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp) THEN
WRITE (unit_nr, '(T2,A4,T7,A,T71,F10.3)') 'BSE|', 'Threshold for oscillator strength f^n', &
mp2_env%bse%eps_nto_osc_str
ELSE
WRITE (unit_nr, '(T2,A4,T7,A,T71,A10)') 'BSE|', 'Threshold for oscillator strength f^n', &
ADJUSTL("---")
END IF
WRITE (unit_nr, '(T2,A4,T7,A,T72,F10.3)') 'BSE|', 'Threshold for NTO weights (λ_I)^2', &
mp2_env%bse%eps_nto_eigval
END IF
! Write the header of NTO info table
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
IF (.NOT. PRESENT(fm_Y)) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'NTOs from solving the BSE within the TDA:'
ELSE
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
'NTOs from solving the BSE without the TDA:'
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T33,A14,T62,A)') 'BSE|', &
'Excitation n', "TDA/ABBA", "Index of NTO I", 'NTO weights (λ_I)^2'
END IF
DO j = 1, mp2_env%bse%num_print_exc_ntos
n_exc = mp2_env%bse%bse_nto_state_list_final(j)
! Takes care of unallocated oscill_str array in case of Triplet
IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp) THEN
! Check actual values
IF (oscill_str(n_exc) < mp2_env%bse%eps_nto_osc_str) THEN
! Print skipped levels to table
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T42,A39)') 'BSE|', &
n_exc, info_approximation, "Skipped (Oscillator strength too small)"
END IF
CYCLE
END IF
END IF
CALL cp_fm_create(fm_m, matrix_struct=fm_struct_m, &
name="single_part_trans_dm")
CALL cp_fm_set_all(fm_m, 0.0_dp)
CALL cp_fm_create(fm_nto_coeff_holes, matrix_struct=fm_struct_nto_holes, &
name="nto_coeffs_holes")
CALL cp_fm_set_all(fm_nto_coeff_holes, 0.0_dp)
CALL cp_fm_create(fm_nto_coeff_particles, matrix_struct=fm_struct_nto_particles, &
name="nto_coeffs_particles")
CALL cp_fm_set_all(fm_nto_coeff_particles, 0.0_dp)
! Reshuffle from X_ia,n_exc to X_i,a
CALL reshuffle_eigvec(fm_X, fm_X_ia, homo, virtual, n_exc, &
.FALSE., unit_nr, mp2_env)
! Copy X to upper block in M, i.e. starting from column homo+1
CALL cp_fm_to_fm_submat(fm_X_ia, fm_m, &
homo, virtual, &
1, 1, &
1, homo + 1)
CALL cp_fm_release(fm_X_ia)
! Copy Y if present
IF (PRESENT(fm_Y)) THEN
! Reshuffle from Y_ia,n_exc to Y_a,i
CALL reshuffle_eigvec(fm_Y, fm_Y_ai, homo, virtual, n_exc, &
.TRUE., unit_nr, mp2_env)
! Copy Y^T to lower block in M, i.e. starting from row homo+1
CALL cp_fm_to_fm_submat(fm_Y_ai, fm_m, &
virtual, homo, &
1, 1, &
homo + 1, 1)
CALL cp_fm_release(fm_Y_ai)
END IF
! Now we compute the SVD of M_{occ+virt,occ+virt}, which yields
! M = U * Lambda * V^T
! Initialize matrices and arrays to store left/right eigenvectors and singular values
CALL cp_fm_create(matrix=fm_eigvl, &
matrix_struct=fm_m%matrix_struct, &
name="LEFT_SINGULAR_MATRIX")
CALL cp_fm_set_all(fm_eigvl, alpha=0.0_dp)
CALL cp_fm_create(matrix=fm_eigvr_t, &
matrix_struct=fm_m%matrix_struct, &
name="RIGHT_SINGULAR_MATRIX")
CALL cp_fm_set_all(fm_eigvr_t, alpha=0.0_dp)
ALLOCATE (eigval_svd(nao_trunc))
eigval_svd(:) = 0.0_dp
info_svd = 0
CALL cp_fm_svd(fm_m, fm_eigvl, fm_eigvr_t, eigval_svd, info_svd)
IF (info_svd /= 0) THEN
IF (unit_nr > 0) THEN
CALL cp_warn(__LOCATION__, &
"SVD for computation of NTOs not successful. "// &
"Skipping print of NTOs.")
IF (info_svd > 0) THEN
CALL cp_warn(__LOCATION__, &
"PDGESVD detected heterogeneity. "// &
"Decreasing number of MPI ranks might solve this issue.")
END IF
END IF
! Release matrices to avoid memory leaks
CALL cp_fm_release(fm_m)
CALL cp_fm_release(fm_nto_coeff_holes)
CALL cp_fm_release(fm_nto_coeff_particles)
ELSE
! Rescale singular values as done in Martin2003 (10.1063/1.1558471)
ALLOCATE (eigval_svd_squ(nao_trunc))
eigval_svd_squ(:) = eigval_svd(:)**2
! Sanity check for TDA: In case of TDA, the sum should be \sum_ia |X_ia|^2 = 1
IF (.NOT. PRESENT(fm_Y)) THEN
IF (ABS(SUM(eigval_svd_squ) - 1) >= 1e-5) THEN
CPWARN("Sum of NTO coefficients deviates from 1!")
END IF
END IF
! Create NTO coefficients for later print to grid via TDDFT routine
! Apply U = fm_eigvl to MO coeffs, which yields hole states
CALL parallel_gemm("N", "N", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvl, 0.0_dp, &
fm_nto_coeff_holes)
! Apply V^T = fm_eigvr_t to MO coeffs, which yields particle states
CALL parallel_gemm("N", "T", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvr_t, 0.0_dp, &
fm_nto_coeff_particles)
!Release intermediary work matrices
CALL cp_fm_release(fm_m)
CALL cp_fm_release(fm_eigvl)
CALL cp_fm_release(fm_eigvr_t)
! Transfer NTO coefficients to sets
nto_name(1) = 'Hole_coord'
nto_name(2) = 'Particle_coord'
ALLOCATE (nto_set(2))
! Extract number of significant NTOs
n_nto = 0
DO i_nto = 1, nao_trunc
IF (eigval_svd_squ(i_nto) > mp2_env%bse%eps_nto_eigval) THEN
n_nto = n_nto + 1
ELSE
! Since svd orders in descending order, we can exit the loop if smaller
EXIT
END IF
END DO
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A4)') 'BSE|'
DO i_nto = 1, n_nto
WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T41,I6,T71,F10.5)') 'BSE|', &
n_exc, info_approximation, i_nto, eigval_svd_squ(i_nto)
END DO
END IF
CALL cp_fm_struct_create(fm_struct_nto_set, template_fmstruct=fm_struct_nto_holes, &
ncol_global=n_nto)
CALL cp_fm_create(fm_nto_set, fm_struct_nto_set)
DO i = 1, 2
CALL allocate_mo_set(nto_set(i), nao_trunc, n_nto, 0, 0.0_dp, 2.0_dp, 0.0_dp)
CALL init_mo_set(nto_set(i), fm_ref=fm_nto_set, name=nto_name(i))
END DO
CALL cp_fm_release(fm_nto_set)
CALL cp_fm_struct_release(fm_struct_nto_set)
! Fill NTO sets
CALL cp_fm_to_fm(fm_nto_coeff_holes, nto_set(1)%mo_coeff, ncol=n_nto)
CALL cp_fm_to_fm(fm_nto_coeff_particles, nto_set(2)%mo_coeff, ncol=n_nto)
! Cube files
nto_section => section_vals_get_subs_vals(bse_section, "NTO_ANALYSIS")
CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
IF (cube_file) THEN
CALL print_bse_nto_cubes(qs_env, nto_set, n_exc, info_approximation, &
stride, append_cube, nto_section)
END IF
CALL cp_fm_release(fm_nto_coeff_holes)
CALL cp_fm_release(fm_nto_coeff_particles)
DEALLOCATE (eigval_svd)
DEALLOCATE (eigval_svd_squ)
DO i = 1, 2
CALL deallocate_mo_set(nto_set(i))
END DO
DEALLOCATE (nto_set)
END IF
END DO
CALL cp_fm_release(fm_mo_coeff)
CALL cp_fm_struct_release(fm_struct_m)
CALL cp_fm_struct_release(fm_struct_nto_holes)
CALL cp_fm_struct_release(fm_struct_nto_particles)
CALL cp_fm_struct_release(fm_struct_mo_coeff)
CALL timestop(handle)
END SUBROUTINE calculate_NTOs
! **************************************************************************************************
!> \brief ...
!> \param exc_descr ...
!> \param fm_eigvec_X ...
!> \param fm_dipole_ij_trunc ...
!> \param fm_dipole_ab_trunc ...
!> \param fm_dipole_ai_trunc ...
!> \param fm_quadpole_ij_trunc ...
!> \param fm_quadpole_ab_trunc ...
!> \param homo ...
!> \param virtual ...
!> \param unit_nr ...
!> \param mp2_env ...
!> \param qs_env ...
!> \param fm_eigvec_Y ...
! **************************************************************************************************
SUBROUTINE get_exciton_descriptors(exc_descr, fm_eigvec_X, &
fm_dipole_ij_trunc, fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
homo, virtual, unit_nr, &
mp2_env, qs_env, &
fm_eigvec_Y)
TYPE(exciton_descr_type), ALLOCATABLE, &
DIMENSION(:) :: exc_descr
TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
INTENT(IN) :: fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
fm_dipole_ai_trunc, &
fm_quadpole_ij_trunc, &
fm_quadpole_ab_trunc
INTEGER, INTENT(IN) :: homo, virtual, unit_nr
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_exciton_descriptors'
INTEGER :: handle, i_dir, i_exc, n_exc
INTEGER, DIMENSION(3) :: mask_quadrupole
LOGICAL :: flag_TDA
REAL(KIND=dp) :: norm_X, norm_XpY, norm_Y
REAL(KIND=dp), DIMENSION(3) :: r_e_h_XX, r_e_h_XY, r_e_h_YX, r_e_h_YY, &
r_e_sq_X, r_e_sq_Y, r_e_X, r_e_Y, &
r_h_sq_X, r_h_sq_Y, r_h_X, r_h_Y
TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia
TYPE(cp_fm_type) :: fm_work_ba, fm_work_ia, fm_work_ia_2, &
fm_X_ia, fm_Y_ia
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(section_vals_type), POINTER :: input, subsys_section
CALL timeset(routineN, handle)
n_exc = homo*virtual
IF (PRESENT(fm_eigvec_Y)) THEN
flag_TDA = .FALSE.
ELSE
flag_TDA = .TRUE.
END IF
! translates 1,2,3 to diagonal entries of quadrupoles xx, yy, zz
! Ordering in quadrupole moments is x, y, z, xx, xy, xz, yy, yz, zz
mask_quadrupole = (/4, 7, 9/)
ALLOCATE (exc_descr(mp2_env%bse%num_print_exc_descr))
CALL cp_fm_struct_create(fm_struct_ia, &
fm_eigvec_X%matrix_struct%para_env, fm_eigvec_X%matrix_struct%context, &
homo, virtual)
CALL cp_fm_struct_create(fm_struct_ab, &
fm_eigvec_X%matrix_struct%para_env, fm_eigvec_X%matrix_struct%context, &
virtual, virtual)
! print actual coords (might be centered and differing from input xyz) for debug runs
IF (mp2_env%bse%bse_debug_print) THEN
CALL get_qs_env(qs_env, particle_set=particle_set, &
qs_kind_set=qs_kind_set, input=input)
subsys_section => section_vals_get_subs_vals(input, "SUBSYS")
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A10,T13,A)') 'BSE|DEBUG|', &
'Printing internal geometry for debug purposes'
END IF
CALL write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label="BSE")
END IF
DO i_exc = 1, mp2_env%bse%num_print_exc_descr
r_e_X(:) = 0.0_dp
r_e_Y(:) = 0.0_dp
r_h_X(:) = 0.0_dp
r_h_Y(:) = 0.0_dp
r_e_sq_X(:) = 0.0_dp
r_h_sq_X(:) = 0.0_dp
r_e_sq_Y(:) = 0.0_dp
r_h_sq_Y(:) = 0.0_dp
r_e_h_XX(:) = 0.0_dp
r_e_h_XY(:) = 0.0_dp
r_e_h_YX(:) = 0.0_dp
r_e_h_YY(:) = 0.0_dp
norm_X = 0.0_dp
norm_Y = 0.0_dp
norm_XpY = 0.0_dp
! Initialize values of exciton descriptors
exc_descr(i_exc)%r_e(:) = 0.0_dp
exc_descr(i_exc)%r_h(:) = 0.0_dp
exc_descr(i_exc)%r_e_sq(:) = 0.0_dp
exc_descr(i_exc)%r_h_sq(:) = 0.0_dp
exc_descr(i_exc)%r_e_h(:) = 0.0_dp
exc_descr(i_exc)%flag_TDA = flag_TDA
exc_descr(i_exc)%norm_XpY = 0.0_dp
CALL reshuffle_eigvec(fm_eigvec_X, fm_X_ia, homo, virtual, i_exc, &
.FALSE., unit_nr, mp2_env)
! Norm of X
CALL cp_fm_trace(fm_X_ia, fm_X_ia, norm_X)
norm_XpY = norm_X
IF (.NOT. flag_TDA) THEN
CALL reshuffle_eigvec(fm_eigvec_Y, fm_Y_ia, homo, virtual, i_exc, &
.FALSE., unit_nr, mp2_env)
! Norm of Y
CALL cp_fm_trace(fm_Y_ia, fm_Y_ia, norm_Y)
norm_XpY = norm_XpY + norm_Y
END IF
exc_descr(i_exc)%norm_XpY = norm_XpY
! <r_h>_X = Tr[ X^T µ_ij X + Y µ_ab Y^T ] = X_ai µ_ij X_ja + Y_ia µ_ab Y_bi
DO i_dir = 1, 3
! <r_h>_X = X_ai µ_ij X_ja + ...
CALL trace_exciton_descr(fm_X_ia, fm_dipole_ij_trunc(i_dir), fm_X_ia, r_h_X(i_dir))
r_h_X(i_dir) = r_h_X(i_dir)/norm_XpY
IF (.NOT. flag_TDA) THEN
! <r_h>_X = ... + Y_ia µ_ab Y_bi
CALL trace_exciton_descr(fm_Y_ia, fm_Y_ia, fm_dipole_ab_trunc(i_dir), r_h_Y(i_dir))
r_h_Y(i_dir) = r_h_Y(i_dir)/norm_XpY
END IF
END DO
exc_descr(i_exc)%r_h(:) = r_h_X(:) + r_h_Y(:)
! <r_e>_X = Tr[ X µ_ab X^T + Y^T µ_ij Y ] = X_ia µ_ab X_bi + Y_ai µ_ij Y_ja
DO i_dir = 1, 3
! <r_e>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
CALL trace_exciton_descr(fm_X_ia, fm_X_ia, fm_dipole_ab_trunc(i_dir), r_e_X(i_dir))
r_e_X(i_dir) = r_e_X(i_dir)/norm_XpY
IF (.NOT. flag_TDA) THEN
! <r_e>_X = ... + Y_ai µ_ij Y_ja
CALL trace_exciton_descr(fm_Y_ia, fm_dipole_ij_trunc(i_dir), fm_Y_ia, r_e_Y(i_dir))
r_e_Y(i_dir) = r_e_Y(i_dir)/norm_XpY
END IF
END DO
exc_descr(i_exc)%r_e(:) = r_e_X(:) + r_e_Y(:)
! <r_h^2>_X = Tr[ X^T M_ij X + Y M_ab Y^T ] = X_ai M_ij X_ja + Y_ia M_ab Y_bi
DO i_dir = 1, 3
! <r_h^2>_X = X_ai M_ij X_ja + ...
CALL trace_exciton_descr(fm_X_ia, fm_quadpole_ij_trunc(mask_quadrupole(i_dir)), &
fm_X_ia, r_h_sq_X(i_dir))
r_h_sq_X(i_dir) = r_h_sq_X(i_dir)/norm_XpY
IF (.NOT. flag_TDA) THEN
! <r_h^2>_X = ... + Y_ia M_ab Y_bi
CALL trace_exciton_descr(fm_Y_ia, fm_Y_ia, &
fm_quadpole_ab_trunc(mask_quadrupole(i_dir)), r_h_sq_Y(i_dir))
r_h_sq_Y(i_dir) = r_h_sq_Y(i_dir)/norm_XpY
END IF
END DO
exc_descr(i_exc)%r_h_sq(:) = r_h_sq_X(:) + r_h_sq_Y(:)
! <r_e^2>_X = Tr[ X M_ab X^T + Y^T M_ij Y ] = X_ia M_ab X_bi + Y_ai M_ij Y_ja
DO i_dir = 1, 3
! <r_e^2>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
CALL trace_exciton_descr(fm_X_ia, fm_X_ia, &
fm_quadpole_ab_trunc(mask_quadrupole(i_dir)), r_e_sq_X(i_dir))
r_e_sq_X(i_dir) = r_e_sq_X(i_dir)/norm_XpY
IF (.NOT. flag_TDA) THEN
! <r_e^2>_X = ... + Y_ai M_ij Y_ja
CALL trace_exciton_descr(fm_Y_ia, fm_quadpole_ij_trunc(mask_quadrupole(i_dir)), &
fm_Y_ia, r_e_sq_Y(i_dir))
r_e_sq_Y(i_dir) = r_e_sq_Y(i_dir)/norm_XpY
END IF
END DO
exc_descr(i_exc)%r_e_sq(:) = r_e_sq_X(:) + r_e_sq_Y(:)
! <r_h r_e>_X
! = Tr[ X^T µ_ij X µ_ab + Y^T µ_ij Y µ_ab + X µ_ai Y µ_ai + Y µ_ai X µ_ai]
! = X_bj µ_ji X_ia µ_ab + Y_bj µ_ji Y_ia µ_ab + X_ia µ_aj Y_jb µ_bi + Y_ia µ_aj X_jb µ_bi
CALL cp_fm_create(fm_work_ia, fm_struct_ia)
CALL cp_fm_create(fm_work_ia_2, fm_struct_ia)
CALL cp_fm_create(fm_work_ba, fm_struct_ab)
DO i_dir = 1, 3
! First term - X^T µ_ij X µ_ab
CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
CALL cp_fm_set_all(fm_work_ia_2, 0.0_dp)
! work_ib = X_ia µ_ab
CALL parallel_gemm("N", "N", homo, virtual, virtual, 1.0_dp, &
fm_X_ia, fm_dipole_ab_trunc(i_dir), 0.0_dp, fm_work_ia)
! work_ja_2 = µ_ji work_ia
CALL parallel_gemm("N", "N", homo, virtual, homo, 1.0_dp, &
fm_dipole_ij_trunc(i_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
! <r_h r_e>_X = work_ia_2 X_bj + ... = X^T work_ia_2 + ...
CALL cp_fm_trace(fm_X_ia, fm_work_ia_2, r_e_h_XX(i_dir))
r_e_h_XX(i_dir) = r_e_h_XX(i_dir)/norm_XpY
IF (.NOT. flag_TDA) THEN
! Second term - Y^T µ_ij Y µ_ab
CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
CALL cp_fm_set_all(fm_work_ia_2, 0.0_dp)
! work_ib = Y_ia µ_ab
CALL parallel_gemm("N", "N", homo, virtual, virtual, 1.0_dp, &
fm_Y_ia, fm_dipole_ab_trunc(i_dir), 0.0_dp, fm_work_ia)
! work_ja_2 = µ_ji work_ia
CALL parallel_gemm("N", "N", homo, virtual, homo, 1.0_dp, &
fm_dipole_ij_trunc(i_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
! <r_h r_e>_X = work_ia_2 Y_bj + ... = Y^T work_ia_2 + ...
CALL cp_fm_trace(fm_Y_ia, fm_work_ia_2, r_e_h_YY(i_dir))
r_e_h_YY(i_dir) = r_e_h_YY(i_dir)/norm_XpY
! Third term - X µ_ai Y µ_ai = X_ia µ_aj Y_jb µ_bi
! Reshuffle for usage of trace (where first argument is transposed)
! = µ_aj Y_jb µ_bi X_ia =
! \__________/
! fm_work_ai
! fm_work_ai = µ_aj Y_jb µ_bi
! fm_work_ia = µ_ib Y_bj µ_ja
! \_____/
! fm_work_ba
CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
! fm_work_ba = Y_bj µ_ja
CALL parallel_gemm("T", "T", virtual, virtual, homo, 1.0_dp, &
fm_Y_ia, fm_dipole_ai_trunc(i_dir), 0.0_dp, fm_work_ba)
! fm_work_ia = µ_ib fm_work_ba
CALL parallel_gemm("T", "N", homo, virtual, virtual, 1.0_dp, &
fm_dipole_ai_trunc(i_dir), fm_work_ba, 0.0_dp, fm_work_ia)
! <r_h r_e>_X = ... + X_ia µ_aj Y_jb µ_bi
CALL cp_fm_trace(fm_work_ia, fm_X_ia, r_e_h_XY(i_dir))
r_e_h_XY(i_dir) = r_e_h_XY(i_dir)/norm_XpY
! Fourth term - Y µ_ai X µ_ai = Y_ia µ_aj X_jb µ_bi
! Reshuffle for usage of trace (where first argument is transposed)
! = µ_aj X_jb µ_bi Y_ia =
! \__________/
! fm_work_ai
! fm_work_ai = µ_aj X_jb µ_bi
! fm_work_ia = µ_ib X_bj µ_ja
! \_____/
! fm_work_ba
CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
! fm_work_ba = Y_bj µ_ja
CALL parallel_gemm("T", "T", virtual, virtual, homo, 1.0_dp, &
fm_X_ia, fm_dipole_ai_trunc(i_dir), 0.0_dp, fm_work_ba)
! fm_work_ia = µ_ib fm_work_ba
CALL parallel_gemm("T", "N", homo, virtual, virtual, 1.0_dp, &
fm_dipole_ai_trunc(i_dir), fm_work_ba, 0.0_dp, fm_work_ia)
! <r_h r_e>_X = ... + X_ia µ_aj Y_jb µ_bi
CALL cp_fm_trace(fm_work_ia, fm_Y_ia, r_e_h_YX(i_dir))
r_e_h_YX(i_dir) = r_e_h_YX(i_dir)/norm_XpY
END IF
END DO
exc_descr(i_exc)%r_e_h(:) = r_e_h_XX(:) + r_e_h_XY(:) + r_e_h_YX(:) + r_e_h_YY(:)
CALL cp_fm_release(fm_work_ia)
CALL cp_fm_release(fm_work_ia_2)
CALL cp_fm_release(fm_work_ba)
! diff_r_abs = |<r_h>_X - <r_e>_X|
exc_descr(i_exc)%diff_r_abs = SQRT(SUM((exc_descr(i_exc)%r_h(:) - exc_descr(i_exc)%r_e(:))**2))
! σ_e = sqrt( <r_e^2>_X - <r_e>_X^2 )
exc_descr(i_exc)%sigma_e = SQRT(SUM(exc_descr(i_exc)%r_e_sq(:)) - SUM(exc_descr(i_exc)%r_e(:)**2))
! σ_h = sqrt( <r_h^2>_X - <r_h>_X^2 )
exc_descr(i_exc)%sigma_h = SQRT(SUM(exc_descr(i_exc)%r_h_sq(:)) - SUM(exc_descr(i_exc)%r_h(:)**2))
! COV(r_e, r_h) = < r_e r_h >_X - < r_e >_X < r_h >_X
exc_descr(i_exc)%cov_e_h_sum = 0.0_dp
exc_descr(i_exc)%cov_e_h(:) = 0.0_dp
DO i_dir = 1, 3
exc_descr(i_exc)%cov_e_h(i_dir) = exc_descr(i_exc)%r_e_h(i_dir) &
- exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(i_dir)
exc_descr(i_exc)%cov_e_h_sum = exc_descr(i_exc)%cov_e_h_sum + &
exc_descr(i_exc)%r_e_h(i_dir) - exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(i_dir)
END DO
! root-mean-square e-h separation
exc_descr(i_exc)%diff_r_sqr = SQRT(exc_descr(i_exc)%diff_r_abs**2 + &
exc_descr(i_exc)%sigma_e**2 + exc_descr(i_exc)%sigma_h**2 &
- 2*exc_descr(i_exc)%cov_e_h_sum)
! e-h-correlation coefficient R_eh = COV(r_e, r_h) / ( σ_e σ_h )
exc_descr(i_exc)%corr_e_h = exc_descr(i_exc)%cov_e_h_sum/(exc_descr(i_exc)%sigma_e*exc_descr(i_exc)%sigma_h)
! Expectation values of r_e and r_h
exc_descr(i_exc)%r_e_shift(:) = exc_descr(i_exc)%r_e(:)
exc_descr(i_exc)%r_h_shift(:) = exc_descr(i_exc)%r_h(:)
CALL cp_fm_release(fm_X_ia)
IF (.NOT. flag_TDA) THEN
CALL cp_fm_release(fm_Y_ia)
END IF
END DO
CALL cp_fm_struct_release(fm_struct_ia)
CALL cp_fm_struct_release(fm_struct_ab)
CALL timestop(handle)
END SUBROUTINE get_exciton_descriptors
END MODULE bse_properties

986
src/bse_util.F
View file

File diff suppressed because it is too large Load diff

23
src/common/bibliography.F
View file

@ -88,7 +88,8 @@ MODULE bibliography
Caldeweyher2017, Caldeweyher2019, Caldeweyher2020, Freeman1977, Gruneis2009, &
Stein2022, Stein2024, Pracht2019, &
Blase2018, Blase2020, Bruneval2015, Golze2019, Gui2018, Jacquemin2017, Liu2020, &
Sander2015, Schreiber2008, vanSetten2015, Setyawan2010, Ahart2024
Sander2015, Schreiber2008, vanSetten2015, Setyawan2010, Ahart2024, Knysh2024, &
Schambeck2024, Mewes2018
CONTAINS
@ -1884,12 +1885,32 @@ CONTAINS
source="J. Chem. Theory Comput.", volume="20", pages="6772-6780", &
year=2024, doi="10.1021/acs.jctc.4c00371")
CALL add_reference(key=Knysh2024, &
authors=s2a("I. Knysh", "F. Lipparini", "A. Blondel", "I. Duchemin", "X. Blase", &
"P.-F. Loos", "D. Jacquemin"), &
title="Reference CC3 Excitation Energies for Organic Chromophores: "// &
"Benchmarking TD-DFT, BSE/GW, and Wave Function Methods", &
source="J. Chem. Theory Comput.", volume="20", pages="8152-8174", &
year=2024, doi="10.1021/acs.jctc.4c00906")
CALL add_reference(key=Schambeck2024, &
authors=s2a("M. Schambeck", "D. Golze", "J. Wilhelm"), &
title="Solving multipole challenges in the GW 100 benchmark enables precise low-scaling GW calculations", &
source="Phys. Rev. B", volume="110", pages="125146", &
year=2024, doi="10.1103/PhysRevB.110.125146")
CALL add_reference(key=Pracht2019, &
authors=s2a("P. Pracht", "E. Caldeweyher", "S. Ehlert", "S. Grimme"), &
title="A Robust Non-Self-Consistent Tight-Binding Quantum Chemistry Method for large Molecules", &
source="ChemRxiv", volume="", pages="", &
year=2019, doi="10.26434/chemrxiv.8326202.v1")
CALL add_reference(key=Mewes2018, &
authors=s2a("S. Mewes", "F. Plasser", "A. Krylov", "A. Dreuw"), &
title="Benchmarking Excited-State Calculations Using Exciton Properties", &
source="J. Chem. Theory Comput.", volume="14", pages="710-725", &
year=2018, doi="10.1021/acs.jctc.7b01145 ")
END SUBROUTINE add_all_references
END MODULE bibliography

138
src/fm/cp_fm_diag.F
View file

@ -41,7 +41,23 @@ MODULE cp_fm_diag
USE cp_fm_types, ONLY: cp_fm_get_info, &
cp_fm_set_element, &
cp_fm_to_fm, &
cp_fm_type, &
cp_fm_create, &
cp_fm_get_info, &
cp_fm_release, &
cp_fm_set_all, &
cp_fm_to_fm, &
cp_fm_to_fm_submat, &
cp_fm_type
USE cp_fm_struct, ONLY: cp_fm_struct_equivalent, &
cp_fm_struct_create, &
cp_fm_struct_release, &
cp_fm_struct_type
USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr, &
cp_get_default_logger, &
cp_logger_get_default_io_unit, &
cp_logger_type, &
cp_to_string
USE cp_log_handling, ONLY: cp_get_default_logger, &
cp_logger_get_default_unit_nr, &
cp_logger_get_unit_nr, &
@ -104,6 +120,7 @@ MODULE cp_fm_diag
cp_fm_power, &
cp_fm_syevd, &
cp_fm_syevx, &
cp_fm_svd, &
cp_fm_geeig, &
cp_fm_geeig_canon, &
diag_init, &
@ -882,6 +899,127 @@ CONTAINS
END SUBROUTINE cp_fm_syevx
! **************************************************************************************************
!> \brief decomposes a quadratic matrix into its singular value decomposition
!> \param matrix_a ...
!> \param matrix_eigvl ...
!> \param matrix_eigvr_t ...
!> \param eigval ...
!> \param info ...
!> \author Maximilian Graml
! **************************************************************************************************
SUBROUTINE cp_fm_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
TYPE(cp_fm_type), INTENT(IN) :: matrix_a
TYPE(cp_fm_type), INTENT(INOUT) :: matrix_eigvl, matrix_eigvr_t
REAL(KIND=dp), DIMENSION(:), POINTER, &
INTENT(INOUT) :: eigval
INTEGER, INTENT(OUT), OPTIONAL :: info
CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_fm_svd'
CHARACTER(LEN=2*default_string_length) :: message
INTEGER :: handle, n, m, myinfo, lwork
REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
TYPE(cp_fm_type) :: matrix_lu
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
#if defined(__parallel)
INTEGER, DIMENSION(9) :: desca, descu, descvt
#endif
CALL timeset(routineN, handle)
CALL cp_fm_create(matrix=matrix_lu, &
matrix_struct=matrix_a%matrix_struct, &
name="A_lu"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
CALL cp_fm_to_fm(matrix_a, matrix_lu)
a => matrix_lu%local_data
m = matrix_lu%matrix_struct%nrow_global
n = matrix_lu%matrix_struct%ncol_global
! Assert that incoming matrix is quadratic
CPASSERT(m == n)
! Prepare for workspace queries
myinfo = 0
lwork = -1
ALLOCATE (work(1))
work(:) = 0.0_dp
#if defined(__parallel)
! To do: That might need a redistribution check as in cp_fm_syevd
desca(:) = matrix_lu%matrix_struct%descriptor(:)
descu(:) = matrix_eigvl%matrix_struct%descriptor(:)
descvt(:) = matrix_eigvr_t%matrix_struct%descriptor(:)
! Workspace query
CALL pdgesvd('V', 'V', m, m, matrix_lu%local_data, 1, 1, desca, eigval, matrix_eigvl%local_data, &
1, 1, descu, matrix_eigvr_t%local_data, 1, 1, descvt, work, lwork, myinfo)
IF (matrix_lu%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
WRITE (message, "(A,I0,A)") "ERROR in PDGESVD: Work space query failed (INFO = ", myinfo, ")"
IF (PRESENT(info)) THEN
CPWARN(TRIM(message))
ELSE
CPABORT(TRIM(message))
END IF
END IF
lwork = INT(work(1))
DEALLOCATE (work)
ALLOCATE (work(lwork))
! SVD
CALL pdgesvd('V', 'V', m, m, matrix_lu%local_data, 1, 1, desca, eigval, matrix_eigvl%local_data, &
1, 1, descu, matrix_eigvr_t%local_data, 1, 1, descvt, work, lwork, myinfo)
IF (matrix_lu%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
WRITE (message, "(A,I0,A)") "ERROR in PDGESVD: Matrix diagonalization failed (INFO = ", myinfo, ")"
IF (PRESENT(info)) THEN
CPWARN(TRIM(message))
ELSE
CPABORT(TRIM(message))
END IF
END IF
#else
! Workspace query
CALL dgesvd('S', 'S', m, m, matrix_lu%local_data, m, eigval, matrix_eigvl%local_data, &
m, matrix_eigvr_t%local_data, m, work, lwork, myinfo)
IF (myinfo /= 0) THEN
WRITE (message, "(A,I0,A)") "ERROR in DGESVD: Work space query failed (INFO = ", myinfo, ")"
IF (PRESENT(info)) THEN
CPWARN(TRIM(message))
ELSE
CPABORT(TRIM(message))
END IF
END IF
! SVD
lwork = INT(work(1))
DEALLOCATE (work)
ALLOCATE (work(lwork))
work(:) = 0.0_dp
CALL dgesvd('S', 'S', m, m, matrix_lu%local_data, m, eigval, matrix_eigvl%local_data, &
m, matrix_eigvr_t%local_data, m, work, lwork, myinfo)
IF (myinfo /= 0) THEN
WRITE (message, "(A,I0,A)") "ERROR in DGESVD: Matrix diagonalization failed (INFO = ", myinfo, ")"
IF (PRESENT(info)) THEN
CPWARN(TRIM(message))
ELSE
CPABORT(TRIM(message))
END IF
END IF
#endif
! Release intermediary matrices
DEALLOCATE (work)
CALL cp_fm_release(matrix_lu)
IF (PRESENT(info)) info = myinfo
CALL timestop(handle)
END SUBROUTINE cp_fm_svd
! **************************************************************************************************
!> \brief ...
!> \param matrix ...

145
src/input_cp2k_mp2.F
View file

@ -49,7 +49,10 @@ MODULE input_cp2k_mp2
section_create, &
section_release, &
section_type
USE input_val_types, ONLY: integer_t, &
USE input_val_types, ONLY: char_t, &
integer_t, &
logical_t, &
real_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: newline, &
@ -1109,7 +1112,7 @@ CONTAINS
CALL keyword_create(keyword, __LOCATION__, name="LOWER_BOUND", &
description="Lower bound for GW-DOS in eV.", &
usage="GW_DOS_LOWER_BOUND -20.0", &
usage="LOWER_BOUND -20.0", &
default_r_val=cp_unit_to_cp2k(value=-20.0_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(gw_dos_section, keyword)
@ -1117,7 +1120,7 @@ CONTAINS
CALL keyword_create(keyword, __LOCATION__, name="UPPER_BOUND", &
description="Upper bound for GW-DOS in eV.", &
usage="GW_DOS_UPPER_BOUND 5.0", &
usage="UPPER_BOUND 5.0", &
default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(gw_dos_section, keyword)
@ -1125,7 +1128,7 @@ CONTAINS
CALL keyword_create(keyword, __LOCATION__, name="STEP", &
description="Difference of two consecutive energy levels for GW-DOS.", &
usage="GW_DOS_PRECISION 0.1", &
usage="STEP 0.1", &
default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(gw_dos_section, keyword)
@ -1133,14 +1136,14 @@ CONTAINS
CALL keyword_create(keyword, __LOCATION__, name="MIN_LEVEL_SPECTRAL", &
description="Lowest energy level to print the self energy to files.", &
usage="MIN_LEVEL_SELF_ENERGY 3", &
usage="MIN_LEVEL_SPECTRAL 3", &
default_i_val=1)
CALL section_add_keyword(gw_dos_section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAX_LEVEL_SPECTRAL", &
description="Highest energy level to print the self energy to files.", &
usage="MAX_LEVEL_SELF_ENERGY 6", &
usage="MAX_LEVEL_SPECTRAL 6", &
default_i_val=0)
CALL section_add_keyword(gw_dos_section, keyword)
CALL keyword_release(keyword)
@ -1298,7 +1301,8 @@ CONTAINS
"$\left( \begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array} \right) = "// &
"\Omega^{(n)}\left(\begin{array}{cc}1&0\\0&-1\end{array}\right)$ "// &
"$\left(\begin{array}{cc}\mathbf{X}^{(n)}\\\mathbf{Y}^{(n)}\end{array}\right)$ "// &
"enables of electronic excitation energies $\Omega^{(n)}$. The BSE can be solved by diagonalizing "// &
"enables, for example, the computation of electronic excitation energies $\Omega^{(n)}$ "// &
"as well as optical properties. The BSE can be solved by diagonalizing "// &
"the full ABBA-matrix or by setting B=0, i.e. within the Tamm-Dancoff approximation (TDA). "// &
"Preliminary reference: Eq. (35) in PRB 92, 045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209", &
n_keywords=8, n_subsections=3, repeats=.FALSE.)
@ -1376,12 +1380,29 @@ CONTAINS
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NUM_PRINT_EXC", &
description="Number of printed excitation energies. Does not affect computation time.", &
description="Number of printed excitation levels with respective "// &
"energies and oscillator strengths. Does not affect computation time.", &
usage="NUM_PRINT_EXC 25", &
default_i_val=25)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NUM_PRINT_EXC_DESCR", &
description="Number of excitation levels for which the exciton "// &
"descriptors are computed. Negative or too large "// &
"NUM_PRINT_EXC_DESCR defaults to NUM_PRINT_EXC.", &
usage="NUM_PRINT_EXC_DESCR 5", &
default_i_val=0)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PRINT_DIRECTIONAL_EXC_DESCR", &
description="Activates printing of exciton descriptors per direction.", &
usage="&PRINT_DIRECTIONAL_EXC_DESCR .TRUE.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_X", &
description="Threshold for printing contributions of singleparticle "// &
"transitions, i.e. elements of the eigenvectors $X_{ia}^{(n)}$ and $Y_{ia}^{(n)}$.", &
@ -1400,8 +1421,96 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
NULLIFY (subsection)
CALL create_bse_nto_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE
SUBROUTINE create_bse_nto_section(print_key)
TYPE(section_type), POINTER :: print_key
TYPE(keyword_type), POINTER :: keyword
CPASSERT(.NOT. ASSOCIATED(print_key))
CALL cp_print_key_section_create(print_key, __LOCATION__, name="NTO_ANALYSIS", &
description="Perform a natural transition orbital analysis, i.e. the transition density matrix "// &
"$T^{(n)}=\left( \begin{array}{cc}0& {X}^{(n)}\\ \left({Y}^{(n)} \right)^T & 0\end{array} \right)$ "// &
"is decomposed into its singular values "// &
"$T^{(n)} = {U}^{(n)} {\Lambda^{(n)}} \left({V}^{(n)}\right)^T$ "// &
"in order to compute the NTO pairs "// &
"$\phi_I^{(n)}(\mathbf{r}_e) = \sum_{p=1}^{N_b} \varphi_p(\mathbf{r}_e) V_{p,I}^{(n)}$ for the electron and "// &
"$\chi_I^{(n)}(\mathbf{r}_h) = \sum_{q=1}^{N_b} \varphi_q(\mathbf{r}_h) U_{q,I}^{(n)}$ for the hole.", &
print_level=debug_print_level + 1, & ! Corresponds to "off" as default behavior
filename="BSE-NTO_ANALYSIS") ! All other print levels will trigger the analysis
! cf. input/cp_output_handling.F:cp_print_key_section_create
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_NTO_EIGVAL", &
description="Threshold for NTO eigenvalues, i.e. only "// &
"${\left(\lambda_I^{(n)}\right)}^2$ > EPS_NTO_EIGVAL are considered.", &
usage="EPS_NTO_EIGVAL 0.01", &
n_var=1, &
type_of_var=real_t, &
default_r_val=0.01_dp)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_OSC_STR", &
description="Threshold of oscillator strengths $f^{(n)}$ for an excitation level. "// &
"In case, the excitation level n has a smaller oscillator strength, the "// &
"NTOs for this excitation level are not printed.", &
usage="EPS_OSC_STR 0.01", &
n_var=1, &
type_of_var=real_t, &
default_r_val=-1.0_dp)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NUM_PRINT_EXC_NTOS", &
description="Number of excitation level $n$ up to which NTOs are printed. "// &
"By default, this is set to NUM_PRINT_EXC. Negative or too large "// &
"NUM_PRINT_EXC_NTOS defaults to NUM_PRINT_EXC.", &
usage="NUM_PRINT_EXC_NTOS 5", &
n_var=1, &
type_of_var=integer_t, &
default_i_val=-1)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="STATE_LIST", &
description="Specifies a list of excitation levels $n$ for which NTOs are printed. "// &
"Overrides NUM_PRINT_EXC_NTOS.", &
usage="STATE_LIST {integer} {integer} .. {integer}", &
n_var=-1, type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CUBE_FILES", &
description="Print NTOs on Cube Files", &
usage="CUBE_FILES {logical}", repeats=.FALSE., n_var=1, &
default_l_val=.TRUE., lone_keyword_l_val=.TRUE., type_of_var=logical_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="STRIDE", &
description="The stride (X,Y,Z) used to write the cube file "// &
"(larger values result in smaller cube files). Provide 3 numbers (for X,Y,Z) or"// &
" 1 number valid for all components.", &
usage="STRIDE 2 2 2", n_var=-1, default_i_vals=(/2, 2, 2/), type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="APPEND", &
description="append the cube files when they already exist", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_bse_nto_section
SUBROUTINE create_bse_spectrum_section(section)
TYPE(section_type), POINTER :: section
@ -1409,21 +1518,27 @@ CONTAINS
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="BSE_SPECTRUM", &
description="Parameters influencing the output of the optical spectrum "// &
"obtained from a BSE calculation.", &
description="Parameters influencing the output of the optical absorption spectrum, i.e. "// &
"the dynamical dipole polarizability tensor $\alpha_{\mu,\mu'}(\omega)$ "// &
"($\mu,\mu'\in\{x,y,z\}$), obtained from a BSE calculation, which is defined as "// &
"$ \alpha_{\mu,\mu'}(\omega) = \sum_n \frac{2 E^{(n)} d^{(n)}_{\mu} d^{(n)}_{\mu'}} "// &
"{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2} $. "// &
"The printed file will contain the specified frequency range $\omega$ and the "// &
"corresponding imaginary part of the average $\bar{\alpha}(\omega)=\frac{1}{3}\mathrm{Tr} "// &
"\left[ \alpha_{\mu,\mu'}(\omega)\right]$ as well as of the elements of $\alpha_{\mu,\mu'}(\omega)$.", &
n_keywords=9, n_subsections=0, repeats=.FALSE.)
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
description="Activates printing of optical spectrum after a BSE calculation.", &
description="Activates printing of optical absorption spectrum from the BSE calculation.", &
usage="&BSE .TRUE.", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY_STEP_SIZE", &
description="Step size of frequencies for the optical spectrum.", &
description="Step size of frequencies for the optical absorption spectrum.", &
usage="FREQUENCY_STEP_SIZE 0.1", unit_str="eV", &
type_of_var=real_t, &
default_r_val=cp_unit_to_cp2k(value=0.1_dp, unit_str="eV"))
@ -1431,14 +1546,14 @@ CONTAINS
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY_STARTING_POINT", &
description="First frequency to print in the optical spectrum.", &
description="First frequency to print in the optical absorption spectrum.", &
usage="FREQUENCY_STARTING_POINT 0", unit_str="eV", &
type_of_var=real_t, default_r_val=0.0_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FREQUENCY_END_POINT", &
description="Last frequency to print in the optical spectrum.", &
description="Last frequency to print in the optical absorption spectrum.", &
usage="FREQUENCY_END_POINT 0", unit_str="eV", &
type_of_var=real_t, &
default_r_val=cp_unit_to_cp2k(value=100.0_dp, unit_str="eV"))
@ -1446,7 +1561,7 @@ CONTAINS
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ETA_LIST", &
description="List of broadening of the peaks in the optical spectrum.", &
description="List of broadening of the peaks in the optical absorption spectrum.", &
usage="ETA_LIST 0.01 ...", unit_str="eV", &
default_r_vals=(/cp_unit_to_cp2k(value=0.01_dp, unit_str="eV")/), &
type_of_var=real_t, n_var=-1)

10
src/mp2_gpw.F
View file

@ -212,7 +212,7 @@ CONTAINS
! check if we want to do imaginary time
do_im_time = mp2_env%do_im_time
do_bse = qs_env%mp2_env%ri_g0w0%do_bse
do_bse = qs_env%mp2_env%bse%do_bse
do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
IF (do_kpoints_cubic_RPA .AND. mp2_env%ri_rpa%do_ri_g0w0) THEN
@ -301,11 +301,13 @@ CONTAINS
CALL m_flush(unit_nr)
END IF
IF (.NOT. mp2_env%ri_g0w0%print_local_bandgap) THEN
! MG: Disable logger layer for BSE, misses some key information to print cube files properly
IF (.NOT. (mp2_env%ri_g0w0%print_local_bandgap .OR. mp2_env%bse%do_nto_analysis)) THEN
! a logger
NULLIFY (logger_sub)
CALL cp_logger_create(logger_sub, para_env=para_env_sub, &
default_global_unit_nr=local_unit_nr, close_global_unit_on_dealloc=.FALSE.)
default_global_unit_nr=local_unit_nr, &
close_global_unit_on_dealloc=.FALSE.)
CALL cp_logger_set(logger_sub, local_filename="MP2_localLog")
! set to a custom print level (we could also have a different print level for para_env%source)
logger_sub%iter_info%print_level = mp2_env%mp2_gpw%print_level
@ -639,7 +641,7 @@ CONTAINS
CALL cp_blacs_env_release(blacs_env_sub)
IF (.NOT. mp2_env%ri_g0w0%print_local_bandgap) THEN
IF (.NOT. (mp2_env%ri_g0w0%print_local_bandgap .OR. mp2_env%bse%do_nto_analysis)) THEN
CALL cp_rm_default_logger()
CALL cp_logger_release(logger_sub)
END IF

76
src/mp2_setup.F
View file

@ -18,7 +18,8 @@ MODULE mp2_setup
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_should_output,&
cp_print_key_unit_nr
cp_print_key_unit_nr,&
debug_print_level
USE cp_parser_methods, ONLY: read_float_object
USE input_constants, ONLY: &
do_eri_mme, do_potential_short, evGW, evGW0, mp2_method_direct, mp2_method_gpw, &
@ -67,7 +68,7 @@ CONTAINS
DIMENSION(:), POINTER :: string_pointer
CHARACTER(LEN=max_line_length) :: error_message
INTEGER :: gw_sc, handle, i, i_special_kp, ival, &
unit_nr
print_level_NTO_analysis, unit_nr
INTEGER, DIMENSION(:), POINTER :: tmplist
LOGICAL :: do_mp2, do_opt_ri_basis, do_ri_mp2, &
do_ri_sos_mp2, do_rpa
@ -190,51 +191,76 @@ CONTAINS
NULLIFY (bse_section)
bse_section => section_vals_get_subs_vals(mp2_section, "RI_RPA%GW%BSE")
CALL section_vals_val_get(bse_section, "_SECTION_PARAMETERS_", &
l_val=mp2_env%ri_g0w0%do_bse)
l_val=mp2_env%bse%do_bse)
CALL section_vals_val_get(bse_section, "SPIN_CONFIG", &
i_val=mp2_env%ri_g0w0%bse_spin_config)
i_val=mp2_env%bse%bse_spin_config)
CALL section_vals_val_get(bse_section, "ENERGY_CUTOFF_OCC", &
r_val=mp2_env%ri_g0w0%bse_cutoff_occ)
r_val=mp2_env%bse%bse_cutoff_occ)
CALL section_vals_val_get(bse_section, "ENERGY_CUTOFF_EMPTY", &
r_val=mp2_env%ri_g0w0%bse_cutoff_empty)
r_val=mp2_env%bse%bse_cutoff_empty)
CALL section_vals_val_get(bse_section, "NUM_PRINT_EXC", &
i_val=mp2_env%ri_g0w0%num_print_exc)
i_val=mp2_env%bse%num_print_exc)
CALL section_vals_val_get(bse_section, "NUM_PRINT_EXC_DESCR", &
i_val=mp2_env%bse%num_print_exc_descr)
CALL section_vals_val_get(bse_section, "PRINT_DIRECTIONAL_EXC_DESCR", &
l_val=mp2_env%bse%print_directional_exc_descr)
CALL section_vals_val_get(bse_section, "EPS_X", &
r_val=mp2_env%ri_g0w0%eps_x)
r_val=mp2_env%bse%eps_x)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%DAVIDSON_ABORT_COND", &
i_val=mp2_env%ri_g0w0%davidson_abort_cond)
i_val=mp2_env%bse%davidson_abort_cond)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%NUM_EXC_EN", &
i_val=mp2_env%ri_g0w0%num_exc_en)
i_val=mp2_env%bse%num_exc_en)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%NUM_ADD_START_Z_SPACE", &
i_val=mp2_env%ri_g0w0%num_add_start_z_space)
i_val=mp2_env%bse%num_add_start_z_space)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%FAC_MAX_Z_SPACE", &
i_val=mp2_env%ri_g0w0%fac_max_z_space)
i_val=mp2_env%bse%fac_max_z_space)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%NUM_NEW_T", &
i_val=mp2_env%ri_g0w0%num_new_t)
i_val=mp2_env%bse%num_new_t)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%EPS_RES", &
r_val=mp2_env%ri_g0w0%eps_res)
r_val=mp2_env%bse%eps_res)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%EPS_EXC_EN", &
r_val=mp2_env%ri_g0w0%eps_exc_en)
r_val=mp2_env%bse%eps_exc_en)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%NUM_DAVIDSON_ITER", &
i_val=mp2_env%ri_g0w0%num_davidson_iter)
i_val=mp2_env%bse%num_davidson_iter)
CALL section_vals_val_get(bse_section, "BSE_ITERAT%Z_SPACE_ENERGY_CUTOFF", &
r_val=mp2_env%ri_g0w0%z_space_energy_cutoff)
r_val=mp2_env%bse%z_space_energy_cutoff)
CALL section_vals_val_get(bse_section, "BSE_DIAG_METHOD", &
i_val=mp2_env%ri_g0w0%bse_diag_method)
i_val=mp2_env%bse%bse_diag_method)
CALL section_vals_val_get(bse_section, "TDA", &
i_val=mp2_env%ri_g0w0%flag_tda)
i_val=mp2_env%bse%flag_tda)
CALL section_vals_val_get(bse_section, "BSE_DEBUG_PRINT", &
l_val=mp2_env%ri_g0w0%bse_debug_print)
l_val=mp2_env%bse%bse_debug_print)
CALL section_vals_val_get(bse_section, "BSE_SPECTRUM%_SECTION_PARAMETERS_", &
l_val=mp2_env%ri_g0w0%bse_print_spectrum)
l_val=mp2_env%bse%bse_print_spectrum)
CALL section_vals_val_get(bse_section, "BSE_SPECTRUM%FREQUENCY_STEP_SIZE", &
r_val=mp2_env%ri_g0w0%bse_spectrum_freq_step_size)
r_val=mp2_env%bse%bse_spectrum_freq_step_size)
CALL section_vals_val_get(bse_section, "BSE_SPECTRUM%FREQUENCY_STARTING_POINT", &
r_val=mp2_env%ri_g0w0%bse_spectrum_freq_start)
r_val=mp2_env%bse%bse_spectrum_freq_start)
CALL section_vals_val_get(bse_section, "BSE_SPECTRUM%FREQUENCY_END_POINT", &
r_val=mp2_env%ri_g0w0%bse_spectrum_freq_end)
r_val=mp2_env%bse%bse_spectrum_freq_end)
CALL section_vals_val_get(bse_section, "BSE_SPECTRUM%ETA_LIST", &
r_vals=mp2_env%ri_g0w0%bse_eta_spectrum_list)
r_vals=mp2_env%bse%bse_eta_spectrum_list)
CALL section_vals_val_get(bse_section, "NTO_ANALYSIS%_SECTION_PARAMETERS_", &
i_val=print_level_NTO_analysis)
! All other print levels besides "off" should trigger the analysis
! cf. input/cp_output_handling.F:cp_print_key_section_create
IF (print_level_NTO_analysis /= debug_print_level + 1) THEN
mp2_env%bse%do_nto_analysis = .TRUE.
ELSE
mp2_env%bse%do_nto_analysis = .FALSE.
END IF
CALL section_vals_val_get(bse_section, "NTO_ANALYSIS%EPS_NTO_EIGVAL", &
r_val=mp2_env%bse%eps_nto_eigval)
CALL section_vals_val_get(bse_section, "NTO_ANALYSIS%EPS_OSC_STR", &
r_val=mp2_env%bse%eps_nto_osc_str)
CALL section_vals_val_get(bse_section, "NTO_ANALYSIS%NUM_PRINT_EXC_NTOS", &
i_val=mp2_env%bse%num_print_exc_ntos)
CALL section_vals_val_get(bse_section, "NTO_ANALYSIS%STATE_LIST", &
EXPLICIT=mp2_env%bse%explicit_nto_list)
IF (mp2_env%bse%explicit_nto_list) THEN
CALL section_vals_val_get(bse_section, "NTO_ANALYSIS%STATE_LIST", &
i_vals=mp2_env%bse%bse_nto_state_list)
END IF
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%IMAGE_CHARGE_MODEL", &
l_val=mp2_env%ri_g0w0%do_ic_model)

67
src/mp2_types.F
View file

@ -226,25 +226,7 @@ MODULE mp2_types
do_aux_bas_gw = .FALSE.
REAL(KIND=dp) :: frac_aux_mos = 0.0_dp
INTEGER :: num_omega_points = 0
INTEGER :: bse_spin_config = 0, &
bse_diag_method = 0, &
flag_tda = 0, &
num_exc_en = 0, &
num_print_exc = 0, &
num_add_start_z_space = 0, &
fac_max_z_space = 0, &
num_new_t = 0, &
num_davidson_iter = 0, &
davidson_abort_cond = 0
REAL(KIND=dp) :: eps_res = 0.0_dp, &
eps_exc_en = 0.0_dp, &
eps_x = 0.0_dp, &
bse_cutoff_occ = 0.0_dp, &
bse_cutoff_empty = 0.0_dp, &
z_space_energy_cutoff = 0.0_dp
LOGICAL :: do_bse = .FALSE., &
bse_debug_print = .FALSE., &
do_ic_model = .FALSE., &
LOGICAL :: do_ic_model = .FALSE., &
print_ic_values = .FALSE.
REAL(KIND=dp) :: eps_dist = 0.0_dp
TYPE(one_dim_real_array), DIMENSION(2) :: ic_corr_list = one_dim_real_array(NULL())
@ -270,13 +252,6 @@ MODULE mp2_types
soc_energy_window = 0.0_dp
INTEGER, DIMENSION(:), POINTER :: stride_loc_bandgap => NULL()
!BSE optical spectrum
REAL(KIND=dp) :: bse_spectrum_freq_step_size = 0.0_dp, &
bse_spectrum_freq_start = 0.0_dp, &
bse_spectrum_freq_end = 0.0_dp
LOGICAL :: bse_print_spectrum = .FALSE.
REAL(KIND=dp), DIMENSION(:), POINTER :: bse_eta_spectrum_list => NULL()
!GW_DOS
REAL(KIND=dp) :: dos_upper = 0.0_dp, &
dos_lower = 0.0_dp, &
@ -326,6 +301,44 @@ MODULE mp2_types
INTEGER :: max_parallel_comm = 0
END TYPE grad_util
TYPE bse_type
INTEGER :: bse_spin_config = 0, &
bse_diag_method = 0, &
flag_tda = 0, &
num_exc_en = 0, &
num_print_exc = 0, &
num_print_exc_descr = 0, &
num_add_start_z_space = 0, &
fac_max_z_space = 0, &
num_new_t = 0, &
num_davidson_iter = 0, &
davidson_abort_cond = 0
REAL(KIND=dp) :: eps_res = 0.0_dp, &
eps_exc_en = 0.0_dp, &
eps_x = 0.0_dp, &
bse_cutoff_occ = 0.0_dp, &
bse_cutoff_empty = 0.0_dp, &
z_space_energy_cutoff = 0.0_dp
LOGICAL :: do_bse = .FALSE., &
bse_debug_print = .FALSE., &
print_directional_exc_descr = .FALSE.
!BSE optical spectrum
REAL(KIND=dp) :: bse_spectrum_freq_step_size = 0.0_dp, &
bse_spectrum_freq_start = 0.0_dp, &
bse_spectrum_freq_end = 0.0_dp
LOGICAL :: bse_print_spectrum = .FALSE.
REAL(KIND=dp), DIMENSION(:), POINTER :: bse_eta_spectrum_list => NULL()
!BSE NTO prints
LOGICAL :: do_nto_analysis = .FALSE., &
explicit_nto_list = .FALSE.
REAL(KIND=dp) :: eps_nto_eigval = 0.0_dp, &
eps_nto_osc_str = 0.0_dp
INTEGER :: num_print_exc_ntos = 0
INTEGER, DIMENSION(:), POINTER :: bse_nto_state_list => NULL()
INTEGER, DIMENSION(:), ALLOCATABLE :: bse_nto_state_list_final
END TYPE bse_type
TYPE mp2_type
INTEGER :: method = mp2_method_none
TYPE(mp2_laplace_type) :: ri_laplace = mp2_laplace_type()
@ -358,11 +371,13 @@ MODULE mp2_types
mo_coeff_o=NULL(), &
mo_coeff_v=NULL(), &
P_mo=NULL(), W_mo=NULL(), L_jb=NULL())
TYPE(bse_type) :: bse = bse_type(bse_nto_state_list_final=NULL())
#else
TYPE(ri_rpa_im_time_type) :: ri_rpa_im_time = ri_rpa_im_time_type()
TYPE(ri_g0w0_type) :: ri_g0w0 = ri_g0w0_type()
TYPE(ri_basis_opt) :: ri_opt_param = ri_basis_opt()
TYPE(grad_util) :: ri_grad = grad_util()
TYPE(bse_type) :: bse = bse_type()
#endif
REAL(KIND=dp) :: mp2_memory = 0.0_dp, &
scale_S = 0.0_dp, &

30
src/rpa_gw.F
View file

@ -1298,6 +1298,7 @@ CONTAINS
E_CBM_SCF, E_CBM_SCF_beta, E_VBM_GW, E_VBM_GW_beta, E_VBM_SCF, E_VBM_SCF_beta, stop_crit
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vec_gw_dos
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: m_value, vec_gw_energ, z_value
TYPE(cp_logger_type), POINTER :: logger
TYPE(kpoint_type), POINTER :: kpoints_Sigma
CALL timeset(routineN, handle)
@ -1560,7 +1561,12 @@ CONTAINS
homo, unit_nr, do_soc_gw=.TRUE., do_soc_scf=.FALSE.)
END IF
iunit = cp_logger_get_default_unit_nr()
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iunit = cp_logger_get_default_unit_nr()
ELSE
iunit = -1
END IF
IF (dos_precision /= 0.0_dp) THEN
IF (iunit > 0) THEN
@ -4341,6 +4347,7 @@ CONTAINS
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vec_omega_fit_gw_sign, &
vec_omega_fit_gw_sign_reorder, &
vec_sigma_imag, vec_sigma_real
TYPE(cp_logger_type), POINTER :: logger
CALL timeset(routineN, handle)
@ -4437,7 +4444,12 @@ CONTAINS
IF (ndos /= 0) THEN
! Hedin shift not implemented
CPASSERT(.NOT. do_hedin_shift)
iunit = cp_logger_get_default_unit_nr()
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iunit = cp_logger_get_default_unit_nr()
ELSE
iunit = -1
END IF
DO idos = 1, ndos
omega_dos = dos_lower_bound + REAL(idos - 1, KIND=dp)*dos_precision
omega_dos_pade_eval = omega_dos - e_fermi
@ -4469,7 +4481,12 @@ CONTAINS
END IF
IF (PRESENT(min_level_self_energy) .AND. PRESENT(max_level_self_energy)) THEN
iunit = cp_logger_get_default_unit_nr()
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iunit = cp_logger_get_default_unit_nr()
ELSE
iunit = -1
END IF
IF (n_level_gw_ref >= min_level_self_energy .AND. &
n_level_gw_ref <= max_level_self_energy .AND. iunit > 0) THEN
@ -4561,7 +4578,12 @@ CONTAINS
END IF
iunit = cp_logger_get_default_unit_nr()
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iunit = cp_logger_get_default_unit_nr()
ELSE
iunit = -1
END IF
CALL open_file(TRIM(filename), unit_number=iunit, file_status="UNKNOWN", file_action="WRITE")
num_omega = 10000

6
src/rpa_gw_sigma_x.F
View file

@ -555,12 +555,12 @@ CONTAINS
gw_corr_lev_virt = mp2_env%ri_g0w0%corr_mos_virt
! If BSE is invoked, manipulate corrected MO number
IF (mp2_env%ri_g0w0%do_bse) THEN
IF (mp2_env%bse%do_bse) THEN
! Logic: If cutoff is negative, all MOs are included in BSE, i.e. we need to correct them all
! If cutoff is positive, we can reduce the number of MOs to be corrected and force gw_corr_lev_...
! to a sufficiently large number by setting it to -2 and read indices afterwards
! Handling for occupied levels
IF (mp2_env%ri_g0w0%bse_cutoff_occ < 0) THEN
IF (mp2_env%bse%bse_cutoff_occ < 0) THEN
gw_corr_lev_occ = -1
ELSE
IF (gw_corr_lev_occ > 0) THEN
@ -568,7 +568,7 @@ CONTAINS
END IF
END IF
! Handling for virtual levels
IF (mp2_env%ri_g0w0%bse_cutoff_empty < 0) THEN
IF (mp2_env%bse%bse_cutoff_empty < 0) THEN
gw_corr_lev_virt = -1
ELSE
IF (gw_corr_lev_virt > 0) THEN

6
src/rpa_main.F
View file

@ -979,7 +979,11 @@ CONTAINS
END IF
! release blacs_env
IF (.NOT. my_blacs_S_ext) CALL cp_blacs_env_release(blacs_env)
IF (.NOT. my_blacs_S_ext) THEN
CALL cp_blacs_env_release(blacs_env)
ELSE
NULLIFY (blacs_env)
END IF
CALL timestop(handle)

83
tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_EXC_DESCR.inp Normal file
View file

@ -0,0 +1,83 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT BSE_H2O_PBE
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END TIMINGS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER WAVELET
&END POISSON
&QS
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 100
SCF_GUESS ATOMIC
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
&WF_CORRELATION
&RI_RPA
QUADRATURE_POINTS 100
&GW
CORR_MOS_OCC 1000
CORR_MOS_VIRT 1000
RI_SIGMA_X
SELF_CONSISTENCY G0W0
&BSE
BSE_DEBUG_PRINT
BSE_DIAG_METHOD FULLDIAG
ENERGY_CUTOFF_EMPTY 50.0
ENERGY_CUTOFF_OCC 60.0
NUM_PRINT_EXC_DESCR -1
TDA OFF
&END BSE
&END GW
&END RI_RPA
&END WF_CORRELATION
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.000 6.000 6.000
PERIODIC NONE
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&TOPOLOGY
COORD_FILE_FORMAT xyz
COORD_FILE_NAME H2O_gas.xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

84
tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_NTOs.inp Normal file
View file

@ -0,0 +1,84 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT BSE_H2O_PBE
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END TIMINGS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER WAVELET
&END POISSON
&QS
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 100
SCF_GUESS ATOMIC
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
&WF_CORRELATION
&RI_RPA
QUADRATURE_POINTS 100
&GW
CORR_MOS_OCC 1000
CORR_MOS_VIRT 1000
RI_SIGMA_X
SELF_CONSISTENCY G0W0
&BSE
BSE_DIAG_METHOD FULLDIAG
ENERGY_CUTOFF_EMPTY 50.0
ENERGY_CUTOFF_OCC 60.0
TDA OFF
&NTO_ANALYSIS
CUBE_FILES F
&END NTO_ANALYSIS
&END BSE
&END GW
&END RI_RPA
&END WF_CORRELATION
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.000 6.000 6.000
PERIODIC NONE
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&TOPOLOGY
COORD_FILE_FORMAT xyz
COORD_FILE_NAME H2O_gas.xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

83
tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_EXC_DESCR.inp Normal file
View file

@ -0,0 +1,83 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT BSE_H2O_PBE
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END TIMINGS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER WAVELET
&END POISSON
&QS
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 100
SCF_GUESS ATOMIC
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
&WF_CORRELATION
&RI_RPA
QUADRATURE_POINTS 100
&GW
CORR_MOS_OCC 1000
CORR_MOS_VIRT 1000
RI_SIGMA_X
SELF_CONSISTENCY G0W0
&BSE
BSE_DEBUG_PRINT
BSE_DIAG_METHOD FULLDIAG
ENERGY_CUTOFF_EMPTY 50.0
ENERGY_CUTOFF_OCC 60.0
NUM_PRINT_EXC_DESCR -1
TDA ON
&END BSE
&END GW
&END RI_RPA
&END WF_CORRELATION
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.000 6.000 6.000
PERIODIC NONE
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&TOPOLOGY
COORD_FILE_FORMAT xyz
COORD_FILE_NAME H2O_gas.xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

84
tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_NTOs.inp Normal file
View file

@ -0,0 +1,84 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT BSE_H2O_PBE
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END TIMINGS
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER WAVELET
&END POISSON
&QS
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 100
SCF_GUESS ATOMIC
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
&WF_CORRELATION
&RI_RPA
QUADRATURE_POINTS 100
&GW
CORR_MOS_OCC 1000
CORR_MOS_VIRT 1000
RI_SIGMA_X
SELF_CONSISTENCY G0W0
&BSE
BSE_DIAG_METHOD FULLDIAG
ENERGY_CUTOFF_EMPTY 50.0
ENERGY_CUTOFF_OCC 60.0
TDA ON
&NTO_ANALYSIS
CUBE_FILES F
&END NTO_ANALYSIS
&END BSE
&END GW
&END RI_RPA
&END WF_CORRELATION
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.000 6.000 6.000
PERIODIC NONE
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&TOPOLOGY
COORD_FILE_FORMAT xyz
COORD_FILE_NAME H2O_gas.xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

6
tests/QS/regtest-bse/TEST_FILES
View file

@ -10,9 +10,15 @@ TDA_H2O_PBE_evGW.inp 113 7e-02
BSE_H2O_PBE_G0W0.inp 114 2e-05 0.014
BSE_H2O_PBE_evGW0.inp 114 2e-05 0.015
BSE_H2O_PBE_evGW.inp 114 2e-04 0.016
TDA_H2O_PBE_G0W0_NTOs.inp 121 2e-03 0.01983
BSE_H2O_PBE_G0W0_NTOs.inp 122 2e-03 0.01164
TDA_H2O_PBE_G0W0_EXC_DESCR.inp 123 5e-04 1.8290
BSE_H2O_PBE_G0W0_EXC_DESCR.inp 124 5e-04 1.8258
#EOF
#Low convergence criteria for evGW make a higher threshold necessary for these tests
#Logic: Absolute errors (for G0W0/evGW0) should be within 1e-4, i.e. relative errors are
# in the range of 2e-5 for energies.
# For oscillator strengths, they should be within 1e-3, so relative errors
# are in the range of 7e-2
# First 6 check excitation energies, then 6 checking oscillator strengths, two check NTOs
# and the last one checks exciton descriptors in TDA

6
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
120
124
Total energy:!3
MD| Potential energy!5
Total energy \[eV\]:!4
@ -119,6 +119,10 @@ MOMENTS_TRACE_RE| 0.20000000E+000!3
POLARIZABILITY| 0.32436607E+002!4
DEBUG:: Total Force !6
SMEAGOL| Number of electrons: !5
BSE| 3 -TDA- 2 !5
BSE| 3 -ABBA- 2 !5
BSE|DEBUG| Exciton descriptor d_exc with TDA for n=1 is !10
BSE|DEBUG| Exciton descriptor d_exc without TDA for n=1 is !10
#
# these are the tests the can be selected for regtesting.
# do regtest will grep for test_grep (first column) and look if the numeric value