Sign learning kink-based (SiLK) quantum Monte Carlo for molecular systems
Xiaoyao Ma, Randall W. Hall, Frank Loffler, Karol Kowalski, Kiran Bhaskaran-Nair, Mark Jarrell, Juana Moreno
aa r X i v : . [ phy s i c s . c o m p - ph ] D ec Sign learning kink-based (SiLK) quantum Monte Carlo for molecular systems
Xiaoyao Ma, Randall W Hall,
2, 3
Frank Löffler, Karol Kowalski, KiranBhaskaran-Nair,
1, 4
Mark Jarrell,
1, 4 and Juana Moreno
1, 4 Department of Physics and Astronomy, Louisiana State University, Baton Rouge,LA 70803, USA Department of Natural Sciences and Mathematics, Dominican University ofCalifornia, San Rafael, CA 94901, USA Department of Chemistry, Louisiana State University, Baton Rouge, LA 70803,USA Center for Computation & Technology, Louisiana State University, Baton Rouge,LA 70803, USA William R. Wiley Environmental Molecular Sciences Laboratory, Battelle,Pacific Northwest National Laboratory, Richland, Washington 99352,USA (Dated: 10 September 2018)
The Sign Learning Kink (SiLK) based Quantum Monte Carlo (QMC) method is usedto calculate the ab initio ground state energies for multiple geometries of the H O,N , and F molecules. The method is based on Feynman’s path integral formulationof quantum mechanics and has two stages. The first stage is called the learning stageand reduces the well-known QMC minus sign problem by optimizing the linear com-binations of Slater determinants which are used in the second stage, a conventionalQMC simulation. The method is tested using different vector spaces and comparedto the results of other quantum chemical methods and to exact diagonalization. Ourfindings demonstrate that the SiLK method is accurate and reduces or eliminates theminus sign problem.PACS numbers: 05,31,82 1 . INTRODUCTION The development of accurate and computationally tractable ab initio methods for study-ing correlated electronic systems ranging from single molecules to bulk materials is an areaof wide interest. Feynman’s path integral formulation of quantum mechanics has long at-tracted 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 degreesof freedom. The path integral formulation is one of a number of methods commonly referredto as Quantum Monte Carlo (QMC)-based algorithms .In general, the use of QMC-based algorithms are hindered by the so-called minus signproblem in which the fluctuating sign of the fermionic density matrix leads to statisticalerrors that scale exponentially with inverse temperature and system size. The minus signproblem remains a great challenge in condensed matter physics and quantum chemistry.In quantum chemistry, there are a number of methods used to include electron corre-lation. Commonly used methods include density functional theory (DFT) , configurationinteraction (CI) , many body perturbation theory(MBPT) and coupled cluster(CC) .The CC method has been regarded as the “gold” standard . These approaches, while veryuseful, 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 importantto investigate alternative approaches.There are three major numerical methods used to study strongly correlated many bodysystems. These are exact diagonalization, density matrix renormalization group (DMRG) ,and QMC. Exact diagonalization is only feasible for small systems since it scales exponen-tially with the system size. DMRG has become useful for certain classes of molecules withan order of 50 strongly correlated electrons .The Monte Carlo method was first introduced and developed by Fermi, Teller, andMetropolis . QMC, unlike exact diagonalization and DMRG, is a scalable method thatcan be applied to multi-dimensional lattice systems. However, QMC does have the minussign 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 ,2nd fixed node diffusion Monte Carlo . More recently, a resummation path integral ap-proach , which is similar to the SiLK method, phaseless auxiliary-field QMC , a fi-nite temperature version of diffusion Monte Carlo , and full configuration interactionQMC have been developed.The Sign-Learning Kink (SiLK) QMC algorithm originally developed by Hall can beused to overcome the minus sign problem. SiLK has previously been used to study the3 × . An approximate version of themethod has been used to study small molecules . This method uses a novel learning processto overcome the minus sign problem.The goal of this work is to investigate the ability of the SiLK method to reduce the signproblem and accurately calculate potential energy surfaces in model systems with relativelysmall basis sets. Investigation of the scalability of the method is left for future work. There-fore, SiLK QMC calculations are performed on H O, N , and F at a number of differentgeometries. The results of the calculations are compared to the results to exact diagonaliza-tion and a variety of quantum chemistry methods and demonstrate that the SiLK methodis accurate and that it reduces the minus sign problem for all geometries. II. SILK FORMALISM AND ALGORITHMA. SiLK Formalism
Assume there are a finite set of states composed of Slater determinants { α i } formed fromorthogonal, one electron spin orbitals. With Hamiltonian, H , and β = 1 /k B T , the canonicalpartition function Q can be written as Q = T r n e − βH o = X j D α j | e − βH | α j E . (1)Using e − βH = ( e − βH/P ) P , (2)and the identity 1 = X j i | α j i ih α j i | , (3)3he partition function becomes Q = X j ,j ,...,j P * α j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) exp − βP H ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + * α j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) exp − βP H ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + · · · * α j P (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) exp − βP H ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + . (4) P is introduced as a discretization variable that allows for the evaluation of the matrixelements by expanding the exponential, vide infra . For a given set of { α j i } , some of thematrix elements in Eqn. 4 may be diagonal. Thus, terms appearing in the summand may beclassified by the number of off-diagonal matrix elements. In the SiLK formalism, off-diagonalmatrix elements are referred to as kinks. By analytically summing over the diagonal matrixelements in Eqn. 4, we obtain a kink-based version of the partition function. Defining x j ≡ * α j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) exp − βP H ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + ≈ * α j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − βP H (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + + O (cid:18) P (cid:19) ≈ exp − βP h α j | H | α j i ! + O (cid:18) P (cid:19) , (5) t ij ≡ * α i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) exp − βP H ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + ≈ * α i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − βP H (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α j + + O (cid:18) P (cid:19) ≈ exp − βP h α i | H | α j i ! − O (cid:18) P (cid:19) , when i = j, (6)the result of this analytical summation is : Q = lim P →∞ Q ( P ) , (7) Q ( P ) = X j x Pj + P X n =2 Pn n Y i =1 X j i n Y k =1 t j k ,j k +1 ! S ( { x j } , n, m, { s j } ) , (8)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 S ( { x j } , n, m, { s j } ) = m Y j =1 s j − d s j − dx s j − j x s j − j m X l =1 x P − n + m − l Q i = l ( x l − x i ) , (9)4here S may be evaluated recursively. Due to the derivatives in Eqn. 9, it is possible for S to be negative. In addition, the off-diagonal matrix elements t j k ,j k +1 can also be negative.Fig. 1 depicts the types of kink configurations that appear in this sum over states. Thetop figure without kinks corresponds to the case where only diagonal matrix elements occursuch that j = j , · · · . The second case contains two “kinks” where two identical off-diagonalmatrix elements are introduced. This so-called kink expansion was used in condensed matterphysics by Anderson and later in the chemical physics literature by Wolynes . FIG. 1. Examples of different kink configurations that occur in Eqn. 4 when P = 5. Horizontal linescorrespond to diagonal matrix elements and slanted lines correspond to off-diagonal matrix elementsthat are referred to as kinks. The lines at the beginning and end of a kink configuration wrap arounddue to the
T race operation required by the partition function. The zero kink configuration containsmatrix elements for a single state, the two kink configuration contains matrix elements for just twostates, etc. Note that the number of kinks and the number of states are not necessarily equal toeach other as it is seen in the four kink configuration.
The first term is non-negative. The n = 2 is also non-negative since the off-diagonalmatrix elements appear as | t j ,j | and with s and s = 1, S >
0. Therefore the sign problemis due to terms with n ≥
3. The SiLK method uses a learning algorithm to construct newstates as linear combinations of the initial { α i } states that minimize the magnitude of thecontributions from terms with n ≥ Q ( P ) = P X n =0 X { α i } ρ ( n, { α i } ) (10)a Monte Carlo simulation will involve sampling different states and inserting and removingkinks. The average energy of the system can be evaluated using < E > = − ddβ ln Q , < E > = − P Pn =0 P { α i } ddβ ρ ( n, { α i } ) P Pn =0 P { α i } ρ ( n, { α i } )= − P Pn =0 P { α i } | ρ ( n, { α i } ) || ρ ( n, { α i } ) | ddβ ρ ( n, { α i } ) P Pn =0 P { α i } | ρ ( n, { α i } ) || ρ ( n, { α i } ) | ρ ( n, { α i } )= − < | ρ | ddβ ρ > | ρ | < sign ( ρ ) > | ρ | (11)where kink configurations are sampled from | ρ | . B. SiLK algorithm
Simulations are performed in two stages. The first stage is a “learning” period and isused to construct an improved description of the states of the system. We choose the lowestenergy 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 haveappeared. At fixed intervals (30 iterations in our calculations) or when the number of kinkspresent at the end of a Monte Carlo pass exceeds a specified number (9 in our case), theHamiltonian is diagonalized in the sub-space of the states that have appeared since the lastdiagonalization (or the start of the simulation). The state is then set to the lowest energystate (a zero-kink configuration) and the simulation is continued. The learning period endswhen there are only zero and two kinks configurations present for an extended number ofMonte Carlo passes. At this point, the expectation is that the partition function will bedominated by kink configurations with a small number of kinks (dominated by configurationswith 0 or 2 kinks) as the current set of states will better approximate the ground state of thesystem than the initial ones. As currently implemented, the learning stage can be thoughtof as using the simulation to construct configuration interaction (CI) states. In the presentwork, the learning period ranged from 8,000 to 119,000 passes.6he second stage in the simulation is the data acquisition during which the states are notmodified and the grand canonical simulation proceeds in the standard way. If the numberof kinks increases dramatically during this stage (perhaps due to the simulation exploring apreviously unexplored region of phase space) a diagonalization is performed and the secondstage is restarted. In the calculations where additional diagonalizations were performedthe diagonalization made an insignificant change in the ground and excited state energies.Between 1,000 and 2,000 Monte Carlo passes were used in the second stage.The Monte Carlo algorithm consists of two types of moves: change of state and inser-tion/removal of states. The former is performed in the standard way using the Metropolisalgorithm. The latter uses the Metropolis algorithm as follows. A potential new kink con-figuration c ′ is sampled based on the current kink configuration c using the normalizedconditional probability T ( c ′ | c ) and accepted with probability A ( c ′ | c ) = min h , ρ ( c ′ ) T ( c | c ′ ) ρ ( c ) T ( c ′ | c ) i . Ifthere are n states in the current kink configuration, there are n + 1 places to insert a newstate into the kink configuration (as state 1, state 2, ..., state n + 1). There are n ways toremove a state. We set T ( c ′ | c ) = T remove ( c ′ | c ) + T add ( c ′ | c ) , (12)with the probability of removing the state at location k in the list of states T remove ( c ′ = { n − , k }| c ) = | ρ (1 , , ..., k − , k + 1 , ..., n − | D ( c ′ | c ) , (13)and with the probability of adding α j at location k : T add ( c ′ = { n + 1 , k, α j }| c ) = | ρ (1 , , ..., k − , α j , k, k + 1 , ..., n + 1) | D ( c ′ | c ) , (14)with D ( c ′ | c ) the normalization for the probability: D ( c ′ | c ) = X k | ρ (1 , , ..., k − , k + 1 , ..., n − | + X k X α j | ρ (1 , , ..., k − , α j , k, k + 1 , ...n + 1) | . (15)7he acceptance probability is then A ( c ′ | c ) = min " , D ( c ′ | c ) D ( c | c ′ ) . (16) III. RESULTS AND DISCUSSION
We use the SiLK algorithm to calculate the energies of H O, N , and F at selected bondlengths and bond angles for H O. The Cartesian Gaussian DZ basis set is used in allcalculations. Computer memory constraints imposed by the current SiLK implementationdictates the number of determinants that can be used in the calculations. The reason isthat at present the CI coefficients for the ground and excited states must be stored. Futurework will focus on alleviating the memory issues. We therefore use either the full vectorspace of determinants generated by all possible excitations of the Hartree-Fock determinant(Full Configuration Interaction, FCI) or the restricted vector spaces generated by the HFdeterminant 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 ofthe SiLK method to the exact result within the vector space is made to assess, as mentionedin the introduction, the ability of SiLK to provide accurate results and alleviate the signproblem. 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 arecreated. Determinants corresponding to excited states are generated by excitations of allmolecular orbitals except the core orbitals (the frozen core approximation). Symmetry isused to restrict the determinants to those with the same symmetry as the ground Hartree-Fock state. For the calculations presented here, we use C v spatial symmetry for H O and D h spatial symmetry for N and F respectively. We use T = 1 K ( β = 3 × 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 inthe energy was obtained. The values of P chosen for the reporting of data ranged from2 × to 2 × . The FCI calculations for H O is performed using Molpro .The ability of SiLK to address the sign problem is evaluated by following the evolutionof the sign (for clarity averaged over every 20 Monte Carlo steps) during the course of thelearning period. Representative of the results from the different molecules is the average8ign for water at the minimum energy FCI geometry . In this calculation, the maximumnumber of states included in a diagonalization is limited to 50 (the entire vector space hada dimension of 128,829). The coarse-grained sign is shown in Fig. 2. The sign fluctuatessignificantly for roughly the first 1,500 diagonalizations, but after approximately 1,600 di-agonalizations it remains 1.0. In the upper panel of Fig. 2 the number of states involvedin each diagonalization is shown. After 1,600 diagonalizations, the coarse-grained sign re-mains at 1.0 even though the number of states involved in subsequent diagonalizations isapproximately 20. This indicates that kinks are being introduced during the Monte Carloprocess but these are not affecting the sign. The number of kinks averaged over the kinkconfigurations between 2 successive diagonalizations ranged from roughly 5 at the beginningof the learning period to roughly 2.5 at the end of the learning period. An examination ofthe kink configurations after the learning period found that the configurations contain eitherzero and two kinks and therefore the average sign is 1.0. N u m b e r o f S t a t e s A v e r a g e S i gn FIG. 2. The evolution of the sign during the SiLK learning period for H O using the DZ basis setat the FCI minimum energy geometry, P = 2 × . The O-H bond length is 1.84345 Bohr andthe HOH angle is 110.565 degrees. The upper plot shows the number of states involved in eachdiagonalization, which was constrained to be less than or equal to 50. The lower plot shows theaverage sign evolution, averaged over every 20 diagonalizations, during the learning process. Accurate calculations of potential energy surfaces are important in understanding reac-tion energetics and rates. The ability of SiLK and other quantum chemical methods tocalculate potential energy surfaces is assessed for H O, F , and N . The goal of a successfulmethod is to achieve the accuracy required to describe energetic differences encountered in9hemical processes such as bond-breaking/bond-forming reactions and reaction activationenergies. This so-called “chemical accuracy” is approximately 0.1 – 1 kcal/mol ≈ − − − Hartrees/mol .Several versions of truncated CC are used in this work. The CCSD method uses single anddouble excitations . The CCSDT method uses single, double, and triple excitations . TheCCSD(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 MR-CCSD(T)(2,2) and MRCCSD(T)(4,4) methods are multi reference CC (MRCC) methodswith single, double, and non-iterative inclusion of perturbative triples . A (2,2) calculationuses 2 electrons and 2 orbitals (one occupied and one virtual) to generate the model space forthe MRCC calculation and a (4,4) calculation uses 4 electrons and 4 orbitals (two occupiedand two virtual) to generate the model space for the MRCC calculation. We also use secondorder many body perturbation theory (MBPT(2)). The NWChem software package is usedto perform all standard ab initio calculations . A. Water
The H O molecule is used to assess the ability of SiLK to describe the variation ofenergy with bond length and bond angle in two separate calculations, one in which the bondlength is varied and another in which the bond angle is varied. Fig. 3 displays the energyand its absolute error as a function of bond length at fixed bond angle of 110.565 degrees ascalculated by different methods. The SiLK method has an absolute error of 10 − Hartree overthe range of bond lengths studied, which is well below the desired chemical accuracy. At theminimum energy geometry (bond length = 1.8434 Bohr), the exact energy is -76.14455299Hartrees and the energy of the lowest energy SiLK state is -76.14454690 Hartrees, an errorof ≈ × − Hartrees. Therefore the SiLK procedure found an excellent approximationto the exact ground state. SiLK is approximately one order of magnitude more accuratethan the most accurate of the other methods in the comparison. Notably, SiLK is accurateat the longer bond lengths (roughly two orders of magnitude more accurate than any othermethod) 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 theentire range of bond lengths studied. 10 E [ H a r t r ee ] MBPT(2)CCSDCCSD(T)CCSDTMRCCSD(T)(2,2)MRCCSD(T)(4,4)SiLKExact -6 -5 -4 -3 -2 -1 ∆ E [ H a r t r ee ] FIG. 3. Potential energy curve of H O molecule as a function of OH bond length. A com-parison of results obtained using MBPT(2), CCSD, CCSD(T), CCSDT, MRCCSD(T)(2,2), MR-CCSD(T)(4,4), and SiLK formalisms. The bottom plot displays the absolute error in the calculatedenergy for the different methods. P = 2 × is used for bond lengths in the range [1.34-3.64]Bohr and P = 2 × is used for bond lengths in the range [3.74-4.34] Bohr. Then we use these methods to calculate the energy as a function of bond angle for theH O molecule with bond length 1.84345 Bohr. Fig. 4 displays the energies and absoluteerrors for bond angles ranging from 95 to 125 degrees. The SiLK method is approximatelytwo orders of magnitude more accurate than the most accurate of the other methods. Allmethods except MBPT(2) and CCSD achieve chemical accuracy.It is also instructive to consider calculations restricted to just single and double excitationsas sometimes computations based on such restricted vector spaces can yield useful resultsusing significantly fewer computational resources. Therefore, the SiLK algorithm is used tocalculate the energies and absolute errors of the H O molecule as a function of bond lengthand bond angle. Figs. 5 and 6 show their comparison with exact results. The SiLK methodis able to reproduce the exact results to 10 − Hartree, well within chemical accuracy.
B. Nitrogen N has a triple bond, which provides a challenging test for ab initio methods due toits large electronic correlation . Due to memory limitations, the SiLK calculations wererestricted to the Hartree-Fock determinant plus either the SD and SDT vector spaces. Fig. 711 E [ H a r t r ee ] MBPT(2)CCSDCCSD(T)CCSDTMRCCSD(T)(2,2)MRCCSD(T)(4,4)SiLKExact
100 110 120Angle [Degree]10 -6 -5 -4 -3 -2 ∆ E [ H a r t r ee ] FIG. 4. Potential energy curves of H O FCI vector space for the DZ basis . A comparison of resultsobtained by SiLK with results from MBPT(2), CCSD, CCSD(T), CCSDT, MRCCSD(T)(2,2), andMRCCSD(T)(4,4). Bottom plot displays the absolute error of energy. P = 2 × is used for allangles in SiLK QMC. -76.1-76-75.9-75.8 E [ H a r t r ee ] Exact_SDSiLK_SDExact_SDTSiLK_SDT2 3 4Bond Length [Bohr]10 -9 -8 -7 -6 ∆ E [ H a r t r ee ] FIG. 5. Potential energy curves for H O molecule using the DZ basis and the SD and SDT vectorspaces. The exact results obtained from exact diagonalization and the SiLK results are shown. P = 2 × is used for all the bond lengths in the SD and the SDT vector spaces. shows that the SiLK QMC results converge to the exact result over a wide range of bondlengths. At the minimum energy geometry (bond length = 2.168 Bohr), the exact energyis -109.0858095 Hartrees and the energy of the lowest energy SiLK state is -109.0858094Hartrees, an error of ≈ × − Hartrees. Therefore the SiLK procedure found an excellent12 E [ H a r t r ee ] Exact_SDSiLK_SDExact_SDTSiLK_SDT100 110 120Angle [Degree]10 -8 -7 -6 -5 ∆ E [ H a r t r ee ] FIG. 6. Potential energy curves for the H O molecule within the SD and SDT vector spaces andthe DZ basis. Exact results obtained from exact diagonalization and SiLK results are shown. P = 2 × is used for all the angles in the SD vector space. Within the SDT space, P = 2 × is used for angles in the range [95.565,120.565] and P = 2 × is used for the 125.565 degreecalculation in SiLK QMC. -109-108.5 E [ H a r t r ee ] Exact_SDSiLK_SDExact_SDTSiLK_SDT2 3 4 5 6 7Bond Length [Bohr]10 -9 -8 -7 -6 ∆ E [ H a r t r ee ] FIG. 7. Potential energy curves for N molecule as a function of bond length using the SD andSDT vector spaces. The exact results obtained from exact diagonalization and SiLK results areshown. P = 2 × is used for all bond lengths. approximation to the exact ground state. The results demonstrate that the SiLK methodis suitable for multi reference systems such as N where more than a single determinant isstrongly coupled in the ground electronic state.13 . Fluoride -198.9-198.8 E [ H a r t r ee ] Exact_SDSiLK_SDExact_SDTSiLK_SDT2 3 4 5 6 7Bond Length [Bohr]10 -9 -8 -7 ∆ E [ H a r t r ee ] FIG. 8. Potential energy curves for F molecule as a function of bond length and using the Hartree-Fock determinant plus either the SD and SDT vector spaces. The exact results obtained from exactdiagonalization and the SiLK results are shown. P = 2 × is used. Electron correlations are difficult to include in the simulation of F as many determinantscontribute small but important contributions to the total energy . This phenomenon is of-ten 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 tothe SD and SDT vector spaces. Fig. 8 shows that the SiLK QMC results converge to theexact results and demonstrate that SiLK is capable of accurately including dynamic corre-lation. At the minimum energy geometry (bond length = 2.86816 Bohr), the exact energyis -198.9494169 Hartrees and the energy of the lowest energy SiLK state is -198.9494169Hartrees, an error of < × − Hartrees. Therefore the SiLK procedure found an excellentapproximation to the exact ground state.
D. Scaling Analysis
It is beyond the scope of this work to make a thorough analysis of the scaling of the SiLKmethod as the memory requirements of the current implementation of the SiLK methodprohibits the use of a wide range of basis set size. However, it is important to assessthe scaling of the SiLK algorithm with the size of the basis set and vector space. No14runcation 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 analysisis presented for the current work using the relatively limited size of the vector spaces. If theSiLK algorithm increases too quickly with the size of the vector space, the computationalrequirements for the SiLK method will make its use in the present form intractable. Thedependence of the length of the learning period on the size of the vector space (numberof determinants) is presented in Fig. 9. As expected, there is an increase in the size of thelearning period. However, a wider range of vector space sizes is necessary to fully understandand quantify the scaling behavior. The number of Slater determinants10 L ea r n i ng M C s t e p s H2ON2F2
FIG. 9. The length of the learning period as a function of the number of Slater determinants. Theresults for the SD, SDT, and FCI vector spaces for all molecules and geometries are shown.
IV. CONCLUSIONS
The minus sign problem in Quantum Monte Carlo simulations of frustrated or correlatedelectronic systems is a challenging problem. It has even been suggested that a generalsolution of this problem is NP-complete . Therefore, one should not expect an effectivesolution for all the Monte Carlo simulations which have the minus sign problem. In thispaper, we demonstrate that SiLK QMC can reduce the minus sign problem by using alearning stage that includes a diagonalization procedure. In this paper, we demonstratethat the energies obtained by the SiLK QMC match the results from exact diagonalization15nd surpass the accuracy obtained using other quantum chemistry methods, particularly forgeometries relatively far from equilibrium. In addition, SiLK can be applied to systems thatrequire a multiple reference state approach. An intriguing possibility for future work is touse the SiLK learning procedure in combination with other QMC algorithms to reduce theminus sign problem.As the learning stage progresses, the states become more complicated linear combinationsof determinants so that more evaluations of matrix elements are required, thereby increasingthe computational expense. However, at the same time the number of non-zero matrixelements between these states decreases. So, further optimization is possible by storingoften-needed matrix elements in memory. For example, storing the off-diagonal matrixelements between the ground and excited states yields a large speed up, since these are theonly matrix elements required once the learning stage reaches the point where mostly zeroand two kink configurations appear in the simulation. It is also possible to halt the learningstage at an earlier point, when the ground SiLK state is not as accurate an approximationto the exact ground state, but when the sign problem is alleviated but not eliminated andrely on the Monte Carlo sampling to provide the exact energy. This would reduce thecomputational effort required to evaluate the matrix elements since fewer diagonalizationswill have occurred. An investigation of the efficacy of a shorter learning period is left forfuture work.The SiLK method requires the knowledge of the off-diagonal matrix elements of theHamiltonian. As the size of the system increases, the number of off-diagonal matrix elementsincreases factorially and it is not possible to store the matrix elements or the CI coefficientsfor the ground and excited states that would allow for on-the-fly evaluation of the matrixelements. As such, without a procedure to accurately truncate the number of determinantsused to describe the ground and excited state wavefunctions, the use of the all possibledeterminants 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 anddouble excitations. Such truncated sets of determinants are often sufficient for the studyof chemical systems. In cases where restrictions to single and double excitations are notsufficient, 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 molecular16ystems. Since the path integral formulation uses the canonical partition function it ispossible 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. ACKNOWLEDGMENTS
REFERENCES A. Dreuw and M. Head-Gordon, Chemical Reviews , 4009 (2005). R. P. Feynman, A. R. Hibbs, and D. F. Styer,
Quantum mechanics and path integrals:Emended edition (Dover Publications, 2005). W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Reviews of Modern Physics , 33(2001). M. Troyer and U.-J. Wiese, Physical Review Letters , 170201 (2005). E. Loh Jr, J. Gubernatis, R. Scalettar, S. White, D. Scalapino, and R. Sugar, PhysicalReview B , 9301 (1990). R. O. Jones and O. Gunnarsson, Reviews of Modern Physics , 689 (1989). R. Harrison and N. Handy, Chemical Physics Letters , 386 (1983). C. Møller and M. S. Plesset, Physical Review , 618 (1934). N. Handy, P. Knowles, and K. Somasundram, Theoretica Chimica Acta , 87 (1985). P. Pendergast and E. F. Hayes, Journal of Computational Physics , 236 (1978). R. J. Bartlett and M. Musiał, Reviews of Modern Physics , 291 (2007). F. Coester, Nuclear Physics , 421 (1958).17 F. Coester and H. Kümmel, Nuclear Physics , 477 (1960). J. Čížek, The Journal of Chemical Physics , 4256 (1966). J. Čížek and J. Paldus, International Journal of Quantum Chemistry , 359 (1971). J. Paldus, J. Čížek, and I. Shavitt, Physical Review A , 50 (1972). G. D. Purvis III and R. J. Bartlett, Journal of Chemical Physics , 1910 (1982). K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chemical PhysicsLetters , 479 (1989). S. R. White, Physical Review B , 10345 (1993). G. K.-L. Chan and S. Sharma, Annual Review of Physical Chemistry , 465 (2011). K. H. Marti and M. Reiher, Zeitschrift für Physikalische Chemie International Journal ofResearch in Physical Chemistry and Chemical Physics , 583 (2010). S. Wouters and D. Van Neck, The European Physical Journal D , 1 (2014). N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Journalof Chemical Physics , 1087 (1953). N. Metropolis and S. Ulam, Journal of the American Statistical Association , 335 (1949). E. Fermi and R. Richtmyer, “Note on census-taking in Monte-Carlo calculations,” Tech.Rep. (Los Alamos Scientific Lab., 1948). N. Rom, E. Fattal, A. K. Gupta, E. A. Carter, and D. Neuhauser, Journal of ChemicalPhysics , 8241 (1998). R. Baer, M. Head-Gordon, and D. Neuhauser, Journal of Chemical Physics , 6219(1998). D. M. Ceperley and B. Alder, Physical Review Letters , 566 (1980). A. J. Thom and A. Alavi, The Journal of chemical physics , 204106 (2005). A. J. Thom, G. H. Booth, and A. Alavi, Physical Chemistry Chemical Physics , 652(2008). M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer, and E. J. Walter, Physical ReviewB , 245123 (2007). E. W. Brown, B. K. Clark, J. L. DuBois, and D. M. Ceperley, Physical review letters , 146405 (2013). E. W. Brown, J. L. DuBois, M. Holzmann, and D. M. Ceperley, Physical Review B ,081102 (2013). G. H. Booth, A. J. Thom, and A. Alavi, Journal of Chemical Physics , 054106 (2009).18 G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Nature , 365 (2013). R. W. Hall, Journal of Chemical Physics , 1 (2002). R. W. Hall, Chemical Physics Letters , 549 (2002). R. W. Hall, Journal of Chemical Physics , 164112 (2005). P. W. Anderson, G. Yuval, and D. Hamann, Physical Review B , 4464 (1970). R. Chiles, G. Jongeward, M. Bolton, and P. Wolynes, Journal of Chemical Physics ,2039 (1984). D. Feller, Journal of Computational Chemistry , 1571 (1996). K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li,and T. L. Windus, Journal of Chemical Information and Modeling , 1045 (2007). T. H. Dunning Jr, Journal of Chemical Physics , 2823 (1970). T. H. Dunning, P. J. Hay, and H. Schaefer, in
Modern Theoretical Chemistry , Vol. 3(Plenum Press New York, 1977) p. 1. M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. Van Dam,D. Wang, J. Nieplocha, E. Apra, T. L. Windus, et al. , Computer Physics Communications , 1477 (2010). H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, Wiley Interdisci-plinary Reviews: Computational Molecular Science , 242 (2012). molpro is a package of ab initio programs written by H.-J. Werner and P. J. Knowles, withcontributions from J. Almlöf, R. D. Amos, A. Berning, D. L. Cooper, M. J. O. Deegan,A. J. Dobbyn, F. Eckert, S. T. Elbert, C. Hampel, R. Lindh, A. W. Lloyd, W. Meyer,A. Nicklass, K. Peterson, R. Pitzer, A. J. Stone, P. R. Taylor, M. E. Mura, P. Pulay, M.Schütz, H. Stoll, and T. Thorsteinsson. P. Saxe, H. F. Shaefer, and N. C. Handy, Chemical Physics Letters , 202 (1981). P. A. Bash, L. L. Ho, A. D. MacKerell, D. Levine, and P. Hallstrom, Proceedings of theNational Academy of Sciences , 3698 (1996). J. Noga and R. J. Bartlett, Journal of Chemical Physics , 7041 (1987). K. Bhaskaran-Nair, J. Brabec, E. Aprà, H. J. van Dam, J. Pittner, and K. Kowalski,Journal of Chemical Physics , 094112 (2012). K. Kowalski and P. Piecuch, Journal of Chemical Physics , 5644 (2000). K. Kowalski and P. Piecuch, Chemical Physics Letters , 165 (2001). D. K. Mok, R. Neumann, and N. C. Handy, The Journal of Physical Chemistry , 6225191996). Y. Lu, M. Höppner, O. Gunnarsson, and M. W. Haverkort,Phys. Rev. B90