Many-body methods related to the GW approximation have proven to
be a highly effective tool for computing core-level excitations [1]. In particular,
the contour deformation (CD) is an efficient, scalable and numerically stable
approach that has enabled core-level calculations on systems up to 100 atoms [2].
In this work, we reduce the scaling of CD applied to core-levels from \(O(N^5)\)
to \(O(N^4)\), using an analytic continuation of the screened Coulomb interaction (W).
Analytic continuation techniques using Padé approximants have been used in GW for
computing the self-energy \(\Re\Sigma_n(\omega)\), but until very recently [3] this idea
was not fully exploited for the W matrices.
The new method (CD-WAC) has been implemented in
FHI-aims. CD-WAC has been extensively
tested on well established benchmark sets like the GW100 and the CORE65 [4],
reporting MAEs of less than 5 meV with respect to CD.
The theoretical scaling has been confirmed by performing scaling experiments
on large acene chains and amorphous carbon. Speedups of 5 times have been
attained with CD-WAC for the largest systems.
This work won a Poster Commendation Prize at the
Psi-k Conference 2022
in Lausanne, Switzerland.
The contour deformation method
A clever application of Jordan's Lemma and
Cauchy's residue theorem.
Integral self-energy expression is a consequence of
the Convolution theorem.
\[
\displaystyle
\small
\begin{align*}
\Sigma_c&(\mathbf{r},\mathbf{r}',\omega) = \oint-\int_{\Im}-\int_{D_{\Gamma^+}}-\int_{D_{\Gamma^-}} =\\
& -\frac{1}{2\pi}\int_{-\infty}^{+\infty}d\omega'G_{0}(\mathbf{r},\mathbf{r}',\omega+i\omega')
W_{0}(\mathbf{r},\mathbf{r}',i\omega')\\[4pt]
& - \sum_i \phi_i(\mathbf{r})\phi_i(\mathbf{r'})
W_{0}(\mathbf{r},\mathbf{r}',|\epsilon_i-\omega|+i\eta)\theta(\epsilon_i-\omega)\\
& + \sum_a \phi_a(\mathbf{r})\phi_a(\mathbf{r'})
W_{0}(\mathbf{r},\mathbf{r}',|\epsilon_a-\omega|+i\eta)\theta(\omega-\epsilon_a)\\
& = R(\omega) - I(i\omega)
\end{align*}
\]
Constructing W using RI-V: the CD way
Many quantities involving 4-center integrals can be efficiently computed in FHI-aims
using the resolution of the identity with coulomb metric (RI-V). The \(W_{mn}(\omega)\)
matrices of CD are no exception:
\[
W_{mn}(\omega) = \sum_{PQ}O_{P}^{nm}[\mathbf{1} -
\mathbf{\Pi}(|\epsilon_m - \omega| + i\eta)]^{-1}_{PQ}O_{Q}^{mn}
\]
where the polarizability matrices are given by:
\[
\Pi_{PQ}(i\omega) = 2 \sum_{ia}O_P^{ia}
\frac{\epsilon_i-\epsilon_a}{\omega^2 + (\epsilon_i-\epsilon_a)^2}O_Q^{ia}
\]
Analytic continuation of W: the rationale for CD-WAC
CD-WAC reduces the formal scaling of CD by computing the W matrices
that appear in the \(R(\omega)\) term using an analytic continuation \(P_i(\omega_j)=W_{mn}(\omega_j)\). The
latter is computed using a modified version of the Thiele's reciprocal
differences algorithm developed in our group:
\[
\small
P_{N}(z) = \cfrac{a_1}{1+\cfrac{a_{2}(z^2-\omega^2_{1})}
{1+\cfrac{a_{3}(z^2-\omega^2_{2})}{1+\cdots+
\cfrac{a_p(z^2-\omega^2_{p-1})}{1 + (z^2-\omega^2_p)g_{p+1}(z)}}}}
\]
The asymptotic behavior of the algorithms for computing the \(R(\omega)\) and \(I(i\omega)\)
terms is:
\[R \sim N_{res}N_{occ}N_{virt}N^2_{aux} \qquad I \sim N_{\omega}N_{occ}N_{virt}N^2_{aux}\]
We are interested in core level spectroscopy. In this case, the integration grid for imaginary frequencies is
independent of the system size, whereas the number of residues is proportional to number of occupied states.
This implies that the global scaling of CD increases to \(O(N^5)\) for the deepest states.
CD-WAC implementation details
The CD-WAC algorithm builds upon the CD one, whose details can be found on Reference
[2]. CD-WAC relays on a series of heuristic optimizations,
presented in the flowchart below (left panel). Depending on the type of atom and orbital considered,
three different sampling strategies for the analytic continuation are possible (right panel).
The full monty example: H2O molecule
The following example illustrates the predictive power of CD-WAC. We have
represented the screened coulomb interaction matrix element \(W_{1s,1s}\) and the self-energy
matrix element \(\Sigma_{1s}\). The level of theory of these calculation is \(G_0W_0\)@PBE
using NAOs of the second tier. The CD-WAC calculation is ~ 30 times faster than CD, at the expense
of only 50 extra frequency points for the interpolation.
Validation of the CD-WAC method with respect to CD using the
GW100 and CORE65 tests suites. The level of theory is \(G_0W_0\)@PBE/def2-QVPZ and
\(G_0W_0\)@PBEh(\(\alpha=0.45\))/cc-pVTZ respectively. For the valence states of the GW100 set only imaginary frequencies
were used, for the CORE65 50 additional frequencies on two regions were employed.
Scaling plot comparing CD and CD-WAC applied to large acene chains C4n+2H2n+4
with n ranging from 1 to 15. We were interested on the C1s core level of the molecules, the level of
theory of the calculations is \(G_0W_0\)@PBEh(\(\alpha=0.5\))/def2-QVPZ. The expected scale reduction
is thus verified by the numerical experiment.
References and web version
-
D. Golze, M. Dvorak, and P. Rinke. Front. Chem. 7 (2019): 377.
- D. Golze, J. Wilhelm, M.J. Van Setten, and P. Rinke.
J. Chem. Theory Comput., 14(9):4856–4869, 2018.
- I. Duchemin and X. Blase. J. Chem. Theory Comput., 16(3):1742–1756, 2020.
- D. Golze, L. Keller, and P. Rinke. J. Phys. Chem. Lett.
11.5 (2020): 1840-1847.
