•17 correlation factors are tested for their performance in MP2-F12 calculations.•Good correlation factors have the near-universal concave shape with a radius of 1.5 a.u.•Long-range behavior and pointwise satisfaction of cusp condition are less important.We analyze the performance of 17 different correlation factors in explicitly correlated second-order many-body perturbation calculations for correlation energies. Highly performing correlation factors are found to have near-universal shape and size in the short range of electron-electron distance (01.5 a.u.) is insignificant insofar as the factor becomes near constant, leaving an orbital expansion to describe decoupled electrons. An analysis based on a low-rank Taylor expansion of the correlation factor seems limited, except that a negative second derivative with the value of around −1.3 a.u. correlates with high performance.Download high-res image (59KB)Download full-size image

In the Monte Carlo second-order many-body perturbation (MC-MP2) method, the long sum-of-product matrix expression of the MP2 energy, whose literal evaluation may be poorly scalable, is recast into a single high-dimensional integral of functions of electron pair coordinates, which is evaluated by the scalable method of Monte Carlo integration. The sampling efficiency is further accelerated by the redundant-walker algorithm, which allows a maximal reuse of electron pairs. Here, a multitude of graphical processing units (GPUs) offers a uniquely ideal platform to expose multilevel parallelism: fine-grain data-parallelism for the redundant-walker algorithm in which millions of threads compute and share orbital amplitudes on each GPU; coarse-grain instruction-parallelism for near-independent Monte Carlo integrations on many GPUs with few and infrequent interprocessor communications. While the efficiency boost by the redundant-walker algorithm on central processing units (CPUs) grows linearly with the number of electron pairs and tends to saturate when the latter exceeds the number of orbitals, on a GPU it grows quadratically before it increases linearly and then eventually saturates at a much larger number of pairs. This is because the orbital constructions are nearly perfectly parallelized on a GPU and thus completed in a near-constant time regardless of the number of pairs. In consequence, an MC-MP2/cc-pVDZ calculation of a benzene dimer is 2700 times faster on 256 GPUs (using 2048 electron pairs) than on two CPUs, each with 8 cores (which can use only up to 256 pairs effectively). We also numerically determine that the cost to achieve a given relative statistical uncertainty in an MC-MP2 energy increases as O(n3) or better with system size n, which may be compared with the O(n5) scaling of the conventional implementation of deterministic MP2. We thus establish the scalability of MC-MP2 with both system and computer sizes.

Electron binding energies are evaluated as differences in total energy between the N- and (N ± 1)-electron systems calculated by the nth-order Møller–Plesset perturbation (MPn) theory using the same set of orbitals. The MPn energies up to n = 30 are, in turn, obtained by the determinant-based method of Knowles et al. (Chem. Phys. Lett. 1985, 113, 8–12). The zeroth- through third-order electron binding energies thus determined agree with those obtained by solving the Dyson equation in the diagonal and frequency-independent approximations of the self-energy. However, as n → ∞, they converge at the exact basis-set solutions from the Dyson equation with the exact self-energy, which is nondiagonal and frequency-dependent. This suggests that the MPn energy differences define an alternative diagrammatic expansion of Koopmans-like electron binding energies, which takes into account the perturbation corrections from the off-diagonal elements and frequency dependence of the irreducible self-energy. Our analysis shows that these corrections are included as semireducible and linked-disconnected diagrams, respectively, which are also found in a perturbation expansion of the electron binding energies of the equation-of-motion coupled-cluster methods. The rate of convergence of the electron binding energies with respect to n and its acceleration by Padé approximants are also discussed.

An embedded-fragment ab initio second-order many-body perturbation (MP2) method is applied to an infinite three-dimensional crystal of carbon dioxide phase I (CO2-I), using the aug-cc-pVDZ and aug-cc-pVTZ basis sets, the latter in conjunction with a counterpoise correction for the basis-set superposition error. The equation of state, phonon frequencies, bulk modulus, heat capacity, Grüneisen parameter (including mode Grüneisen parameters for acoustic modes), thermal expansion coefficient (α), and thermal pressure coefficient (β) are computed. Of the factors that enter the expression of α, MP2 reproduces the experimental values of the heat capacity, Grüneisen parameter, and molar volume accurately. However, it proves to be exceedingly difficult to determine the remaining factor, the bulk modulus (B0), the computed value of which deviates from the observed value by 50–100%. As a result, α calculated by MP2 is systematically too low, while having the correct temperature dependence. The thermal pressure coefficient, β = αB0, which is independent of B0, is more accurately reproduced by theory up to 100 K.

Molecular crystals are chemists’ solids in the sense that their structures and properties can be understood in terms of those of the constituent molecules merely perturbed by a crystalline environment. They form a large and important class of solids including ices of atmospheric species, drugs, explosives, and even some organic optoelectronic materials and supramolecular assemblies. Recently, surprisingly simple yet extremely efficient, versatile, easily implemented, and systematically accurate electronic structure methods for molecular crystals have been developed. The methods, collectively referred to as the embedded-fragment scheme, divide a crystal into monomers and overlapping dimers and apply modern molecular electronic structure methods and software to these fragments of the crystal that are embedded in a self-consistently determined crystalline electrostatic field. They enable facile applications of accurate but otherwise prohibitively expensive ab initio molecular orbital theories such as Møller–Plesset perturbation and coupled-cluster theories to a broad range of properties of solids such as internal energies, enthalpies, structures, equation of state, phonon dispersion curves and density of states, infrared and Raman spectra (including band intensities and sometimes anharmonic effects), inelastic neutron scattering spectra, heat capacities, Gibbs energies, and phase diagrams, while accounting for many-body electrostatic (namely, induction or polarization) effects as well as two-body exchange and dispersion interactions from first principles. They can fundamentally alter the role of computing in the studies of molecular crystals in the same way ab initio molecular orbital theories have transformed research practices in gas-phase physical chemistry and synthetic chemistry in the last half century.In this Account, after a brief summary of formalisms and algorithms, we discuss applications of these methods performed in our group as compelling illustrations of their unprecedented power in addressing some of the outstanding problems of solid-state chemistry, high-pressure chemistry, or geochemistry. They are the structure and spectra of ice Ih, in particular, the origin of two peaks in the hydrogen-bond-stretching region of its inelastic neutron scattering spectra, a solid–solid phase transition from CO2-I to elusive, metastable CO2-III, pressure tuning of Fermi resonance in solid CO2, and the structure and spectra of solid formic acid, all at the level of second-order Møller–Plesset perturbation theory or higher.

Second-order many-body perturbation theory [MBPT(2)] is the lowest-ranked member of a systematic series of approximations convergent at the exact solutions of the Schrödinger equations. It has served and continues to serve as the testing ground for new approximations, algorithms, and even theories. This article introduces this basic theory from a variety of viewpoints including the Rayleigh–Schrödinger perturbation theory, the many-body Green’s function theory based on the Dyson equation, and the related Feynman–Goldstone diagrams. It also explains the important properties of MBPT(2) such as size consistency, its ability to describe dispersion interactions, and divergence in metals. On this basis, this article surveys three major advances made recently by the authors to this theory. They are a finite-temperature extension of MBPT(2) and the resolution of the Kohn–Luttinger conundrum, a stochastic evaluation of the correlation and self-energies of MBPT(2) using the Monte Carlo integration of their Laplace-transformed expressions, and an extension to anharmonic vibrational zero-point energies and transition frequencies based on the Dyson equation.

A Monte Carlo (MC) integration of the second-order many-body perturbation (MP2) corrections to energies and self-energies eliminates the usual computational bottleneck of the MP2 algorithm (i.e., the basis transformation of two-electron integrals), thereby achieving near-linear size dependence of its operation cost, a negligible core and disk memory cost, and a naturally parallel computational kernel. In this method, the correlation correction expressions are recast into high-dimensional integrals, which are approximated as the sums of integrands evaluated at coordinates of four electron random walkers guided by a Metropolis algorithm for importance sampling. The gravest drawback of this method, however, is the inevitable statistical uncertainties in its results, which decay slowly as the inverse square-root of the number of MC steps. We propose an algorithm that can increase the number of MC sampling points in each MC step by many orders of magnitude by having 2m electron walkers (m > 2) and then using m(m – 1)/2 permutations of their coordinates in evaluating the integrands. Hence, this algorithm brings an O(m2)-fold increase in the number of MC sampling points at a mere O(m) additional cost of propagating redundant walkers, which is a net O(m)-fold enhancement in sampling efficiency. We have demonstrated a large performance increase in the Monte Carlo MP2 calculations for the correlation energies of benzene and benzene dimer as well as for the correlation corrections to the energy, ionization potential, and electron affinity of C60. The calculation for C60 has been performed with a parallel implementation of this method running on up to 400 processors of a supercomputer, yielding an accurate prediction of its large electron affinity, which makes its derivative useful as an electron acceptor in bulk heterojunction organic solar cells.

The mathematical constructs of the Dyson coordinates and geometry are introduced. The former are a unitary transformation of the normal coordinates and the anharmonic vibrational counterpart of the Dyson orbitals in electronic structure theory. The first-order Dyson coordinates bring the sums of the harmonic force constants and their first-order diagrammatic perturbation corrections (the first-order Dyson self-energy) to a diagonal form. The first-order Dyson geometry has no counterpart in electronic structure theory. It is the point on the potential energy surface at which the sums of the energy gradients and their first-order diagrammatic perturbation corrections vanish. It agrees with the vibrationally averaged geometry of vibrational self-consistent field (VSCF) theory in the bulk limit. These constructs provide a unified view of the relationship of VSCF and its diagrammatically size-consistent modifications as well as the self-consistent phonon method widely used in solid-state physics.

A pedagogical, semi-rigorous proof is presented for the existence of the thermodynamic (infinite-volume) limit of the energy per volume for an electrically neutral, metallic or nonmetallic crystal. The proof is based on the demonstration of the same for individual energy components, namely, the kinetic, Coulomb, and exchange contributions to the Hartree–Fock energy as well as the correlation contribution obtained by many-body perturbation or coupled-cluster theory.

A graphical abstract is available for this content

A linear-scaling, embedded-fragment, second-order many-body perturbation (MP2) method with basis sets up to aug-cc-pVTZ is applied to the antiparallel structure of solid hydrogen fluoride and deuterium fluoride under 0–20 GPa of ambient pressure. The optimized structures, including the lattice parameters and molar volume, and phonon dispersion as well as phonon density of states (DOS), are determined as a function of pressure. The basis-set superposition errors are removed by the counterpoise correction. The structural parameters at 0 GPa calculated by MP2 agree accurately with the observed, making the predicted values at higher pressures a useful pilot for future experiments. The corresponding values obtained by the Hartree–Fock method have large, systematic errors. The MP2/aug-cc-pVDZ frequencies of the infrared- and Raman-active vibrations of the three-dimensional solids are in good agreement with the observed and also justify previous vibrational analyses based on one-dimensional chain models; the non-coincidence of the infrared and Raman mode pairs can be explained as factor-group (Davydov) splitting. The exceptions are one pair of modes in the librational region, for which band assignments based on a one-dimensional chain model need to be revised, as well as the five pseudo-translational modes that exist only in a three-dimensional treatment. The observed pressure dependence of Raman bands in the stretching region, which red-shift with pressure, is accounted for by theory only qualitatively, while that in the pseudo-translational region is reproduced with quantitative accuracy. The present calculation proves to be limited in explaining the complex pressure dependence of the librational modes. The hydrogen-amplitude-weighted phonon DOS at 0 GPa is much less structured than the DOS obtained from one-dimensional models and may be more realistic in view of the also broad, structureless observed inelastic neutron scattering spectra. All major observed peaks can be straightforwardly assigned to the calculated peaks in the DOS. With increasing pressure, MP2 predicts further broadening of bands and breach of the demarcation between the pseudo-translational and librational bands.

This article aims to dispel confusions about the definition of size consistency as well as some incompatibility that exists between different criteria for judging whether an electronic structure theory is size consistent and thus yields energies and other quantities having correct asymptotic size dependence. It introduces extensive and intensive diagram theorems, which provide unambiguous sufficient conditions for size consistency for extensive and intensive quantities, respectively, stipulated in terms of diagrammatic topology and vertex makeup. The underlying algebraic size-consistency criterion is described, which relies on the polynomial dependence of terms in the formalism on the number of wave vector sampling points in Brillouin-zone integrations. The physical meanings of two types of normalization of excitation amplitudes in electron-correlation theories, namely, the intermediate and standard normalization, are revealed. The amplitudes of the operator that describes an extensive quantity (the extensive operator) are subject to the intermediate normalization, while those of the operator that corresponds to an intensive quantity (the intensive operator) must be normalized. The article also introduces the extensive-intensive consistency theorem which specifies the relationship between the spaces of determinants reached by the extensive and intensive operators in a size-consistent method for intensive quantities. Furthermore, a more fundamental question is addressed as to what makes energies extensive and thus an application of thermodynamics to chemistry valid. It is shown that the energy of an electrically neutral, periodic, non-metallic solid is extensive. On this basis, a strictly size-consistent redefinition of the Hartree--Fock theory is proposed.

A pedagogical, semi-rigorous proof is presented for the existence of the thermodynamic (infinite-volume) limit of the energy per volume for an electrically neutral, metallic or nonmetallic crystal. The proof is based on the demonstration of the same for individual energy components, namely, the kinetic, Coulomb, and exchange contributions to the Hartree–Fock energy as well as the correlation contribution obtained by many-body perturbation or coupled-cluster theory.

A linear-scaling, embedded-fragment, second-order many-body perturbation (MP2) method with basis sets up to aug-cc-pVTZ is applied to the antiparallel structure of solid hydrogen fluoride and deuterium fluoride under 0–20 GPa of ambient pressure. The optimized structures, including the lattice parameters and molar volume, and phonon dispersion as well as phonon density of states (DOS), are determined as a function of pressure. The basis-set superposition errors are removed by the counterpoise correction. The structural parameters at 0 GPa calculated by MP2 agree accurately with the observed, making the predicted values at higher pressures a useful pilot for future experiments. The corresponding values obtained by the Hartree–Fock method have large, systematic errors. The MP2/aug-cc-pVDZ frequencies of the infrared- and Raman-active vibrations of the three-dimensional solids are in good agreement with the observed and also justify previous vibrational analyses based on one-dimensional chain models; the non-coincidence of the infrared and Raman mode pairs can be explained as factor-group (Davydov) splitting. The exceptions are one pair of modes in the librational region, for which band assignments based on a one-dimensional chain model need to be revised, as well as the five pseudo-translational modes that exist only in a three-dimensional treatment. The observed pressure dependence of Raman bands in the stretching region, which red-shift with pressure, is accounted for by theory only qualitatively, while that in the pseudo-translational region is reproduced with quantitative accuracy. The present calculation proves to be limited in explaining the complex pressure dependence of the librational modes. The hydrogen-amplitude-weighted phonon DOS at 0 GPa is much less structured than the DOS obtained from one-dimensional models and may be more realistic in view of the also broad, structureless observed inelastic neutron scattering spectra. All major observed peaks can be straightforwardly assigned to the calculated peaks in the DOS. With increasing pressure, MP2 predicts further broadening of bands and breach of the demarcation between the pseudo-translational and librational bands.