From 4dc10b34915a2b21df1dca61771b8f2b68eee459 Mon Sep 17 00:00:00 2001 From: Max Graml <78024843+MGraml@users.noreply.github.com> Date: Wed, 20 Nov 2024 16:34:36 +0100 Subject: [PATCH] BSE: new features and update manual (#3781) --- .../properties/optical/bethe-salpeter.md | 513 +++++++-- src/CMakeLists.txt | 2 + src/bse_full_diag.F | 313 +++--- src/bse_iterative.F | 24 +- src/bse_main.F | 17 +- src/bse_print.F | 793 ++++++++++++++ src/bse_properties.F | 981 +++++++++++++++++ src/bse_util.F | 986 +++++++++--------- src/common/bibliography.F | 23 +- src/fm/cp_fm_diag.F | 138 +++ src/input_cp2k_mp2.F | 145 ++- src/mp2_gpw.F | 10 +- src/mp2_setup.F | 76 +- src/mp2_types.F | 67 +- src/rpa_gw.F | 30 +- src/rpa_gw_sigma_x.F | 6 +- src/rpa_main.F | 6 +- .../BSE_H2O_PBE_G0W0_EXC_DESCR.inp | 83 ++ .../QS/regtest-bse/BSE_H2O_PBE_G0W0_NTOs.inp | 84 ++ .../TDA_H2O_PBE_G0W0_EXC_DESCR.inp | 83 ++ .../QS/regtest-bse/TDA_H2O_PBE_G0W0_NTOs.inp | 84 ++ tests/QS/regtest-bse/TEST_FILES | 6 + tests/TEST_TYPES | 6 +- 23 files changed, 3600 insertions(+), 876 deletions(-) create mode 100644 src/bse_print.F create mode 100644 src/bse_properties.F create mode 100644 tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_EXC_DESCR.inp create mode 100644 tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_NTOs.inp create mode 100644 tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_EXC_DESCR.inp create mode 100644 tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_NTOs.inp diff --git a/docs/methods/properties/optical/bethe-salpeter.md b/docs/methods/properties/optical/bethe-salpeter.md index 7f67f70e12..ec600ca573 100644 --- a/docs/methods/properties/optical/bethe-salpeter.md +++ b/docs/methods/properties/optical/bethe-salpeter.md @@ -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*0*W*0, ev*GW*0 or ev*GW* eigenvalues, see -details in [GW] and in Ref. \[[](#Golze2019)\], i.e. we perform -BSE@*G*0*W*0/ev*GW*0/ev*GW*@DFT. Thus, also input parameters for a -DFT and *GW* calculation influence the BSE calcualation (see full input in Sec. -[3.1](#header-input-file)). We obtain optical properties from BSE solving the following generalized +It is possible to use *G*0*W*0 , ev*GW*0 or ev*GW* eigenvalues, see +details in [GW] and in Ref. \[[](#Golze2019)\], i.e. we perform BSE@*G* 0 +*W*0/ev*GW*0/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*0*W*0, ev*GW*0 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*0 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*0@PBE or the PBE0 functional, - when using BSE@*G*0*W*0@PBE0 (see also + \[[](#Knysh2024)\]. Motivated by the discussion in \[[](#Schambeck2024)\], we strongly recommend + to use BSE@ev*GW*0@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 `aug-cc-pVDZ` basis set should be sufficient for most calculations, but needs to be checked regarding convergence, e.g. using - `aug-cc-pVTZ`. + `aug-cc-pVTZ` + . 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*0@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 H2. 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: -aug-cc-pVDZ, -aug-cc-pVDZ-RIFIT. +`aug-cc-pVDZ`, +`aug-cc-pVDZ-RIFIT`. The geometry for H2 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 H2, 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), diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 05470ff3d9..dd6c95777a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/bse_full_diag.F b/src/bse_full_diag.F index 89aef6b1f0..5facff40ae 100644 --- a/src/bse_full_diag.F +++ b/src/bse_full_diag.F @@ -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 diff --git a/src/bse_iterative.F b/src/bse_iterative.F index 01b1a98a3f..3f4612d08e 100644 --- a/src/bse_iterative.F +++ b/src/bse_iterative.F @@ -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 diff --git a/src/bse_main.F b/src/bse_main.F index 768da2e8f7..7c42cc70b9 100644 --- a/src/bse_main.F +++ b/src/bse_main.F @@ -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) diff --git a/src/bse_print.F b/src/bse_print.F new file mode 100644 index 0000000000..0aa569b690 --- /dev/null +++ b/src/bse_print.F @@ -0,0 +1,793 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! 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 = | _exc |' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', & + 'σ_e = sqrt( _exc - _exc^2 )' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', & + 'σ_h = sqrt( _exc - _exc^2 )' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', & + 'COV_eh = _exc - _exc _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 = | _exc |' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', & + 'σ_e^x = sqrt( _exc - _exc^2 )' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', & + 'σ_h^x = sqrt( _exc - _exc^2 )' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', & + 'COV_eh^x = _exc - _exc _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 diff --git a/src/bse_properties.F b/src/bse_properties.F new file mode 100644 index 0000000000..a4e74ca683 --- /dev/null +++ b/src/bse_properties.F @@ -0,0 +1,981 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! 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 + + ! _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 + ! _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 + ! _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(:) + + ! _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 + ! _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 + ! _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(:) + + ! _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 + ! _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 + ! _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(:) + + ! _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 + ! _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 + ! _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(:) + + ! _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) + ! _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) + ! _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) + ! _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) + ! _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 = |_X - _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( _X - _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( _X - _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 diff --git a/src/bse_util.F b/src/bse_util.F index 403bb4bcfa..0dc59265e9 100644 --- a/src/bse_util.F +++ b/src/bse_util.F @@ -11,7 +11,10 @@ !> 11.2023 created [Maximilian Graml] ! ************************************************************************************************** MODULE bse_util + USE atomic_kind_types, ONLY: atomic_kind_type + USE cell_types, ONLY: cell_type USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_api, ONLY: dbcsr_create,& dbcsr_init_p,& dbcsr_p_type,& @@ -21,9 +24,7 @@ MODULE bse_util USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,& dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set - USE cp_files, ONLY: close_file,& - open_file - USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& + USE cp_fm_basic_linalg, ONLY: cp_fm_trace,& cp_fm_upper_to_full USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& cp_fm_cholesky_invert @@ -32,15 +33,20 @@ MODULE bse_util cp_fm_struct_type USE cp_fm_types, ONLY: cp_fm_create,& cp_fm_get_info,& - 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_logger_get_default_unit_nr + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_type + USE cp_output_handling, ONLY: cp_print_key_finished_output,& + cp_print_key_unit_nr + USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube USE input_constants, ONLY: use_mom_ref_coac - USE kinds, ONLY: dp,& + USE input_section_types, ONLY: section_vals_type + USE kinds, ONLY: default_path_length,& + dp,& int_8 USE message_passing, ONLY: mp_para_env_type,& mp_request_type @@ -48,13 +54,25 @@ MODULE bse_util USE mp2_types, ONLY: integ_mat_buffer_type,& mp2_type USE parallel_gemm_api, ONLY: parallel_gemm + USE particle_list_types, ONLY: particle_list_type + USE particle_types, ONLY: particle_type USE physcon, ONLY: evolt + USE pw_env_types, ONLY: pw_env_get,& + pw_env_type + USE pw_pool_types, ONLY: pw_pool_p_type,& + pw_pool_type + USE pw_types, ONLY: pw_c1d_gs_type,& + pw_r3d_rs_type + USE qs_collocate_density, ONLY: calculate_wavefunction USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type + USE qs_kind_types, ONLY: qs_kind_type USE qs_mo_types, ONLY: get_mo_set,& mo_set_type USE qs_moments, ONLY: build_local_moment_matrix USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_subsys_types, ONLY: qs_subsys_get,& + qs_subsys_type USE rpa_communication, ONLY: communicate_buffer USE util, ONLY: sort,& sort_unique @@ -66,11 +84,11 @@ MODULE bse_util CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_util' - PUBLIC :: mult_B_with_W, fm_general_add_bse, truncate_fm, fm_write_thresh, print_BSE_start_flag, & + PUBLIC :: mult_B_with_W, fm_general_add_bse, truncate_fm, & deallocate_matrices_bse, comp_eigvec_coeff_BSE, sort_excitations, & estimate_BSE_resources, filter_eigvec_contrib, truncate_BSE_matrices, & - get_oscillator_strengths, compute_absorption_spectrum, determine_cutoff_indices, & - print_transition_probabilities + determine_cutoff_indices, adapt_BSE_input_params, get_multipoles_mo, & + reshuffle_eigvec, print_bse_nto_cubes, trace_exciton_descr CONTAINS @@ -158,7 +176,7 @@ CONTAINS ! ************************************************************************************************** !> \brief Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,ab -!> to A_ia, which needs MPI communication. +!> to A_ia, jb which needs MPI communication. !> \param fm_out ... !> \param fm_in ... !> \param beta ... @@ -173,7 +191,8 @@ CONTAINS SUBROUTINE fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, & nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env) - TYPE(cp_fm_type), INTENT(INOUT) :: fm_out, fm_in + TYPE(cp_fm_type), INTENT(INOUT) :: fm_out + TYPE(cp_fm_type), INTENT(IN) :: fm_in REAL(kind=dp) :: beta INTEGER, INTENT(IN) :: nrow_secidx_in, ncol_secidx_in, & nrow_secidx_out, ncol_secidx_out @@ -227,7 +246,7 @@ CONTAINS nrow_block=nrow_block_in, & ncol_block=ncol_block_in) - 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,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, & fm_out%matrix_struct%nrow_global WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, & @@ -271,7 +290,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_2_comm_entry_nums", handle2) - 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,A27)') 'BSE|DEBUG|', 'Communicating entry numbers' END IF @@ -280,7 +299,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_3_alloc_buffer", handle2) - 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,A18)') 'BSE|DEBUG|', 'Allocating buffers' END IF @@ -320,7 +339,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_4_buf_from_fmin_"//fm_out%name, handle2) - 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,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers' END IF @@ -358,7 +377,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_5_comm_buffer", handle2) - 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,A21)') 'BSE|DEBUG|', 'Communicating buffers' END IF @@ -369,7 +388,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_6_buffer_to_fmout"//fm_out%name, handle2) - 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,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name END IF @@ -392,7 +411,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_7_cleanup", handle2) - 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,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers' END IF @@ -493,7 +512,7 @@ CONTAINS nrow_block=nrow_block_in, & ncol_block=ncol_block_in) - 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,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, & fm_out%matrix_struct%nrow_global WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, & @@ -551,7 +570,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_2_comm_entry_nums", handle2) - 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,A27)') 'BSE|DEBUG|', 'Communicating entry numbers' END IF @@ -560,7 +579,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_3_alloc_buffer", handle2) - 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,A18)') 'BSE|DEBUG|', 'Allocating buffers' END IF @@ -600,7 +619,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_4_buf_from_fmin_"//fm_out%name, handle2) - 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,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers' END IF @@ -656,7 +675,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_5_comm_buffer", handle2) - 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,A21)') 'BSE|DEBUG|', 'Communicating buffers' END IF @@ -667,7 +686,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_6_buffer_to_fmout"//fm_out%name, handle2) - 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,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name END IF @@ -690,7 +709,7 @@ CONTAINS CALL timestop(handle2) CALL timeset(routineN//"_7_cleanup", handle2) - 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,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers' END IF @@ -711,113 +730,6 @@ CONTAINS END SUBROUTINE truncate_fm -! ************************************************************************************************** -!> \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 ... -!> \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 fm_mat_S_bar_ia_bse ... @@ -826,13 +738,15 @@ CONTAINS !> \param fm_mat_S_ij_trunc ... !> \param fm_mat_S_ab_trunc ... !> \param fm_mat_Q_static_bse_gemm ... +!> \param mp2_env ... ! ************************************************************************************************** SUBROUTINE deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, & fm_mat_S_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) TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_trunc, & fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, fm_mat_Q_static_bse_gemm + TYPE(mp2_type) :: mp2_env CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_matrices_bse' @@ -846,6 +760,9 @@ CONTAINS CALL cp_fm_release(fm_mat_S_ij_trunc) CALL cp_fm_release(fm_mat_S_ab_trunc) CALL cp_fm_release(fm_mat_Q_static_bse_gemm) + IF (mp2_env%bse%do_nto_analysis) THEN + DEALLOCATE (mp2_env%bse%bse_nto_state_list_final) + END IF CALL timestop(handle) @@ -1108,7 +1025,7 @@ CONTAINS DO ii = 1, nrow_local eigvec_idx = row_indices(ii) eigvec_entry = fm_eigvec%local_data(ii, jj)/SQRT(2.0_dp) - IF (ABS(eigvec_entry) > mp2_env%ri_g0w0%eps_x) THEN + IF (ABS(eigvec_entry) > mp2_env%bse%eps_x) THEN num_entries_to_comm(para_env%mepos) = num_entries_to_comm(para_env%mepos) + 1 END IF END DO @@ -1137,7 +1054,7 @@ CONTAINS DO ii = 1, nrow_local eigvec_idx = row_indices(ii) eigvec_entry = fm_eigvec%local_data(ii, jj)/SQRT(2.0_dp) - IF (ABS(eigvec_entry) > mp2_env%ri_g0w0%eps_x) THEN + IF (ABS(eigvec_entry) > mp2_env%bse%eps_x) THEN buffer_entries(para_env%mepos)%indx(kk, 1) = (eigvec_idx - 1)/virtual + 1 buffer_entries(para_env%mepos)%indx(kk, 2) = MOD(eigvec_idx - 1, virtual) + 1 buffer_entries(para_env%mepos)%msg(kk) = eigvec_entry @@ -1216,16 +1133,16 @@ CONTAINS CALL timeset(routineN, handle) ! Determine index in homo and virtual for truncation - ! Uses indices of outermost orbitals within energy range (-mp2_env%ri_g0w0%bse_cutoff_occ,mp2_env%ri_g0w0%bse_cutoff_empty) - IF (mp2_env%ri_g0w0%bse_cutoff_occ > 0 .OR. mp2_env%ri_g0w0%bse_cutoff_empty > 0) THEN - IF (-mp2_env%ri_g0w0%bse_cutoff_occ .LT. Eigenval(1) - Eigenval(homo) & - .OR. mp2_env%ri_g0w0%bse_cutoff_occ < 0) THEN + ! Uses indices of outermost orbitals within energy range (-mp2_env%bse%bse_cutoff_occ,mp2_env%bse%bse_cutoff_empty) + IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN + IF (-mp2_env%bse%bse_cutoff_occ .LT. Eigenval(1) - Eigenval(homo) & + .OR. mp2_env%bse%bse_cutoff_occ < 0) THEN homo_red = homo homo_incl = 1 ELSE homo_incl = 1 DO i_homo = 1, homo - IF (Eigenval(i_homo) - Eigenval(homo) .GT. -mp2_env%ri_g0w0%bse_cutoff_occ) THEN + IF (Eigenval(i_homo) - Eigenval(homo) .GT. -mp2_env%bse%bse_cutoff_occ) THEN homo_incl = i_homo EXIT END IF @@ -1233,14 +1150,14 @@ CONTAINS homo_red = homo - homo_incl + 1 END IF - IF (mp2_env%ri_g0w0%bse_cutoff_empty .GT. Eigenval(homo + virtual) - Eigenval(homo + 1) & - .OR. mp2_env%ri_g0w0%bse_cutoff_empty < 0) THEN + IF (mp2_env%bse%bse_cutoff_empty .GT. Eigenval(homo + virtual) - Eigenval(homo + 1) & + .OR. mp2_env%bse%bse_cutoff_empty < 0) THEN virt_red = virtual virt_incl = virtual ELSE virt_incl = homo + 1 DO j_virt = 1, virtual - IF (Eigenval(homo + j_virt) - Eigenval(homo + 1) .GT. mp2_env%ri_g0w0%bse_cutoff_empty) THEN + IF (Eigenval(homo + j_virt) - Eigenval(homo + 1) .GT. mp2_env%bse%bse_cutoff_empty) THEN virt_incl = j_virt - 1 EXIT END IF @@ -1307,7 +1224,7 @@ CONTAINS CALL timeset(routineN, handle) ! Determine index in homo and virtual for truncation - ! Uses indices of outermost orbitals within energy range (-mp2_env%ri_g0w0%bse_cutoff_occ,mp2_env%ri_g0w0%bse_cutoff_empty) + ! Uses indices of outermost orbitals within energy range (-mp2_env%bse%bse_cutoff_occ,mp2_env%bse%bse_cutoff_empty) CALL determine_cutoff_indices(Eigenval_scf, & homo, virtual, & @@ -1316,15 +1233,15 @@ CONTAINS mp2_env) IF (unit_nr > 0) THEN - IF (mp2_env%ri_g0w0%bse_cutoff_occ > 0) THEN + IF (mp2_env%bse%bse_cutoff_occ > 0) THEN WRITE (unit_nr, '(T2,A4,T7,A29,T71,F10.3)') 'BSE|', 'Cutoff occupied orbitals [eV]', & - mp2_env%ri_g0w0%bse_cutoff_occ*evolt + mp2_env%bse%bse_cutoff_occ*evolt ELSE WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No cutoff given for occupied orbitals' END IF - IF (mp2_env%ri_g0w0%bse_cutoff_empty > 0) THEN + IF (mp2_env%bse%bse_cutoff_empty > 0) THEN WRITE (unit_nr, '(T2,A4,T7,A26,T71,F10.3)') 'BSE|', 'Cutoff empty orbitals [eV]', & - mp2_env%ri_g0w0%bse_cutoff_empty*evolt + mp2_env%bse%bse_cutoff_empty*evolt ELSE WRITE (unit_nr, '(T2,A4,T7,A34)') 'BSE|', 'No cutoff given for empty orbitals' END IF @@ -1362,7 +1279,7 @@ CONTAINS CALL cp_fm_create(fm_mat_S_ab_trunc, fm_struct_ab, "fm_S_ab_trunc") !Copy parts of original matrices to truncated ones - IF (mp2_env%ri_g0w0%bse_cutoff_occ > 0 .OR. mp2_env%ri_g0w0%bse_cutoff_empty > 0) THEN + IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN !Truncate eigenvals ALLOCATE (Eigenval_reduced(homo_red + virt_red)) Eigenval_reduced(:) = Eigenval(homo_incl:homo + virt_incl) @@ -1403,181 +1320,330 @@ CONTAINS END SUBROUTINE truncate_BSE_matrices ! ************************************************************************************************** -!> \brief Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | µ | ψ_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 dipole_exc BSE dipole vectors in real space per excitation level -!> \param oscill_str Oscillator strength per excitation level -!> \param mp2_env ... -!> \param qs_env ... -!> \param mo_coeff ... -!> \param homo_red ... -!> \param virtual_red ... +!> \brief ... +!> \param fm_eigvec ... +!> \param fm_eigvec_reshuffled ... +!> \param homo ... +!> \param virtual ... +!> \param n_exc ... +!> \param do_transpose ... !> \param unit_nr ... -!> \param fm_eigvec_Y ... +!> \param mp2_env ... ! ************************************************************************************************** - SUBROUTINE get_oscillator_strengths(fm_eigvec_X, Exc_ens, & - dipole_exc, oscill_str, & - mp2_env, qs_env, mo_coeff, & - homo_red, virtual_red, unit_nr, & - fm_eigvec_Y) + SUBROUTINE reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, & + unit_nr, mp2_env) - TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(IN) :: Exc_ens - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & - INTENT(OUT) :: dipole_exc - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(OUT) :: oscill_str - TYPE(mp2_type), INTENT(IN) :: mp2_env - TYPE(qs_environment_type), POINTER :: qs_env - TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff - INTEGER, INTENT(IN) :: homo_red, virtual_red, unit_nr - TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y + TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec + TYPE(cp_fm_type), INTENT(INOUT) :: fm_eigvec_reshuffled + INTEGER, INTENT(IN) :: homo, virtual, n_exc + LOGICAL, INTENT(IN) :: do_transpose + INTEGER, INTENT(IN) :: unit_nr + TYPE(mp2_type), INTENT(INOUT) :: mp2_env - CHARACTER(LEN=*), PARAMETER :: routineN = 'get_oscillator_strengths' + CHARACTER(LEN=*), PARAMETER :: routineN = 'reshuffle_eigvec' - INTEGER :: handle, idir, n, n_occ, n_virt, nao - TYPE(cp_fm_struct_type), POINTER :: fm_struct_dipole_exc, fm_struct_dipoles, & - fm_struct_dipoles_reordered, & - fm_struct_dipoles_trunc - TYPE(cp_fm_type) :: fm_eigvec_XYsum - TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_exc, fm_dipole_per_dir, & - fm_dipole_per_dir_reordered, & - fm_dipoles_per_dir_trunc + INTEGER :: handle, my_m_col, my_n_row + INTEGER, DIMENSION(4) :: reordering + TYPE(cp_fm_struct_type), POINTER :: fm_struct_eigvec_col, & + fm_struct_eigvec_reshuffled + TYPE(cp_fm_type) :: fm_eigvec_col CALL timeset(routineN, handle) - ! Extract dimension from qs_env - nao = qs_env%mos(1)%nao - ! Writing homo to n_occ instead if nmo, - ! takes care of ADDED_MOS, which would overwrite nmo, if invoked - n_occ = qs_env%mos(1)%homo - n_virt = nao - n_occ - ! Original fm_struct_dipoles is invoked with a different blacs env to ensure correct - ! multiplication in get_dipoles_mo using cp_dbcsr_sm_fm_multiply - ! We use general copy afterwards, such that fm_dipoles_per_dir_trunc is already - ! back to BSE blacs envs - CALL cp_fm_struct_create(fm_struct_dipoles, mo_coeff(1)%matrix_struct%para_env, & - mo_coeff(1)%matrix_struct%context, n_virt, n_occ) - CALL cp_fm_struct_create(fm_struct_dipoles_trunc, fm_eigvec_X%matrix_struct%para_env, & - fm_eigvec_X%matrix_struct%context, virtual_red, homo_red) - CALL cp_fm_struct_create(fm_struct_dipoles_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_dipole_exc, fm_eigvec_X%matrix_struct%para_env, & - fm_eigvec_X%matrix_struct%context, 1, homo_red*virtual_red) - - DO idir = 1, 3 - !Create final dipole matrix - CALL cp_fm_create(fm_dipole_per_dir(idir), matrix_struct=fm_struct_dipoles, name="dipoles_mo") - CALL cp_fm_set_all(fm_dipole_per_dir(idir), 0.0_dp) - END DO - - ! Obtain dipoles in MO basis with index structure n_virt x n_occ: \vec{D}_ai = - CALL get_dipoles_mo(fm_dipole_per_dir, qs_env, mo_coeff) - - ! 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 ovver i and a by parallel gemm - DO idir = 1, 3 - CALL cp_fm_create(fm_dipoles_per_dir_trunc(idir), matrix_struct=fm_struct_dipoles_trunc, & - name="dipoles_mo_trunc") - - CALL cp_fm_to_fm_submat_general(fm_dipole_per_dir(idir), & - fm_dipoles_per_dir_trunc(idir), & - virtual_red, & - homo_red, & - 1, & - n_occ - homo_red + 1, & - 1, & - 1, & - fm_struct_dipoles_trunc%context) - CALL cp_fm_create(fm_dipole_per_dir_reordered(idir), matrix_struct=fm_struct_dipoles_reordered, & - name="dipoles_mo_reordered") - CALL cp_fm_set_all(fm_dipole_per_dir_reordered(idir), 0.0_dp) - CALL fm_general_add_bse(fm_dipole_per_dir_reordered(idir), fm_dipoles_per_dir_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_dipole_exc(idir), matrix_struct=fm_struct_dipole_exc, & - name="excitonic_dipoles") - CALL cp_fm_set_all(fm_dipole_exc(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) + ! Define reordering: + ! (ia,11) to (a1,i1) for transposition + ! (ia,11) to (i1,a1) for default + IF (do_transpose) THEN + reordering = (/2, 3, 1, 4/) + my_n_row = virtual + my_m_col = homo + ELSE + reordering = (/1, 3, 2, 4/) + my_n_row = homo + my_m_col = virtual 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_per_dir_reordered(idir), fm_eigvec_XYsum, 0.0_dp, fm_dipole_exc(idir)) - END DO - ! Get oscillator strengths themselves - ALLOCATE (oscill_str(homo_red*virtual_red)) - ALLOCATE (dipole_exc(3, 1, homo_red*virtual_red)) - dipole_exc(:, :, :) = 0.0_dp + CALL cp_fm_struct_create(fm_struct_eigvec_col, & + fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, & + homo*virtual, 1) + CALL cp_fm_struct_create(fm_struct_eigvec_reshuffled, & + fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, & + my_n_row, my_m_col) - ! Sum over all directions - DO idir = 1, 3 - CALL cp_fm_get_submatrix(fm_dipole_exc(idir), dipole_exc(idir, :, :)) - END DO + ! Resort indices + CALL cp_fm_create(fm_eigvec_col, fm_struct_eigvec_col, name="BSE_column_vector") + CALL cp_fm_set_all(fm_eigvec_col, 0.0_dp) + CALL cp_fm_create(fm_eigvec_reshuffled, fm_struct_eigvec_reshuffled, name="BSE_reshuffled_eigenvector") + CALL cp_fm_set_all(fm_eigvec_reshuffled, 0.0_dp) + ! Fill matrix + CALL cp_fm_to_fm_submat(fm_eigvec, fm_eigvec_col, & + homo*virtual, 1, & + 1, n_exc, & + 1, 1) + ! Reshuffle + CALL fm_general_add_bse(fm_eigvec_reshuffled, fm_eigvec_col, 1.0_dp, & + virtual, 1, & + 1, 1, & + unit_nr, reordering, mp2_env) - DO n = 1, homo_red*virtual_red - oscill_str(n) = 2.0_dp/3.0_dp*Exc_ens(n)*SUM(ABS(dipole_exc(:, :, n))**2) - END DO - - CALL cp_fm_struct_release(fm_struct_dipoles) - CALL cp_fm_struct_release(fm_struct_dipoles_reordered) - CALL cp_fm_struct_release(fm_struct_dipole_exc) - CALL cp_fm_struct_release(fm_struct_dipoles_trunc) - DO idir = 1, 3 - CALL cp_fm_release(fm_dipole_per_dir_reordered(idir)) - CALL cp_fm_release(fm_dipole_exc(idir)) - CALL cp_fm_release(fm_dipoles_per_dir_trunc(idir)) - CALL cp_fm_release(fm_dipole_per_dir(idir)) - END DO - CALL cp_fm_release(fm_eigvec_XYsum) + CALL cp_fm_release(fm_eigvec_col) + CALL cp_fm_struct_release(fm_struct_eigvec_col) + CALL cp_fm_struct_release(fm_struct_eigvec_reshuffled) CALL timestop(handle) - END SUBROUTINE get_oscillator_strengths + END SUBROUTINE reshuffle_eigvec ! ************************************************************************************************** -!> \brief Computes and returns dipole in MO basis, i.e. < ψ_i | µ | ψ_a > -!> \param fm_dipole_per_dir Dipoles in MO basis with blacs environment from DFT +!> \brief Borrowed from the tddfpt module with slight adaptions !> \param qs_env ... -!> \param mo_coeff mo_coeff in layout nao x nao from mp2, i.e. all occ and virt indices in both axes +!> \param mos ... +!> \param istate ... +!> \param info_approximation ... +!> \param stride ... +!> \param append_cube ... +!> \param print_section ... ! ************************************************************************************************** - SUBROUTINE get_dipoles_mo(fm_dipole_per_dir, qs_env, mo_coeff) + SUBROUTINE print_bse_nto_cubes(qs_env, mos, istate, info_approximation, & + stride, append_cube, print_section) - TYPE(cp_fm_type), DIMENSION(3), INTENT(INOUT) :: fm_dipole_per_dir + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos + INTEGER, INTENT(IN) :: istate + CHARACTER(LEN=10) :: info_approximation + INTEGER, DIMENSION(:), POINTER :: stride + LOGICAL, INTENT(IN) :: append_cube + TYPE(section_vals_type), POINTER :: print_section + + CHARACTER(LEN=*), PARAMETER :: routineN = 'print_bse_nto_cubes' + + CHARACTER(LEN=default_path_length) :: filename, info_approx_trunc, & + my_pos_cube, title + INTEGER :: handle, i, iset, nmo, unit_nr_cube + LOGICAL :: mpi_io + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(cp_logger_type), POINTER :: logger + TYPE(dft_control_type), POINTER :: dft_control + TYPE(particle_list_type), POINTER :: particles + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(pw_c1d_gs_type) :: wf_g + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools + TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + TYPE(pw_r3d_rs_type) :: wf_r + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(qs_subsys_type), POINTER :: subsys + + logger => cp_get_default_logger() + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env) + CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools) + CALL auxbas_pw_pool%create_pw(wf_r) + CALL auxbas_pw_pool%create_pw(wf_g) + + CALL get_qs_env(qs_env, subsys=subsys) + CALL qs_subsys_get(subsys, particles=particles) + + my_pos_cube = "REWIND" + IF (append_cube) THEN + my_pos_cube = "APPEND" + END IF + + CALL get_qs_env(qs_env=qs_env, & + atomic_kind_set=atomic_kind_set, & + qs_kind_set=qs_kind_set, & + cell=cell, & + particle_set=particle_set) + + DO iset = 1, 2 + CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo) + DO i = 1, nmo + CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, & + cell, dft_control, particle_set, pw_env) + IF (iset == 1) THEN + WRITE (filename, '(A6,I3.3,A5,I2.2,a11)') "_NEXC_", istate, "_NTO_", i, "_Hole_State" + ELSEIF (iset == 2) THEN + WRITE (filename, '(A6,I3.3,A5,I2.2,a15)') "_NEXC_", istate, "_NTO_", i, "_Particle_State" + END IF + info_approx_trunc = TRIM(ADJUSTL(info_approximation)) + info_approx_trunc = info_approx_trunc(2:LEN_TRIM(info_approx_trunc) - 1) + filename = TRIM(info_approx_trunc)//TRIM(filename) + mpi_io = .TRUE. + unit_nr_cube = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", & + middle_name=TRIM(filename), file_position=my_pos_cube, & + log_filename=.FALSE., ignore_should_output=.TRUE., mpi_io=mpi_io) + IF (iset == 1) THEN + WRITE (title, *) "Natural Transition Orbital Hole State", i + ELSEIF (iset == 2) THEN + WRITE (title, *) "Natural Transition Orbital Particle State", i + END IF + CALL cp_pw_to_cube(wf_r, unit_nr_cube, title, particles=particles, stride=stride, mpi_io=mpi_io) + CALL cp_print_key_finished_output(unit_nr_cube, logger, print_section, '', & + ignore_should_output=.TRUE., mpi_io=mpi_io) + END DO + END DO + + CALL auxbas_pw_pool%give_back_pw(wf_g) + CALL auxbas_pw_pool%give_back_pw(wf_r) + + CALL timestop(handle) + END SUBROUTINE print_bse_nto_cubes + +! ************************************************************************************************** +!> \brief Borrowed from the tddfpt module with slight adaptions +!> \param homo ... +!> \param virtual ... +!> \param unit_nr ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE adapt_BSE_input_params(homo, virtual, unit_nr, mp2_env) + + INTEGER, INTENT(IN) :: homo, virtual, unit_nr + TYPE(mp2_type) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'adapt_BSE_input_params' + + INTEGER :: handle, i, j, n, & + num_state_list_exceptions + + CALL timeset(routineN, handle) + + ! Handle negative NUM_PRINT_EXC + IF (mp2_env%bse%num_print_exc < 0 .OR. & + mp2_env%bse%num_print_exc > homo*virtual) THEN + mp2_env%bse%num_print_exc = homo*virtual + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Keyword NUM_PRINT_EXC is either negative or too large. "// & + "Printing all computed excitations.") + END IF + END IF + + ! Default to NUM_PRINT_EXC if too large or negative, + ! but only if NTOs are called - would be confusing for the user otherwise + ! Prepare and adapt user inputs for NTO analysis + ! Logic: Explicit state list overrides NUM_PRINT_EXC_NTOS + ! If only NUM_PRINT_EXC_NTOS is given, we write the array 1,...,NUM_PRINT_EXC_NTOS to + ! bse_nto_state_list + IF (mp2_env%bse%do_nto_analysis) THEN + IF (mp2_env%bse%explicit_nto_list) THEN + IF (mp2_env%bse%num_print_exc_ntos > 0) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Keywords NUM_PRINT_EXC_NTOS and STATE_LIST are both given in input. "// & + "Overriding NUM_PRINT_EXC_NTOS.") + END IF + END IF + ! Check if all states are within the range + ! Count them and initialize new array afterwards + num_state_list_exceptions = 0 + DO i = 1, SIZE(mp2_env%bse%bse_nto_state_list) + IF (mp2_env%bse%bse_nto_state_list(i) < 1 .OR. & + mp2_env%bse%bse_nto_state_list(i) > mp2_env%bse%num_print_exc) THEN + num_state_list_exceptions = num_state_list_exceptions + 1 + END IF + END DO + IF (num_state_list_exceptions > 0) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "STATE_LIST contains indices outside the range of included excitation levels. "// & + "Ignoring these states.") + END IF + END IF + n = SIZE(mp2_env%bse%bse_nto_state_list) - num_state_list_exceptions + ALLOCATE (mp2_env%bse%bse_nto_state_list_final(n)) + mp2_env%bse%bse_nto_state_list_final(:) = 0 + i = 1 + DO j = 1, SIZE(mp2_env%bse%bse_nto_state_list) + IF (mp2_env%bse%bse_nto_state_list(j) >= 1 .AND. & + mp2_env%bse%bse_nto_state_list(j) <= mp2_env%bse%num_print_exc) THEN + mp2_env%bse%bse_nto_state_list_final(i) = mp2_env%bse%bse_nto_state_list(j) + i = i + 1 + END IF + END DO + + mp2_env%bse%num_print_exc_ntos = SIZE(mp2_env%bse%bse_nto_state_list_final) + ELSE + IF (mp2_env%bse%num_print_exc_ntos > mp2_env%bse%num_print_exc .OR. & + mp2_env%bse%num_print_exc_ntos < 0) THEN + mp2_env%bse%num_print_exc_ntos = mp2_env%bse%num_print_exc + END IF + ALLOCATE (mp2_env%bse%bse_nto_state_list_final(mp2_env%bse%num_print_exc_ntos)) + DO i = 1, mp2_env%bse%num_print_exc_ntos + mp2_env%bse%bse_nto_state_list_final(i) = i + END DO + END IF + END IF + + ! Takes care of triplet states, when oscillator strengths are 0 + IF (mp2_env%bse%bse_spin_config /= 0 .AND. & + mp2_env%bse%eps_nto_osc_str > 0) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Cannot apply EPS_OSC_STR for Triplet excitations. "// & + "Resetting EPS_OSC_STR to default.") + END IF + mp2_env%bse%eps_nto_osc_str = -1.0_dp + END IF + + ! Take care of number for computed exciton descriptors + IF (mp2_env%bse%num_print_exc_descr < 0 .OR. & + mp2_env%bse%num_print_exc_descr > mp2_env%bse%num_print_exc) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Keyword NUM_PRINT_EXC_DESCR is either negative or too large. "// & + "Printing exciton descriptors up to NUM_PRINT_EXC.") + END IF + mp2_env%bse%num_print_exc_descr = mp2_env%bse%num_print_exc + END IF + + CALL timestop(handle) + END SUBROUTINE adapt_BSE_input_params + +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... +!> \param fm_multipole_ai_trunc ... +!> \param fm_multipole_ij_trunc ... +!> \param fm_multipole_ab_trunc ... +!> \param qs_env ... +!> \param mo_coeff ... +!> \param rpoint ... +!> \param n_moments ... +!> \param homo_red ... +!> \param virtual_red ... +!> \param context_BSE ... +! ************************************************************************************************** + SUBROUTINE get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, & + qs_env, mo_coeff, rpoint, n_moments, & + homo_red, virtual_red, context_BSE) + + TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), & + INTENT(INOUT) :: fm_multipole_ai_trunc, & + fm_multipole_ij_trunc, & + fm_multipole_ab_trunc TYPE(qs_environment_type), POINTER :: qs_env TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff + REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: rpoint + INTEGER, INTENT(IN) :: n_moments, homo_red, virtual_red + TYPE(cp_blacs_env_type), POINTER :: context_BSE - CHARACTER(LEN=*), PARAMETER :: routineN = 'get_dipoles_mo' + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_multipoles_mo' - INTEGER :: handle, idir, n_occ, n_virt, nao - REAL(dp), ALLOCATABLE, DIMENSION(:) :: rpoint + INTEGER :: handle, idir, n_multipole, n_occ, & + n_virt, nao REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point - TYPE(cp_fm_struct_type), POINTER :: fm_struct_coeff_occ, & - fm_struct_coeff_virt, & - fm_struct_dipoles_ao, fm_struct_nao_nao - TYPE(cp_fm_type) :: fm_coeff_occ, fm_coeff_virt, & - fm_dip_intermed - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_r, matrix_s + TYPE(cp_fm_struct_type), POINTER :: fm_struct_mp_ab_trunc, fm_struct_mp_ai_trunc, & + fm_struct_mp_ij_trunc, fm_struct_multipoles_ao, fm_struct_nao_nao + TYPE(cp_fm_type) :: fm_work + TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_multipole_per_dir + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_multipole, matrix_s TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(mp_para_env_type), POINTER :: para_env_BSE TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb @@ -1591,23 +1657,26 @@ CONTAINS sab_orb=sab_orb) ! Use the same blacs environment as for the MO coefficients to ensure correct multiplication dbcsr x fm later on - fm_struct_dipoles_ao => mos(1)%mo_coeff%matrix_struct - - NULLIFY (matrix_r) - CALL dbcsr_allocate_matrix_set(matrix_r, 3) - DO idir = 1, 3 - CALL dbcsr_init_p(matrix_r(idir)%matrix) - CALL dbcsr_create(matrix_r(idir)%matrix, name="DIPOLE MATRIX X/Y/Z", & + fm_struct_multipoles_ao => mos(1)%mo_coeff%matrix_struct + ! BSE has different contexts and blacsenvs + para_env_BSE => context_BSE%para_env + ! Get size of multipole tensor + n_multipole = (6 + 11*n_moments + 6*n_moments**2 + n_moments**3)/6 - 1 + NULLIFY (matrix_multipole) + CALL dbcsr_allocate_matrix_set(matrix_multipole, n_multipole) + ALLOCATE (fm_multipole_per_dir(n_multipole)) + DO idir = 1, n_multipole + CALL dbcsr_init_p(matrix_multipole(idir)%matrix) + CALL dbcsr_create(matrix_multipole(idir)%matrix, name="ao_multipole", & template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric) - CALL cp_dbcsr_alloc_block_from_nbl(matrix_r(idir)%matrix, sab_orb) - CALL dbcsr_set(matrix_r(idir)%matrix, 0._dp) + CALL cp_dbcsr_alloc_block_from_nbl(matrix_multipole(idir)%matrix, sab_orb) + CALL dbcsr_set(matrix_multipole(idir)%matrix, 0._dp) END DO - ALLOCATE (rpoint(3)) CALL get_reference_point(rpoint, qs_env=qs_env, reference=use_mom_ref_coac, ref_point=ref_point) - CALL build_local_moment_matrix(qs_env, matrix_r, 1, ref_point=rpoint) - DEALLOCATE (rpoint) + CALL build_local_moment_matrix(qs_env, matrix_multipole, n_moments, ref_point=rpoint) + NULLIFY (sab_orb) ! Now we transform them to MO @@ -1617,218 +1686,143 @@ CONTAINS CALL get_mo_set(mo_set=mos(1), homo=n_occ, nao=nao) n_virt = nao - n_occ - ! At the end, we need four different layouts of matrices in this multiplication: - ! D = dipoles - ! Final result: D = C_{mu a} <\mu|\vec{r}|\nu> C_{\nu i} EQ.I + ! At the end, we need four different layouts of matrices in this multiplication, e.g. for a dipole: + ! D_pq = full multipole matrix for occupied and unoccupied + ! Final result:D_pq= C_{mu p} <\mu|\vec{r}|\nu> C_{\nu q} EQ.I ! \_______/ \___________/ \______/ - ! fm_coeff_virt matrix_r fm_coeff_occ + ! fm_coeff matrix_multipole fm_coeff ! (EQ.Ia) (EQ.Ib) (EQ.Ia) ! Intermediate work matrices: - ! fm_dip_intermed = <\mu|\vec{r}|\nu> C_{\nu i} EQ.II + ! fm_work = <\mu|\vec{r}|\nu> C_{\nu q} EQ.II + + ! Struct for the full multipole matrix CALL cp_fm_struct_create(fm_struct_nao_nao, & - fm_struct_dipoles_ao%para_env, fm_struct_dipoles_ao%context, & + fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, & nao, nao) - CALL cp_fm_struct_create(fm_struct_coeff_virt, fm_struct_dipoles_ao%para_env, fm_struct_dipoles_ao%context, & - nao, n_virt) - CALL cp_fm_struct_create(fm_struct_coeff_occ, fm_struct_dipoles_ao%para_env, fm_struct_dipoles_ao%context, & - nao, n_occ) - ! Split between unoccupied and occupied MO coefficients - ! First, the unoccupied coefficients "C_\mu a" - CALL cp_fm_create(fm_coeff_virt, matrix_struct=fm_struct_coeff_virt, name="mo_coeff_virt") - CALL cp_fm_set_all(fm_coeff_virt, 0.0_dp) - ! Coefficients from Diagonalization of the KS matrix, starting from the first unoccupied index, i.e. n_occ+1 - ! We are writing between different contexts here, so we need submat_general - CALL cp_fm_to_fm_submat_general(mo_coeff(1), fm_coeff_virt, & - nao, n_virt, & - 1, n_occ + 1, & - 1, 1, & - mo_coeff(1)%matrix_struct%context) - - ! Now occupied ones "C_\mu i" - CALL cp_fm_create(fm_coeff_occ, matrix_struct=fm_struct_coeff_occ, name="mo_coeff_occ") - CALL cp_fm_set_all(fm_coeff_occ, 0.0_dp) - - ! Coefficients from Diagonalization of the KS matrix, starting from the first occupied index - CALL cp_fm_to_fm_submat_general(mo_coeff(1), fm_coeff_occ, & - nao, n_occ, & - 1, 1, & - 1, 1, & - mo_coeff(1)%matrix_struct%context) + ! At the very end, we copy the multipoles corresponding to truncated BSE indices in i and a + CALL cp_fm_struct_create(fm_struct_mp_ai_trunc, para_env_BSE, & + context_BSE, virtual_red, homo_red) + CALL cp_fm_struct_create(fm_struct_mp_ij_trunc, para_env_BSE, & + context_BSE, homo_red, homo_red) + CALL cp_fm_struct_create(fm_struct_mp_ab_trunc, para_env_BSE, & + context_BSE, virtual_red, virtual_red) + DO idir = 1, n_multipole + CALL cp_fm_create(fm_multipole_ai_trunc(idir), matrix_struct=fm_struct_mp_ai_trunc, & + name="dipoles_mo_ai_trunc") + CALL cp_fm_set_all(fm_multipole_ai_trunc(idir), 0.0_dp) + CALL cp_fm_create(fm_multipole_ij_trunc(idir), matrix_struct=fm_struct_mp_ij_trunc, & + name="dipoles_mo_ij_trunc") + CALL cp_fm_set_all(fm_multipole_ij_trunc(idir), 0.0_dp) + CALL cp_fm_create(fm_multipole_ab_trunc(idir), matrix_struct=fm_struct_mp_ab_trunc, & + name="dipoles_mo_ab_trunc") + CALL cp_fm_set_all(fm_multipole_ab_trunc(idir), 0.0_dp) + END DO ! Need another temporary matrix to store intermediate result from right multiplication ! D = C_{mu a} <\mu|\vec{r}|\nu> C_{\nu i} - CALL cp_fm_create(fm_dip_intermed, matrix_struct=fm_struct_coeff_occ, name="dip_intermed") - CALL cp_fm_set_all(fm_dip_intermed, 0.0_dp) + CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nao, name="multipole_work") + CALL cp_fm_set_all(fm_work, 0.0_dp) - DO idir = 1, 3 - ! Fill final (MO) dipole matrix - ! Specifying beta=1.0_dp would sum up all directions - CALL cp_dbcsr_sm_fm_multiply(matrix_r(idir)%matrix, fm_coeff_occ, & - fm_dip_intermed, ncol=n_occ) - ! Now obtain the dipoles by the final multiplication; - ! We do that inside the loop to obtain dipoles per axis for print - CALL parallel_gemm('T', 'N', n_virt, n_occ, nao, 1.0_dp, fm_coeff_virt, fm_dip_intermed, 0.0_dp, fm_dipole_per_dir(idir)) + DO idir = 1, n_multipole + ! Create the full multipole matrix per direction + CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nao_nao, name="multipoles_mo") + CALL cp_fm_set_all(fm_multipole_per_dir(idir), 0.0_dp) + ! Fill final (MO) multipole matrix + CALL cp_dbcsr_sm_fm_multiply(matrix_multipole(idir)%matrix, mo_coeff(1), & + fm_work, ncol=nao) + ! Now obtain the multipoles by the final multiplication; + ! We do that inside the loop to obtain multipoles per axis for print + CALL parallel_gemm('T', 'N', nao, nao, nao, 1.0_dp, mo_coeff(1), fm_work, 0.0_dp, fm_multipole_per_dir(idir)) + ! Truncate full matrix to the BSE indices + ! D_ai + CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), & + fm_multipole_ai_trunc(idir), & + virtual_red, & + homo_red, & + n_occ + 1, & + n_occ - homo_red + 1, & + 1, & + 1, & + fm_multipole_per_dir(idir)%matrix_struct%context) + ! D_ij + CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), & + fm_multipole_ij_trunc(idir), & + homo_red, & + homo_red, & + n_occ - homo_red + 1, & + n_occ - homo_red + 1, & + 1, & + 1, & + fm_multipole_per_dir(idir)%matrix_struct%context) + ! D_ab + CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), & + fm_multipole_ab_trunc(idir), & + virtual_red, & + virtual_red, & + n_occ + 1, & + n_occ + 1, & + 1, & + 1, & + fm_multipole_per_dir(idir)%matrix_struct%context) END DO !Release matrices and structs - NULLIFY (fm_struct_dipoles_ao) + NULLIFY (fm_struct_multipoles_ao) + CALL cp_fm_struct_release(fm_struct_mp_ai_trunc) + CALL cp_fm_struct_release(fm_struct_mp_ij_trunc) + CALL cp_fm_struct_release(fm_struct_mp_ab_trunc) CALL cp_fm_struct_release(fm_struct_nao_nao) - CALL cp_fm_struct_release(fm_struct_coeff_virt) - CALL cp_fm_struct_release(fm_struct_coeff_occ) - CALL cp_fm_release(fm_coeff_virt) - CALL cp_fm_release(fm_coeff_occ) - CALL cp_fm_release(fm_dip_intermed) - CALL dbcsr_deallocate_matrix_set(matrix_r) - - CALL timestop(handle) - - END SUBROUTINE get_dipoles_mo - -! ************************************************************************************************** -!> \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 Exc_ens ... -!> \param info_approximation ... -!> \param unit_nr ... -!> \param mp2_env ... -! ************************************************************************************************** - SUBROUTINE compute_absorption_spectrum(oscill_str, Exc_ens, & - info_approximation, unit_nr, mp2_env) - - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(IN) :: oscill_str, 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 - CHARACTER(LEN=28) :: file_name_spectrum - INTEGER :: handle, i, j, k, num_steps, unit_nr_file - 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 - - CALL timeset(routineN, handle) - - freq_step = mp2_env%ri_g0w0%bse_spectrum_freq_step_size - freq_start = mp2_env%ri_g0w0%bse_spectrum_freq_start - freq_end = mp2_env%ri_g0w0%bse_spectrum_freq_end - eta_list => mp2_env%ri_g0w0%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) - ! Convert to str - WRITE (eta_str, '(F6.3)') eta - file_name_spectrum = 'BSE'//TRIM(ADJUSTL(info_approximation))//'eta='//TRIM(eta_str)//'.spectrum' - ALLOCATE (abs_spectrum(num_steps, 2)) - 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 - 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)) - END DO - END DO - ! Print it to file - unit_nr_file = cp_logger_get_default_unit_nr() - 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, '(A75,A6)') "# Imaginary part of polarizability from Bethe Salpeter equation for method ", & - TRIM(ADJUSTL(info_approximation)) - DO i = 1, num_steps - WRITE (unit_nr_file, '(F17.10, F17.10)') abs_spectrum(i, 1)*evolt, abs_spectrum(i, 2) - END DO - CALL close_file(unit_nr_file) - END IF - DEALLOCATE (abs_spectrum) + DO idir = 1, n_multipole + CALL cp_fm_release(fm_multipole_per_dir(idir)) 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,A59)') & - 'BSE|', "with oscillator strengths f_n and excitation energies Ω^n." - WRITE (unit_nr, '(T2,A4)') 'BSE|' - END IF + DEALLOCATE (fm_multipole_per_dir) + CALL cp_fm_release(fm_work) + CALL dbcsr_deallocate_matrix_set(matrix_multipole) CALL timestop(handle) - END SUBROUTINE + END SUBROUTINE get_multipoles_mo - ! ************************************************************************************************** - ! ************************************************************************************************** ! ************************************************************************************************** -!> \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 ... +!> \brief Computes trace of form Tr{A^T B C} for exciton descriptors +!> \param fm_A Full Matrix, typically X or Y, in format homo x virtual +!> \param fm_B ... +!> \param fm_C ... +!> \param alpha ... ! ************************************************************************************************** - SUBROUTINE print_transition_probabilities(fm_eigvec, direction_excitation, info_approximation, & - i_exc, virtual, homo, homo_irred, & - unit_nr, mp2_env) + SUBROUTINE trace_exciton_descr(fm_A, fm_B, fm_C, alpha) - 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 + TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B, fm_C + REAL(KIND=dp), INTENT(OUT) :: alpha - CHARACTER(LEN=*), PARAMETER :: routineN = 'print_transition_probabilities' + CHARACTER(LEN=*), PARAMETER :: routineN = 'trace_exciton_descr' - INTEGER :: handle, k, num_entries - INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries + INTEGER :: handle, ncol_A, ncol_B, ncol_C, nrow_A, & + nrow_B, nrow_C + TYPE(cp_fm_type) :: fm_work_ia 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 cp_fm_create(fm_work_ia, fm_A%matrix_struct) + CALL cp_fm_get_info(fm_A, nrow_global=nrow_A, ncol_global=ncol_A) + CALL cp_fm_get_info(fm_B, nrow_global=nrow_B, ncol_global=ncol_B) + CALL cp_fm_get_info(fm_C, nrow_global=nrow_C, ncol_global=ncol_C) + + ! Check matrix sizes + CPASSERT(nrow_A == nrow_B .AND. ncol_A == ncol_C .AND. ncol_B == nrow_C) + + CALL cp_fm_set_all(fm_work_ia, 0.0_dp) + + CALL parallel_gemm("N", "N", nrow_A, ncol_A, nrow_C, 1.0_dp, & + fm_B, fm_C, 0.0_dp, fm_work_ia) + + CALL cp_fm_trace(fm_A, fm_work_ia, alpha) + + CALL cp_fm_release(fm_work_ia) + CALL timestop(handle) - END SUBROUTINE + END SUBROUTINE trace_exciton_descr + END MODULE diff --git a/src/common/bibliography.F b/src/common/bibliography.F index c98793772e..09b6b60eaf 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -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 diff --git a/src/fm/cp_fm_diag.F b/src/fm/cp_fm_diag.F index f4e9b6dbe8..b9e0417257 100644 --- a/src/fm/cp_fm_diag.F +++ b/src/fm/cp_fm_diag.F @@ -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 ... diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index 2398c55ce6..fc7b0e9c53 100644 --- a/src/input_cp2k_mp2.F +++ b/src/input_cp2k_mp2.F @@ -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) diff --git a/src/mp2_gpw.F b/src/mp2_gpw.F index e52f24001e..272acba4bb 100644 --- a/src/mp2_gpw.F +++ b/src/mp2_gpw.F @@ -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 diff --git a/src/mp2_setup.F b/src/mp2_setup.F index 9cb85265ac..b3a708453f 100644 --- a/src/mp2_setup.F +++ b/src/mp2_setup.F @@ -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) diff --git a/src/mp2_types.F b/src/mp2_types.F index 1550ca1782..f355eede2e 100644 --- a/src/mp2_types.F +++ b/src/mp2_types.F @@ -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, & diff --git a/src/rpa_gw.F b/src/rpa_gw.F index 27f49d5faf..8b75280373 100644 --- a/src/rpa_gw.F +++ b/src/rpa_gw.F @@ -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 diff --git a/src/rpa_gw_sigma_x.F b/src/rpa_gw_sigma_x.F index 6934775c70..a59bf9aacb 100644 --- a/src/rpa_gw_sigma_x.F +++ b/src/rpa_gw_sigma_x.F @@ -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 diff --git a/src/rpa_main.F b/src/rpa_main.F index a09a614bba..d2b51594d6 100644 --- a/src/rpa_main.F +++ b/src/rpa_main.F @@ -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) diff --git a/tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_EXC_DESCR.inp b/tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_EXC_DESCR.inp new file mode 100644 index 0000000000..119539c580 --- /dev/null +++ b/tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_EXC_DESCR.inp @@ -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 diff --git a/tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_NTOs.inp b/tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_NTOs.inp new file mode 100644 index 0000000000..3e6d056bcc --- /dev/null +++ b/tests/QS/regtest-bse/BSE_H2O_PBE_G0W0_NTOs.inp @@ -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 diff --git a/tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_EXC_DESCR.inp b/tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_EXC_DESCR.inp new file mode 100644 index 0000000000..e4101bc2a7 --- /dev/null +++ b/tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_EXC_DESCR.inp @@ -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 diff --git a/tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_NTOs.inp b/tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_NTOs.inp new file mode 100644 index 0000000000..9a61675145 --- /dev/null +++ b/tests/QS/regtest-bse/TDA_H2O_PBE_G0W0_NTOs.inp @@ -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 diff --git a/tests/QS/regtest-bse/TEST_FILES b/tests/QS/regtest-bse/TEST_FILES index abc0333233..2ccc430922 100644 --- a/tests/QS/regtest-bse/TEST_FILES +++ b/tests/QS/regtest-bse/TEST_FILES @@ -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 diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 5cff120bfb..f1d7d6a972 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -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