{'fraction_non_alphanumeric': 0.041120080726538846, 'fraction_numerical': 0.0093970736629667, 'mean_word_length': 4.849133161195131, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of interatomic forces and optimization of molecular geometry with auxiliary-field quantum Monte Carlo

We propose an algorithm for accurate, systematic and scalable computation of interatomic forces within the auxiliary-field Quantum Monte Carlo (AFQMC) method. The algorithm relies on the Hellman-Fenyman theorem, and incorporates Pulay corrections in the presence of atomic orbital basis sets. We benchmark the method for small molecules by comparing the computed forces with the derivatives of the AFQMC potential energy surface, and by direct comparison with other quantum chemistry methods. We then perform geometry optimizations using the steepest descent algorithm in larger molecules. With realistic basis sets, we obtain equilibrium geometries in agreement, within statistical error bars, with experimental values. The increase in computational cost for computing forces in this approach is only a small prefactor over that of calculating the total energy. This paves the way for a general and efficient approach for geometry optimization and molecular dynamics within AFQMC.

We propose an algorithm for accurate, systematic and scalable computation of interatomic forces within the auxiliary-field Quantum Monte Carlo (AFQMC) method. The algorithm relies on the Hellman-Fenyman theorem, and incorporates Pulay corrections in the presence of atomic orbital basis sets. We benchmark the method for small molecules by comparing the computed forces with the derivatives of the AFQMC potential energy surface, and by direct comparison with other quantum chemistry methods. We then perform geometry optimizations using the steepest descent algorithm in larger molecules. With realistic basis sets, we obtain equilibrium geometries in agreement, within statistical error bars, with experimental values. The increase in computational cost for computing forces in this approach is only a small prefactor over that of calculating the total energy. This paves the way for a general and efficient approach for geometry optimization and molecular dynamics within AFQMC.
Calculating interatomic forces in molecules is important for a quantitative understanding of their physical properties. Forces are the basic ingredient in a variety of fundamental studies including molecular dynamics simulations [1], optimization of molecular geometries [2][3][4], computation of vibrational properties [4] and reaction path following [5].
Despite the incredible success within the framework of density functional theory (DFT) [6], much effort has been devoted to developing many-body schemes that can describe more accurately physical and chemical situations where DFT is less reliable, for instance in cases of large electronic correlations and long-range dispersive forces [7][8][9]. Quantum Monte Carlo (QMC) methods are one of the promising alternative many-body approaches, which treat electronic correlations by stochastically sampling of correlated wavefunctions.
The auxiliary-field quantum Monte Carlo (AFQMC) method [10], in particular, has seen rapid development and applications to a wide variety of quantum chemistry and condensed-matter systems [11][12][13][14][15][16][17][18][19][20], yielding state-ofthe-art, benchmark-quality results [21] for the energies of ground and excited states. It is intrinsically parallel, having tremendous capacity to take advantage of petascale (and bebyond) computing resources [22]. The computational cost scales as the third or fourth power of the system size, offering the potential to treat large manyelectron systems.
In this paper, we present an algorithm to directly calculate interatomic forces within AFQMC. As the AFQMC algorithm works by sampling a non-orthogonal Slater determinant space, its formalism allows the adaptation of several key ingredients from other electronic structure or quantum chemistry methods. Combining this with recent advances in the back-propagation technique [32], we achieve an efficient approach for computing atomic forces within AFQMC. In this paper we describe the algorithm, and demonstrate it by full QMC geometry optimization in molecules using a simple steepest descent approach. Internal consistency of the algorithm is verified by measuring the agreement between computed forces and the gradients of the potential energy surface. We benchmark the accuracy of the algorithm by comparison with available high-level quantum chemistry and/or experimental results.
The AFQMC method [10,33,34] reaches the manybody ground state Ψ 0 of a HamiltonianĤ by an iterative process, |Ψ 0 ∝ lim n→∞ exp(−n ∆τĤ) |Ψ T , where ∆τ is a small parameter and, for convenience, we take the initial state Ψ T , which must be non-orthogonal to Ψ 0 , to be a single Slater determinant. The many-body propagator is written as whereB(x) is an independent-particle propagator that depends on the multi-dimensional vector x, and p(x) is a probability distribution function [10]. AFQMC thus represents the many-body wave function in the iteration as an ensemble of Slater determinants, where dΨ denotes integration over the manifold of Slater determinants Ψ [35].
in the limit of large n and m, where ρ m,n (Φ, Ψ) = µ * m (Φ)λ n (Ψ). Formally µ m in the distribution ρ should be the same as λ m ; however, we have used a different symbol to emphasize that in AFQMC the paths leading to it are obtained by the back-propagation algorithm [32,36]. Operationally this can be thought of as generating N w Monte Carlo (MC) path configurations of auxiliary fields, {x 0 , x 1 , · · · , x n−1 ; x n , · · · , x n+m−1 }, in the stochastic sampling of Ψ (m) |Ψ (n) = Ψ T |Ψ (n+m) , and then back-propagate Ψ T | for m steps along the path to obtain Ψ (m) |, (4) In Eq. (3), |Ψ and Φ| in the overlap W (Φ, Ψ) = Φ|Ψ and the local expectation A loc (Φ, Ψ) = Φ|Â|Ψ / Φ|Ψ are defined on the path at "time-slices" (n − 1) and n, respectively. The expectation is evaluated over the MC samples (labeled by w) as Because the propagatorB(x) contains stochastically fluctuating fields, the MC sampling will lead to negative (indeed complex) overlaps W , which will cause the variance of this estimator to grow exponentially with the projection times, n and m. Control of this phase problem is achieved by the introduction of an importance sampling transformation and a generalized gauge condition to constrain the random walks [10].
This framework eliminates the phase problem, at the cost of modifying the distribution ρ m,n . The constraint introduces a bias in λ n . For µ m an additional subtlety arises because the constraint on the path is imposed in the time-reversed direction [34]. Previous studies in a variety of systems have shown that the bias from the constraint tends to be small, in both models [37] and realistic materials [21], making AFQMC of the most accurate many-body approaches for general interacting fermion systems.
To address the problem of computing atomic forces, let us consider a molecule comprising N n ions, whose spatial positions R = (R 1 . . . R Nn ) define a molecular geometry. Given a molecular geometry R, ground-state expectation values of physical observables are given by ForÂ(R) =Ĥ(R), Eq. (6) gives the potential energy surface (PES), E(R). The gradient of (6) is in general given by In the case of the PES, the first two terms vanish for sufficiently large m and n. However, the expression does not reduce to Hellmann-Feynman theorem [38,39], because Ψ and Φ are represented by incomplete atom-centered basis sets which depend on R. Performing the partial derivative in Eq. (3), one obtains with The terms in Eq. (9) are detailed in the Appendix, and encompass a Pulay's correction [2] to take into account the explicit dependence of the one-electron Green's function on atomic coordinates R.
Ground-state expectation values of the form in Eq. (8) can be computed accurately with the back-propagation algorithm, originally formulated for lattice models of correlated electrons [40], generalized to weakly correlated systems with complexB(x) [36], and recently adapted to molecular systems with a "path restoration" (BP-PRes) technique [32]. The computational cost of evaluating F α loc (Φ, Ψ; R) is similar to that of computing the local energy. Additional speedup can be achieved by exploiting the local nature of the atomic forces and the sparsity of the Hamiltonian matrix elements. In larger molecules and solids, the ability to evaluate all forces simultaneously (as opposed to separate total energy calculations for each force component) is crucial, and paves the way for geometry optimizations.
We note that the formalism and estimator detailed above hold for any other method operating in an overcomplete manifold of non-orthogonal Slater determinants. Differences arise in the generation of the distributions µ and λ which, as discussed above, need not be differentiated to compute the atomic forces.
In Fig. 1 we assess the internal consistency and accuracy of our method using the symmetric stretching of the CH bond in CH 4 as a test case. The computed gradients, across the entire bondlength range, are in agreement with results from explicit (numerical) derivative of the PES within statistical uncertainties. As the inset in the upper panel illustrates, it is crucial to have BP in order for Eq. (8) to achieve an accurate representation of Eq. (7). Performing calculations with increasingly large basis sets, we observe convergence of the equilibrium bondlength  As a further test, we perform a geometry optimization in Fig. 2 for the H 2 O molecule using the steepestdescent algorithm. The computed equilibrium geometry is in agreement with the global minimum of the PES, obtained from AFQMC calculations of the ground-state energy on a dense mesh of points. In Table I we list the optimized geometries for different basis sets, together with those from the coupled-cluster CCSD(T) for comparison and benchmark. The statistical error bars on the QMC geometry are obtained by averaging over configurations of {R i } with converged and statistically compatible energies, for example, at the cc-pVDZ level, i = 6 . . . 11 (see bottom panel in Fig. 2). From Table I, convergence to the CBS limit is seen at the cc-pVQZ level, where the AFQMC results are in good agreement with experiment  [41].
In Fig. 3 we move to the more challenging case of ethane. Geometries respecting the D 3 h symmetry of the molecule can be expressed in terms of 3 parameters, as sketched in the inset. As the number of geometric parameters grows, it quickly becomes impractical to compute the entire PES with AFQMC as was done in Fig. 2, and forced-based methods become essential. The optimized parameters, with the modest but realistic cc-pVDZ basis, are within 0.02 a.u. (or 1%) of the experimental equilibrium geometry, as seen in Table II. The excellent agreement between the AFQMC geometry and those from high-level QC methods at this basis level suggests that the discrepancy with experimental data is likely a result of basis set incompleteness, which is quantitatively confirmed by a cc-pVTZ calculation. The ground-state energies computed at the optimized geometry and at the experimental equilibrium geometry are in agreement with each other: E(R op ) − E(R eq ) = −0.5(4) mHa. As a last example, we optimize the molecular geometry of nitric acid. We use increasingly larger basis sets to reach the continuum limit near the optimal geometry. As shown in Fig. 4, there is a large residual error in the STO-6G basis set, and the optimized geometry is significantly different from experiment. As more realistic basis sets are employed, systematically improved results are obtained, and the computed geometry approaches the experimental equilibrium geometry. At the cc-pVTZ level, AFQMC results are in agreement with experiment to within 0.01 a.u. Results with the cc-pVDZ basis set are compared with several other QC methods. Results from cc-pVTZ indicate reasonable convergence with respect to basis set, and good agreement with experiment. DFT and CCSD(T) data were computed using the NWChem software [42]. Geometries are expressed in terms of the parameters p1, p2, p3 sketched in Fig. 3. We have demonstrated the direct computation of interatomic forces and molecular geometry optimization within AFQMC. We proposed an internally consistent, numerically stable and computationally efficient algorithm based on the Hellman-Feynman theorem and Pulay's corrections. Results from a first application are presented. Accurate forces are obtained using simple RHF Slater determinant as constraining trial wave function. These results pave the way for systematic geometry optimization and potentially molecular dynamics using one of the most accurate many-body methods with low-power computational scaling.
A variety of future directions are possible, including both applications and further generalization and improvement of the algorithm. More refined Ansätze [14][15][16] for the trial wave function can be adopted straightforwardly for more challenging molecules. Interfacing with better optimization strategies and improving the efficiency of the force computation algorithm with AFQMC itself can both lead to major increases in the capability of the method. Alternative representations [19,43] of the Hamiltonian operator should be explored. Applications to a variety of molecular systems are within reach including those containing post-second-row elements. Geometry optimization in crystalline solids is being investigated. We . . M } are the basis functions and C Ψ,σ is a M × N σ matrix that contains the orbital coefficients. The Slater determinant |Φ is similarly defined in terms of orbitals |φ lσ and matrix C Φ,σ . In Eq. (A1), the overlap is given by where S(R) ij = χ i |χ j is the overlap matrix of the AO basis set. The local energy is given by where H 0 (R) is a constant from the internuclear repulsion, the matrix h(R) is the one-body part of the Hamiltonian in the AO basis (h(R) ij = χ i |h 1 (R)|χ j ), the matrices L γ (R) arise from decomposing the twobody part of the Hamiltonian (e.g., via the Cholesky or density-fitting decomposition [43,44], (ik|jl) = γ L γ (R) ik L γ (R) jl ), and the (spin-dependent) oneelectron Green's function is which obeys the idempotence relation The gradient of W can be obtained from Jacobi's formula for the derivative of the determinant, The gradient of E loc (Φ, Ψ; R) is readily computed by recalling the cyclicity of the trace, and that differentiation of (A4) yields complying with the idempotence relation (A5) [2]. The gradient ∂ Rα G σ (R) incorporates Pulay's corrections in the framework of AFQMC, which are important for an accurate evaluation of forces when atom-centered Gaussian orbitals are used as basis for the one-electron Hilbert space. The presence of Pulay's corrections is reminiscent of the covariant derivatives that appear in fiber bundle theory and differential geometry [45]. This reminiscence is not accidental: the purpose of corrective terms in covariant derivatives is to preserve some properties of the space on which they act. For example, in gauge theories corrections preserve the gradient of the wavefunction under gauge transformations and, on Riemann manifolds, connections preserve parallelism. In this framework, Pulay's correction preserve the idempotency (A5) of the one-particle reduced Green's function under changes in molecular geometries.
When m = n = 0, ρ m,n (Φ, Ψ) = δ(Φ − Ψ T )δ(Ψ − Ψ T ), so that E(R) = E loc (Ψ T , Ψ T ; R) and the estimator reduces to the familiar Hartree-Fock expression [2]. In this work, the derivatives ∂ Rα S, ∂ Rα h and ∂ Rα L γ were all obtained by numerical differentiation. A finite difference step size of 10 −4 a.u. was used, and it was verified that using a smaller step does not change the results within statistical error bars.