Beyond Single-Reference Fixed-Node Approximation in Ab Initio Diffusion Monte Carlo Using Antisymmetrized Geminal Power Applied to Systems with Hundreds of Electrons

Diffusion Monte Carlo (DMC) is an exact technique to project out the ground state (GS) of a Hamiltonian. Since the GS is always bosonic, in Fermionic systems, the projection needs to be carried out while imposing antisymmetric constraints, which is a nondeterministic polynomial hard problem. In practice, therefore, the application of DMC on electronic structure problems is made by employing the fixed-node (FN) approximation, consisting of performing DMC with the constraint of having a fixed, predefined nodal surface. How do we get the nodal surface? The typical approach, applied in systems having up to hundreds or even thousands of electrons, is to obtain the nodal surface from a preliminary mean-field approach (typically, a density functional theory calculation) used to obtain a single Slater determinant. This is known as single reference. In this paper, we propose a new approach, applicable to systems as large as the C60 fullerene, which improves the nodes by going beyond the single reference. In practice, we employ an implicitly multireference ansatz (antisymmetrized geminal power wave function constraint with molecular orbitals), initialized on the preliminary mean-field approach, which is relaxed by optimizing a few parameters of the wave function determining the nodal surface by minimizing the FN-DMC energy. We highlight the improvements of the proposed approach over the standard single-reference method on several examples and, where feasible, the computational gain over the standard multireference ansatz, which makes the methods applicable to large systems. We also show that physical properties relying on relative energies, such as binding energies, are affordable and reliable within the proposed scheme.


Introduction
Ab initio electronic structure calculations, which compute the electronic structure of materials non-empirically, have become an essential methodology in the materials science and condensed matter physics communities.Density functional theory (DFT), a mean field approach which was originally proposed by Kohn and Hohenberg1 , is the most widely used methodology for ab initio electronic structure calculations.DFT has enjoyed widespread success, despite its reliance on the so-called exchange-correlation (XC) functional, whose exact form is yet to be discovered.Although many XCs have been proposed, no functional that performs universally well for all materials is established.largest source of error: the fixed-node (FN) approximation.DMC yields the exact ground state in bosonic systems.In fermionic systems (for instance, in electronic structure calculations) DMC suffers from the so-called negative sign problem, arising from the fact that the fermionic ground state has positive and negative regions.The negative sign problem in the DMC method for fermions has been proven to be a nondeterministic polynomial hard problem 14 ; thus, it seems unrealistic to find a general solution at present.This problem is avoided, in practice, by modifying the DMC projection algorithm with the introduction of the fixed-node (FN) approximation, where the projected wave function Φ FN is constrained to have the nodes of a predetermined guiding function Ψ T .
The FN approximation keeps the projected wave function Φ FN antisymmetric, but Φ FN is the exact ground state Φ 0 only if its nodes are exact.A general property of Φ FN is that it is always the closest function to Φ 0 within the FN constraint.For trial functions obtained from mean field approaches, such as Hartree-Fock(HF) or DFT, it is generally believed that the error associated to the FN approximation is small and benefits from a large error cancellation in the evaluation of binding energies 8 .However, the FN error is typically not accessible, as Φ 0 is unknown, and this yields an uncontrollable error in FN-DMC.
In standard FN-DMC simulations, the nodal surface is given by an approximate wave function, which is typically obtained starting from a mean-field approach, such as HF or DFT.The variational principle can still apply to FN-DMC, 1 and so to go beyond the meanfield solution, one should optimize the given nodal surface by minimizing the the FN-DMC energy E FN (which is the expectation value of Φ FN ), going in the direction of the exact wavefunction Φ 0 and the exact energy E 0 .This procedure is seldomly followed in DMC simulations, especially on large systems (say, with hundreds or thousands of electrons), as it is hardly affordable computationally and the uncertainty on the optimization of the FN surface could be easily comparable, if not larger, than the binding energy under consideration.Thus, the standard approach is to just keep the nodal surface of the Slater determinant built with the Kohn-Sham orbitals obtained from a DFT calculation.Whilst the FN surface from DFT might be suboptimal, this approach typically yields quite reliable results, especially in the evaluations of non-covalent interactions, due to very favourable error cancellations. 5,80][41][42][43][44][45] The standard approach here is to optimize the wave function parameters at the level of theory of variational Monte Carlo (VMC) 19,[46][47][48][49][50] .Indeed, optimization at the FN-DMC (FN-opt) level implies further difficulties, as we will discuss below.However, optimization at the VMC (VMC-opt) level has some flaws.In VMC-opt the object that is optimized is the variational wave function Ψ T , which is obtained from the product of one of the ansatze discussed above and the Jastrow factor 2 .The closer Ψ T gets to the ground state Φ 0 , the smaller its VMC energy (variational principle) and its VMC variance (zero-variance property) are.VMC-opt explores the parameters variational space, seeking the set which minimises the VMC energy or the VMC variance, and it is often done by employing the VMC gradient.It is not guaranteed that the parameters obtained from VMC-opt are those giving the best possible nodal surface allowed by the employed ansatz (unless we are in the limit case where Ψ T yields VMC with zero variance, such that we know that Ψ T is an eigenstate of the Hamiltonian).Although this approach, in practice, gives a better nodal surface than the DFT one, it sometimes gives unreasonable outcomes, e.g., it overestimates binding energies, as revealed in this work.It would be desirable, instead, to implement an optimization at the FN-DMC level of theory, where the parameters of the function Ψ T giving the nodal surface are optimized so as to minimize the FN energy.This would guarantee to find the best nodal surface allowed by the adopted wave function ansatz.
To the best of our knowledge, the first attempt to directly optimize the variational parameters included in a trial WF at the FN-DMC level was done by Reboredo et al. 51 in the ab initio framework.They proposed a way to iteratively generate new trial wavefunctions to get a better nodal surface.They generalized the method to excited states 52 and finite temperatures 53 and also applied for large systems such as C 20 54 .Very recently, McFarland and Monousakis 55 reported successful energy minimizations with approximated and exact FN gradients.They proposed to optimize nodes using a combination of FN gradients and the projected gradient descent (PGD) method.The PGD method works for Be, Li 2 , and Ne using all-electron DMC calculations 55 , while it has been successful only for small molecules.
When it comes to optimizing the nodal surface of a large system, the main problem is that the number of variational parameters determining the nodal surface often scales more than linearly with the size of the system.For instance, the number of variational parameters in the determinant part of the AGP ansatz scales with O(L 2 ), where L is the number of basis functions in a system.It makes the parameter space to be optimized so complex that the optimization is easily trapped in local minima and one cannot find the true ground state.Moreover, since the optimization algorithms are stochastic, there is always an additional uncertainty on the optimized parameters, which are not going to be exact and the corresponding QMC energy has therefore an optimization bias.The optimization bias increases with the system size and with the number of variational parameters, and can be reduced only at the cost of increasing the statistical sampling (and the computational cost).
The evaluation of a binding energy implies the difference between two or more DMC energies, and it is often a tiny fraction of the total energy.Therefore the optimization uncertainty can often be comparable to the binding energy, making the evaluation of the interaction energy unreliable.Moreover, we need to verify that the adopted approach satisfies basic physical properties, such as being size-consistent 3 .At the VMC level, the size-consistency is a property of the wave function ansatz employed, and it depends on the optimization procedure.At the FN-DMC level, size consistency might depend on some choices on the algorithm, 56 on the ansatz of the wave function providing the FN constraint, and on the optimization.
In this paper, we propose a scheme which aims to address these issues.In particular, our scheme satisfies the following points: (i) it is systematically more accurate than the standard approach of employing a single Slater determinant, (ii) it is size consistent, and (iii) it is applicable also to large systems.The idea underlying the present work is the combination of the AGP wavefunction consisting of molecular orbitals (AGPn) 21 , the use of natural orbitals, and the optimization of its nodal surface using FN gradients on a selected subset of the AGPn parameters.In particular, we initialize the orbitals in the AGPn wave function using natural orbitals, which are kept fixed afterwords, such that only the coefficients combining them are optimized to relax the nodal surface.We call this scheme the fixed node antisymmetrized geminal power active-space (FNAGPAS).Since the orbitals are fixed, this results in a much smaller number of variational parameters in the ansatz; thus, one can apply it for larger systems, such as C 60 fullerene.We show that our scheme gives a better nodal surface (i.e., a lower energy in the FN-DMC calculation) compared to the typical Slater-Jastrow ansatz, and it reliably describes also strongly correlated systems (such as diradicals).We show that the use of FN-opt is important to fulfill the size-consistency property.

The FNAGPAS scheme
We describe here the scheme that we suggest to improve the accuracy of FN-DMC over the traditional single determinant Slater-Jastrow ansatz.The key idea is the combination of the AGP wavefunction constraint with molecular orbitals (AGPn) 21 and the optimization of the 3 An approach is size consistent if the energy of a system constituted by two or more non-interacting subsystems (e.g., two molecules far away) is the same of the sum of the energies of the subsystems.ansatz using approximated FN gradients. 55We describe the ansatz in the following section, assuming an unpolarized system for simplicity.The schematic illustration explaining the key concept and its workflow is shown in figure 1.

Computation of Gradients
Computing gradients of the MO or NO weights using FN-DMC.

Update of nodal surface
Update the weights of the MOs ( ) according to the gradients λ μ → λ′ μ We perform a preliminary mean field calculation to obtain molecular orbitals (MOs), followed by a correlated calculation yielding natural orbitals (NOs).The AGPn ansatz correspond to a multideterminant expansion built on the NOs and depending on the coefficients λ i associated to each orbital i and optimized in order to minimise the FN energy.Panel b: Flowchart illustrating the FNAGPAS scheme workflow.

Final FN-DMC evaluation
The real-space quantum Monte Carlo typically employs a many-body WF ansatz Ψ written as the product of two terms, Ψ QMC = Φ AS × exp J.The term exp J, conventionally dubbed Jastrow factor, is symmetric, and the term Φ AS is antisymmetric.The Jastrow factor is explicitly dependent on electron-electron distance, and often includes electron-nucleus and electron-electron-nucleus terms. 4The nodal surface of a WF is determined by the antisymmetric part Φ AS (because exp J ≥ 0).Thus, in FN-DMC the accuracy of the results depends crucially on the quality of the nodes of Φ AS .
The antisymmetric part of a trial WF is initially constructed from a mean-field selfconsistent-field (SCF) approach, such as DFT or HF.The standard QMC setup in large systems is to define Φ AS as the single Slater determinant (SD) obtained from such preliminary SCF calculation.The corresponding Ψ QMC is dubbed JSD.Therefore, the nodes of JSD are predefined before any QMC calculation and unrelaxed.Initializing the SD using different setups for the SCF calculations (e.g., different exchange correlation functionals) leads to slightly different total energies, but most of the times, the interaction energies (which are evaluated from energy differences between two or more systems) are almost unaffected by the details of the preliminary SCF calculation, expecially for weak non-covalent interactions.
This is an indication that there is an almost perfect cancellation of the error induced by the FN approximation within the JSD ansatz, provided that the SD is initialized consistently in all systems.
However, changing the setup of the SCF calculation only allows the nodes to move within the variational freedom of a single Slater determinant.By contrast, giving Φ AS the variational freedom to relax the nodes beyond the JSD ansatz leads to an improvement of the FN-DMC total energy of the system, 58 and possibly also the interaction energies could change.The challenge that we take here is to generalize the ansatz in a way that large systems are still doable.
Here we suggest to use the AGP ansatz as Φ AS .AGP is an implicitly multideterminant ansatz, 56,59 which corresponds to a constrained zero-seniority expansion, as illustrated schematically in figure 1.The evaluation of an AGP function can be reduced to the computation of a determinant, therefore the AGP ansatz is computationally comparable to a Slater determinant (differently from explicitly multideterminant functions), thus ensuring the cubic scaling with the system size of both the variational and FN algorithms 5 .The AGP ansatz for a system of N el electrons is (we are assuming for simplicity an unpolarized system with even numbers of electrons, but the ansatz can be generalized as discussed in Ref. 19), where Â is the antisymmetrization operator and the function g is the geminal function , which is a paring function between two electrons with coordinates x 1 and x 2 forming a spin singlet.
The spatial part f (r 1 , r 2 ) is symmetric, and it can be written in terms of a basis set {χ µ } for the single electron orbital space as follows: where µ and ν runs over all the L basis orbitals, and c µν are variational parameters.Notice that in general L ≫ N , and the number of variational parameters c µν is equal to L 2 .The parameters define a L × L symmetric matrix C (the symmetry of f implies c µν = c νµ ), so there is an orthogonal transformation U which diagonalizes C and allows to rewrite f as: where ϕ µ = ν U µν χ ν .With no loss of generality we can assume that λ's are ranked in decreasing order of their absolute value.Notice that if only the first N el /2 λ's are different from zero then Ψ AGP corresponds to a single Slater determinant built on the orbitals ϕ 1 , . . ., ϕ N el /2 occupied with both spin up and spin down electrons.Since such Slater determinant built on orbitals from an SCF calculation is the standard QMC setup, and it typically delivers good results, we tried to relax the nodes by considering a subset n orb (larger than N el /2 but ≪ L) of the orbitals obtained from the SCF calculation.This is what we call the AGPn ansatz.
For an efficient and effective use in QMC the AGP and AGPn functions shall be multiplied by a Jastrow factor, yielding the so-called JAGP and JAGPn functions.The Jastrow factor can have the same variational form used also in JSD, which allows for the JSD, JAGP and JAGPn functions to satisfy the cusp conditions and to effectively recover the dynamical correlations.Indeed, the main improvement of JAGP and JAGPn over JSD is their ability to capture static correlations, yielding to qualitatively different results on systems with an underlying multireference character, both at the variational and at the FN level of theory. 56,59e optimization of the parameters in the Jastrow is usually quite a feasible problem also on large systems, as their number does not grow uncontrollably with the size of the system.
In practice, every QMC code implements a slightly different functional form of the Jastrow, but their share the general features mentioned above.The QMC code used in this work is TurboRVB, 57 an open-source package.In TurboRVB the implemented the Jastrow factor (described in Ref. 57) has a number of parameters growing linearly with the size of the system (as shown in the results section).
In this work, we keep the orbital frozen and optimize the coefficients λ 1 , . . ., λ n orb of the JAGPn ansatz using FN-DMC gradients.A similar idea, but at the variational level, was also mentioned in a seminal work by Casula and Sorella to decribe the BCS paring function in iron-based superconductors 60 .JAGPn dramatically reduces the number of variational parameters with respect to the JAGP ansatz, such that the optimization of the JAGPn function is doable even in pretty large systems, in contrast to JAGP which is affordable only on relatively small systems.Nevertheless, employing JAGPn significantly improves the FN-DMC energy (as well as the variational QMC energy) over the results within the traditional JSD function, as we will show on the results section.Of course, the JAGP ansatz has higher variational freedom than JAGPn, so JAGP can in principle improve further over JAGPn.
However, in practice, we observe that FN-DMC energies obtained from the JAGP ansatz are comparable to those obtained from JAGPn on small systems (and both JAGP and JAGPn are significantly better that JSD), while, in large systems, JAGP is unaffordable because the optimization can be stuck at local minima at the variational level and can become unstable at the FN level.The latter instability is probably due to insufficient signal-to-noise ratios 61 that the QMC optimization always suffers from, but the origin of the instability is yet unclear.
On intermediate systems, we notice that JAGP FN-DMC energy is worse than the JAGPn FN-DMC energy, as a clear indication that despite the higher variational freedom on JAGP, the optimization of that many parameters is not converging and there is too much noise on the parameters.
The main problem of the AGP ansatz (and AGPn) is that it is not size-consistent at the variational level of theory, but JAGP (JAGPn) is size consistent if we employ a very flexible Jastrow factor. 62,63Since the FN-DMC corresponds to applying an infinitely flexible Jastrow factor to the determinant part, optimizing the AGPn parameters at the FN level ensures the size-consistency of our approach.
A crucial point to make JAGPn almost as accurate as JAGP, despite employing only a small number n orb of parameters λ's, is to carefully choose the orbitals.We notice that the virtual orbitals obtained from SCF calculations are typically not optimal, as we need a large number of them (of the order of L) to converge to the best JAGPn FN energy.
Moreover, if we cannot afford a systematic test of the convergence of n orb for each system of interest, it is difficult to define a sensible criterion to decide which n orb to pick.We solved both the problems by employing Natural Orbitals (NOs) for expanding the paring function, instead of using MOs.NOs were constructed from second-order Møller-Plesset (MP2) calculations.This is because the MP2 unoccupied orbitals incorporate perturbation effects and are physically better than those obtained with HF or DFT, 64 as shown in the Supplemental Information.More specifically, we constructed natural orbitals by diagonalizing the density matrix obtained by MP2 calculations.We also notice that a method to construct NOs should be affordable also for large systems.This is also a reason why we chose MP2 for constructing NOs in this study.In practice, from the weight of the NOs we can easily define a cutoff value to select n orb on each system, and we notice that we get to converged results already with a value of n that is not much larger than N el /2 (n orb = N el /2 would correspond to a single SD).

Computational details
We applied our scheme to planar and twisted ethylenes, eight hydrocarbons (CH files.We employed the cc-pVQZ basis set accompanied with the ccECP pseudopotentials 69 for the eight hydrocarbons and C 60 fullerene, while the cc-pVTZ basis set accompanied with the ccECP pseudopotentials 69 for the water-methane and for the torsion calculation of ethylene.We employed [3s], [3s1p], and [3s1p] primitive Jastrow basis for H, C, and O atoms, respectively.The Jastrow factor and the weights of the natural orbitals in the paring function (i.e., the nodal surface of a WF) were optimized using the stochastic reconfiguration method 70 implemented in TurboRVB 57 with an adaptive hyperparameter 71 .The Jastrow factor was optimized only with VMC gradients, and it was held fixed during optimization with FN gradients.The FN gradients were computed from a standard walker distribution using mixed estimators, which corresponds to Method A in Ref. 55.The lattice discretized version of the FN-DMC calculations (LRDMC) 72,73 was used in this study.The singleshot LRDMC calculations were performed by the single-grid scheme 72 with lattice spaces a = 0.30, 0.25, 0.20, and 0.10 Bohr, and the energies were extrapolated to a → 0 using The LRDMC calculations for computing those gradients were performed by the single-grid scheme 72 with lattice spaces a = 0.20 Bohr.The Determinant Locality approximation (DLA) 18 was employed for the LRDMC calculations 6 .We notice that the LRDMC framework guarantees the variational principle even with the presence of non-local pseudopotentials, as proven in the Appendix.The molecular structures are depicted using VESTA 74 .
4 Results and Discussion

The FNAGPAS captures strong correlation
We show that the proposed FNAGPAS is able to incorporate the correlation effect that the JSD ansatz cannot do at all.We apply our scheme for the torsion energy estimation of ethylene (C 2 H 4 ).The torsion energy is defined as the energy difference between the ground state ethylene structure (denoted as planar ethylene) and the orthogonally rotated ethylene structure (denoted as twisted ethylene), which are both shown in the inset of Fig. 2. Here, we consider only the singlet states for both configurations.It was shown 59 that the JSD ansatz cannot describe the torsion energy correctly since the ansatz cannot consider the static electronic correlation of the twisted ethylene, which has a diradical character.This is true both at the variational and at the FN level of theory 59 .The lack of reliability in the FN results based on a JSD ansatz indicates that projection schemes cannot recover strong correlation if the FN constraints are given from a wave function with qualitative issues, due to the constraint on the projection coming from the trial wave function.Thus, the way to improve the quality of the FN results is to adopt a more general ansatz, able to improve the nodes of the trial wave function and enhance the reliability of FN estimations.
The planar ethylene has an electronic structure characterized by a highest occupied molecular orbital (HOMO) of type π and a lowest unoccupied molecular orbital (LUMO) of type π * , and the HOMO-LUMO gap is finite.A single Slater determinant having two electrons of unlike spin on the HOMO and no electrons on the LUMO captures qualitatively well the nature of the wave function and there is no static correlation.However, when the molecule is twisted, the HOMO-LUMO gap decreases because the overlap between the p orbitals (or-thogonal to the plane of the −CH 2 atoms) of the two carbons decrease.At a torsional angle on 90 degrees (i.e., twisted ethylene) the two p orbitals become orthogonal and the frontier orbitals become degenerate, forming two singly occupied molecular orbitals (SUMOs).We can define three independent (orthogonal) wave functions having two electrons on two degenerate orbitals forming a spin singlet, a diradical and two zwitterionic states. 75Their wave functions imply the use of more than one Slater determinant, i.e., their electronic structure shows strong correlation.Thus, a multireference ansatz is needed to correctly describe the diradical character of the orthogonally twisted ethylene 59 .
Figure 2 shows the torsion energies of ethylene computed with the JSD ansatz with a HF nodal surface, and the same energies computed with the JAGPn ansatz with HF molecular orbitals 7 , whose weights are optimized using DMC gradients.As a comparison, we also show results obtained with the full JAGP ansatz optimized using VMC gradients, which was taken from Ref. 59.The reference value in Fig. 2 is taken from Ref. 76, and it is computed using MR-CISD+Q 8 .The JSD ansatz gives 133.1(4) kcal/mol for the torsion energy, which is far from the reference value obtained by MR-CISD+Q (i.e., 69.2 kcal/mol 76 ).Our JAGPn ansatz gives a FN energy of 73.0(4) kcal/mol for the torsion energy, which is close to the reference values.This result demonstrates that the JAGPn ansatz optimized using FN gradients correctly describes the diradical character of the orthogonally twisted ethylene, something that the JSD ansatz cannot do.

Application of FNAGPAS to small and large systems
We now show that the FNAGPAS scheme leads to a systematic improvement over the traditional JSD ansatz in molecular systems of increasing size, showing an accuracy in line with the full JAGP ansatz (and better on systems where the optimization error for the JAGP ansatz is large), while being affordable on much larger systems.We consider the eight hydrocarbons and the C 60 fullerene, represented in Figure 3.
Figure 4 (top panel) shows the energy gain in the LRDMC total energies (a → 0) by the nodal-surface optimizations of JAGP and JAGPn over the traditional JSD ansatz (with the nodal surface taken from the DFT LDA calculations).Our proposed FNAGPAS scheme (JAGPn ansatz optimized using FN gradients) shows positive gains for all molecules, indicating that the nodal surface optimizations improve the nodes of the Slater determinant obtained from DFT.Therefore, there is a systematic improvement in the description of the correlation energy.The energy gain scales linearly with the number of electrons in the system.The traditional JAGP ansatz (optimized using VMC gradients) was computationally affordable only on the four smallest systems, due to the rapid increase of the number of variational parameters (see the bottom panel in Figure 4), which makes the optimization unstable or not converging.In addition, we could only use VMC-opt for the JAGP ansatz, because FN-opt is not stable.This highlights an additional crucial advantage of FNAGPAS Figure 3: Molecular systems considered in this work, whose FN energy has been computed with the traditional JSD ansatz and with the JAGPm ansatz (within the FNAGPAS scheme) discussed in this work.The energy gain (i.e., the improvement of the FNAGPAS scheme over the traditional scheme which employs the JSD ansatz) and the number of variational parameters in the wave function for each system are shown in figure 4.
over the traditional JAGP approach.In the four systems where we have both the traditional JAGP and the FNAGPAS results, the latter is equivalent to the former on ethane, and it recovers more correlation energy in methane, ethylene and benzene.Larger systems were computationally unaffordable with JAGP, while JAGPn optimization remains feasible both at the variational and at the FN level.In fact, FNAGPAS has been successfully performed up to C 60 fullerene.The gain in C 60 is ∼ 2 meV/valence electron, as shown in the inset of Fig. 4.This is a reasonable value, considering a previous study by Marchi et al. reporting ∼ 3meV/valence electron for the finite-size graphene calculations with the same atoms as the C 60 77 .
Let us consider more closely the medium-size molecules.Figure 4 shows that the gains of JAGPn (optimized with FN gradients) are larger than those of JAGP (optimized using VMC gradients) in spite of the compactness of the AGPn ansatz.In fact, the number of variational parameters in the benzene molecule is 86 for the JAGPn ansatz, and is 17,629 for the JAGP ansatz.Moreover, JAGP is a generalization of JAGPn.Therefore, one could naively expect that the larger the number of variational parameters, the lower the energy.Here, we observe an exception to this expectation.For this point, we recall that the calculations reported in figure 4 are obtained with a quite small Jastrow factor, employing a [3s1p] basis set for C atoms and a [3s] for H atoms.This is because we target large systems with FNAGPAS, for which the use of large Jastrow factors is unaffordable.It has been reported that an incomplete Jastrow factor leads to misdirection of the nodal surface within the variational optimization of the JAGP ansatz in the square H 4 78 .To confirm if this is the case in the present calculations, we performed additional calculations with a larger Jastrow factor in the JAGP ansatz calculations (i.e., a basis set of [4s3p1d] and [3s1p] for C and H atoms, respectively) and obtained that the larger Jastrow factor leads to a much larger energy gain than that obtained with the JAGP ansatz with a small Jastrow (see results in the SI (Table S-I and Fig. S-I).The result indicates that the small Jastrow factor leads to misdirection of the nodal surface of the JAGP ansatz also in this study.On the other hand, figure 4 demonstrates that the FNAGPAS scheme works even with a small Jastrow factor and a minimal number of parameters in the antisymmetric part, making the approach applicable to larger systems.
As mentioned in the method part, see Section 2, the two main features over which FNAGPAS is built are: 1) the AGPn ansatz, and 2) the optimization of its nodal surface using FN gradients.To reveal which of the two is more crucial for the success of the method, i.e., the ansatz or the gradient, we tried the following combinations: (i) JAGPn with VMCopt; (ii) JAGPn with FN-opt, (iii) JAGP with VMC-opt; (iv) JAGP with FN-opt.Note that (ii) corresponds to FNAGPAS.The scheme (iv), unfortunately, is not possible as the JAGP has too many parameters and the FN optimization becomes unstable.Results obtained with schemes (i) to (iii) are reported in the SI (Table S-I and Fig. S-I).We observe that scheme (ii) gives the best gains.Scheme (i) gives gains close to (ii), and they both are much better than (iii).Thus, it emerges that freezing the orbitals to those obtained by a mean-field approach plays a crucial role in avoiding misdirection of the nodes optimization.

The FNAGPAS scheme is size-consistent
We have shown that the AGPn ansatz is able to gain correlation energies at the FN level using very few variational parameters.In addition to their role in improving the nodal surface, FN gradients also appear to be crucial when calculating binding energies of molecules, preserving size consistency.As shown in Table 1 and discussed hereafter for the particular case of the water-methane dimer, this is not the case when VMC gradients are used.Therefore, when calculating binding energies of molecules, the use of VMC gradients in the JAGPn ansatz gives incorrect results, while the use of FN gradients plays a crucial role on it.program 9 .We chose the CCSD(T) value as a reference because the bounded water-methane dimer is not a strongly-correlated system, thus CCSD(T) should describe the binding energy correctly.In this system the JSD ansatz gives a binding energy of −27(2) meV, which is in good agreement with the CCSD(T) values of -27.0 meV.Thus, a new DMC approach with nodal surface optimization should lower the value of the total energies but should not affect the energy differences.The FNAGPAS scheme, which optimizes the JAGPn parameters with the FN gradients, behaves as expected, yielding a binding energy of -29(2) meV, still in good agreement with the reference value.However, this is not the case for the JAGPn We can interpret the deterioration of the binding energy as follows: Binding energies are computed from relative energies among two or more molecules; thus, the accuracy relies on its error cancellation.The error cancellation in DMC was reviewed and discussed by Dubecký in 2016 8 .Their conclusion is that one can rely on error cancellation as long as one keeps the constructions and optimizations of the corresponding wave functions as systematic as possible.Indeed, this cancellation works when the nodes are kept at the same systematic accuracy at every step of the trial wave function constructions.In fact, for the water-methane dimer calculations in this study, our JSD ansatz fully satisfies the size consistency and gives satisfactory binding energy, which means that the error cancellation works with the DFT nodal surfaces.In this study, we found that error cancellation was deteriorated by the nodal surface optimizations using the VMC gradients while recovered by those using the FN-DMC gradients.When one computes the binding energy of a complex system, one usually uses the same Jastrow basis sets for each element in the complex and the isolated systems.The use of the same Jastrow basis sets does not guarantee the same contribution to the total energy of both the complex and the isolated systems at the VMC level.Indeed, during the nodal surface optimization at the VMC level, the incomplete Jastrow factor affects the nodal surface differently between the complex and isolated systems; thus, the resultant nodal sets, cc-pCVnZ, with quadruple-zeta (n = Q) and quintuple-zeta (n = 5) basis set.We performed both estimations with and without counterpoise correction, both yielding a binding energy of -27.2 meV.
surface gives the incorrect binding energy.The recovery should be because FN-DMC is a projection method to relax the amplitude of the AGPn ansatz, which corresponds to adding an unlimited flexible Jastrow factor to a given ansatz.
The Jastrow incompleteness is also related to the deterioration of the size consistency for JAGPn and JAGP with VMC optimization.The size consistency is a property that guarantees the consistency of the energy behavior when the interaction between the involved molecular system is nullified (e.g., by a long distance).If the size consistency is fulfilled, the energy of the far-away system should be equal to the sum of the energies of the two isolated molecules.The last column in Table 1 shows the difference in energies of the faraway watermethane complex (at a distance of ∼ 11 Å) and the sum of the isolated molecules, which can be considered the size-consistency error and is here dubbed E SCE .JSD ansatz is size consistent, as expected. 81The table clearly shows that the size consistency is deteriorated by the optimization using VMC gradients, i.e., the difference between the isolated and far-away energies is finite.In contrast, the size consistency is perfectly retrieved by the optimization using FN gradients.Neuscamman 63 pointed out that the deterioration of the size consistency comes from an incomplete Jastrow factor.More specifically, the real-space three/four-body Jastrow factor, which was employed in the present study, cannot completely remove the size consistency error unless we use unlimited flexibility in the Jastrow 63 .To solve the problem, Goetz and Neuscamman proposed the so-called number-counting Jastrow factors that can suppress the unfavorable ionic terms, and is able to solve the size-consistency problem 82,83 within the VMC framework.In this regard, our proposed scheme can be interpreted as an alternative approach because, again, FN-DMC is a projection method to relax the amplitude of the AGPn ansatz, which corresponds to adding an unlimited flexible Jastrow factor to a given ansatz.

Discussion
First, we compare our approach with others that also target to go beyond the single-reference fixed-node approximation.A well-established strategy is to use the multi-determinant ansatz, which has witnessed numerous successes so far [84][85][86][87][88][89][90][91][92] .The multi-determinant approach offers the advantage of systematic improvement by increasing the number of SDs.Nonetheless, the number of SDs for a comprehensive representation exponentially scales with system size, imposing substantial computational demands for large systems.5][86] However, there have been successful efforts to reduce the number of required determinants by neglecting less important ones 87,88 using, for instance, the configuration interaction using a perturbative selection made iteratively (CIPSI), which mitigates the exponential character of the multi-determinant approach in practice 90,92 .Recently, Benali et al. successfully applied the multi-determinant approach for solids with more than a hundred electrons by combining the CIPSI technique with a restricted active space built using natural orbitals 91 , which is a similar idea as we present in this study.Indeed, they demonstrated that one can go beyond the single-reference nodal surface in large systems by the multi-determinant approach in practice, though its naive asymptotic scaling is exponential.The multi-determinant approach is becoming as practical and promising as the single-determinant approach.
Concerning the actual computational costs of our proposed methods, the choice of ansatz (i.e., JSD or AGPn) does not significantly affect the cost of WF optimization, while the choice of gradients does.For instance, for C 60 , Jastrow optimization with the JSD ansatz and Jastrow+nodal surface (i.e., weights of Natural Orbitals) optimization with the JAGPn ansatz using VMC gradients require 11.9 and 43.6 cores • hours per optimization step with ∼ 7 mHa accuracy on the total energy evaluation at each optimization step, respectively 10 .
However, if one uses FN gradients for WF optimization, one needs more computational time. 10We measured the computational times on the Numerical Materials Simulator at National Institute for Materials Science (NIMS) using 1536 cores (32 nodes × Intel Xeon Platinum 8268 (2.9GHz, 24cores) × 2 per node).
For instance, for C 60 , the nodal surface (i.e., weights of Natural Orbitals) optimization with the JAGPn ansatz using FN gradients with a = 0.20 a.u.requires 195.3 cores • hours per optimization step with ∼ 7 mHa accuracy on the total energy evaluation at each optimization step.Thus, our FNAGPAS scheme using FN gradients shows the same scaling of the number of variational parameters as the single-reference FN DMC with JSD ansatz while it increases the prefactor of computational cost.
Based on the results obtained in this work so far, we finally discuss how to improve a fermionic ansatz in ab initio QMC calculations, in general.Recently, there have been many successful reports about machine-learning-inspired ansatz with a huge degree of freedom in describing electronic and spin states, such as deep neural networks 93 , restricted Boltzmann machines [94][95][96] and transformers 97 , which are utilized as ansatz of wave functions to solve the Schrödinger equation with lattice Hamiltonians.Also, in the ab initio community, ansatz using deep neural networks have been successfully applied for realistic problems, such as PauliNet 42 , FermiNet 39 , and others 40,41,[43][44][45] .In light of the present results, let us consider exploiting an ansatz with a huge degree of freedom (i.e., many variational parameters) in ab initio QMC calculations to pursue an exact fermionic ground state.If we stop at the VMC level, we may apply such a flexible ansatz to Jastrow factors, determinant part, or both parts, and it is expected that the larger the degree of freedom an ansatz has, the larger the energy gain should be.However, improvements at the VMC level do not necessarily lead to improvements at the FN level, especially if the determinant part is optimized at the variational level.A variational optimization improves the overall shape of the trial wave function Ψ T , whilst the nodal surface might not be as optimized as the Ψ T .In this work, indeed, we have seen how the JAGPn ansatz optimized at the FN level leads to much better results than the JAGP ansatz optimized at the VMC level, despite the latter having many more variational parameters and it is much better at the VMC level.Moreover, we also observed how the JAGPn (and JAGP, for that matter) ansatz itself yields a sizeconsistency error at the FN level if the parameters are optimized at the VMC level, while the same ansatz with parameters optimized at the FN level is not affected by this issue.
Thus, caution should be used when employing these new highly flexible machine-learningbased wave function parametrizations, as it is not guaranteed that improvements in the VMC energy are reflected in improvements in the FN energy in a consistent way.Basic physical properties, which were present in the most standard wave functions (such as the JSD), might not appear in the more fancy approaches, similar to the mentioned problem of size-inconsistency in JAGP and JAGPn.

Conclusions and perspectives
In this study, we propose a method for variational optimization of the AGP wavefunction expressed in terms of natural orbitals, with pairing coefficients optimized using FN gradients.
Within our scheme, the variational parameter space increases only linearly with the system size, as opposed to the quadratic scaling of the standard parametrization of AGP, with the result that our proposed method allows the optimization of the nodal surfaces for large systems, which has been difficult to achieve with conventional approaches.In addition to demonstrating that our scheme can be applied to systems as large as C 60 , we showed that our scheme also achieves better (i.e., lower) DMC energies than the single-reference fixed-node DMC calculations.Moreover, we have shown that our approach is size-consistent and can be used to estimate binding energies.
We showed that the Jastrow incompleteness affecting nodal surface optimizations can be mitigated by using FN gradients combined with the JAGPn ansatz.However, in this study, we did not investigate the effect of the basis set incompleteness on the determinant part (i.e., nodal surface).The basis set incompleteness is believed to be less severe in QMC calculations than in quantum chemistry methods because the Jastrow factor (at the variational level) or the projection (at the FN level) mitigates its error.However, to the best of our knowledge, no one has seriously investigated the error so far.Considering binding energy calculations done by DMC reported so far, 8 the basis set incompleteness should have a small effect on small molecules, but it should be carefully considered when studying large molecules using DMC done with localized basis sets.This is one of the intriguing future works for applying the single-reference-DMC and our proposed methods to large systems. (EP/F036884/1).
We dedicate this paper to one of the authors, Prof. Sandro Sorella (SISSA), who passed away during the collaboration.He initially devised an idea to use FN gradients for WF optimizations, which is one of the keys for the success of the present work.The QMC community will remember that he is one of the most influential contributors of the past and of the beginning of the present century to the community, and in particular for deeply inspiring this work with his developing of the ab initio QMC package, TurboRVB 57 .
Appendix: Proof for the variational principle of the LRDMC optimzation with DLA As pointed out in seminal works by Casula et al. 72,98 , the use of a pseudo potential that has the so-called non-local term induces an additional sign problem in the standard DMC approach with the locality approximation (LA); thus the variatioanl principle, which justifies the energy minimization strategy, is deteriorated.Instead, one of the advantages of the LRDMC is that the use of pseudo potentials does not deteriorate the variational princple 72 ; thus, the energy-minimization is justified.Recently, we implemented the Determinant Locality approximation (DLA) 18 into the TurboRVB package.In this study, we combine the DLA with the LRDMC framework implemented in the TurboRVB package.We prove here that the variational principle holds also in the LRDMC with the DLA.This proof is inspired by the proof by Haaf et al. 99 that the lattice Green's function Monte Carlo method is variational.
In LRDMC calculations with the DLA, the effective Hamiltonian (i.e., the fixed-node Hamiltonian) reads: where , by which the original term in the LRDMC approach, V sf x,x = for sf x ′ ̸ =x Ψ T (x ′ )H x,x ′ /Ψ T (x), is replaced, and sf means that all x ′ (̸ = x) satisfying Ψ T (x ′ )H x,x ′ Ψ T (x) > 0. Here, we omit the lattice-space dependency of the Hamiltonian (i.e., H ≡ H a ) because one can extrapolate energies to the a → 0 limit.Notice that, we assume that a trial WF can be decomposed into the Jastrow and determinant parts, i.e., Ψ T = J T D T .We also notice that since the Jastrow factor does not affect the sign of a wavefunction.We define the following notations: where |Φ FN ⟩ is the fixed-node ground state of ĤFN .In the following, we will show the following equations hold: The first equal (E MA = E FN ) holds because |Φ FN ⟩ is the exact ground state of H FN (i.e. ĤFN |Φ FN ⟩ = E FN |Φ FN ⟩).This is also true with the non-local terms of pseudo potentials.Now we define the difference between the effective fixed-node energy obtained with the effective Hamiltonian ĤFN and that obtained with the true Hamiltonian Ĥ: We want to prove that ∆E ≥ 0 for the fixed-node state and the equal holds for Φ FN = Ψ T = Ψ 0 , where we denote Ψ 0 as the exact WF of the original Hamiltonian Ĥ, i.e., Ĥ |Ψ 0 ⟩ = E 0 |Ψ 0 ⟩.Hereafter, we will do the same exercise written in Ref. 99.We define the difference between the effective fixed-node energy obtained with the effective Hamiltonian H FN and that obtained with the true Hamiltonian H: where, we define a truncated Hamiltonian H tr and a spin-flip Hamiltonian H sf , by and Indeed, the matrix elements are: and ∆E can be written explicitly in terms of the matrix elements: and rewriten as: where, sf means that all x ′ (̸ = x) satisfying Ψ T (x ′ )H x,x ′ Ψ T (x) > 0. In this double summation, each pair of configurations (x, x ′ ) appear twice.Therefore, we can combine these terms and rewrite it as a summation over the pairs: Notice that the hamiltonian is hermitian: H x,x ′ = H x ′ ,x .Since all the pair satisfies D T (x ′ ) D T (x) H x,x ′ > 0 , then, Then, − sgn(x, x ′ )Φ * FN (x)Φ FN (x ′ ) − sgn(x, x ′ )Φ * FN (x ′ )Φ FN (x) (20)   where sgn(x, x ′ ) denotes the sign of H x,x ′ .Finally, we get: indicating that ∆E is positive for any wave function Φ FN .Thus, the ground-state energy of H FN is an upper bound for the ground-state energy of the original Hamiltonian H (i.e., E FN ≥ E).Hereafter, we consider the case that one uses the true ground-state Ψ 0 for the determinant of the trial wave-function (i.e., Ψ T = J T • Ψ 0 ), to prove that E FN = E holds with Ψ T = J T • Ψ 0 (i.e., D T = Ψ 0 ): For all the pairs (x, x ′ ), Ψ 0 H x,x ′ Ψ 0 > 0 is satisfied, meaning sgn(x, x ′ ) → + and Ψ 0 (x) Ψ 0 (x ′ ) → +, or sgn(x, x ′ ) → − and Ψ 0 (x) Ψ 0 (x ′ ) → −.Thus, the above condition is fulfilled when the following condition is satisfied: In the DLA approach, the spin-flip term is composed only of the determinant of the trial WF.
Therefore, the fixed-node outcome with the DLA approach is not affected by the presence of the Jastrow factor in the trial WF (in the a → 0 limit).Therefore, one gets Φ FN = Ψ 0 with Ψ T = J T • Ψ 0 .Thus, ∆E = 0 is fulfilled with Ψ T = J T • Ψ 0 , and the following relations hold: meaning that the effective Hamiltonian ĤFN and the true Hamiltonian Ĥ has the same ground-state energy E 0 and the same ground state Φ FN = Ψ 0 with Ψ T = J T • Ψ 0 , where the final equal E = E 0 comes from the usual variational principle.
In the DLA approach, we can update the trial WF Ψ T such that E MA goes down according to the gradient ∂ α E MA or using a more sophisticated optimization scheme.As written above, the equals E FN = E = E 0 are met when Ψ T = J • Ψ 0 .It implies that one can look for the true ground-state energy and wavefunction by variation of the determinant part of the trial wavefunction.Indeed, in the LRDMC calculations with the DLA, one can access the mixed-average energy E MA , and its derivative ∂ ⃗ α E MA , where ⃗ α is a set of the variational parameters.Since E MA satisfies the variational principle, i.e., E MA ≥ E 0 , the equal holds when Ψ T = J T • Ψ 0 , as proven above, one can update the determinant part of the trial WF, D T , such that E MA goes down, then, it is expected that D T finally reaches D T → Ψ 0 , and TOC Graphic The torsion energy is defined as the energy difference between the ground state ethylene (denoted as planer ethylene) and the orthogonally rotated ethylene (denoted as twisted ethylene).Table S-V shows the total energies of the ethylenes computed with the JSD and JAGPn ansatz.Table S-VI summarizes the obtained torsion energies and reference values obtained in previous works.

FNFigure 1 :
Figure1: Panel a: Schematic illustration of the FNAGPAS scheme.We perform a preliminary mean field calculation to obtain molecular orbitals (MOs), followed by a correlated calculation yielding natural orbitals (NOs).The AGPn ansatz correspond to a multideterminant expansion built on the NOs and depending on the coefficients λ i associated to each orbital i and optimized in order to minimise the FN energy.Panel b: Flowchart illustrating the FNAGPAS scheme workflow.

Figure 2 :
Figure 2: Torsion energies between the planar and twisted ethylene.The values of JAGP and MR-CISD+Q (horizontal broken line) are taken from Ref. 59 and Ref. 76, respectively.

Figure 4 :
Figure 4: The top panel shows the improvement, dubbed energy gain, of the JAGP (red) and JAGPn (blue) ansatz with respect to the traditional JSD ansatz for each of the considered systems, as a function of the number of valence electrons.The energy gain is difference between FN energy of the JSD ansatz and the JAGP (or JAGPn) ansatz.The bottom panel shows the number of parameters in the Jastrow, in the AGP and in the AGPn parts of the wave function.The dashed lines show the linear (gray for JSD, cyan for JAGPn) and quadratic (orange for JAGP) fitting curves.

Figure. S- 2
Figure.S-2 shows the schematic figure of the ethylene torsion.The torsion energy is defined as the energy difference between the ground state ethylene (denoted as planer ethylene) and the orthogonally rotated ethylene (denoted as twisted ethylene).TableS-V shows the total energies of the ethylenes computed with the JSD and JAGPn ansatz.TableS-VI summarizes the obtained torsion energies and reference values obtained in previous works.

S4FIG. S- 2 .
FIG. S-2.The schematic figure of the torsion between the planar and twisted ethylene.

Table 1 :
FN binding energy E b and size consistency energy error E SCE , computed with LRDMC a → 0, as obtained with the JSD, JAGPn and JAGP wave functions.For JAGPn we consider both the case of using VMC and FN gradients to optimize the nodal surface.The latter is the scheme dubbed FNAGPAS in this work.

Table 1
79,80ins the binding energies of the methane-water dimer computed with the JSD ansatz, with the JAGPn ansatz optimized using either VMC or FN gradients (the FNAGPAS approach), and with the JAGP ansatz optimized with VMC gradients.The binding energy is evaluated as the energy difference between the dimer and the sum of the energies of the two molecules: E b = E water-methane − E water − E methane .The reference value for the binding energy of the water-methane dimer, -27 meV, was computed by CCSD(T) implemented in Orca79,80

TABLE S -
V. Comparison of LRDMC energies (a → 0) of the planar and twisted ethylene.

TABLE S -
VI.The torsion energies between the planar and twisted ethylene, obtained with various approaches.