Auxiliary second-order Green’s function perturbation theory (AGF2)#

Modules: pyscf.agf2

Introduction#

Auxiliary second-order Green’s function perturbation theory (AGF2) [42, 43] is an iterative, $O [N^{5}]$ scaling post-Hartree–Fock method primarily intended for the calculation of charged excitation spectra, ionisation potentials (IPs) and electron affinities (EAs), with energetics and single-particle static properties also available. It can also be considered loosely as an iterative approximate ADC(2) method. One advantage is that the entire spectrum of charged excitations is found simultaneously, as the full Green’s function is self-consistently optimised. AGF2 calculations can be performed with or without density fitting, and for restricted or unrestricted Hartree–Fock references, and is available with hybrid parallelism capabilities.

Basic usage of AGF2 as given in examples/agf2/00-simple_agf2.py, as

#!/usr/bin/env python
#
# Author: Oliver J. Backhouse <olbackhouse@gmail.com>
#         George H. Booth <george.booth@kcl.ac.uk>
#

'''
A simple example of restricted AGF2. AGF2 will compute correlation energies, one-particle
properties and charged excitations / energy levels via an iterated, renormalized perturbation
theory.

Default AGF2 corresponds to the AGF2(1,0) method outlined in the papers:
  - O. J. Backhouse, M. Nusspickel and G. H. Booth, J. Chem. Theory Comput., 16, 1090 (2020).
  - O. J. Backhouse and G. H. Booth, J. Chem. Theory Comput., 16, 6294 (2020).
'''

from pyscf import gto, scf, agf2

mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='cc-pvdz')

mf = scf.RHF(mol)
mf.conv_tol = 1e-12
mf.run()

# Run an AGF2 calculation
gf2 = agf2.AGF2(mf)
gf2.conv_tol = 1e-7
gf2.run(verbose=4)

# Print the first 3 ionization potentials
# Note that there is no additional cost to write out larger numbers of excitations.
gf2.ipagf2(nroots=3)

# Print the first 3 electron affinities
gf2.eaagf2(nroots=3)

Theory#

In AGF2, one iteratively solves the Dyson equation

G (ω) = G_{0} (ω) + G_{0} (ω) Σ (ω) G (ω),

using the self-energy correct through second-order perturbation theory, which can be written as

\begin{aligned} Σ_{x y} (ω) & = \sum_{i j a} \frac{(x i | j a) [2 (y i | j a) - (y j | i a)]}{ω - (E_{i} + E_{j} - E_{a})} \\ + \sum_{a b i} \frac{(x a | b i) [2 (y a | b i) - (y b | a i)]}{ω - (E_{a} + E_{b} - E_{i})} . \end{aligned}

In AGF2, the first two spectral moments of the MP2 self-energy,

\begin{aligned} T_{x y}^{(0)} & = \sum_{i j a} (x i | j a) [2 (y i | j a) - (y j | i a)] \\ + \sum_{a b i} (x a | b i) [2 (y a | b i) - (y b | a i)] \\ T_{x y}^{(1)} & = \sum_{i j a} (x i | j a) [2 (y i | j a) - (y j | i a)] (E_{i} + E_{j} - E_{a}) \\ + \sum_{a b i} (x a | b i) [2 (y a | b i) - (y b | a i)] (E_{a} + E_{b} - E_{i}), \end{aligned}

are used as a conserved quantity in order to coarse-grain this frequency dependence as a set of compressed and renormalized $2 p 1 h$ and $1 p 2 h$ states. These are used to form an effective single-particle Hamiltonian whose size now only scales with $O [N]$ , rather than $O [N^{3}]$ if the full frequency-dependence was maintained. This allows Dyson equation to be solved via the hermitian eigenvalue problem

\begin{array}{r} (\begin{array}{c} F & v \\ v^{†} & ϵ \end{array}) ϕ_{n} = λ_{n} ϕ_{n}, \end{array}

where $F$ is the Fock matrix, $ϵ$ the compressed states coupling to $F$ with coupling strength $v$ , and $ϕ_{n}$ are termed the quasi-molecular orbitals (QMOs), which represent the Dyson orbitals at each iteration, and span both the original MOs and these auxiliary functions. The QMOs can be reinserted into the self-energy in place of the $i, j, a$ and $a, b, i$ indices, permitting self-consistency to convergence in the Green’s function with this coarse-grained description of the frequency-dependence of the self-energy.

Additionally, once projected into the physical space, the QMOs define a correlated (non-idempotent) density matrix, permitting a further self-consistency in the Fock matrix in a similar fashion to a Hartree–Fock calculation. The correct number of electrons remain in the physical space throughout this process via optimisation of a chemical potential $μ$ . Finally, we note that self-consistency based on higher-order spectral moments of either the self-energy or resulting Green’s function at each iteration is possible, with a limiting behaviour to the traditional second-order Green’s function method (GF2). However, numerical investigations have shown that this higher-order self-consistency leads to a deterioration of results, and therefore the consistency based on the first two self-energy spectral moments is the default behaviour.

Photoemission spectra#

The compression of the effective dynamics performed in AGF2 permits the calculation of the full spectrum of charged excitations, as no additional computational effort is required for this (in contrast to iterative eigensolvers of effective hamiltonians), and the computational effort is therefore independent of the number of excitations requested. An example of a calculation which provides the full spectral function can be found in examples/agf2/03-photoemission_spectra.py:

#!/usr/bin/env python
#
# Author: Oliver J. Backhouse <olbackhouse@gmail.com>
#         George H. Booth <george.booth@kcl.ac.uk>
#

'''
Use a converged AGF2 calculation to build the full photoemission (quasiparticle) spectrum

Default AGF2 corresponds to the AGF2(1,0) method outlined in the papers:
  - O. J. Backhouse, M. Nusspickel and G. H. Booth, J. Chem. Theory Comput., 16, 1090 (2020).
  - O. J. Backhouse and G. H. Booth, J. Chem. Theory Comput., 16, 6294 (2020).
'''

import numpy
from pyscf import gto, scf, agf2

mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='cc-pvdz')

mf = scf.RHF(mol)
mf.conv_tol = 1e-12
mf.run()

# Run an AGF2 calculation
gf2 = agf2.AGF2(mf)
gf2.conv_tol = 1e-7
gf2.run()

# Access the GreensFunction object and compute the spectrum
gf = gf2.gf
grid = numpy.linspace(-10.0, 10.0, 1000)
eta = 0.02
spectrum = gf.real_freq_spectrum(grid, eta=eta)
# The array `spectrum` is now a length nfreq array of the
# spectral function -1/pi * Tr[Im[G(\omega + i\eta)]].
# Note that for UHF AGF2 calculations, the individual gf 
# elements can be passed to real_freq_spectrum in order 
# to obtain the spin-resolved spectra.

# We can also build the self-energy on the real-frequency axis
# by accessing the renormalized auxiliary states:
e = gf2.se.energy - gf2.se.chempot
v = gf2.se.coupling
denom = grid[:,None] - (e + numpy.sign(e)*eta*1.0j)[None]
se = numpy.einsum('xk,yk,wk->wxy', v, v.conj(), 1./denom)

Additionally, the Dyson orbitals corresponding to individual excitations are accessible directly from the method, which can be seen in examples/agf2/10-dyson_orbitals.py.

Restart a calculation#

The contents of an AGF2 calculation can be dumped to the disk via the familiar PySCF chkfile utilities. By default, an AGF2 method will inherit the chkfile attribute of the underlying Hartree–Fock reference object. An example of restoring an AGF2 calculation from a checkpoint file can be found in examples/agf2/04-restart.py:

#!/usr/bin/env python
#
# Author: Oliver J. Backhouse <olbackhouse@gmail.com>
#         George H. Booth <george.booth@kcl.ac.uk>
#

'''
An example of restarting an AGF2 calculation.

The agf2.chkfile module provides similar functionality to the existing
chkfile utilities in pyscf, but prevents failure during MPI runs.
'''

import numpy
from pyscf import gto, scf, agf2, lib

mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='cc-pvdz', verbose=5)

mf = scf.RHF(mol)
mf.conv_tol = 1e-12
# if we are using MPI, we only want pyscf to save a chkfile on the root process:
if agf2.mpi_helper.rank == 0:
    mf.chkfile = 'agf2.chk'
mf.run()

# Run an AGF2 calculation:
gf2 = agf2.AGF2(mf)
gf2.conv_tol = 1e-7
gf2.max_cycle = 3
gf2.run()

# Restore the Mole and SCF first
mol = agf2.chkfile.load_mol('agf2.chk')
mf = scf.RHF(mol)
mf.__dict__.update(agf2.chkfile.load('agf2.chk', 'scf'))

# Restore the AGF2 calculation
gf2a = agf2.AGF2(mf)
gf2a.__dict__.update(agf2.chkfile.load_agf2('agf2.chk')[1])
gf2a.max_cycle = 50
gf2a.run()

Parallelisation#

The AGF2 module supports both MPI an OpenMP parallelisation in an aim to provide a scalable method applicable to interesting problems in quantum chemistry. Furthermore, the dominant scaling step is embarrassingly parallel. Distribution of computational load is handled by the optional dependency mpi4py, and will run without MPI using OpenMP threads if an installation cannot be found.