Low-scaling GW: The space-time formalism: Difference between revisions

From VASP Wiki
No edit summary
No edit summary
 
(3 intermediate revisions by 2 users not shown)
Line 1: Line 1:
Available as of VASP.6 are low-scaling algorithms for [[RPA/ACFDT: Correlation energy in the Random Phase Approximation|ACFDT/RPA]].<ref name="kaltak2"/> This page describes the formalism of the corresponding low-scaling GW approach.<ref name="liu"/>
Available as of VASP.6 are low-scaling algorithms for [[RPA/ACFDT: Correlation energy in the Random Phase Approximation|ACFDT/RPA]].{{cite|kaltak:prb:2014}} This page describes the formalism of the corresponding low-scaling GW approach.{{cite|liu:prb:2016}}
A theoretical description of the ACFDT/RPA total energies is found [[ACFDT/RPA calculations#ACFDTR/RPAR|here]]. A brief summary regarding GW theory is given below, while a practical guide can be found [[GW calculations#LowGW|here]].  
A theoretical description of the ACFDT/RPA total energies is found [[ACFDT/RPA calculations#ACFDTR/RPAR|here]]. A brief summary regarding GW theory is given below, while a practical guide can be found [[GW calculations#LowGW|here]].  


== Theory ==  
== Theory ==  


The GW implementations in VASP described in the papers of Shishkin ''et al.''<ref name="shishkin-PRB74"/><ref name="shishkin-PRB75"/> avoid storage of the Green's function <math>G</math> as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of <math>\chi</math> and <math>\Sigma</math> in reciprocal space and results in a relatively high computational cost that scales with <math>N^4</math> (number of electrons).
The GW implementations in VASP described in the papers of Shishkin ''et al.''{{cite|shishkin:prb:2006}}{{cite|shishkin:prb:2007}} avoid storage of the Green's function <math>G</math> as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of <math>\chi</math> and <math>\Sigma</math> in reciprocal space and results in a relatively high computational cost that scales with <math>N^4</math> (number of electrons).


The scaling with system size can, however, be reduced to <math>N^3</math> by performing a so-called Wick-rotation to imaginary time <math>t\to i\tau</math>.<ref name="rojas"/>
The scaling with system size can, however, be reduced to <math>N^3</math> by performing a so-called Wick-rotation to imaginary time <math>t\to i\tau</math>.{{cite|rojas:prl:1995}}


Following the [[Groundstate in the Random Phase Approximation#ACFDTR/RPAR| low scaling ACFDT/RPA algorithms]] the space-time implementation determines first, the non-interacting Green's function on the imaginary time axis in real space  
Following the [[Groundstate in the Random Phase Approximation#ACFDTR/RPAR| low scaling ACFDT/RPA algorithms]] the space-time implementation determines first, the non-interacting Green's function on the imaginary time axis in real space  
Line 12: Line 12:
<math>G({\bf r},{\bf r}',i\tau)=-\sum_{n{\bf k}}\phi_{n{\bf k}}^{(0)}({\bf r}) \phi_{n{\bf k}}^{*(0)}({\bf r}') e^{-(\epsilon_{n{\bf k}}-\mu)\tau}\left[\Theta(\tau)(1-f_{n{\bf k}})-\Theta(-\tau)f_{n{\bf k}}\right]</math>  
<math>G({\bf r},{\bf r}',i\tau)=-\sum_{n{\bf k}}\phi_{n{\bf k}}^{(0)}({\bf r}) \phi_{n{\bf k}}^{*(0)}({\bf r}') e^{-(\epsilon_{n{\bf k}}-\mu)\tau}\left[\Theta(\tau)(1-f_{n{\bf k}})-\Theta(-\tau)f_{n{\bf k}}\right]</math>  


Here <math>\Theta</math> is the step function and <math>f_{n{\bf k}}</math> the occupation number of the state <math>\phi_{n{\bf k}}^{(0)}</math>. Because the Green's function is non-oscillatory on the imaginary time axis it can be represented on a coarse grid <math>\tau_{m}</math>, where the number of time points can be selected in VASP via the {{TAG|NOMEGA}} tag. Usually 12 to 16 points are sufficient for insulators and small band gap systems.<ref name="kaltak"/>
Here <math>\Theta</math> is the step function and <math>f_{n{\bf k}}</math> the occupation number of the state <math>\phi_{n{\bf k}}^{(0)}</math>. Because the Green's function is non-oscillatory on the imaginary time axis it can be represented on a coarse grid <math>\tau_{m}</math>, where the number of time points can be selected in VASP via the {{TAG|NOMEGA}} tag. Usually 12 to 16 points are sufficient for insulators and small band gap systems.{{cite|kaltak:2014}}


Subsequently, the irreducible polarizability is calculated from a contraction of two imaginary time Green's functions   
Subsequently, the irreducible polarizability is calculated from a contraction of two imaginary time Green's functions   
Line 20: Line 20:
</math>
</math>


Afterwards, the same compressed Fourier transformation as for the [[Groundstate in the Random Phase Approximation#ACFDTR/RPAR| low scaling ACFDT/RPA algorithms]] is employed to obtain the irreducible polarizability in reciprocal space on the imaginary frequency axis <math>\chi({\bf r},{\bf r}',i\tau_m) \to \chi_{{\bf G}{\bf G}'}({\bf q},i \omega_n) </math>.<ref name="kaltak"/><ref name="liu"/>
Afterwards, the same compressed Fourier transformation as for the [[Groundstate in the Random Phase Approximation#ACFDTR/RPAR| low scaling ACFDT/RPA algorithms]] is employed to obtain the irreducible polarizability in reciprocal space on the imaginary frequency axis <math>\chi({\bf r},{\bf r}',i\tau_m) \to \chi_{{\bf G}{\bf G}'}({\bf q},i \omega_n) </math>.{{cite|kaltak:2014}}{{cite|liu:prb:2016}}


The next step is the computation of the screened potential
The next step is the computation of the screened potential
Line 46: Line 46:
on the real frequency axis <math>z=\omega</math>.  
on the real frequency axis <math>z=\omega</math>.  


Because preceding Fourier transformations have been carried out with exponentially suppressed errors, the analytical continuation <math>\Sigma(z)</math> of the self-energy can be determined with high accuracy. The analytical continuation typically yields energies that differ less than 20 meV from quasi-particle energies obtained from the real-frequency calculation.<ref name="liu"/>
Because preceding Fourier transformations have been carried out with exponentially suppressed errors, the analytical continuation <math>\Sigma(z)</math> of the self-energy can be determined with high accuracy. The analytical continuation typically yields energies that differ less than 20 meV from quasi-particle energies obtained from the real-frequency calculation.{{cite|liu:prb:2016}}


In addition, the space-time formulation allows to solve the full Dyson equation for <math>G({\bf r,r'},i\tau)</math> with decent computational cost.<ref name="grumet"/> This approach is known as the self-consistent GW approach (scGW) and is available as of VASP6.<noinclude>
In addition, the space-time formulation allows to solve the full Dyson equation for <math>G({\bf r,r'},i\tau)</math> with decent computational cost.{{cite|grumet:prb:2018}} This approach is known as the self-consistent GW approach (scGW) and is available as of VASP6.<noinclude>


=== Finite temperature formalism===
=== Finite temperature formalism===
Line 54: Line 54:


== References ==
== References ==
<references>
<ref name="kaltak2">[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.90.054115 M. Kaltak, J. Klimeš, and G. Kresse, Phys. Rev. B 90, 054115 (2014)]</ref>
<ref name="liu">[http://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.165109 P. Liu, M. Kaltak, J. Klimes and G. Kresse, Phys. Rev. B 94, 165109 (2016).]</ref>
<ref name="rojas">[https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.74.1827 H. N. Rojas, R. W. Godby and R. J. Needs, Phys. Rev. Lett. 74, 1827 (1995)]</ref>
<ref name="shishkin-PRB74">[http://link.aps.org/doi/10.1103/PhysRevB.74.035101 M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006).]</ref>
<ref name="shishkin-PRB75">[http://link.aps.org/doi/10.1103/PhysRevB.75.235102 M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).]</ref>
<ref name="liu">[http://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.165109 P. Liu, M. Kaltak, J. Klimes and G. Kresse, Phys. Rev. B 94, 165109 (2016).]</ref>
<ref name="kaltak">[http://pubs.acs.org/doi/abs/10.1021/ct5001268 M. Kaltak, J. Klimeš, and G. Kresse, Journal of Chemical Theory and Computation 10, 2498-2507 (2014). ]</ref>
<ref name="grumet">[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.155143 M. Grumet, P. Liu, M. Kaltak, J. Klimeš and Georg Kresse, Phys. Ref. B 98, 155143 (2018). ]</ref>
</references>
----


[[Category:VASP]][[Category:Many-body perturbation theory]][[Category:VASP6]][[Category:Low-scaling GW and RPA]][[Category:Theory]]
 
[[Category:VASP]][[Category:Many-body perturbation theory]][[Category:Low-scaling GW and RPA]][[Category:Theory]]

Latest revision as of 10:18, 22 September 2025

Available as of VASP.6 are low-scaling algorithms for ACFDT/RPA.[1] This page describes the formalism of the corresponding low-scaling GW approach.[2] A theoretical description of the ACFDT/RPA total energies is found here. A brief summary regarding GW theory is given below, while a practical guide can be found here.

Theory

The GW implementations in VASP described in the papers of Shishkin et al.[3][4] avoid storage of the Green's function [math]\displaystyle{ G }[/math] as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of [math]\displaystyle{ \chi }[/math] and [math]\displaystyle{ \Sigma }[/math] in reciprocal space and results in a relatively high computational cost that scales with [math]\displaystyle{ N^4 }[/math] (number of electrons).

The scaling with system size can, however, be reduced to [math]\displaystyle{ N^3 }[/math] by performing a so-called Wick-rotation to imaginary time [math]\displaystyle{ t\to i\tau }[/math].[5]

Following the low scaling ACFDT/RPA algorithms the space-time implementation determines first, the non-interacting Green's function on the imaginary time axis in real space

[math]\displaystyle{ G({\bf r},{\bf r}',i\tau)=-\sum_{n{\bf k}}\phi_{n{\bf k}}^{(0)}({\bf r}) \phi_{n{\bf k}}^{*(0)}({\bf r}') e^{-(\epsilon_{n{\bf k}}-\mu)\tau}\left[\Theta(\tau)(1-f_{n{\bf k}})-\Theta(-\tau)f_{n{\bf k}}\right] }[/math]

Here [math]\displaystyle{ \Theta }[/math] is the step function and [math]\displaystyle{ f_{n{\bf k}} }[/math] the occupation number of the state [math]\displaystyle{ \phi_{n{\bf k}}^{(0)} }[/math]. Because the Green's function is non-oscillatory on the imaginary time axis it can be represented on a coarse grid [math]\displaystyle{ \tau_{m} }[/math], where the number of time points can be selected in VASP via the NOMEGA tag. Usually 12 to 16 points are sufficient for insulators and small band gap systems.[6]

Subsequently, the irreducible polarizability is calculated from a contraction of two imaginary time Green's functions

[math]\displaystyle{ \chi({\bf r},{\bf r}',i\tau_m) = -G({\bf r},{\bf r}',i\tau_m)G({\bf r}',{\bf r},-i\tau_m) }[/math]

Afterwards, the same compressed Fourier transformation as for the low scaling ACFDT/RPA algorithms is employed to obtain the irreducible polarizability in reciprocal space on the imaginary frequency axis [math]\displaystyle{ \chi({\bf r},{\bf r}',i\tau_m) \to \chi_{{\bf G}{\bf G}'}({\bf q},i \omega_n) }[/math].[6][2]

The next step is the computation of the screened potential

[math]\displaystyle{ W_{{\bf G}{\bf G}'}({\bf q},i\omega_m)=\left[\delta_{{\bf G}{\bf G}'}-\chi_{{\bf G}{\bf G}'}({\bf q},i\omega_m)V_{{\bf G}{\bf G}'}({\bf q})\right]^{-1}V_{{\bf G}{\bf G}'}({\bf q}) }[/math]


followed by the inverse Fourier transform [math]\displaystyle{ W_{{\bf G}{\bf G}'}({\bf q},i \omega_n) \to \chi({\bf r},{\bf r}',i\tau_m) }[/math] and the calculation of the self-energy

[math]\displaystyle{ \Sigma({\bf r},{\bf r}',i\tau_m) = -G({\bf r},{\bf r}',i\tau_m)W({\bf r}',{\bf r},i\tau_m) }[/math]

From here, several routes are possible including all approximations mentioned above, that is the single-shot, EVG0 and QPEVG0 approximation. All approximations have one point in common.

In contrast to the real-frequency implementation, the low-scaling GW algorithms require an analytical continuation of the self-energy from the imaginary frequency axis to the real axis. In general, this is an ill-defined problem and usually prone to errors, since the self-energy is known on a finite set of points. VASP determines internally a Padé approximation of the self-energy [math]\displaystyle{ \Sigma(z) }[/math] from the calculated set of NOMEGA points [math]\displaystyle{ \Sigma(i\omega_n) }[/math] and solves the non-linear eigenvalue problem

[math]\displaystyle{ \left[ T+V_{ext}+V_h+\Sigma(z) \right]\left|\phi_{n\bf k}\right\rangle = z\left| \phi_{n\bf k} \right\rangle }[/math]

on the real frequency axis [math]\displaystyle{ z=\omega }[/math].

Because preceding Fourier transformations have been carried out with exponentially suppressed errors, the analytical continuation [math]\displaystyle{ \Sigma(z) }[/math] of the self-energy can be determined with high accuracy. The analytical continuation typically yields energies that differ less than 20 meV from quasi-particle energies obtained from the real-frequency calculation.[2]

In addition, the space-time formulation allows to solve the full Dyson equation for [math]\displaystyle{ G({\bf r,r'},i\tau) }[/math] with decent computational cost.[7] This approach is known as the self-consistent GW approach (scGW) and is available as of VASP6.

Finite temperature formalism

The zero-temperature formalism of many-body perturbation theory breaks down for metals (systems with zero energy band-gap) as pointed out by Kohn and Luttinger.[8] This conundrum is lifted by considering diagrammatic perturbation theory at finite temperature [math]\displaystyle{ T\gt 0 }[/math], which may be understood by an analytical continuation of the real-time [math]\displaystyle{ t }[/math] to the imaginary time axis [math]\displaystyle{ -i\tau }[/math]. Matsubara has shown that this Wick rotation in time [math]\displaystyle{ t\to-i\tau }[/math] reveals an intriguing connection to the inverse temperature [math]\displaystyle{ \beta=1/T }[/math] of the system.[9] More precisely, Matsubara has shown that all terms in perturbation theory at finite temperature can be expressed as integrals of imaginary time quantities (such as the polarizability [math]\displaystyle{ \chi(-i\tau) }[/math]) over the fundamental interval [math]\displaystyle{ -\beta\le\tau\le\beta }[/math].

As a consequence, one decomposes imaginary time quantities into a Fourier series with period [math]\displaystyle{ \beta }[/math] that determines the spacing of the Fourier modes. For instance the imaginary polarizability can be written as

[math]\displaystyle{ \chi(-i\tau)=\frac1\beta\sum_{m=-\infty}^\infty \tilde \chi(i\nu_m)e^{-i\nu_m\tau},\quad \nu_m=\frac{2m}\beta\pi }[/math]

and the corresponding random-phase approximation of the correlation energy at finite temperature becomes a series over (in this case, bosonic) Matsubara frequencies

[math]\displaystyle{ \Omega_c^{\rm RPA}=\frac12\frac1\beta \sum_{m=-\infty}^\infty {\rm Tr}\left\lbrace \ln\left[ 1 -\tilde \chi(i\nu_m) V \right] -\tilde \chi(i\nu_m) V \right\rbrace,\quad \nu_m=\frac{2m}\beta\pi }[/math]

The Matsubara formalism has the advantage that all contributions to the Green's function and the polarizability are mathematically well-defined, including contributions from states close to the chemical potential [math]\displaystyle{ \epsilon_{n{\bf k}}\approx \mu }[/math], such that Matsubara series also converge for metallic systems.

Although formally convenient, the Matsubara series converges poorly with the number of considered terms in practice. VASP, therefore, uses a compressed representation of the Fourier modes by employing the Minimax-Isometry method.[10] This approach converges exponentially with the number of considered frequency points.

References