Diffusion Monte Carlo Calculations Research Articles

Finite-Size Error Cancellation in Diffusion Monte Carlo Calculations of Surface Chemistry.

The accurate prediction of reaction mechanisms in heterogeneous (surface) catalysis is one of the central challenges in computational chemistry. Quantum Monte Carlo methods─Diffusion Monte Carlo (DMC) in particular─are being recognized as higher-accuracy, albeit more computationally expensive, alternatives to Density Functional Theory (DFT) for energy predictions of catalytic systems. A major computational bottleneck in the broader adoption of DMC for catalysis is the need to perform finite-size extrapolations by simulating increasingly large periodic cells (supercells) to eliminate many-body finite-size effects and obtain energies in the thermodynamic limit. Here, we show that it is possible to significantly reduce this computational cost by leveraging the cancellation of many-body finite-size errors that accompanies the evaluation of energy differences when calculating quantities like adsorption (binding) energies and mapping potential energy surfaces. We analyze the cancellation and convergence of many-body finite-size errors in two well-known adsorbate/slab systems, H2O/LiH(001) and CO/Pt(111). Based on this analysis, we identify strategies for obtaining binding energies in the thermodynamic limit that optimally utilize error cancellation to balance accuracy and computational efficiency. Using one such strategy, we then predict the correct order of adsorption site preference on CO/Pt(111), a challenging problem for a wide range of density functionals. Our accurate and inexpensive DMC calculations are found to unambiguously recover the top > bridge > hollow site order, in agreement with experimental observations. We proceed to use this DMC method to map the potential energy surface of CO hopping between Pt(111) adsorption sites. This reveals the existence of an L-shaped top-bridge-hollow diffusion trajectory characterized by energy barriers that provide an additional kinetic justification for experimental observations of CO/Pt(111) adsorption. Overall, this work demonstrates that it is routinely possible to achieve order-of-magnitude speedups and memory savings in DMC calculations by taking advantage of error cancellation in the calculation of energy differences that are ubiquitous in heterogeneous catalysis and surface chemistry more broadly.

Fast Near Ab Initio Potential Energy Surfaces Using Machine Learning.

A machine-learning based approach for evaluating potential energies for quantum mechanical studies of properties of the ground and excited vibrational states of small molecules is developed. This approach uses the molecular-orbital-based machine learning (MOB-ML) method to generate electronic energies with the accuracy of CCSD(T) calculations at the same cost as a Hartree-Fock calculation. To further reduce the computational cost of the potential energy evaluations without sacrificing the CCSD(T) level accuracy, GPU-accelerated Neural Network Potential Energy Surfaces (NN-PES) are trained to geometries and energies that are collected from small-scale Diffusion Monte Carlo (DMC) simulations, which are run using energies evaluated using the MOB-ML model. The combined NN+(MOB-ML) approach is used in variational calculations of the ground and low-lying vibrational excited states of water and in DMC calculations of the ground states of water, CH5+, and its deuterated analogues. For both of these molecules, comparisons are made to the results obtained using potentials that were fit to much larger sets of electronic energies than were required to train the MOB-ML models. The NN+(MOB-ML) approach is also used to obtain a potential surface for C2H5+, which is a carbocation with a nonclassical equilibrium structure for which there is currently no available potential surface. This potential is used to explore the CH stretching vibrations, focusing on those of the bridging hydrogen atom. For both CH5+ and C2H5+ the MOB-ML model is trained using geometries that were sampled from an AIMD trajectory, which was run at 350 K. By comparison, the structures sampled in the ground state calculations can have energies that are as much as ten times larger than those used to train the MOB-ML model. For water a higher temperature AIMD trajectory is needed to obtain accurate results due to the smaller thermal energy. A second MOB-ML model for C2H5+ was developed with additional higher energy structures in the training set. The two models are found to provide nearly identical descriptions of the ground state of C2H5+.

Q-AQUA: A Many-Body CCSD(T) Water Potential, Including Four-Body Interactions, Demonstrates the Quantum Nature of Water from Clusters to the Liquid Phase.

Many model potential energy surfaces (PESs) have been reported for water; however, none are strictly from "first-principles". Here we report such a potential, based on a many-body representation at the CCSD(T) level of theory up to the four-body interaction. The new PES is benchmarked for the isomers of the water hexamer for dissociation energies, harmonic frequencies, and unrestricted diffusion Monte Carlo (DMC) calculations of the zero-point energies of the Prism, Book, and Cage isomers. Dissociation energies of several isomers of the 20-mer agree well with recent benchmark energies. Exploratory DMC calculations on this cluster verify the robustness of the new PES for quantum simulations. The accuracy and speed of the new PES are demonstrated for standard condensed phase properties, i.e., the radial distribution function and the self-diffusion constant. Quantum effects are shown to be substantial for these observables and also needed to bring theory into excellent agreement with experiment.

Weighted nodal domain averages of eigenstates for quantum Monte Carlo and beyond

We study the nodal properties of many-body eigenstates of stationary Schr\"odinger equation that affect the accuracy of real-space quantum Monte Carlo calculations. In particular, we introduce weighted nodal domain averages that provide a new probe of nodal surfaces beyond the usual expectations. Particular choices for the weight function reveal, for example, that the difference between two arbitrary fermionic eigenvalues is given by the nodal hypersurface integrals normalized by overlaps with the bosonic ground state of the given Hamiltonian. Noninteracting and fully interacting Be atom with corresponding almost exact and approximate wave functions are used to illustrate several aspects of these concepts. Variational formulations that employ different weights are proposed for prospective improvement of nodes in variational and fixed-node diffusion Monte Carlo calculations.

The binding of atomic hydrogen on graphene from density functional theory and diffusion Monte Carlo calculations.

In this work, density functional theory (DFT) and diffusion Monte Carlo (DMC) methods are used to calculate the binding energy of a H atom chemisorbed on the graphene surface. The DMC value of the binding energy is about 16% smaller in magnitude than the Perdew-Burke-Ernzerhof (PBE) result. The inclusion of exact exchange through the use of the Heyd-Scuseria-Ernzerhof functional brings the DFT value of the binding energy closer in line with the DMC result. It is also found that there are significant differences in the charge distributions determined using PBE and DMC approaches.

A quantum Monte Carlo study of systems with effective core potentials and node nonlinearities

We study beryllium dihydride (BeH$_2$) and acetylene (C$_2$H$_2$) molecules using real-space diffusion Monte Carlo (DMC) method. The molecules serve as perhaps the simplest prototypes that illustrate the difficulties with biases in the fixed-node DMC calculations that might appear with the use of effective core potentials (ECPs) or other nonlocal operators. This is especially relevant for the recently introduced correlation consistent ECPs (ccECPs) for $2s2p$ elements. Corresponding ccECPs exhibit deeper potential functions due to higher fidelity to all-electron counterparts, which could lead to larger local energy fluctuations. We point out that the difficulties stem from issues that are straightforward to address by upgrades of basis sets, use of T-moves for nonlocal terms, inclusion of a few configurations into the trial function and similar. The resulting accuracy corresponds to the ccECP target (chemical accuracy) and it is in consistent agreement with independent correlated calculations. Further possibilities for upgrading the reliability of the DMC algorithm and considerations for better adapted and more robust Jastrow factors are discussed as well.

Permutationally invariant polynomial regression for energies and gradients, using reverse differentiation, achieves orders of magnitude speed-up with high precision compared to other machine learning methods.

Permutationally invariant polynomial (PIP) regression has been used to obtain machine-learned potential energy surfaces, including analytical gradients, for many molecules and chemical reactions. Recently, the approach has been extended to moderate size molecules with up to 15 atoms. The algorithm, including "purification of the basis," is computationally efficient for energies; however, we found that the recent extension to obtain analytical gradients, despite being a remarkable advance over previous methods, could be further improved. Here, we report developments to further compact a purified basis and, more significantly, to use the reverse differentiation approach to greatly speed up gradient evaluation. We demonstrate this for our recent four-body water interaction potential. Comparisons of training and testing precision on the MD17 database of energies and gradients (forces) for ethanol against numerous machine-learning methods, which were recently assessed by Dral and co-workers, are given. The PIP fits are as precise as those using these methods, but the PIP computation time for energy and force evaluation is shown to be 10-1000 times faster. Finally, a new PIP potential energy surface (PES) is reported for ethanol based on a more extensive dataset of energies and gradients than in the MD17 database. Diffusion Monte Carlo calculations that fail on MD17-based PESs are successful using the new PES.

The three-center two-positron bond.

Computational studies have shown that one or more positrons can stabilize two repelling atomic anions through the formation of two-center positronic bonds. In the present work, we study the energetic stability of a system containing two positrons and three hydride anions, namely 2e+[H3 3-]. To this aim, we performed a preliminary scan of the potential energy surface of the system with both electrons and positrons in a spin singlet state, with a multi-component MP2 method, that was further refined with variational and diffusion Monte Carlo calculations, and confirmed an equilibrium geometry with D 3h symmetry. The local stability of 2e+[H3 3-] is demonstrated by analyzing the vertical detachment and adiabatic energy dissociation channels. Bonding properties of the positronic compound, such as the equilibrium interatomic distances, force constants, dissociation energies, and bonding densities are compared with those of the purely electronic H3 + and Li3 + systems. Through this analysis, we find compelling similarities between the 2e+[H3 3-] compound and the trilithium cation. Our results strongly point out the formation of a non-electronic three-center two-positron bond, analogous to the well-known three-center two-electron counterparts, which is fundamentally distinct from the two-center two-positron bond [D. Bressanini, J. Chem. Phys., 2021, 155, 054306], thus extending the concept of positron bonded molecules.

Energy Derivatives in Real-Space Diffusion Monte Carlo.

We present unbiased, finite-variance estimators of energy derivatives for real-space diffusion Monte Carlo calculations within the fixed-node approximation. The derivative dλE is fully consistent with the dependence E(λ) of the energy computed with the same time step. We address the issue of the divergent variance of derivatives related to variations of the nodes of the wave function both by using a regularization for wave function parameter gradients recently proposed in variational Monte Carlo and by introducing a regularization based on a coordinate transformation. The essence of the divergent variance problem is distilled into a particle-in-a-box toy model, where we demonstrate the algorithm.

Machine-learning approach to finite-size effects in systems with strongly interacting fermions

We investigate the applicability of machine learning techniques in studying the finite-size effects associated with many-body physics. These techniques have an emerging presence in many-body theory as they have been used for interpolations, extrapolations, and in modeling wavefunctions. We will resolve several issues associated with machine learning and many-body calculations such as small datasets, outliers, and discontinuities, for the purpose of extrapolating finite calculations to macroscopic scales. We carry out a systematic investigation of two related systems by developing metrics that aim to avoid spurious effects and capture desired features. This work uses neural networks to extrapolate the Unitary Gas to the thermodynamic limit at zero-range, which is otherwise difficult to reach. The effective mass of strongly interacting neutron matter is also studied and makes use of the non-interacting problem to resolve discontinuous predictions. For this investigation, we also carried out new Auxiliary Field Diffusion Monte Carlo (AFDMC) calculations for a variety of densities and particle numbers. Ultimately, we demonstrate an effective utility for neural networks in this context.

Metastable Metallic Phase of a Bilayer Blue Phosphorene Induced by Interlayer Bonding and Intralayer Charge Redistributions.

We have carried out diffusion Monte Carlo calculations for an A1B-1-stacked bilayer blue phosphorene to find that it undergoes a semiconductor-metal transition as the interlayer distance decreases. While the most stable bilayer structure is a semiconducting one with two monolayers coupled through a weak van der Waals interaction, the metallic bilayer at a shorter interlayer distance is found to be only metastable. This is in contrast to a recent theoretical prediction based on a random phase approximation that the metallic phase would be the most stable bilayer configuration of blue phosphorene. Our analysis of charge density distributions reveals that the metastable metallic phase is induced by interlayer chemical bonding and intralayer charge redistributions. This study enriches our understanding of interlayer binding of a blue phosphorene and contributes to the establishment of correct energetic order between its different phases, which will be essential in devising an experimental pathway for a metallic phosphorene.

GPU-Accelerated Neural Network Potential Energy Surfaces for Diffusion Monte Carlo.

Diffusion Monte Carlo (DMC) provides a powerful method for understanding the vibrational landscape of molecules that are not well-described by conventional methods. The most computationally demanding step of these calculations is the evaluation of the potential energy. In this work, a general approach is developed in which a neural network potential energy surface is trained by using data generated from a small-scale DMC calculation. Once trained, the neural network can be evaluated by using highly parallelizable calls to a graphics processing unit (GPU). The power of this approach is demonstrated for DMC simulations on H2O, CH5+, and (H2O)2. The need to include permutation symmetry in the neural network potentials is explored and incorporated into the molecular descriptors of CH5+ and (H2O)2. It is shown that the zero-point energies and wave functions obtained by using the neural network potentials are nearly identical to the results obtained when using the potential energy surfaces that were used to train the neural networks at a substantial savings in the computational requirements of the simulations.

Growth of rare gases on coronene

We study the growth of rare gases (Rg) Ne, Ar and Kr on coronene (Cor) molecules, Cor-(Rg)\(_N\). The study is taken up to sizes N=60, and we show the global energy minima for those clusters. Improved Lennard-Jones and atom-bond potentials are used to represent the Rg-Rg and Rg-Cor interactions, respectively. The Basin-Hopping (BH) global optimization technique is employed to locate the putative global energy minimum structures. Results are presented for Cor-(Ne), Cor-(Ar) and Cor-(Kr) systems. Both classical Cor-(Ar)\(_N\) and Cor-(Kr)\(_N\) clusters present the same “magic numbers” for N=3, 6, 10, 14, 19 and 38. This last number marks the first layer of solvation. For Cor-(Ne)\(_N\) clusters, we consider quantum effects adding the zero-point energy (ZPE) into the classical potential energy minima (BH-ZPE). This semiclassical approximation is compared with quantum diffusion Monte Carlo (DMC) calculations up to N=20, obtaining a good agreement and showing that this approximation works reasonably well. In this BH-ZPE approach, high stability configurations for Cor-(Ne)\(_N\) clusters are obtained at N=6, 14, 42, 51 and 57. However, we have not found any evidence of a first solvation layer using either the semiclassical BH-ZPE or the classical BH one.

Origin of the ν=1/2 fractional quantum Hall effect in wide quantum wells

The nature of the fractional quantum Hall effect at $\ensuremath{\nu}=1/2$, observed in wide quantum wells almost three decades ago, is still under debate. Previous studies have investigated it using the variational Monte Carlo method, which assumes that the transverse wave function and the gap between the symmetric and antisymmetric subbands obtained in a local density approximation at zero magnetic field remain valid even at high perpendicular magnetic fields; this method also ignores the effect of Landau level mixing. We develop in this work a three-dimensional fixed-phase diffusion Monte Carlo method, which gives, in a single framework, the total energies of various candidate states in a finite width quantum well, including Landau level mixing, directly in a large magnetic field. This method can be applied to one-component states and also to two-component states in the limit where the symmetric and antisymmetric bands are nearly degenerate. Our three-dimensional fixed-phase diffusion Monte Carlo calculations find that the one-component composite-fermion Fermi sea and the one-component Pfaffian states are very close in energy for a range of quantum-well widths and densities, suggesting that the observed 1/2 fractional quantum Hall state in wide quantum wells is likely to be the one-component Pfaffian state. We hope that this will motivate further experimental studies of this state.

Adsorption of a single Pt atom on graphene: spin crossing between physisorbed triplet and chemisorbed singlet states.

Diffusion Monte Carlo (DMC) calculations have been performed to study the adsorption of a single Pt atom on pristine graphene. We obtain the adsorption energy curves of a single Pt atom adsorbed at three different adsorption sites (bridge, on-top, hollow) as functions of the vertical distance from a graphene surface for both spin singlet and triplet states. The bridge-site adsorption in a singlet spin state is found to be energetically most stable, which is consistent with previous theoretical predictions. As the Pt atom moves away from a graphene surface, spin triplet states are favored over spin singlet states for all three adsorption sites, reflecting that the ground state of an isolated Pt atom is in a spin triplet state. Furthermore, our DMC calculations reveal local-minimum features in the triplet region which is understood to be due to van der Waals interaction between the Pt atom and graphene. This provides a comprehensive understanding for a spin crossing from a physisorbed triplet state to a chemisorbed singlet state in the adsorption process of a single Pt atom on graphene.

The role of high-order electron correlation effects in a model system for non-valence correlation-bound anions.

The diffusion Monte Carlo (DMC), auxiliary field quantum Monte Carlo (AFQMC), and equation-of-motion coupled cluster (EOM-CC) methods are used to calculate the electron binding energy (EBE) of the non-valence anion state of a model (H2O)4 cluster. Two geometries are considered, one at which the anion is unbound and the other at which it is bound in the Hartree-Fock (HF) approximation. It is demonstrated that DMC calculations can recover from the use of a HF trial wave function that has collapsed onto a discretized continuum solution, although larger EBEs are obtained when using a trial wave function for the anion that provides a more realistic description of the charge distribution and, hence, of the nodal surface. For the geometry at which the cluster has a non-valence correlation-bound anion, both the inclusion of triples in the EOM-CC method and the inclusion of supplemental diffuse d functions in the basis set are important. DMC calculations with suitable trial wave functions give EBE values in good agreement with our best estimate EOM-CC result. AFQMC using a trial wave function for the anion with a realistic electron density gives a value of the EBE nearly identical to the EOM-CC result when using the same basis set. For the geometry at which the anion is bound in the HF approximation, the inclusion of triple excitations in the EOM-CC calculations is much less important. The best estimate EOM-CC EBE value is in good agreement with the results of DMC calculations with appropriate trial wave functions.

First-principles quantum Monte Carlo studies for prediction of double minima for positronic hydrogen molecular dianion.

We studied the positron (e+) interaction with the hydrogen molecular dianion H2 2- to form the positronic bound state of [H-; e+; H-] using the first-principles quantum Monte Carlo method combined with the multi-component molecular orbital one. H2 2- itself is unstable, but it was shown that such an unbound H2 2- may become stable by intermediating a positron and forming the positronic covalent bond of the [H-; e+; H-] system [J. Charry et al., Angew. Chem., Int. Ed. 57, 8859-8864 (2018)]. We newly found that [H-; e+; H-] has double minima containing another positronic bound state of [H2; Ps-]-like configuration with the positronium negative ion Ps- at the bond distance approximately equal to the equilibrium H2 molecule. Our multi-component variational Monte Carlo calculation and the multi-component configuration interaction one resulted in the positronic covalent bonded structure being the global minimum, whereas a more sophisticated multi-component diffusion Monte Carlo calculation clearly showed that the [H2; Ps-]-like structure at the short bond distance is energetically more stable than the positronic covalent bonded one. The relaxation due to interparticle correlation effects pertinent to Ps- (or Ps) formation is crucial for the formation of the Ps-A2-like structure for binding a positron to the non-polar negatively charged dihydrogen.

Spin Susceptibility in Neutron Matter from Quantum Monte Carlo Calculations

The spin susceptibility in pure neutron matter is computed from auxiliary field diffusion Monte Carlo calculations over a wide range of densities. The calculations are performed for different spin asymmetries, while using twist-averaged boundary conditions to reduce finite-size effects. The employed nuclear interactions include both the phenomenological Argonne AV8′ + UIX potential and local interactions that are derived from chiral effective field theory up to next-to-next-to-leading order.

Toward a systematic improvement of the fixed-node approximation in diffusion Monte Carlo for solids-A case study in diamond.

While Diffusion Monte Carlo (DMC) is in principle an exact stochastic method for ab initio electronic structure calculations, in practice, the fermionic sign problem necessitates the use of the fixed-node approximation and trial wavefunctions with approximate nodes (or zeros). This approximation introduces a variational error in the energy that potentially can be tested and systematically improved. Here, we present a computational method that produces trial wavefunctions with systematically improvable nodes for DMC calculations of periodic solids. These trial wavefunctions are efficiently generated with the configuration interaction using a perturbative selection made iteratively (CIPSI) method. A simple protocol in which both exact and approximate results for finite supercells are used to extrapolate to the thermodynamic limit is introduced. This approach is illustrated in the case of the carbon diamond using Slater-Jastrow trial wavefunctions including up to one million Slater determinants. Fixed-node DMC energies obtained with such large expansions are much improved, and the fixed-node error is found to decrease monotonically and smoothly as a function of the number of determinants in the trial wavefunction, a property opening the way to a better control of this error. The cohesive energy extrapolated to the thermodynamic limit is in close agreement with the estimated experimental value. Interestingly, this is also the case at the single-determinant level, thus, indicating a very good error cancellation in carbon diamond between the bulk and atomic total fixed-node energies when using single-determinant nodes.

Pairing in SU(6) × SU(2) one-dimensional fermionic clusters

We present diffusion Monte Carlo calculations on the behavior of a mixture of $^{173}\mathrm{Yb}$ and $^{171}\mathrm{Yb}$ fermionic isotopes in a one-dimensional environment. The interaction parameters between species were modeled by $\ensuremath{\delta}$ potentials, whose strengths were taken from their experimental scattering lengths and varied by changing the transversal confinement. This implies a repulsive interaction for the $^{173}\mathrm{Yb}\ensuremath{-}^{173}\mathrm{Yb}$ pair, a strong attractive one for the $^{173}\mathrm{Yb}\ensuremath{-}^{171}\mathrm{Yb}$ set, and a weak attraction between $^{171}\mathrm{Yb}\ensuremath{-}^{171}\mathrm{Yb}$ atoms. Those arrangements were described by a corrected geminal multiplied by the appropriate set of Jastrow functions. We found that, for the same number of strong attractive pairs, the width of the cluster decreases when the number of fermionic species increases, the narrower cluster being the one that includes a mixture of six $^{173}\mathrm{Yb}$ species and two $^{171}\mathrm{Yb}$ spin types.