Expeditious Stochastic Calculation of Random-Phase Approximation Energies for Thousands of Electrons in Three Dimensions

Reference Energies for Valence Ionizations and Satellite Transitions

Upon ionization of an atom or a molecule, another electron (or more) can be simultaneously excited. These concurrently generated states are called "satellites" (or shakeup transitions) as they appear in ionization spectra as higher-energy peaks with weaker intensity and larger width than the main peaks associated with single-particle ionizations. Satellites, which correspond to electronically excited states of the cationic species, are notoriously challenging to model using conventional single-reference methods due to their high excitation degree compared to the neutral reference state. This work reports 42 satellite transition energies and 58 valence ionization potentials (IPs) of full configuration interaction quality computed in small molecular systems. Following the protocol developed for the quest database [Véril, M.; Scemama, A.; Caffarel, M.; Lipparini, F.; Boggio-Pasqua, M.; Jacquemin, D.; and Loos, P.-F. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2021, 11, e1517], these reference energies are computed using the configuration interaction using a perturbative selection made iteratively (CIPSI) method. In addition, the accuracy of the well-known coupled-cluster (CC) hierarchy (CC2, CCSD, CC3, CCSDT, CC4, and CCSDTQ) is gauged against these new accurate references. The performances of various approximations based on many-body Green's functions (GW, GF2, and T-matrix) for IPs are also analyzed. Their limitations in correctly modeling satellite transitions are discussed.

Massively parallel implementation of gradients within the random phase approximation: Application to the polymorphs of benzene

The Random-Phase approximation (RPA) provides an appealing framework for semi-local density functional theory. In its Resolution-of-the-Identity (RI) approach, it is a very accurate and more cost-effective method than most other wavefunction-based correlation methods. For widespread applications, efficient implementations of nuclear gradients for structure optimizations and data sampling of machine learning approaches are required. We report a well scaling implementation of RI-RPA nuclear gradients on massively parallel computers. The approach is applied to two polymorphs of the benzene crystal obtaining very good cohesive and relative energies. Different correction and extrapolation schemes are investigated for further improvement of the results and estimations of error bars.

Time-Dependent Second-Order Green's Function Theory for Neutral Excitations

We develop a time-dependent second-order Green's function theory (GF2) for calculating neutral excited states in molecules. The equation of motion for the lesser Green's function (GF) is derived within the adiabatic approximation to the Kadanoff-Baym (KB) equation, using the second-order Born approximation for the self-energy. In the linear response regime, we recast the time-dependent KB equation into a Bethe-Salpeter-like equation (GF2-BSE), with a kernel approximated by the second-order Coulomb self-energy. We then apply our GF2-BSE to a set of molecules and atoms and find that GF2-BSE is superior to configuration interaction with singles (CIS) and/or time-dependent Hartree-Fock (TDHF), particularly for charge-transfer excitations, and is comparable to CIS with perturbative doubles (CIS(D)) in most cases.

Reduced Scaling of Optimal Regional Orbital Localization via Sequential Exhaustion of the Single-Particle Space

Wannier functions have become a powerful tool in the electronic structure calculations of extended systems. The generalized Pipek-Mezey Wannier functions exhibit appealing characteristics (e.g., reaching an optimal localization and the separation of the σ–π orbitals) compared with other schemes. However, when applied to giant nanoscale systems, the orbital localization suffers from a large computational cost overhead if one is interested in localized states in a small fragment of the system. Herein, we present a swift, efficient, and robust approach for obtaining regionally localized orbitals of a subsystem within the generalized Pipek-Mezey scheme. The proposed algorithm introduces a reduced work space and sequentially exhausts the entire orbital space until the convergence of the localization functional. It tackles systems with ∼10000 electrons within 0.5 h with no loss in localization quality compared to the traditional approach. Regionally localized orbitals with a higher extent of localization are obtained via judiciously extending the subsystem’s size. Exemplifying on large bulk and a 4 nm wide slab of diamond with an NV^– center, we demonstrate the methodology and discuss how the choice of the localization region affects the excitation energy of the defect. Furthermore, we show how the sequential algorithm is easily extended to stochastic methodologies that do not provide individual single-particle eigenstates. It is thus a promising tool to obtain regionally localized states for solving the electronic structure problems of a subsystem embedded in giant condensed systems.

Stochastic evaluation of four-component relativistic second-order many-body perturbation energies: A potentially quadratic-scaling correlation method

A second-order many-body perturbation correction to the relativistic Dirac-Hartree-Fock energy is evaluated stochastically by integrating 13-dimensional products of four-component spinors and Coulomb potentials. The integration in the real space of electron coordinates is carried out by the Monte Carlo (MC) method with the Metropolis sampling, whereas the MC integration in the imaginary-time domain is performed by the inverse-CDF (cumulative distribution function) method. The computational cost to reach a given relative statistical error for spatially compact but heavy molecules is observed to be no worse than cubic and possibly quadratic with the number of electrons or basis functions. This is a vast improvement over the quintic scaling of the conventional, deterministic second-order many-body perturbation method. The algorithm is also easily and efficiently parallelized with demonstrated 92% strong scalability going from 64 to 4096 processors for a fixed job size.

O(N) Stochastic Evaluation of Many-Body van der Waals Energies in Large Complex Systems

We propose a new strategy to solve the key equations of the many-body dispersion (MBD) model by Tkatchenko, DiStasio Jr., and Ambrosetti. Our approach overcomes the original O(N3) computational complexity that limits its applicability to large molecular systems within the context of O(N) density functional theory. First, to generate the required frequency-dependent screened polarizabilities, we introduce an efficient solution to the Dyson-like self-consistent screening equations. The scheme reduces the number of variables and, coupled to a direct inversion of the iterative subspace extrapolation, exhibits linear-scaling performances. Second, we apply a stochastic Lanczos trace estimator resolution to the equations evaluating the many-body interaction energy of coupled quantum harmonic oscillators. While scaling linearly, it also enables communication-free pleasingly parallel implementations. As the resulting O(N) stochastic massively parallel MBD approach is found to exhibit minimal memory requirements, it opens up the possibility of computing accurate many-body van der Waals interactions of millions-atoms' complex materials and solvated biosystems with computational times in the range of minutes.

Stochastic Vector Techniques in Ground-State Electronic Structure

We review a suite of stochastic vector computational approaches for studying the electronic structure of extended condensed matter systems. These techniques help reduce algorithmic complexity, facilitate efficient parallelization, simplify computational tasks, accelerate calculations, and diminish memory requirements. While their scope is vast, we limit our study to ground-state and finite temperature density functional theory (DFT) and second-order perturbation theory. More advanced topics, such as quasiparticle (charge) and optical (neutral) excitations and higher-order processes, are covered elsewhere. We start by explaining how to use stochastic vectors in computations, characterizing the associated statistical errors. Next, we show how to estimate the electron density in DFT and discuss highly effective techniques to reduce statistical errors. Finally, we review the use of stochastic vector techniques for calculating correlation energies within the second-order Møller-Plesset perturbation theory and its finite temperature variational form. Example calculation results are presented and used to demonstrate the efficacy of the methods. Expected final online publication date for the Annual Review of Physical Chemistry, Volume 73 is April 2022. Please see http://www.annualreviews.org/page/journal/pubdates for revised estimates.

Monte Carlo MP2-F12 for Noncovalent Interactions: The C60 Dimer

A scalable stochastic algorithm is presented that can evaluate explicitly correlated (F12) second-order many-body perturbation (MP2) energies of weak, noncovalent, intermolecular interactions. It first transforms the formulas of the MP2 and F12 energy differences into a short sum of high-dimensional integrals of Green's functions in real space and imaginary time. These integrals are then evaluated by the Monte Carlo method augmented by parallel execution, redundant-walker convergence acceleration, direct-sampling autocorrelation elimination, and control-variate error reduction. By sharing electron-pair walkers across the supermolecule and its subsystems spanned by the joint basis set, the statistical uncertainty is reduced by one to 2 orders of magnitude in the MP2 binding energy corrected for the basis-set incompleteness and superposition errors. The method predicts the MP2-F12/aug-cc-pVDZ binding energy of 19.1 ± 4.0 kcal mol^-1 for the C₆₀ dimer at the center distance of 9.748 Å.

Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space-Time Approach

We present an implementation of the GW space-time approach that allows cubic-scaling all-electron calculations with standard Gaussian basis sets without exploiting any localization or sparsity considerations. The independent-electron susceptibility is constructed in a time representation over a nonuniform distribution of real-space locations {r_k} optimized within a separable resolution-of-the-identity framework to reproduce standard Coulomb-fitting calculations with meV accuracy. The compactness of the obtained {r_k} distribution leads to a crossover with the standard Coulomb-fitting scheme for system sizes below a few hundred electrons. The needed analytic continuation follows a recent approach that requires the continuation of the screened Coulomb potential rather than the much more structured self-energy. The present scheme is benchmarked over large molecular sets, and scaling properties are demonstrated on a family of defected hexagonal boron-nitride flakes containing up to 6000 electrons.

Stochastic evaluation of fourth-order many-body perturbation energies

A scalable, stochastic algorithm evaluating the fourth-order many-body perturbation (MP4) correction to energy is proposed. Three hundred Goldstone diagrams representing the MP4 correction are computer generated and then converted into algebraic formulas expressed in terms of Green's functions in real space and imaginary time. They are evaluated by the direct (i.e., non-Markov, non-Metropolis) Monte Carlo (MC) integration accelerated by the redundant-walker and control-variate algorithms. The resulting MC-MP4 method is efficiently parallelized and is shown to display O(n^5.3) size-dependence of cost, which is nearly two ranks lower than the O(n⁷) dependence of the deterministic MP4 algorithm. It evaluates the MP4/aug-cc-pVDZ energy for benzene, naphthalene, phenanthrene, and corannulene with the statistical uncertainty of 10 mE_h (1.1% of the total basis-set correlation energy), 38 mE_h (2.6%), 110 mE_h (5.5%), and 280 mE_h (9.0%), respectively, after about 10⁹ MC steps.

Stochastically Realized Observables for Excitonic Molecular Aggregates

We show that a stochastic approach enables calculations of the optical properties of large 2-dimensional and nanotubular excitonic molecular aggregates. Previous studies of such systems relied on numerically diagonalizing the dense and disordered Frenkel Hamiltonian, which scales approximately as O(N3) for N dye molecules. Our approach scales much more efficiently as O(Nlog(N)), enabling quick study of systems with a million of coupled molecules on the micrometer size scale. We calculate several important experimental observables, including the optical absorption spectrum and density of states, and develop a stochastic formalism for the participation ratio. Quantitative agreement with traditional matrix diagonalization methods is demonstrated for both small- and intermediate-size systems. The stochastic methodology enables the study of the effects of spatial-correlation in site energies on the optical signatures of large 2D aggregates. Our results demonstrate that stochastic methods present a path forward for screening structural parameters and validating experiments and theoretical predictions in large excitonic aggregates.

Separable resolution-of-the-identity with all-electron Gaussian bases: Application to cubic-scaling RPA

We explore a separable resolution-of-the-identity (RI) formalism built on quadratures over limited sets of real-space points designed for all-electron calculations. Our implementation preserves, in particular, the use of common atomic orbitals and their related auxiliary basis sets. The setup of the present density fitting scheme, i.e., the calculation of the system specific quadrature weights, scales cubically with respect to the system size. Extensive accuracy tests are presented for the Fock exchange and MP2 correlation energies. We finally demonstrate random phase approximation (RPA) correlation energy calculations with a scaling that is cubic in terms of operations, quadratic in memory, with a small crossover with respect to our standard RI-RPA implementation.

Exchange-correlation energy from pairing matrix fluctuation and the particle-particle random phase approximation

Despite their unmatched success for many applications, commonly used local, semi-local, and hybrid density functionals still face challenges when it comes to describing long-range interactions, static correlation, and electron delocalization. Density functionals of both the occupied and virtual orbitals are able to address these problems. The particle-hole (ph-) Random Phase Approximation (RPA), a functional of occupied and virtual orbitals, has recently known a revival within the density functional theory community. Following up on an idea introduced in our recent communication [H. van Aggelen, Y. Yang, and W. Yang, Phys. Rev. A 88, 030501 (2013)], we formulate more general adiabatic connections for the correlation energy in terms of pairing matrix fluctuations described by the particle-particle (pp-) propagator. With numerical examples of the pp-RPA, the lowest-order approximation to the pp-propagator, we illustrate the potential of density functional approximations based on pairing matrix fluctuations. The pp-RPA is size-extensive, self-interaction free, fully anti-symmetric, describes the strong static correlation limit in H2, and eliminates delocalization errors in H2(+) and other single-bond systems. It gives surprisingly good non-bonded interaction energies--competitive with the ph-RPA--with the correct R(-6) asymptotic decay as a function of the separation R, which we argue is mainly attributable to its correct second-order energy term. While the pp-RPA tends to underestimate absolute correlation energies, it gives good relative energies: much better atomization energies than the ph-RPA, as it has no tendency to underbind, and reaction energies of similar quality. The adiabatic connection in terms of pairing matrix fluctuation paves the way for promising new density functional approximations.

Breaking the theoretical scaling limit for predicting quasiparticle energies: the stochastic GW approach

We develop a formalism to calculate the quasiparticle energy within the GW many-body perturbation correction to the density functional theory. The occupied and virtual orbitals of the Kohn-Sham Hamiltonian are replaced by stochastic orbitals used to evaluate the Green function G, the polarization potential W, and, thereby, the GW self-energy. The stochastic GW (sGW) formalism relies on novel theoretical concepts such as stochastic time-dependent Hartree propagation, stochastic matrix compression, and spatial or temporal stochastic decoupling techniques. Beyond the theoretical interest, the formalism enables linear scaling GW calculations breaking the theoretical scaling limit for GW as well as circumventing the need for energy cutoff approximations. We illustrate the method for silicon nanocrystals of varying sizes with N_{e}>3000 electrons.

A Guided Stochastic Energy-Domain Formulation of the Second Order Møller-Plesset Perturbation Theory

We develop an alternative formulation in the energy-domain to calculate the second order Møller-Plesset (MP2) perturbation energies. The approach is based on repeatedly choosing four random energies using a nonseparable guiding function, filtering four random orbitals at these energies, and averaging the resulting Coulomb matrix elements to obtain a statistical estimate of the MP2 correlation energy. In contrast to our time-domain formulation, the present approach is useful for both quantum chemistry and real-space/plane wave basis sets. The scaling of the MP2 calculation is roughly linear with system size, providing a useful tool to study dispersion energies in large systems. This is demonstrated on a structure of 64 fullerenes within the SZ basis as well as on silicon nanocrystals using real-space grids.

Self-averaging stochastic Kohn-Sham density-functional theory

We formulate the Kohn-Sham density functional theory (KS-DFT) as a statistical theory in which the electron density is determined from an average of correlated stochastic densities in a trace formula. The key idea is that it is sufficient to converge the total energy per electron to within a predefined statistical error in order to obtain reliable estimates of the electronic band structure, the forces on nuclei, the density and its moments, etc. The fluctuations in the total energy per electron are guaranteed to decay to zero as the system size increases. This facilitates "self-averaging" which leads to the first ever report of sublinear scaling KS-DFT electronic structure. The approach sidesteps calculation of the density matrix and thus, is insensitive to its evasive sparseness, as demonstrated here for silicon nanocrystals. The formalism is not only appealing in terms of its promise to far push the limits of application of KS-DFT, but also represents a cognitive change in the way we think of electronic structure calculations as this stochastic theory seamlessly converges to the thermodynamic limit.