Skip to content

Commit

Permalink
refs are converted to abb
Browse files Browse the repository at this point in the history
  • Loading branch information
aziziph committed Aug 22, 2023
1 parent b5f27ee commit 6b14f03
Show file tree
Hide file tree
Showing 4 changed files with 10 additions and 13 deletions.
Binary file modified JOSS/flowchart.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
17 changes: 8 additions & 9 deletions JOSS/paper.md
Original file line number Diff line number Diff line change
Expand Up @@ -65,19 +65,19 @@ The package comprises minimax time and frequency grids [@Takatsuka2008; @kaltak2

# Statement of need

RPA is an accurate approach to compute the electronic correlation energy. It is non-local, includes long-range dispersion interactions and dynamic electronic screening, and is applicable to a wide range of systems from 0 to 3 dimensions [@eshuis2012electron; @ren2012random]. The \textit{GW} method [@hedin1965new] is based on the RPA susceptibility and has become the method of choice for the calculation of direct and inverse photoemission spectra of molecules and solids [@golze2019gw;@reining2018gw;@stankovski2011g;@li2005quasiparticle;@bruneval2008accurate;@nabok2016accurate]. \textit{GW} forms the basis for Bethe-Salpeter Equation (BSE) calculations of optical spectra [@onida2002electronic], where the \textit{GW} results are used as input.
RPA is an accurate approach to compute the electronic correlation energy. It is non-local, includes long-range dispersion interactions and dynamic electronic screening, and is applicable to a wide range of systems from 0 to 3 dimensions [@eshuis2012electron; @ren2012random]. The \textit{GW} method [@hedin1965new] is based on the RPA susceptibility and has become the method of choice for the calculation of direct and inverse photoemission spectra of molecules and solids [@golze2019gw;@reining2018gw]. \textit{GW} forms the basis for Bethe-Salpeter Equation (BSE) calculations of optical spectra [@onida2002electronic], where the \textit{GW} results are used as input.

Despite their wide adoption, RPA and \textit{GW} face computational challenges, especially for large systems. Conventional RPA and \textit{GW} implementations scale with the fourth power of the system size $N$ and are therefore usually limited to systems of a few tens to at most a hundred atoms [@delben2013electron; @wilhelm2016gw; @stuke2020atomic]. To tackle larger and more realistic systems, scaling reductions present a promising strategy to decrease the computational cost. Such low-scaling algorithms utilize real-space representations and time-frequency transformations, such as the real-space/imaginary-time approach [@rojas1995space] that reduces the complexity to $\mathcal{O}(N^3)$. Several such cubic-scaling \textit{GW} algorithms have recently been implemented, e.g. in a plane-wave/projector-augmented-wave (PAW) \textit{GW} code [@kaltak2014low;@liu2016cubic] or with localized basis sets using Gaussian [@wilhelm2018toward;@wilhelm2021low;@duchemin2021cubic;@Graml2023] or Slater-type orbitals [@forster2020low;@foerster2021GW100;@foerster2021loworder;@foerster2023twocomponent]. Similarly, low-scaling RPA algorithms were implemented with different basis sets [@kaltak2014cubic;@kaltak2014low;@wilhelm2016rpa;@luenser2017;@graf2018accurate;@duchemin2019separable;@drontschenko2022efficient;@kutepov2012electronic;@kutepov2017linearized].
Despite their wide adoption, RPA and \textit{GW} face computational challenges, especially for large systems. Conventional RPA and \textit{GW} implementations scale with the fourth power of the system size $N$ and are therefore usually limited to systems of a few tens to at most a hundred atoms [@wilhelm2016gw; @stuke2020atomic]. To tackle larger and more realistic systems, scaling reductions present a promising strategy to decrease the computational cost. Such low-scaling algorithms utilize real-space representations and time-frequency transformations, such as the real-space/imaginary-time approach [@rojas1995space] that reduces the complexity to $\mathcal{O}(N^3)$. Several such cubic-scaling \textit{GW} algorithms have recently been implemented, e.g. in a plane-wave/projector-augmented-wave (PAW) \textit{GW} code [@kaltak2014low;@liu2016cubic,@kutepov2017linearized] or with localized basis sets using Gaussian [@wilhelm2018toward;@wilhelm2021low;@duchemin2021cubic;@Graml2023] or Slater-type orbitals [@forster2020low;@foerster2021GW100;@foerster2021loworder;@foerster2023twocomponent]. Similarly, low-scaling RPA algorithms were implemented with different basis sets [@kaltak2014cubic;@kaltak2014low;@wilhelm2016rpa;@graf2018accurate;@duchemin2019separable;@drontschenko2022efficient].

An important consideration for low-scaling algorithms is the crossover point. Due to their larger pre-factor, low-scaling algorithms are typically more expensive for smaller systems and only become more cost-effective than canonical implementations for larger systems due to their reduced scaling [@wilhelm2018toward]. Furthermore, the numerical precision of low-scaling \textit{GW} algorithms is strongly coupled to the time-frequency treatment [@wilhelm2021low]. Early low-scaling \textit{GW} algorithms did not reach the same precision as canonical implementations [@vlcek2017stochastic; @wilhelm2018toward; @forster2020low]. Although appropriate Fourier transforms and corresponding time-frequency grids have been implemented [@liu2016cubic;@wilhelm2021low; @duchemin2021cubic; @foerster2021GW100], these implementations and grids are tied to particular codes and are often buried deeply inside the code. Furthermore, reuse of such implementations elsewhere is often restricted by license requirements or dependencies on definitions made in the host code.

In this work, we present the GX-TimeFrequency component of the GreenX library, an open-source package distributed under the Apache license (Version 2.0). GX-TimeFrequency provides time and frequency grids and the corresponding integration weights to compute correlation energies for Green's function implementations. It also provides Fourier weights to convert between imaginary time and imaginary frequency. The library can be used for low-scaling RPA and \textit{GW} implementations. BSE codes, which use (low-scaling) \textit{GW} as input, can also utilize our library. The minimax grids are also suitable for RPA implementations with conventional scaling [@delben2015enabling]: The minimax grids are more compact than, e.g., Gauss-Legendre grids, resulting in a reduction of the computational prefactor, while yielding the same accuracy [@delben2015enabling;@azizi_minimax]. We note that the minimax grids are not recommended for conventional imaginary-frequency-only \textit{GW} implementations [@ren2012resolution;@wilhelm2016gw] since they have not been optimized for the frequency integral of the self-energy.

While not being the main target of the library, the minimax time grids can also be utilized to calculate the LT-dMP2 correlation energy [@almloef1991elimination; @jung2004scaled; @Takatsuka2008; @kaltak2014low; @glasbrenner2020efficient]. The dMP2 term is one of two terms of the MP2 correlation energy. In a diagrammatic representation, dMP2 is the lowest order of the RPA correlation energy [@ren2012random]. The dMP2 correlation energy can be reformulated using the Laplace transform to obtain the LT-dMP2 expression which scales cubically in contrast to the $O(N^5)$ scaling of standard MP2.
While not being the main target of the library, the minimax time grids can also be utilized to calculate the LT-dMP2 correlation energy [@almloef1991elimination; @jung2004scaled; @Takatsuka2008; @kaltak2014low; @glasbrenner2020efficient]. The dMP2 term is one of two terms of the MP2 correlation energy and, in a diagrammatic representation, corresponds to the lowest order of the RPA correlation energy [@ren2012random]. The dMP2 correlation energy can be reformulated using the Laplace transform to obtain the LT-dMP2 expression which scales cubically in contrast to the $O(N^5)$ scaling of standard MP2.

# Mathematical framework

The single-particle Green's function $G$ and the non-interacting susceptibility $\chi^0$ are the starting point for a set of many-body perturbation theory (MBPT) methods. In canonical implementations, the non-interacting susceptibility $\chi^0(i\omega)$ is often expressed in the Adler-Wiser form [@Adler1962;@Wiser1963], where the sums over occupied (index $j$) and unoccupied (index $a$) single-particle states $\psi$ are coupled via their corresponding energies $\varepsilon$.
The single-particle Green's function $G$ and the non-interacting susceptibility $\chi^0$ are the starting point for a set of many-body perturbation theory (MBPT) methods. In canonical implementations, $\chi^0(\mathbf{r},\mathbf{r'},i\omega)$ is often expressed in the Adler-Wiser form [@Adler1962;@Wiser1963], where the sums over occupied (index $j$) and unoccupied (index $a$) single-particle states $\psi$ are coupled via their corresponding energies $\varepsilon$.

![Sketch of the methods supported by GX-TimeFrequency which start from $\hat{\chi}^0(i\tau)$. In addition to the discrete time and frequency grids $\{\tau_{j}\}$ and $\{\omega_{k}\}$, the library provides the corresponding weights $\{\sigma_{j}\}$ and $\{\gamma_{k}\}$ for the integration of the correlation energy $E_c$ as well as the Fourier weights $\delta_{kj}$, $\eta_{jk}$ and $\lambda_{kj}$ defined in Eqs. \eqref{ct_st_even}-\eqref{st_odd_t_to_w}. The bare and screened Coulomb interactions are indicated by $v(\mathbf{r},\mathbf{r}')=1/|\mathbf{r}-\mathbf{r}|'$ and $W(i\omega)$, respectively. $\epsilon(i\omega)$ is the dynamical dielectric function, $\Sigma$ the \textit{GW} self-energy, and AC stands for analytic continuation.\label{fig:flowchart}](flowchart.png)

Expand All @@ -90,12 +90,11 @@ Low-scaling RPA and \textit{GW} algorithms include the Fourier transform of $\ha
To convert between imaginary time and frequency grids, we introduce a nonuniform discrete cosine and sine transformation for even and odd functions $F^\text{even/odd}$, respectively [@liu2016cubic]. If the function $F$ is neither odd nor even we split the computation of functions $\hat{F}(i\tau)$ and $F(i\omega)$ into even and an odd parts [@liu2016cubic]

\begin{align}
\hat{F}(i\tau) &= \hat{F}^\text{even}(i\tau) + \hat{F}^\text{odd}(i\tau)
\\
{F}(i\omega) &= {F}^\text{even}(i\omega) + {F}^\text{odd}(i\omega)\label{Fw_split}
\hat{F}(i\tau) = \hat{F}^\text{even}(i\tau) + \hat{F}^\text{odd}(i\tau) \hspace{1.9em} &\text{and} \hspace{1.9em}
{F}(i\omega) = {F}^\text{even}(i\omega) + {F}^\text{odd}(i\omega)\,,\label{Fw_split}
\end{align}
with
\begin{align}
\begin{align}/
\hat{F}^\text{even}(i\tau)=\hat{F}^\text{even}(-i\tau) \hspace{2.4em} &\text{and} \hspace{1.9em} F^\text{even}(i\omega)=F^\text{even}(-i\omega)\,, \\
\hat{F}^\text{odd}(i\tau)=-\hat{F}^\text{odd}(-i\tau) \hspace{2em} &\text{and} \hspace{2em} F^\text{odd}(i\omega)=-F^\text{odd}(-i\omega)\,.
\end{align}
Expand All @@ -112,7 +111,7 @@ The corresponding discrete Fourier transforms read [@liu2016cubic]
\end{align}
where $\{\tau_j\}_{j=1}^n, \tau_j\,{>}\,0$ are again the time grid points, $\{\omega_k\}_{k=1}^n,\omega_k\,{>}\,0$ frequency grid points and $\{\delta_{kj}\}_{k,j=1}^n$, $\{\eta_{jk}\}_{k,j=1}^n$,$\{\lambda_{kj}\}_{k,j=1}^n$, $\{\zeta_{jk}\}_{k,j=1}^n$ the corresponding Fourier integration weights. $\hat{\chi}^0(i\tau)$ is an even function and we need thus the transform defined in Eq. \eqref{ct_st_even} to obtain $\chi^0(i\omega)$. The screened Coulomb interaction is also even and we use expression \eqref{ct_even_w_to_t} to convert $W(i\omega)$ to $\widehat{W}(i\tau)$. The self-energy is neither odd nor even and we use Eq. \eqref{Fw_split} in combination with Eqs. \eqref{ct_st_even} and  \eqref{st_odd_t_to_w} to transform $\widehat{\Sigma}(i\tau)$ to $\Sigma(i\omega)$ [@liu2016cubic]. The transformation defined in Eq. \eqref{ct_st_odd} is not required for the methods summarized in Figure \ref{fig:flowchart}, but is added for the sake of completeness and clarity.

Ideal grid parameters $\tau_j$, $\sigma_j$, $\omega_k$, $\gamma_k$, $\delta_{kj}$, $\eta_{jk}$, $\lambda_{kj}$ feature a vanishing error for the LT-dMP2 and RPA correlation energy integrations and Fourier transforms of $\chi^0,W$ and $\Sigma$ (Figure \ref{fig:flowchart}). We compute minimax grid parameters $\tau_j$, $\sigma_j$, $\omega_k$, $\gamma_k$ that minimize the maximum error of the LT-dMP2 and RPA correlation energy integration (Fig. \ref{fig:flowchart}) over all possible functions $\hat{\chi}^0( \mathbf{r}, \mathbf{r'}, i\tau )$ and $\chi^0( \mathbf{r}, \mathbf{r'}, i\omega )$
Ideal grid parameters $\tau_j$, $\sigma_j$, $\omega_k$, $\gamma_k$, $\delta_{kj}$, $\eta_{jk}$, $\lambda_{kj}$ feature a vanishing error for the LT-dMP2 and RPA correlation energy integrations and Fourier transforms of $\chi^0,W$ and $\Sigma$ (Figure \ref{fig:flowchart}). We compute minimax grid parameters $\tau_j$, $\sigma_j$, $\omega_k$, $\gamma_k$ that minimize the maximum error of the LT-dMP2 and RPA correlation energy integration (Figure \ref{fig:flowchart}) over all possible functions $\hat{\chi}^0( \mathbf{r}, \mathbf{r'}, i\tau )$ and $\chi^0( \mathbf{r}, \mathbf{r'}, i\omega )$
[@Takatsuka2008;@kaltak2014low;@liu2016cubic]. For this minimax grid optimization, we use the Remez algorithm [@kaltak2014low] which is an iterative, numerically ill-conditioned procedure requiring high numerical precision. As the generation of the minimax parameters $\tau_j$, $\sigma_j$, $\omega_k$, $\gamma_k$ is tedious, the most feasible strategy is to tabulate the computed minimax parameters $\{\tau_j\}_{j=1}^n$, $\{\sigma_j\}_{j=1}^n$, $\{\omega_k\}_{k=1}^n$, $\{\gamma_k\}_{k=1}^n$ for their later use in LT-dMP2, RPA, and \textit{GW} calculations.

It has been shown that minimax time and frequency grids $\{\tau_j\}_{j=1}^n$, $\{\omega_k\}_{k=1}^n$ are also suitable for performing Fourier transforms of $\chi^0,W$ and $\Sigma$ [@liu2016cubic]. With the knowledge of the tabulated $\{\tau_j\}_{j=1}^n$ and $\{\omega_k\}_{k=1}^n$ parameters, it is possible to use least-squares optimization to calculate the Fourier integration weights $\delta_{kj}$, $\eta_{jk}$, $\lambda_{kj}$ [@azizi_minimax]. The least-squares optimization can be executed by simple non-iterative linear matrix algebra which is straightforward and is done during the run time of the GreenX library [@azizi_minimax].
Expand Down
4 changes: 1 addition & 3 deletions JOSS/refs.bib
Original file line number Diff line number Diff line change
Expand Up @@ -543,9 +543,7 @@ @article{foerster2023twocomponent
title = {Two-Component GW Calculations: Cubic Scaling Implementation and Comparison of Vertex-Corrected and Partially Self-Consistent GW Variants},
year = {2023},
doi = {10.1021/acs.jctc.3c00512},
eprint = { https://doi.org/10.1021/acs.jctc.3c00512 },
url = {https://doi.org/10.1021/acs.jctc.3c00512
},
url = {https://doi.org/10.1021/acs.jctc.3c00512},
}
@article{glasbrenner2020efficient,
author = {Glasbrenner, Michael and Graf, Daniel and Ochsenfeld, Christian},
Expand Down
2 changes: 1 addition & 1 deletion docs/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,4 +22,4 @@ Currently, our library supports the following electronic structure methods:
# Technical Details
GreenX is written in Fortran 2008. Functionality needed for testing and error handling is written in C and Python. In the following we provide a few more details for current and future developers of GreenX.
- [Structure of the library](structure.md)
- [Coding conventions and regression testing](test.md)
- [Coding conventions and regression testing](tests.md)

0 comments on commit 6b14f03

Please sign in to comment.