Analytic Interatomic Forces in the Random Phase Approximation

We discuss that in the random phase approximation (RPA) the first derivative of the energy with respect to the Green ’ s function is the self-energy in the GW approximation. This relationship allows us to derive compact equations for the RPA interatomic forces. We also show that position dependent overlap operators are elegantly incorporated in the present framework. The RPA force equations have been implemented in the projector augmented wave formalism, and we present illustrative applications, including ab initio molecular dynamics simulations, the calculation of phonon dispersion relations for diamond and graphite, as well as structural relaxations for water on boron nitride. The present derivation establishes a concise framework for forces within perturbative approaches and is also applicable to more involved approximations for the correlation energy.

Density functional theory (DFT) is almost ubiquitous in present day computational chemistry and materials science.Compared to other methods, one of the great advantages of DFT is that forces between the atoms are readily computable, so that relaxed ground state structures, vibrational properties, as well as finite temperature properties of any material are straightforwardly accessible.Although quantum chemists have devised ways to compute analytic forces, e.g., for Møller-Plesset perturbation theory (MP) already decades ago [1,2], in condensed matter simulations forces beyond DFT have been rarely available.This has limited the success of more involved electronic structure methods such as the random phase approximation (RPA) for solids.Only recently, this situation has improved, with first derivatives being computable for second order MP perturbation theory (MP2) for solids [3], as well as for the RPA for molecules [4,5].However, to date all these implementations are limited to real valued Gaussian-type orbitals and even then the algebra is often intricate.
Here we use Green's function theory, to derive compact equations for the forces in the random phase approximation.These can be recapped in two lines.One first calculates a Hermitian first order one-particle density matrix involving the Fock operator F, the correlation part of the self-energy in the GW approximation Σ GW , and the independent particle Green's function GðiνÞ of the Kohn-Sham (KS) Hamiltonian.Then one needs to trace this density matrix over the first derivate of the KS Hamiltonian F i ¼ −Tr½ρ ð1Þ ∂H KS =∂R i , where R i are the ionic coordinates and F i is the corresponding component of the force.We have implemented the corresponding equations in the projector augmented wave [6,7] formalism and show the potential of this method in the present Letter.
The fundamental problem we address here is that perturbation theory (including the RPA) is nonvariational; i.e., the kinetic, Hartree, exchange, and correlation energies are evaluated using a set of orbitals calculated using some other zero-order Hamiltonian.For the purpose of this presentation, we assume that these one-electron orbitals are always determined using a KS Hamiltonian H.We rely on Green's function theory and establish the essentials first and then proceed to an outline of the derivation.Given the generalized eigenvalue equation with a KS Hamiltonian H and an overlap operator S, the time dependent Schrödinger equation reads Then the single particle KS Green's function G can be written as In the present Letter, the Green's functions are represented in imaginary frequency iν and time iτ, where they are smooth functions of their arguments.The Fermi-level is placed at μ ¼ 0; occupied (unoccupied) states have all negative (positive) one-electron energies.
Green's function theory allows us to write compactly the change of the Green's function upon perturbations.This is achieved by the Dyson equation, which relates the perturbed Green's function G þ ΔG corresponding to HþΔH and SþΔS to the unperturbed Green's function G: Neglecting ΔG on the rhs of the second line, one obtains for the change of the Green's function ΔG the first order term It is noted that ΔH is the total change of the KS Hamiltonian including the self-consistent field changes.The Green's function in imaginary time can be obtained by a Fourier transformation: Using standard mathematical relations, it is easy to show that one can also write [8,9] GðiτÞ ¼ 8 > < > : Here the sums are restricted to the occupied (occ) manifold for negative τ and to the unoccupied or virtual (virt) manifold for positive τ, and thus the exponents are always bound.Clearly [see Eq. ( 8)], the one-particle density operator ρ is given by the Green's function for τ → 0 − , i.e., for a slightly negative imaginary time where η ¼ 0 − is a small negative infinitesimal.
We now proceed to calculate first derivatives of the Hartree-Fock energy functional.This functional depends on the density matrix ρ and an external potential V ext : The trace operator usually involves integration over one (for the variations of the external potential) or two spatial coordinates (for variations of the density matrix).If one calculates the HF energy derivatives at the HF density matrix ρ HF , terms related to variations of the density matrix do not need to be calculated, since the Hartree-Fock density matrix fulfills the equation where δρ is subject to the orthonormality constraint of the corresponding basis functions and δE HF =δρ is the Fock operator F. In this case the second term in Eq. ( 10) vanishes, which is just the famous Hellmann-Feynman theorem.
If the random phase approximation is used, the Hartree-Fock energy functional, however, needs to be evaluated at the KS density matrix.The corresponding energy is commonly referred to as exact exchange energy (EXX) [10].The non-Hellmann-Feynman term is then nonzero.Using Green's functions, we can write this term compactly by combining Eqs. ( 6), (7), and (9).For the non-Hellmann-Feynman term in Eq. ( 10) one then obtains where the Fock operator F is constructed from the KS density matrix.Equation ( 12) looks expensive to evaluate, since, seemingly, for any change of the Hamiltonian GðiνÞδHGðiνÞ needs to be calculated.However, since the trace is invariant under cyclic permutations, one can calculate two Hermitian matrices, once and then efficiently calculate the energy change for any variation of the KS Hamiltonian and overlap operator as We note that ρ ð1Þ is a density matrix, whereas γ ð1Þ has the dimension of an energy (first moment iν).Furthermore, since F is frequency independent, the integrals in Eq. ( 13) can be calculated analytically using the residue theorem, where the term involving η ¼ 0 − specifies how the contours need to be closed (see Supplemental Material [11]).This equation can be straightforwardly generalized to correlation energy functionals that depend on the Green's function E c ¼ E c ½G.Since the correlation energy usually does not depend on the external potential, the Hellmann-Feynman term is zero for this contribution.The non-Hellmann-Feynman term yields an additional contribution to the matrices, where the self-energy Σ c is the functional derivative of the correlation energy with respect to the Green's function Σ c ðiνÞ ¼ δE c =δGðiνÞ.Equation ( 15) is obviously analogous to Eq. ( 13), only the Fock operator F has been replaced by the self-energy.Furthermore, the exponential e −iνη can be dropped, since the integrand decays like oð1=ν 2 Þ, and the integral is unconditionally convergent.
For the correlation energy in the random phase approximation, E RPA , the corresponding self-energy is exactly the correlation part of the GW approximation to the self-energy in Hedin's approach [12], which we rationalize now.The RPA correlation energy is defined as [13][14][15][16][17][18][19] where χðiνÞ is the independent particle polarizability calculated using KS orbitals, and v is the Coulomb kernel.
The derivative with respect to χðiνÞ can be calculated by Taylor expanding the logarithm and taking the derivative (see also Ref. [20]).This yields the correlation part of the screened Coulomb potential Using the relation between the independent particle polarizability and the Green's function [12] χðr; r 0 ; iτÞ ¼ Gðr; r 0 ; iτÞGðr 0 ; r; −iτÞ; ð17Þ one formally obtains the desired relation (dropping the position coordinates): Details are given in the Supplemental Material [11].In summary, non-Hellmann-Feynman forces are simple to evaluate.One only needs to determine the one particle density matrix ρ ð1Þ and its first energy moment γ ð1Þ , Eqs. ( 13) and ( 15), and then evaluate Eq. ( 14).The first derivative of the self-consistent density functional theory Hamiltonian ∂H=∂R i and of the overlap operator ∂S=∂R i need to be calculated for each ionic displacement, δR i , yielding one component in the atomic forces.
Concerning the scaling, we note that we have recently implemented an algorithm that calculates the GW selfenergy with linear scaling in the number of k points and cubic scaling in the number of plane waves [21] following the work of Rojas, Godby, and Needs [22,23].The algorithm is designed to use a minimal number of time and frequency points, with 12 points yielding μeV accuracy per atom for gapped systems [24,25].The scaling with the number of plane waves and k points carries over to the calculation of the density matrices, so that the complexity of the entire calculation is, in principle, identical to standard DFT, but with much larger prefactors.In practice, one force calculation requires about 3 orders of magnitude more compute time than for a semilocal functional, about 10-50 times longer than for hybrid functionals (albeit scaling quadratically in the number of k points), and about 10 times longer than for a single RPA total energy calculation.One reason for the steep computational cost is that the calculation of the GW self-energy requires about 6 times more operations than an RPA total energy calculation.Furthermore, the evaluation of ∂H=∂R i for each ionic displacement is currently done using linear response theory and can take about the same time as the calculation of ρ ð1Þ and γ ð1Þ for systems with hundreds of atoms.On the positive side, G 0 W 0 quasiparticle energies for all orbitals are readily available and can be calculated en passant.
We have implemented the RPA forces in the projector augmented wave (PAW) [6,7] based Vienna ab initio simulation package (VASP).Typically we can obtain 5 digits accuracy compared to a numerical differentiation, and the precision increases as the number of time points increases (see Supplemental Material [11]).A rigorous test for the validity of forces are first principles molecular dynamics (MD) simulations, since errors in the forces will result in a drift of the total energy.In Fig. 1 we show results for short simulations of crystalline diamond at 2000 K.The configurations were initially equilibrated using DFT, and then about 250 RPA-MD steps were performed.For the 64 atom supercell, one time step took 20 min on 40 Intel Xeon v3 cores.It is clear that the energy conservation is excellent for both the small 8 atom diamond cell, as well as the 64 atom cell.
Relaxations to the ionic ground state structure are also feasible and fairly efficient using the present code.One advantage of the present approach is that the second ionic derivatives of the density functional E KS , that is used to 106403-3 calculate the orbitals, are calculated when the linear response code is used.This implies that the KS Hessian matrix A is available, and this Hessian can be used as the preconditioner for the update of the ionic positions R R ← R − A −1 F; where F are the RPA forces.If the RPA Hessian is close to the KS Hessian, convergence to the ionic RPA ground state is quadratic and hence, in practice, the forces often decrease by an order of magnitude in each step.
Being able to relax structures adds a new important facet to the modeling of materials beyond density functional theory.Here we have looked at one prototypical problem, H 2 O adsorption on a BN layer [26,27].We have relaxed the H 2 O molecule using the RPA and find after few steps a low symmetry configuration atop the N atom that was previously not considered.The main reason for finding this configuration was the use of the DFT Hessian as preconditioner which results in a rapid rotation of the molecule along a very soft mode.All used functionals confirmed that this position is lower in energy than the previously considered structures (6 meV in the RPA).In Table I, we summarize the geometry and adsorption energy.The equilibrium distance to the surface decreases with increasing adsorption energy, with PBE and HSE yielding too small adsorption energies and too large distances.The van der Waals (vdW) corrected functionals generally yield larger adsorption energies and shorter distances than RPA.We expect the RPA to be accurate here, since RPA parameterized force fields yield excellent contact angles for water droplets [27].However, RPA generally overestimates the Pauli repulsion, with singles increasing the adsorption energy by 20 meV [28].
To demonstrate that the random phase approximation yields also very good vibrational properties, we have calculated the phonon dispersion relation of diamond and graphite (see Fig. 3).We note that the random phase approximation yields excellent structural properties for both, diamond and graphite, with lattice constants in virtually perfect agreement with experiment [35,36].For the phonons, agreement between the RPA and experiment is also excellent for diamond and very good for graphite.Although phonon dispersions from DFT (not shown) can be in very good agreement with experiment, they are often under-(semilocal functionals) or overestimated (hybrid functionals) [37,38].
In summary, we have presented an elegant scheme to calculate interatomic forces using many body perturbation theory.The derivation relies on Green's functions and the self-energy and seamlessly incorporates position dependent overlap operators.Although applied to the RPA here, analogous expressions can easily be derived for, e.g., MP2, with the RPA self-energy Σ GW simply replaced by the self-energy in second order Σ ð2Þ .Our derivation involves a first order one-particle density matrix, and we  note that seemingly related density matrices can be found in the quantum chemistry literature but without reference to a self-energy [1,2,4,5,44].We have put our equations to the test by implementing them in a plane wave code, and we have demonstrated that MD simulations, relaxations, and predictions of phonon frequencies are possible.We are convinced that the present work paves the way towards accurate beyond DFT modeling in condensed matter physics and materials sciences.

FIG. 1 .
FIG.1.Potential and kinetic energy per atom for diamond supercells of 64 and 8 atoms, using the Γ point and 2 × 2 × 2 k points, respectively.Energies are shifted so that the total energies approach zero.The slight fluctuations in the total energy are similar to those observed in DFT for a 1 fs time step using a Verlet algorithm.

TABLE I .
[39][40][41][42][43]ls with hydrogen to surface distance d and angle θ between the H-H axis and the surface normal (see Fig.2) and adsorption energy E. Phonon dispersion relation of diamond and graphite calculated using 128 atom supercells.RPA calculations were performed at the Γ point, with k-point sampling errors corrected for using DFT.Experimental data and lattice constants are taken from Refs.[39][40][41][42][43].