Sign Learning Kink-based (SiLK) Quantum Monte Carlo for molecular systems

The development of accurate and computationally tractable ab initio methods for studying correlated electronic systems ranging from single molecules to bulk materials¹ is an area of wide interest. Feynman’s path integral formulation of quantum mechanics² has long attracted attention due to its ability to include exact correlation and finite temperature effects, as well as providing a method that can simultaneously treat electronic and geometric degrees of freedom. The path integral formulation is one of a number of methods commonly referred to as Quantum Monte Carlo (QMC)-based algorithms.³ 

In general, the use of QMC-based algorithms are hindered by the so-called minus sign problem in which the fluctuating sign of the fermionic density matrix leads to statistical errors that scale exponentially with inverse temperature and system size. The minus sign problem^4,5 remains a great challenge in condensed matter physics and quantum chemistry.

In quantum chemistry, there are a number of methods used to include electron correlation. Commonly used methods include density functional theory (DFT),⁶ configuration interaction (CI),⁷ many body perturbation theory (MBPT),^8–10 and coupled cluster (CC).^11–18 The CC method has been regarded as the “gold” standard.¹¹ These approaches, while very useful, have well-known deficiencies such as the approximate inclusion of correlation (DFT), size inconsistency (truncated CI, such as with single and double excitations or with single, double, and triple excitations), or non-variational energies (CC). Therefore, it is important to investigate alternative approaches.

There are three major numerical methods used to study strongly correlated many body systems. These are exact diagonalization, density matrix renormalization group (DMRG),¹⁹ and QMC. Exact diagonalization is only feasible for small systems since it scales exponentially with the system size. DMRG has become useful for certain classes of molecules with an order of 50 strongly correlated electrons.^20–22 

The Monte Carlo method was first introduced and developed by Fermi, Teller, and Metropolis.^23–25 QMC, unlike exact diagonalization and DMRG, is a scalable method that can be applied to multi-dimensional lattice systems. However, QMC does have the minus sign problem in fermionic and frustrated quantum systems.

A variety of methods have been proposed to alleviate the minus sign problem in QMC. These include auxiliary field Monte Carlo,²⁶ shifted contour auxiliary field Monte Carlo,²⁷ and fixed node diffusion Monte Carlo.^3,28 More recently, a resummation path integral approach,^29,30 which is similar to the SiLK method, phaseless auxiliary-field QMC,³¹ a finite temperature version of diffusion Monte Carlo,^32,33 and full configuration interaction QMC^34,35 have been developed.

The Sign-Learning Kink (SiLK) QMC algorithm originally developed by Hall^36,37 can be used to overcome the minus sign problem. SiLK has previously been used to study the 3 × 3 Hubbard model and atoms using a small basis set.^36,37 An approximate version of the method has been used to study small molecules.³⁸ This method uses a novel learning process to overcome the minus sign problem.

The goal of this work is to investigate the ability of the SiLK method to reduce the sign problem and accurately calculate potential energy surfaces in model systems with relatively small basis sets. Investigation of the scalability of the method is left for future work. Therefore, SiLK QMC calculations are performed on H₂O, N₂, and F₂ at a number of different geometries. The results of the calculations are compared to the results to exact diagonalization and a variety of quantum chemistry methods and demonstrate that the SiLK method is accurate and that it reduces the minus sign problem for all geometries.

Assume there are a finite set of states composed of Slater determinants $α i$ formed from orthogonal, one electron spin orbitals. With Hamiltonian, H, and β = 1/k_BT, the canonical partition function Q can be written as

Using

and the identity

the partition function becomes

P is introduced as a discretization variable that allows for the evaluation of the matrix elements by expanding the exponential, vide infra. For a given set of {α_ji}, some of the matrix elements in Eq. (4) may be diagonal. Thus, terms appearing in the summand may be classified by the number of off-diagonal matrix elements. In the SiLK formalism, off-diagonal matrix elements are referred to as kinks. By analytically summing over the diagonal matrix elements in Eq. (4), we obtain a kink-based version of the partition function. Defining

the result of this analytical summation is³⁶ 

where n is the number of kinks, m is the number of distinct α_j’s in a given set of states with n kinks, s_j is the number of times the state α_j appears in a given set of states with n kinks, and

where S may be evaluated recursively. Due to the derivatives in Eq. (9), it is possible for S to be negative. In addition, the off-diagonal matrix elements t_jk,jk+1 can also be negative.

Fig. 1 depicts the types of kink configurations that appear in this sum over states. The top figure without kinks corresponds to the case where only diagonal matrix elements occur such that j₁ = j₂, …. The second case contains two “kinks” where two identical off-diagonal matrix elements are introduced. This so-called kink expansion was used in condensed matter physics by Anderson³⁹ and later in the chemical physics literature by Wolynes.⁴⁰ 

The first term is non-negative. The n = 2 is also non-negative since the off-diagonal matrix elements appear as |t_j1,j2|² and with s₁ and s₂ = 1, S > 0. Therefore the sign problem is due to terms with n ≥ 3. The SiLK method uses a learning algorithm to construct new states as linear combinations of the initial {α_i} states that minimize the magnitude of the contributions from terms with n ≥ 3 and thereby reduces or eliminates the sign problem.

Eq. (8) has the form of a grand canonical partition function and thus Monte Carlo methods may be used to evaluate the partition function and its properties. Writing this equation as

a Monte Carlo simulation will involve sampling different states and inserting and removing kinks. The average energy of the system can be evaluated using $〈E〉=− d d β lnQ$⁠,

where kink configurations are sampled from |ρ|.

Simulations are performed in two stages. The first stage is a “learning” period and is used to construct an improved description of the states of the system. We choose the lowest energy Hartree-Fock state as the initial state. As the grand canonical simulation proceeds, additional states are inserted and removed and a list is maintained of states that have appeared. At fixed intervals (30 iterations in our calculations) or when the number of kinks present at the end of a Monte Carlo pass exceeds a specified number (9 in our case), the Hamiltonian is diagonalized in the sub-space of the states that have appeared since the last diagonalization (or the start of the simulation). The state is then set to the lowest energy state (a zero-kink configuration) and the simulation is continued. The learning period ends when there are only zero and two kinks configurations present for an extended number of Monte Carlo passes. At this point, the expectation is that the partition function will be dominated by kink configurations with a small number of kinks (dominated by configurations with 0 or 2 kinks) as the current set of states will better approximate the ground state of the system than the initial ones. As currently implemented, the learning stage can be thought of as using the simulation to construct CI states. In the present work, the learning period ranged from 8000 to 119 000 passes.

The second stage in the simulation is the data acquisition during which the states are not modified and the grand canonical simulation proceeds in the standard way. If the number of kinks increases dramatically during this stage (perhaps due to the simulation exploring a previously unexplored region of phase space), a diagonalization is performed and the second stage is restarted. In the calculations, where additional diagonalizations were performed, the diagonalization made an insignificant change in the ground and excited state energies. Between 1000 and 2000 Monte Carlo passes were used in the second stage.

The Monte Carlo algorithm consists of two types of moves: change of state and insertion/removal of states. The former is performed in the standard way using the Metropolis algorithm. The latter uses the Metropolis algorithm as follows. A potential new kink configuration c′ is sampled based on the current kink configuration c using the normalized conditional probability T(c′|c) and accepted with probability $A ( c ′ | c ) =min 1 , ρ ( c ′ ) T ( c | c ′ ) ρ ( c ) T ( c ′ | c )$⁠. If there are n states in the current kink configuration, there are n + 1 places to insert a new state into the kink configuration (as state 1, state 2, …, state n + 1). There are n ways to remove a state. We set

with the probability of removing the state at location k in the list of states

and with the probability of adding α_j at location k,

with D(c′|c) the normalization for the probability,

The acceptance probability is then

We use the SiLK algorithm to calculate the energies of H₂O, N₂, and F₂ at selected bond lengths and bond angles for H₂O. The Cartesian Gaussian DZ basis set^41–44 is used in all calculations. Computer memory constraints imposed by the current SiLK implementation dictates the number of determinants that can be used in the calculations. The reason is that at present the CI coefficients for the ground and excited states must be stored. Future work will focus on alleviating the memory issues. We therefore use either the full vector space of determinants generated by all possible excitations of the Hartree-Fock determinant (Full Configuration Interaction, FCI) or the restricted vector spaces generated by the HF determinant and either all possible single and double excitations (SD) or all possible single, double and triple excitations (SDT) of the HF determinant. In all cases, a comparison of the SiLK method to the exact result within the vector space is made to assess, as mentioned in the Introduction, the ability of SiLK to provide accurate results and alleviate the sign problem. At each geometry, a Hartree-Fock computation using the NWChem ab initio package⁴⁵ is used to generate the initial molecular orbitals from which the determinants are created. Determinants corresponding to excited states are generated by excitations of all molecular orbitals except the core orbitals (the frozen core approximation). Symmetry is used to restrict the determinants to those with the same symmetry as the ground Hartree-Fock state. For the calculations presented here, we use C_2v spatial symmetry for H₂O and D_2h spatial symmetry for N₂ and F₂, respectively. We use T = 1 K (β = 3 × 10⁵ hartree⁻¹). Exact energies are obtained by numerical diagonalization for the SD and SDT vector spaces. A series of calculations with increasing values for P were performed until a convergence in the energy was obtained. The values of P chosen for the reporting of data ranged from 2 × 10⁷ to 2 × 10¹⁰. The FCI calculations for H₂O is performed using Molpro.^46,47

The ability of SiLK to address the sign problem is evaluated by following the evolution of the sign (for clarity averaged over every 20 Monte Carlo steps) during the course of the learning period. Representative of the results from the different molecules is the average sign for water at the minimum energy FCI geometry.⁴⁸ In this calculation, the maximum number of states included in a diagonalization is limited to 50 (the entire vector space had a dimension of 128 829). The coarse-grained sign is shown in Fig. 2. The sign fluctuates significantly for roughly the first 1500 diagonalizations, but after approximately 1600 diagonalizations it remains 1.0. In the upper panel of Fig. 2 the number of states involved in each diagonalization is shown. After 1600 diagonalizations, the coarse-grained sign remains at 1.0 even though the number of states involved in subsequent diagonalizations is approximately 20. This indicates that kinks are being introduced during the Monte Carlo process but these are not affecting the sign. The number of kinks averaged over the kink configurations between 2 successive diagonalizations ranged from roughly 5 at the beginning of the learning period to roughly 2.5 at the end of the learning period. An examination of the kink configurations after the learning period found that the configurations contain either zero and two kinks and therefore the average sign is 1.0.

Accurate calculations of potential energy surfaces are important in understanding reaction energetics and rates. The ability of SiLK and other quantum chemical methods to calculate potential energy surfaces is assessed for H₂O, F₂, and N₂. The goal of a successful method is to achieve the accuracy required to describe energetic differences encountered in chemical processes such as bond-breaking/bond-forming reactions and reaction activation energies. This so-called “chemical accuracy” is approximately 0.1–1 kcal/mol ≈10⁻³ to 10⁻⁴ Hartree.⁴⁹ 

Several versions of truncated CC are used in this work. The CCSD method uses single and double excitations.¹⁷ The CCSDT method uses single, double, and triple excitations.⁵⁰ The CCSD(T) method uses single, double, and non-iterative inclusion of perturbative triples¹⁸ and is considered to be the “gold standard” of ab initio quantum chemistry. The MRCCSD(T)(2,2) and MRCCSD(T)(4,4) methods are multi-reference CC (MRCC) methods with single, double, and non-iterative inclusion of perturbative triples.⁵¹ A (2,2) calculation uses 2 electrons and 2 orbitals (one occupied and one virtual) to generate the model space for the MRCC calculation and a (4,4) calculation uses 4 electrons and 4 orbitals (two occupied and two virtual) to generate the model space for the MRCC calculation. We also use second order MBPT (MBPT(2)). The NWChem software package is used to perform all standard ab initio calculations.

The H₂O molecule is used to assess the ability of SiLK to describe the variation of energy with bond length and bond angle in two separate calculations, one in which the bond length is varied and another in which the bond angle is varied. Fig. 3 displays the energy and its absolute error as a function of bond length at fixed bond angle of 110.565° as calculated by different methods. The SiLK method has an absolute error of 10⁻⁵ hartree over the range of bond lengths studied, which is well below the desired chemical accuracy. At the minimum energy geometry (bond length = 1.8434 bohrs), the exact energy is −76.144 552 99 hartree and the energy of the lowest energy SiLK state is −76.144 546 90 hartree, an error of ≈6 × 10⁻⁶ hartree. Therefore, the SiLK procedure found an excellent approximation to the exact ground state. SiLK is approximately one order of magnitude more accurate than the most accurate of the other methods in the comparison. Notably, SiLK is accurate at the longer bond lengths (roughly two orders of magnitude more accurate than any other method) which is crucial to a description of bond dissociation and bond breaking processes. None of the other methods (except MRCCSD(T)(4,4)) achieves chemical accuracy over the entire range of bond lengths studied.

Then we use these methods to calculate the energy as a function of bond angle for the H₂O molecule with bond length 1.843 45 bohrs. Fig. 4 displays the energies and absolute errors for bond angles ranging from 95° to 125°. The SiLK method is approximately two orders of magnitude more accurate than the most accurate of the other methods. All methods except MBPT(2) and CCSD achieve chemical accuracy.

It is also instructive to consider calculations restricted to just single and double excitations as sometimes computations based on such restricted vector spaces can yield useful results using significantly fewer computational resources. Therefore, the SiLK algorithm is used to calculate the energies and absolute errors of the H₂O molecule as a function of bond length and bond angle. Figs. 5 and 6 show their comparison with exact results. The SiLK method is able to reproduce the exact results to 10⁻⁵ hartree, well within chemical accuracy.

N₂ has a triple bond, which provides a challenging test for ab initio methods due to its large electronic correlation.⁵² Due to memory limitations, the SiLK calculations were restricted to the Hartree-Fock determinant plus either the SD and SDT vector spaces. Fig. 7 shows that the SiLK QMC results converge to the exact result over a wide range of bond lengths. At the minimum energy geometry (bond length = 2.168 bohrs), the exact energy is −109.085 809 5 hartree and the energy of the lowest energy SiLK state is −109.085 809 4 hartree, an error of ≈1 × 10⁻⁷ hartree. Therefore, the SiLK procedure found an excellent approximation to the exact ground state. The results demonstrate that the SiLK method is suitable for multi-reference systems such as N₂ where more than a single determinant is strongly coupled in the ground electronic state.

Electron correlations are difficult to include in the simulation of F₂ as many determinants contribute small but important contributions to the total energy.⁵³ This phenomenon is often referred to as dynamic correlation.⁵⁴ Therefore, the SiLK method is applied to the F₂ molecule. As with Nitrogen, due to memory limitations, the SiLK calculations is limited to the SD and SDT vector spaces. Fig. 8 shows that the SiLK QMC results converge to the exact results and demonstrate that SiLK is capable of accurately including dynamic correlation. At the minimum energy geometry (bond length = 2.868 16 bohrs), the exact energy is −198.949 416 9 hartree and the energy of the lowest energy SiLK state is −198.949 416 9 hartree, an error of <1 × 10⁻⁷ hartree. Therefore, the SiLK procedure found an excellent approximation to the exact ground state.

It is beyond the scope of this work to make a thorough analysis of the scaling of the SiLK method as the memory requirements of the current implementation of the SiLK method prohibit the use of a wide range of basis set size. However, it is important to assess the scaling of the SiLK algorithm with the size of the basis set and vector space. No truncation methods, such as a truncation in the space of the single-particle density matrix,⁵⁵ will improve the efficiency of the algorithm for certain systems. Therefore, a scaling analysis is presented for the current work using the relatively limited size of the vector spaces. If the SiLK algorithm increases too quickly with the size of the vector space, the computational requirements for the SiLK method will make its use in the present form intractable. The dependence of the length of the learning period on the size of the vector space (number of determinants) is presented in Fig. 9. As expected, there is an increase in the size of the learning period. However, a wider range of vector space sizes is necessary to fully understand and quantify the scaling behavior.

The minus sign problem in Quantum Monte Carlo simulations of frustrated or correlated electronic systems is a challenging problem. It has even been suggested that a general solution of this problem is NP-complete.⁴ Therefore, one should not expect an effective solution for all the Monte Carlo simulations which have the minus sign problem. In this paper, we demonstrate that SiLK QMC can reduce the minus sign problem by using a learning stage that includes a diagonalization procedure. In this paper, we demonstrate that the energies obtained by the SiLK QMC match the results from exact diagonalization and surpass the accuracy obtained using other quantum chemistry methods, particularly for geometries relatively far from equilibrium. In addition, SiLK can be applied to systems that require a multiple reference state approach. An intriguing possibility for future work is to use the SiLK learning procedure in combination with other QMC algorithms to reduce the minus sign problem.

As the learning stage progresses, the states become more complicated linear combinations of determinants so that more evaluations of matrix elements are required, thereby increasing the computational expense. However, at the same time, the number of non-zero matrix elements between these states decreases. So, further optimization is possible by storing often-needed matrix elements in memory. For example, storing the off-diagonal matrix elements between the ground and excited states yields a large speed up, since these are the only matrix elements required once the learning stage reaches the point where mostly zero and two kink configurations appear in the simulation. It is also possible to halt the learning stage at an earlier point, when the ground SiLK state is not as accurate an approximation to the exact ground state, but when the sign problem is alleviated but not eliminated and relies on the Monte Carlo sampling to provide the exact energy. This would reduce the computational effort required to evaluate the matrix elements since fewer diagonalizations will have occurred. An investigation of the efficacy of a shorter learning period is left for future work.

The SiLK method requires the knowledge of the off-diagonal matrix elements of the Hamiltonian. As the size of the system increases, the number of off-diagonal matrix elements increases factorially and it is not possible to store the matrix elements or the CI coefficients for the ground and excited states that would allow for on-the-fly evaluation of the matrix elements. As such, without a procedure to accurately truncate the number of determinants used to describe the ground and excited state wavefunctions, the use of the all possible determinants in a SiLK calculation will be limited to relatively small systems. However, SiLK can certainly be used when determinants are restricted to, for example, single and double excitations. Such truncated sets of determinants are often sufficient for the study of chemical systems. In cases where restrictions to single and double excitations are not sufficient, more sophisticated methods of truncation, such as the one developed by Maurits,⁵⁵ will be needed.

The SiLK QMC is a versatile method to calculate the ground state energy of molecular systems. Since the path integral formulation uses the canonical partition function it is possible to use the SiLK method to simulate the motion of the atoms at a finite temperature. Future work will investigate the use of the SiLK method in finite temperature simulations.

This work is funded by the NSF EPSCoR LA-SiGMA project under Award No. EPS-1003897 and the US Department of Energy under Contract No. DE-AC06.76RLO-1830. R.W.H. acknowledges support from Dominican University of California’s Lillian L. Y. Wang Yin, Ph.D. Endowed Professorship of Chemistry. This work used the high performance computing resources provided by Louisiana State University (http://www.hpc.lsu.edu) and the computational resources of the Environmental Molecular Sciences Laboratory at Pacific Northwest National Laboratory (PNNL), which is sponsored by the Department of Energy Office of Biological and Environmental Research.