TurboRVB: a many-body toolkit for {\it ab initio} electronic simulations by quantum Monte Carlo

TurboRVB is a computational package for {\it ab initio} Quantum Monte Carlo (QMC) simulations of both molecular and bulk electronic systems. The code implements two types of well established QMC algorithms: Variational Monte Carlo (VMC), and Diffusion Monte Carlo in its robust and efficient lattice regularized variant. A key feature of the code is the possibility of using strongly correlated many-body wave functions. The electronic wave function (WF) is obtained by applying a Jastrow factor, which takes into account dynamical correlations, to the most general mean-field ground state, written either as an antisymmetrized geminal product with spin-singlet pairing, or as a Pfaffian, including both singlet and triplet correlations. This wave function can be viewed as an efficient implementation of the so-called resonating valence bond (RVB) ansatz, first proposed by L. Pauling and P. W. Anderson in quantum chemistry and condensed matter physics, respectively. The RVB ansatz implemented in TurboRVB has a large variational freedom, including the Jastrow correlated Slater determinant as its simplest, but nontrivial case. Moreover, it has the remarkable advantage of remaining with an affordable computational cost, proportional to the one spent for the evaluation of a single Slater determinant. The code implements the adjoint algorithmic differentiation that enables a very efficient evaluation of energy derivatives, comprising the ionic forces. Thus, one can perform structural optimizations and molecular dynamics in the canonical NVT ensemble at the VMC level. For the electronic part, a full WF optimization is made possible thanks to state-of-the-art stochastic algorithms for energy minimization. The code has been efficiently parallelized by using a hybrid MPI-OpenMP protocol, that is also an ideal environment for exploiting the computational power of modern GPU accelerators.

TurboRVB is a computational package for ab initio Quantum Monte Carlo (QMC) simulations of both molecular and bulk electronic systems. The code implements two types of well established QMC algorithms: Variational Monte Carlo (VMC), and Diffusion Monte Carlo in its robust and efficient lattice regularized variant. A key feature of the code is the possibility of using strongly correlated many-body wave functions, capable of describing several materials with very high accuracy, even when standard mean-field approaches (e.g., density functional theory (DFT)) fail. The electronic wave function (WF) is obtained by applying a Jastrow factor, which takes into account dynamical correlations, to the most general mean-field ground state, written either as an antisymmetrized geminal product with spin-singlet pairing, or as a Pfaffian, including both singlet and triplet correlations. This wave function can be viewed as an efficient implementation of the so-called resonating valence bond (RVB) ansatz, first proposed by L. Pauling and P. W. Anderson in quantum chemistry and condensed matter physics, respectively. The RVB ansatz implemented in TurboRVB has a large variational freedom, including the Jastrow correlated Slater determinant as its simplest, but nontrivial case. Moreover, it has the remarkable advantage of remaining with an affordable computational cost, proportional to the one spent for the evaluation of a single Slater determinant. Therefore, its application to large systems is computationally feasible. The WF is expanded in a localized basis set. Several basis set functions are implemented, such as Gaussian, Slater, and mixed types, with no restriction on the choice of their contraction. The code implements the adjoint algorithmic differentiation that enables a very efficient evaluation of energy derivatives, comprising the ionic forces. Thus, one can perform structural optimizations and molecular dynamics in the canonical NVT ensemble at the VMC level. For the electronic part, a full WF optimization (Jastrow and antisymmetric parts together) is made possible thanks to state-ofthe-art stochastic algorithms for energy minimization. In the optimization procedure, the first guess can be obtained at the mean-field level by a built-in DFT driver. The code has been efficiently parallelized by using a hybrid MPI-OpenMP protocol, that is also an ideal environment for exploiting the computational power of modern GPU accelerators. The solution of the many-body Schrödinger equation, which describes the interaction between electrons and ions at the quantum mechanical level, represents a fundamental challenge in computational chemistry, condensed matter, and materials science. Since about a century, there has been a relentless theoretical and computational effort to find an accurate solution to this problem, which also features many recent interdisciplinary applications in machine learning 1 and materials informatics 2 . While computationally, the scaling of the problem is exponential with the number of electrons, several numerical approximate methods have been put forward in the last decades. Among them, the Density Functional Theory (DFT) method, proposed by Kohn and Sham 3 in 1965, is one of the most successful approaches. In this framework, the original interacting 3N many-body problem (N being the number of electrons in a system) is mapped to a noninteracting electron system, defined by an effective mean-field potential, to be determined self-consistently 4 . While DFT is an exact theory in principle, the exact form for the exchangecorrelation functional, which is an essential part of the DFT mean-field potential, remains unknown. Unfortunately, the progress in generating increasingly successful approximations of this functional is rather slow 5 , partly because there is no established strategy for systematic improvement 6 , while maintaining an efficient scaling with the system size. The commonly adopted approximations for the exchange-correlation DFT functional have well-known limitations especially in describing weak dispersive interactions, strongly correlated materials, and extreme environments (e.g., high pressure 7,8 ).
Alternative strategies are represented by the so-called WF-based approaches popular in quantum chemistry applications 9 .
Among them, coupled-cluster with single, double, and perturbative triple excitations, or CCSD(T), is considered to be the gold standard in quantum chemistry as it typically provides results in good agreement with experiments, and an optimal balance between accuracy and computational affordability (despite its cost grows as the seventh power of the number of electrons). While quantum chemistry methods have the advantage of treating electronic exchange and correlation effects in a systematically improvable fashion, they are also much more computationally demanding compared to DFT, and their applicability to large or periodic systems is often computationally prohibitive.
Another way to tackle the problem of the electron correlation and the huge dimension of a many-body WF, exponentially large in the number of electrons, is by means of stochastic approaches, in this context referenced with the widely used expression of quantum Monte Carlo (QMC) methods 14,15 . Since the invention of Markov Chain Monte Carlo (MCMC) in the 1940s, 16 stochastic approaches to numerical algorithms had a pervasive influence in a wide range of fields, from physics, engineering to finance. In this respect, the QMC framework is qualitatively different from the deterministic methods mentioned above, and it represents, therefore, an original and alternative approach for the solution of the Schrödinger equation, overcoming some of the drawbacks of DFT and the deterministic WF-based approaches of quantum chemistry. In particular, QMC does not rely on uncontrolled approximations; its accuracy can be systematically improved, the scaling with system size is good, and it is straightforwardly applicable to both isolated and periodic systems. While the scaling with the system size is comparable with standard DFT methods, the prefactor is typically much larger. However, methods based on stochastic sampling are well suited for massively parallel computing architectures, as the algorithms can sustain an almost ideal scaling with the number of cores. As a result, the feasibility and popularity of QMC are expected to increase with the foreseeable substantial improvements in high-performance computing (HPC) facilities.
TurboRVB includes two of the most popular ground-state QMC algorithms: variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC) 14 . In VMC, a systematically improvable approximation for the ground state is obtained by direct minimization of the energy, evaluated by a parametrized many-body WF. In an explicitly correlated WF, the electron coordinates are not separable, and the expectation value of the Hamiltonian has to be calculated with Monte Carlo integration, hence the name. In realistic calculations, a faithful and nevertheless compact parametrization for the trial state is essential for a successful energy optimization. The first property allows an unbiased treatment of the electronic correlations across different regimes, such as bond dissociation and electronic phase transitions. The second is needed for a stable optimization of the variational parameters, which may become a too difficult task when the chosen parametrization is redundant or unnecessarily detailed. Therefore, a central effort, in the TurboRVB project, has been devoted to the development of efficient and systematic parametrizations of cor-related WFs.
As we will see in the following, an accurate trial WF is also fundamental in DMC, which is the second main algorithm present in TurboRVB. DMC is an imaginary-time projection technique 15 , that, when combined with the so-called fixed-node (FN) approximation, represents a powerful route to electronic ground-state calculations. In this approximation, the nodal-surface, which is kept fixed during the projection, is determined by an accurate variational optimization.
A typical workflow for a QMC calculation performed with TurboRVB is shown in Fig. 1. It shall be noticed that the choice of the most suited ansatz for the WF, as well as the optimization of its parameters, is a prerequisite of both VMC and DMC evaluations. Ionic forces, if evaluated, can be used to perform QMC-based structural relaxations, Langevin molecular dynamics (MD) or path integral molecular dynamics (PIMD). Moreover, several post-processing tools are available to analyze the outcomes. The TurboRVB code is peculiar because it relies on the assumption that a complex electronic problem, though, until now, cannot be solved exactly with a black-box computer software, can be described very accurately by approriate variational wavefunctions, guided by robust physical and chemical requirements and human ingenuity. In particular TurboRVB is different from other codes in the following features: ( i ) it employs the Resonating Valence Bond (RVB) WF, which includes statical and dynamical correlation effects beyond the commonly used Slater determinant, by keeping the computational cost at the singledeterminant level, thanks to its efficient implementation; (ii) the code implements a VMC algorithm based on localized orbitals (e.g., Gaussians) and state-of-the-art optimization methods, such as the stochastic reconfiguration. Therefore, at the VMC level, one can optimize not only the amplitude of the WF (i.e., the Jastrow factor), but also the nodal surfaces (e.g., the Slater determinant). This leads to a better variational energy in general, and also improves the corresponding FN-DMC energy; (iii) the energy derivatives (e.g., atomic forces) are calculated very efficiently thanks to an implementation based on the Adjoint Algorithmic Differentiation (AAD). As a consequence, one can perform structural optimizations and Langevin molecular dynamics; (iv) The code implements the newly developed Lattice Regularized Diffusion Monte Carlo (LRDMC), a stable DMC algorithm that, very recently, has shown to have a better scaling with the atomic number Z, compared with standard DMC. 35 This review is organized as follows: in Sec. II, we briefly explain the fundamental QMC algorithms implemented in TurboRVB, namely, the VMC and LRDMC; in Sec. III, we describe the RVB WF; in Sec. IV, we extend the WF to treat periodic systems; in Sec. V we introduce the DFT algorithm efficiently implemented in TurboRVB; in Sec. VI, we describe the way energy derivatives are computed (e.g., atomic forces) by means of AAD; in Sec. VII, we list the TurboRVB, steps to optimize a many-body WF; in Sec. VIII, we introduce the first and the second-order Langevin molecular dy-namics implemented in TurboRVB; in Sec. IX, we summarize the typical QMC workflow from the WF generation to the final LRDMC calculation; in Sec. X, we show weak and strong scaling results of TurboRVB, measured on the Marconi/CINECA supercomputer; in Sec XI, we list the physical properties that can be calculated by TurboRVB; in Sec. XII, we review the major TurboRVB applications done so far; in Sec. XIII, we introduce a python-based workflow system, named Turbo-Genius.

II. METHODS
TurboRVB implements two types of well established Quantum Monte Carlo methods: Variational Monte Carlo (VMC) and Lattice Regularized Diffusion Monte Carlo (LRDMC). We summarize these methods in this section. The interested readers should also refer to the comprehensive re-view of QMC 14 for details. TurboRVB also implements an original DFT engine to generate trial WFs, that are more suitable for the VMC and LRDMC calculations, as shown in Sec. V.

A. Variational Monte Carlo
Starting from the variational principle, the expectation value of the energy evaluated for a given WF Ψ can be written Jastrow factor, composed of J 1 , J 2 , and J 3/4 . Ψ T Many-body trial WF. Ψ G Many-body guiding function.
Pairing function between electrons i ↔ (r i , σ i ) and j ↔ (r j , σ j ). f r i , r j Spatial part of the pairing function. ψ i (r) Primitive atomic orbital for the determinant part. φ k (r) k-th molecular orbital for the determinant part.
Primitive atomic orbital for the Jastrow part.
Variance-covariance matrix of the logarithmic derivative of a many-body WF. S QMC (R) Variance-covariance matrix of QMC forces. as: where x = (r 1 σ 1 , r 2 σ 2 , . . . r N σ N ) here and henceforth is a shorthand notation for the N electron coordinates and their spins, whereas are the so-called local energy and the probability of the configuration x, respectively. This multidimensional integration can be evaluated stochastically by generating a set {x i } according to the distribution π (x) using the Markov chain Monte Carlo such as the (accelerated 36,37 ) Metropolis method, and by averaging the obtained local energies e L (x i ): which has an associated statistical error of is the variance of the sampled local energies, andM is the sampling size M divided by the autocorrelation time. This indicates that the precision of the VMC evaluation is inversely proportional to the square root of the number of samplings (i.e., of the computational cost). It worth to notice that, if Ψ(x) is an eigenfunction ofĤ, say with eigenvalue E 0 , than e L (x) = E 0 for each x, implying that the variance of the local energy is zero and E VMC = E 0 with no stochastic uncertainty. This feature is known as the zero-variance property. The probability distribution used for the importance sampling can also differ from π (x). Indeed, one can use an arbitrary probability distribution function π ′ (x) = Ψ 2 G (x) / dxΨ 2 G (x), and estimate a generic local observable O (x) either by usingŌ VMC = O (x) π(x) or: where W (x) = |Ψ (x) /Ψ G (x) | 2 , and the points x ′ i are distributed according to π ′ . This reweighting scheme is very important when evaluating forces with finite variances, as discussed in Sec. VI B. Since the evaluations of the standard deviation is nontrivial in this case, TurboRVB employs the bootstrap and jackknife methods in order to estimate the mean value and the statistical error 15 , which are also used when evaluating those of the local energy, forces, and so on. Indeed, the code outputs the history of the local energies, forces, or other properties in appropriate files (when the corresponding option is true), thus allowing the error bar estimates by simple post-processing. The user can also use the reblocking (binning) technique to remove the autocorrelation bias. 38 One can optimize the WF based on the variational theorem by introducing a set of parameters α 1 , α 2 , · · · , α p to the WF Ψ (x, α): However, the optimization of a many-body WF remains a difficult challenge not only because optimizing a cost function containing several parameters is a complex numerical task, due to the presence of several local minima in the energy landscape, but also because this difficult task is further complicated by the presence of statistical errors in the QMC evaluation of any quantity. Nevertheless, a great improvement in this field has been achieved when the QMC optimization techniques have made use of the explicit evaluation of energy derivatives with finite statistical errors. In particular, in TurboRVB the adjoint algorithmic differentiation (AAD) has been implemented, by allowing the efficient calculation of generalized forces ( f i = − ∂E ∂α i ) 39 , and very efficient optimization methods, the so-called "stochastic reconfiguration" 40,41 and "the linear method" [42][43][44] . These methods are discussed in Sec. VII A and VII B.

B. Lattice regularized Monte Carlo
Lattice regularized diffusion Monte Carlo (LRDMC), which was initially proposed by M. Casula et al. 45 , is a projection technique that allows us to improve a variational ansatz systematically. This method is based on Green's function Monte Carlo (GFMC) [46][47][48] , filtering out the ground state WF |Υ 0 from a given trial WF |Ψ T : since the eigenstates of the Hamiltonian have the completeness property, the trial WF can be expanded as: where a n is the coefficient for the n-th eigenvectors (Υ n ). Therefore, by applying G M = Λ −Ĥ M , one can obtain where Λ is a diagonal matrix with Λ x,x ′ = λδ x,x ′ (λ should be sufficiently large to obtain the ground state), and E n is nth eigenvalue ofĤ. Since Λ−E n Λ−E 0 < 1, the projection filters out the ground state WF Υ 0 from a given trial WF |Ψ T , as long as the trial WF is not orthogonal to the true ground state (i.e., a 0 ≡ Ψ T |Υ 0 0). To apply the GFMC for ab initio electron calculations, the original continuous Hamiltonian is regularized by allowing electron hopping with step size a, in order to mimic the electronic kinetic energy. The corresponding HamiltonianĤ a is then defined such thatĤ a →Ĥ for a → 0. Namely, the kinetic part is approximated by a finite difference form: and the potential term is modified as: The corresponding Green's function matrix elements are: and the single LRDMC iteration step is given by the following equation: The sketch of the LRDMC algorithm, a Markov chain that evolves the many-body WF according to the Eq. 10, is as follows 45 : (STEP 1) Prepare a walker with configuration x and weight w (w 0 = 1). (STEP 2) A new configuration x ′ is generated by the transition probability: where is a normalization factor. By applying the discretized Hamiltonian to a given configuration (Ĥ a |x ), (6N + 1) configurations |x ′ are determined according to the probability p x,x ′ in Eq. 11, where N is the number of electrons in the system 49 . This allows the evaluation of the normalization factor b x in Eq. 12 even in a continuous model. Notice that 6N comes from the diffusion of each electron in two directions (±a) and the remaining 1 stands for the starting configuration before the possible hoppings (all N electrons) x (i.e., x ′ = x). (STEP 3) Finally, update the weight with w n+1 = w n b x , and return to the STEP I. After a sufficiently large number of iterations (the Markov process is equilibrated), one can calculate the ground state energy E 0 : where · · · denotes the statistical average over many independent samples generated by the Markov chains, and e L (x) is called the (bare) local energy that reads: Indeed, the ground state energy can be calculated after many independent n-step calculations. A more efficient computation can be realized by using the so-called "correcting factor" technique: after a single simulation that is much larger than the equilibration time, one can imagine starting a projection of length p from each (n − p) th iteration. The accumulated weight for each projection is: Then, the ground state energy E 0 can be estimated by: This straightforward implementation of the above simple method is not suitable for realistic simulations due to fluctuations of weights, large correlation times, the sign problems, and so on. TurboRVB implements the following state-of-art techniques for real electronic structure calculations. If the potential term (Eq. 8) is unbounded (it is the case in ab initio calculations), the bare weight b x (Eq. 12) and the local energy e L (x) (Eq. 14) significantly fluctuate, making the numerical simulation very unstable and inefficient. To overcome this difficulty, the code employs the importance sampling scheme 15 , in which the original Green's function is modified using the so-called guiding function Ψ G as: and the projection is modified as: In practice, the guiding function is prepared by a VMC calculation. The modified Green's function for importance sam-plingG x ′ ,x has the same eigenvalues as the original one, and this transformation does not change the formalism of LRDMC. The weight is updated by: and the local energy with importance sampling is: Eq. 20 implies that if the guiding function Ψ G is an exact eigenstate of the Hamiltonian, there are no statistical fluctuations, implying the zero variance property, namely the computational efficiency to obtain a given statistical error on the energy improves with the quality of the variational WF. In this respect, it is also important to emphasize that a meaningful reduction of the statistical fluctuations is obtained by satisfying the so-called cusp conditions. As long as they are satisfied, the resulting local energy does not diverge at the coalescence points where two particles are overlapped, despite the singularity of the Coulomb potential term (V(x) in Eq. 8) 49 . In addition, the importance sampling maintains the electrons in a region away from the nodal surface, since the guiding function vanishes there (i.e., the RHS of Eq.17 → 0). This clearly enhances the efficiency of the sampling because the local energy diverges at the nodal surface. The Green's function cannot be made strictly positive for fermions; therefore, the fixed-node (FN) approximation has to be introduced 15 in order to avoid the sign problem. Indeed, the Hamiltonian is modified using the spin-flip term V sf (x) = and γ ≥ 0 is a real parameter. The use of the fixed-node Green's function: can prevent the crossing of regions where the configuration space yields a sign flip of the Green's function; therefore, the walkers are constrained in the same nodal pockets, and avoid the sign problem. TurboRVB also implements the many-walker technique and the branching (denoted as reconfiguration 50 in Tur-boRVB) scheme for a more efficient computation 15 . The code performs the branching as follows: (1) Set the new weights equal to the average of the old ones: (2) Select the new walkers from the old ones with a probability that is proportional to the old walkers' weights: which does not change the statistical average of weights, but suppresses the fluctuations by dropping walkers having small weights. The code performs branching (reconfiguration) after a projection time τ, that can be chosen as a user input parameter. In practice, within the many-walker and branching schemes, the average weights are stored and are set to one for all walkers after each branching. The user can retrieve the accumulated weights at the end of the simulation: and calculate the ground state energy: where e L (x n ) is the mean local energy averaged over the walkers, which reads: and is evaluated just before each reconfiguration. Notice that p is also an input parameter, that has to be carefully chosen by the user to allow energy convergence.
When λ is sufficiently large, the correlation time also becomes large because the diagonal terms of the Green's function become very close to one (i.e., a walker remains in the same configuration), which causes a very large correlation time. In TurboRVB, this difficulty is solved by considering in a different way the diagonal and non-diagonal moves. In a given interval of M iterations, The non-diagonal updates are efficiently calculated using a random number according to the probability that the configuration remains in the same one (diagonal updates). This technique can be generalized to the continuous-time limit, namely, M → ∞, at M Λ = τ fixed. In the M → ∞ limit, the projection Λ −Ĥ M is equal to the imaginary time evolution exp −τĤ , apart for an irrelevant constant Λ M . Thus the user should specify only τ as input parameter. Indeed, the walker weight is updated by w → w exp −τ ξ e L (x) and the imaginary time is updated τ left → τ left − τ ξ at each non-diagonal update until τ left becomes 0, where τ ξ = − log (1 − ξ) /b x is a diagonal move time step determined by a uniform random number 0 ≤ ξ < 1. The branching (reconfiguration) is performed each τ that a user puts in the input file within the many walker and the branching implementation.
In practice, there are three important features in LRDMC. First, there is not a time-step error in LRDMC because the Suzuki-Trotter decomposition is not necessary, unlike the standard DMC algorithm 15 . Instead, there is a finite-size lattice error due to the discretization of the Hamiltonian (a). Therefore, in order to obtain an unbiased FN energy, it is important to extrapolate the LRDMC energy to the a → 0 limit by using several results corresponding to different lattice spaces 45 . This is then consistent with the standard DMC energy estimate (Fig. 2) obtained in the limit of an infinitely small time step. Probably one of the most important advantages of the LRDMC method is that the extrapolation to the a → 0 limit is very smooth and reliable, so that unbiased FN energies are easily obtained with low order polynomial fits. Secondly, LRDMC can straightforwardly handle different length scales of a WF by introducing different mesh sizes (a and a ′ ), so that electrons in the vicinity of the nuclei and those in the valence region can be appropriately diffused 35,45 , which defines the so-called double-grid LRDMC. This scheme saves a substantial computational cost in all-electron calculations, especially for a system including large atomic number atoms 35 , with a typical computational cost scaling with Z ∼6 where Z is the maximum atomic number. TurboRVB makes use of an appropriate ratio of the mesh sizes (i.e., a/a ′ ) the smaller one a used when electrons are close to the nuclei and the larger one a ′ ≫ a adopted in the valence region. By choosing a proper Thomas-Fermi characteristic length around the nuclei, where short hops of lengths a mostly occur, a significant improvement of the scaling (i.e., Z ∼5.6 → Z ∼5 ) has been recently reported. 35 Finally, the inclusion of non-local pseudopotentials in this framework is straightforward by means of an additional spherical grid defined in an appropriate mesh 49 . As described in Ref. 45,51, LRDMC provides an upper bound for the true ground-state energy and allows the estimation of E FN , even in the pres-ence of non-local pseudopotentials. Notice that this variational property has also been extended to the standard DMC framework 52 , with the introduction of the so-called T-moves. Moreover, the recently introduced determinant locality approximation (DLA) 53 to deal with non-local pseudopotentials is also implemented in TurboRVB and can be optionally used in LRDMC.

III. WAVE FUNCTIONS
Both the accuracy and the computational efficiency of QMC approaches crucially depend on the WF ansatz. The optimal ansatz is typically a tradeoff between accuracy and efficiency. On the one side, a very accurate ansatz can be involved and cumbersome, having many parameters and being expensive to evaluate. On the other hand, an efficient ansatz is described only by the most relevant parameters and can be quickly and easily evaluated. In particular, in the previous sections, we have seen that QMC algorithms, both at the variational and fixed-node level, imply several calculations of the local energy e L (x) and the ratio Ψ(x)/Ψ(x ′ ) for different electronic configurations x and x ′ . The computational cost of these operations determines the overall efficiency of QMC and its scaling with systems size.
TurboRVB employs a many-body WF ansatz Ψ which can be written as the product of two terms: where the term exp J, conventionally dubbed Jastrow factor, is symmetric under electron exchange, and the term Φ AS , also referred to as the determinant part of the WF, is antisymmetric. The resulting WF Ψ is antisymmetric, thus fermionic.
In the majority of QMC applications, the chosen Φ AS is a single Slater determinant (SD) Φ SD , i.e., an antisymmetrized product of single-electron WFs. Clearly, SD alone does not include any correlation other than the exchange. However, when a Jastrow factor, explicitly depending on the interelectronic distances, is applied to Φ SD the resulting ansatz Ψ JSD = Φ SD * exp J often provides over 70% of the correlation energy 54 at the variational level. Thus, the Jastrow factor proves very effective in describing the correlation, employing only a relatively small number of parameters, and therefore providing a very efficient way to improve the ansatz. 55 A Jastrow correlated SD (JSD) function yields a computational cost for QMC simulations -both VMC and FN -about ∝ N 3 , namely the same scaling of most DFT codes. Therefore, although QMC has a much larger prefactor, it represents an approach much cheaper than traditional quantum chemistry ones, at least for large enough systems.
While the JSD ansatz is quite satisfactory in several applications, there are situations where very high accuracy is required, and a more advanced ansatz is necessary. The route to improve JSD is not unique, and different approaches have been attempted within the QMC community. First, it should be mentioned that improving the Jastrow factors is not an effective approach to achieve higher accuracy at the FN level, as the Jastrow is positive and cannot change the nodal surface. A popular approach is through the employment of backflow 56 , which is a remapping of the electronic configurations that enters into Φ AS (SD as a special case) where each electron position is appropriately changed depending on nearby electrons and nuclei. Backflow is an effective way to recover correlation energy, both at the variational and FN level. However, it can be used at a price to increase significantly an already large computational cost. 57 Another possibility is to improve Ψ SD similarly to conventional quantum chemistry approaches, namely by considering Φ AS as a linear expansion of several Slater determinants. While this second approach can provide very high accuracy, it may be extremely expensive, as the number of determinants necessary to remain with a predefined accuracy grows combinatorially with the system size.
The vision embraced in TurboRVB, is that the route toward an improved ansatz should not compromise the efficiency and good scaling of QMC. For this reason, neither backflow nor explicit multideterminant expansions are implemented in the code. Within the TurboRVB project, the main goal is instead to consider an ansatz that can be implicitly equivalent to a multideterminant expansion, but remains in practice as efficient as a single determinant. There are five alternatives for the choice of Φ AS , which correspond to ( i ) the Pfaffian (Pf), (ii) the Pfaffian with constrained number of molecular orbitals (Pfn) (iii) the Antisymmetrized Geminal Power (AGP), (iv) the Antisymmetrized Geminal Power with constrained number of molecular orbitals (AGPn), and (v) the single Slater determinant. It is interesting to observe that all the latter four WFs are obtained by introducing specific constraints on the most general Pf ansatz. The hierarchy of the five ansätze is represented in the Venn diagram of Fig. 3. Clearly, a more general ansatz is more accurate in the total energy but not necessarily in the energy differences. Moreover, it is described by more variational parameters, that could imply a more challenging optimization (see next section) and a slightly higher cost. TurboRVB includes several tools to go from one ansatz to another, as represented in Fig. 4. Typically the starting point is SD, which can be obtained from a Hartree-Fock (HF) or a DFT calculation with different exchange-correlation functionals. Both methods are not expected to provide the optimal parameters when the Jastrow factor is included in the WF. Indeed, in QMC, the WFs are always meant to include the Jastrow factor, which proves fundamental to improve the properties of the overall WF. For instance, AGP carries more correlation than SD. However, it is not size-consistent unless it is multiplied by a Jastrow factor. Thus, a fundamental step to take advantage of the WF ansatz is the possibility to perform reliable optimizations of the parameters. Optimization will be discussed in section VII. In this section, we will de-scribe the functional form of the Jastrow factor implemented in TurboRVB (sec. III A), the Pfaffian (sec. III B), the AGP (sec. III C), the AGPn and Pfn (sec. III D), the SD (sec. III E), the multiconfigurational character of the AGP (sec. III F), the basis set (sec. III G, the pseudopotentials (sec. III H), the contractions of the orbitals (sec. III I), and the conversion tool (sec. III J).

A. Jastrow factor (J)
The Jastrow factor (exp J) plays an important role in improving the correlation of the WF and in fulfilling Kato's cusp condition 58 . TurboRVB implements the Jastrow term composed of one-body, two-body, and three/four-body factors (J = J 1 + J 2 + J 3/4 ). The one-body and two-body factors are used to fulfill the electron-ion and electron-electron cusp conditions, respectively, and the three/four-body factors are employed to consider a systematic expansion, in principle converging to the most general independent electron pairs contribution. The one-body Jastrow factor J 1 is the sum of two parts, the homogeneous part (enforcing the electron-ion cusp condition): and the corresponding inhomogeneous part: where r i are the electron positions, R a are the atomic positions with corresponding atomic number Z a , l runs over atomic orbitals χ a,l (e.g., GTO) centered on the atom a, σ i represents the electron spin (↑ or ↓), {M σ i a,l } are variational parameters, and u a (r) is a simple bounded function. In TurboRVB, the most common choice for u a is: depending on a single variational parameter b ea , that may be optimized independently for each atomic species. The two-body Jastrow factor is defined as: where v σ i ,σ j is another simple bounded function. There are several possible choices for v σ i ,σ j implemented in TurboRVB (all listed in the file input.tex in the doc folder), and one of them is, for instance, the following spin-dependent form: where r i, j = r i − r j , and b para ee and b anti ee are variational parameters.
The three/four-body Jastrow factor reads: where the indices l and m again indicate different orbitals centered on corresponding atoms a and b, and {M σ i ,σ j a,l,b,m } are variational parameters. Sometimes it is convenient to set to zero part of the coefficients of the four-body Jastrow factor, namely those corresponding to a b, as they increase the overall variational space significantly and make the optimization more challenging, without being much more effective in improving the variational WF.

B. Pfaffian Wave function (Pf)
The SD is an antisymmetrized product of single-electron WFs. Thus, SD neglects almost entirely the correlation between electrons. A natural way to improve this description is to include explicitly in the ansatz pairwise correlations among electrons. This is precisely what the Pfaffian WF does. [59][60][61] The building block of the Pfaffian WF is the pairing function g(i, j) between any pair of electrons i and j. Henceforth we denote with the generic bolded index i both the space r i coordinates and the spin values σ i : corresponding to the i th electron. For simplicity, let us first consider a system with an even number N of electrons. The WF, written in terms of pairing functions, is: where A is the antisymmetrization operator: S N the permutation group of N elements,P the operator corresponding to the generic permutation P, and ǫ P its sign. Let us define G the N×N matrix with elements G i, j = g(i, j). Notice that as a consequence of the statistics of fermionic particles, thus G is skew-symmetric (i.e., G T = −G, being T the transpose operator), so the diagonal is zero and the number of independent entries is N−1 n=1 n = N(N − 1)/2. The Pfaffian 62 of a N × N skew-symmetrix matrix G is defined as: (1),P(2) · · · G P(N−1),P(N) (39) if N is even, and it is zero if N is odd. 63 Therefore, the WF Φ AS defined in the right-hand side (RHS) of Eq. 36 equals 1 N!! Pf(G) where the semifactorial N!! ≡ N! 2 N/2 (N/2)! is irrelevant in QMC, as it affects only the normalization of the WF. Thus, we can define our electronic WF as: Notice that the Φ Pf here defined allows the description of any system with N u electrons with spin-up and N d electrons with spin-down, provided that N = N u + N d is even. Indeed, with no loss of generality, we can assume that electrons i = 1, . . . , N u have σ i = 1/2 and electrons with i = N u + 1, . . . , N have σ i = −1/2. Thus, the N × N skew-symmetric matrix G is written as: where G uu is a N u × N u skew-symmetric matrix with elements g uu (i, j), G dd is a N d × N d skew-symmetric matrix with elements g dd (i, j), G ud is a N u × N d matrix with elements g ud (i, j), and G du = −G ud T , i.e., g du (i, j) = −g ud (j, i). g uu describes the pairing between a pair of electrons with spin-up: where the function f uu describes the spatial dependence on the coordinates r i , r j for i, j ≤ N u . Notice that f uu (r j , r i ) = − f uu (r i , r j ) as a consequence of the properties of g. The spin part |↑↑ describes a system with unit total spin and spin projection along the z-axis, and will be indicated by |1, +1 . Similarly, g dd describes the pairing between pairs of electrons with spin-down for i, j > N u : with f dd (r j , r i ) = − f dd (r i , r j ), and the spin part |↓↓ describes a system with total unit spin and negative spin projection along the z-axis, indicated with |1, −1 . g ud describes the pairing between pairs of electrons with unlike spins. Since two electrons with unlike spins can form a singlet |0, 0 = |↑↓ −|↓↑ √ 2 or a triplet |1, 0 = |↑↓ +|↓↑ √ 2 , in the general case the pairing function g ud will be a linear combination of the the two components: scribes the spatial dependence of the triplet part. Therefore, the generic pairing function g(i, j) is the sum of all the four components mentioned above, namely : The pairing functions f S , f T , f uu , and f dd are expanded over atomic orbitals (see sec. III G). Say, for a generic pairing function f we have where ψ a,l and ψ b,m are primitive or contracted atomic orbitals, their indices l and m indicate different orbitals centered on atoms a and b, while i and j label the electron coordinates. Symmetries on the system, or properties of the underlying pairing function f imply constraints on the coefficients. For instance, the coefficients of f S are such that Let us consider now the remaining case of a system with an odd number of electrons N. The simplest way to handle this case is to consider a system with an extra fictitious electron N + 1 ↔ (r N+1 , σ N+1 ) that we set at infinity r N+1 → ∞, thus non interacting with all the N physical electrons. The extra matrix elements are easily computed: where φ(i) ≡ φ(r i , σ i ) can be considered an extra spin dependent unpaired orbital vanishing at infinity. The antisymmetric part of the WF of the overall system is then: where the antisymmetrization operator acts on the N electrons in the first like, and on the N + 1 electrons in the second line. This is a perfectly allowed N electron fermionic WF, because by definition: ( i ) it is antisymmetric for all permutations of the N + 1 particles and in particular for the N− physical ones; (ii) it does not depend by definition on the extra N + 1 th coordinates and spin as the fictitious particle is at infinity. On the other hand, it is easy to show that this is a nontrivial and nonvanishing WF, because we can use the basic Pfaffian formula in Eq. 39 for the antisymmetrization of the N + 1 particles and we readily obtain that: Its Pfaffian is quite generally nonvanishing for at least some configuration provided the diagonalization of the Pfaffian (see later) has at least (N + 1)/2 non-zero eigenvalues (likewise for the Slater determinant it is enough that the N−molecular orbitals are linearly independent) and we obtain in this way that the Pfaffian is perfectly defined even in the odd number of electron case. Finally, we would like to emphasize that the above argument can be generalized to arbitrary number k of unpaired orbitals and arbitrary boundary conditions including calculations on periodic supercells. We can indeed assume at the beginning that, when N is odd (even), we add an odd (even) number k of unpaired orbitals, such that there are N/2 pairs of electrons plus k unpaired electrons. The k unpaired electrons are paired with k fictitious electrons (with a pairing function satisfying conditions analogous to those in Eq. 47). This yields a Pfaffian of even linear dimension N + k of the same form as the one in Eq. 50 but with Φ being a block N × k rectangular matrix, determined by k different spin-dependent orbitals. The same argument holds that the antisymmetrized WF is written as the Pfaffian of the corresponding (N + k) × (N + k) skew-symmetric matrix, and does not depend on the fictitious particle coordinates introduced and therefore is an allowed physical WF antisymmetric over the N physical particles.
A special case that is of interest is when all the N electrons are set to be unpaired (i.e., when k = N and Φ is a generic N × N matrix) and it is assumed that there is no pairing among them (i.e., G = 0 in the N × N block matrix in Eq. 50). Thus, using a well-known relation of a Pfaffian we recover the standard Slater determinant. In this way, we can clearly see that the Pfaffian represents a quite remarkable generalization of the single determinant picture, and that a larger variational freedom is exploited by allowing pairing correlations with a non zero matrix G.
Recently we have noticed that, following the above derivation, it is possible to generalize further the Pfaffian WF ansatz by introducing fictitious particles and relaxing the constraints on theG matrix, even in the lower block diagonal part, that can be assumed to be an arbitrary non-zero skew-symmetric matrix. Work is in progress in order to understand whether this extended formulation can further improve the accuracy within a single Pfaffian WF.

C. Antisymmetrized Geminal Power (AGP)
If we consider only the case of a pairing function g(i, j) that is a spin singlet (namely, f uu , f dd f T in Eq. 45 are set to zero, yielding g(i, j) = f S (r i , r j ) |0, 0 ) then we obtain the singlet Antisymmetrized Geminal Power (denoted as AGPs).
Let us consider first an unpolarized system, having an even number N on electrons, and without loss of generality, we can assume that the electrons i = 1, . . . , N/2 have spin up and electrons j = N/2 + 1, . . . , N have spin down. Then the matrices G uu and G dd defined in the previous section are both zero matrices of size N/2 × N/2, and the matrix G ud has only the contribution coming from the singlet, that we dub G S . The antisymmetrization operator implies the computation of where the equality follows from a property of the Pfaffian (see Eq. 51). The overall sign is arbitrary for a WF; thus the antisymmetrized product of singlet pairs (geminals) is indeed equivalent to the computation of the determinant of the matrix G S : It should be noticed that it is not necessary that the matrix G ud is symmetric to reduce the Pfaffian to a single determinant evaluation. As long as the matrices G uu and G dd are zero, the Pfaffian is indeed equivalent to det(G ud ) and describes an antisymmetric WF. However, if G ud is not symmetric the function is not an eigenstate of the spin. In other terms, there is a spin contamination, similarly to the case of unrestricted HF calculations.
The AGP ansatz can be generalized to describe polarized systems, i.e., systems where the number N u of electrons of spin up is different from the number N d of electrons with spin down. With no loss of generality, we can assume that N u > N d , thus the system is constituted by a number p = N d of electron pairs and a number k = N u − N d of unpaired electrons (clearly, N = N u + N d = 2p + k). We aim at evaluating: As above, we assume that the pairing function is zero for likespin pairs. With no loss of generality, we can assume that electrons i = 1, . . . , N u have spin up and electrons j = N u + 1, . . . , N have spin down. Moreover, as it was done at the end of the previous subsection, we can add k fictitious entries to g(i, j), such that g(i, N + l) = φ l (i) for l = 1, . . . , k and i = 1, . . . , N u . Thus, the antisymmetrization implies the use of Eq. 51, i.e., : with the N u ×N u matrixG = G ud Φ , the N u ×N d matrix G ud describing the pairing between the N u spin up electrons and the N d spin down electrons, and the N u ×k matrix Φ describing the k unpaired orbitals. Thus, we need to evaluate: One of the most important advantages of the AGP ansatz is that it is equivalent to a linear combination of Slater determinants (i.e., multi configurations), but the computational cost remains at the level of a single-determinant one [64][65][66] , see Section III F. This is especially important for large systems because the naive multireference approach requires an exponentially large number of Slater determinants, which drastically increases the computational cost. The AGP ansatz was applied to the ab initio calculation by M. Casula and S. Sorella 64 for the first time in 2003; then it has been also implemented in other QMC codes.

D. AGP and Pf with constrained number of molecular orbital (AGPn and Pfn)
A convenient way to impose constraints on the variational parameters defining the AGP or Pf WF is obtained by rewriting the expansion of the geminal in terms of molecular orbitals (MOs). As shown in Eq. 46, a geminal g(i, j) is natively expressed in terms of the atomic orbitals {ψ a,l (r)}, by summing over all the atoms a, the corresponding orbitals l and the spin index σ (as the elements of the basis may in principle also depend on the spin, even though the most common choice is to take the same orbital for each of the two possible values of the spins). In order to simplify the notation here, let us merge the atomic orbital and spin indices in a unique one that is indicated with a greek symbol (e.g., µ ↔ (a, l)) running from 1 to the total dimension 2L of the atomic orbitals used, L for each spin component. Therefore, and clearly the symmetry of g implies that Moreover, the overlaps S µ,ν ≡ ψ µ |ψ ν between atomic orbitals define the overlap matrix S , that in the case of a spin-dependent basis is block diagonal: with S uu and S dd positive definite L×L square matrices (S uu = S dd when orbitals are the same for spin up and spin down). In TurboRVB, the overlap matrix S is computed on a suitable uniform mesh with an efficient and general parallel algorithm.
Then an orthonormal basis is defined: where S −1/2 is well defined since S is strictly positive definite. 67 The matrix elements of g can be recasted in the new orthonormal basis ψ µ (i) , yielding At this point, from the spectral theory of skew-symmetric matrix, it is possible to perform the Youla decomposition ofÃ (see Ref. 68), which can be written in the formÃ = QΣQ T , where Q is unitary (also real ifÃ is real), and the matrix Σ is block diagonal with Σ 2k−1,2k = λ k = −Σ 2k,2k−1 for k = 1, ..., L, and zero everywhere else, with λ k ≥ 0. 69 So, the pairing function g(i, j) can be written as This defines a basis of L MOs {φ k (i)} and corresponding twinned ones φ k (i) , forming together a basis of 2L mutually orthonormal elements for which the original geminal function reads: with λ k ≥ 0. After these transformations these MOs can be finally written in the chosen (hybrid) atomic basis: by appropriate p × 2L rectangular matrices P andP. Then, with no loss of generality, we can assume that the molecular The above expression highlights that the most important MOs are those corresponding to the larger values of λ k . Therefore, it is possible to constrain the variational freedom by neglecting all the orbitals with k > n, yielding the pairing function: where n is conveniently chosen and is ≪ L. This yields the AGPn ansatz and the Φ AGPn WF, which can be useful to improve the stability of the wave-function optimization. The corresponding algorithm (enabled by setting molopt=-1 in the optimization section input of TurboRVB), based on projection operators in the space of the n molecular orbitals considered, has been described extensively in Ref. 15. Moreover, in the original paper 70 introducing the AGPn, a precise recipe was given to improve the evaluation of the binding energies. Indeed, despite a constraint on the variational parameters necessarily increases the variational energy expectation value, energy differences may actually improve by an appropriate choice of n. In the mentioned work 70 , this promising approach was applied with an AGP containing only singlet correlations, but the binding energies were defined without using a rigorous size consistent criterium. This drawback can be now removed, by exploiting the full variational freedom of the Pf WF combined with a general spin-dependent Jastrow factor (see e.g., Fig. 14). Work is in progress in this interesting research direction. The variational optimization of an AGP with a fixed number n of molecular orbitals can be easily generalized to the Pf case, by exploiting that the constrained Pf WF, dubbed Pfn, can be written either in the canonical form with MOs as in Eq. 63 or in the localized basis set expansion, as in Eq. 58, with a corre- to Eq. 63 an arbitrary small variation δg n of the constrained pairing function g n reads: and therefore satisfies the following property, as it will be shown later: whereÎ is the identity operator, L and R are projection operators over the occupied MOs, i.e., L 2 (i, j) = dkL(i, k)L(k, j) = L(i, j), and similarly R 2 = R, where here and henceforth the shorthand integration symbol dk = σ k dr k contains implicitly also the spin summation. These operators are then defined as follows: With the above definitions, Eq. 64 is easily verified because each term of Eq. 64 is annihilated either by the left (Î − L) or the right (Î − R) projection over the unoccupied MOs. Notice that L = R in the real case and L = R * in the most general complex case. In this way, in order to implement a constrained variation δg n of the Pfn WF, corresponding to an appropriate variation of its matrix δA n µ,ν , it is useful to work with a small free variation δg (with corresponding δA µ,ν ). This is then projected onto the chosen restricted ansatz by means of the following equation: Indeed, it is easy to show that the RHS of the above equation vanishes if we applyÎ − L andÎ − R to its left and its right, respectively, just becauseÎ − R andÎ − L are projection operators, being such R and L, yielding (Î − R) 2 = (Î − R) and (Î − L) 2 = (Î − L), from which Eq. 67 fulfills Eq. 65. Eq. 67 represents, therefore, a linear relation applied to the variational parameter matrix change δA µ,ν corresponding to the unconstrained geminal g in Eq. 58, yielding the new constrained variation δA n µ,ν . Indeed, by using the definitions of the projector operators in Eq. 66 and the expansion of the MOs in the atomic (hybrid) basis (see Eq. 62), Eq. 67 turns to a number of matrix-matrix operations acting on δA, P,P and the overlap matrix S that can be easily and efficiently implemented 15 .
This linear relation between A and A n can be therefore easily implemented together with the corresponding derivatives necessary to the optimization of the energy 71 and allows the explicit calculation of the new matrix A n µ,ν , yielding the new constrained geminal g n +δg n . Then the new geminal can be recasted in the form of Eq. 63 by the mentioned diagonalization of skew-symmetric matrices, in this way implicitly neglecting nonlinear contributions that are irrelevant close to convergence, when δg n → 0. After employing several iterations of this type, the lowest energy ansatz of the JPfn type can be obtained in a relatively simple and very efficient way.
Finally, it is also important to emphasize that this constrained optimization algorithm allows the reduction of the number of parameters, by efficiently exploiting locality, namely that variational parameters A µ,ν corresponding to atoms at a distance larger than a reasonable cutoff (an input named rmax in TurboRVB) can be safely disregarded with negligible error 15 .

E. Single Slater determinant (SD)
An important special case of the AGPn and Pfn ansätze discussed in the previous section is when we constrain the pairing function g n to use the minimum possible number of molecular orbitals providing a non zero WF. The minimum number for an unpolarized system with N electrons is equal to the number of electron pairs, that is n = N/2. The WF obtained in this way starting from the AGP is indeed the single Slater determinant (SD) 15,70 , and we dub it Φ SD . In principle also the Pfn WF with n = N/2 corresponds to a single Slater determinant with spin dependent molecular orbitals, a case that has never been considered so far, but it represents an available option within the most recent versions of TurboRVB.
Similarly, in a polarized system having N u spin up electrons and N d spin down electrons (we assume N u > N d ), the SD ansatz is obtained by using N u − N d unpaired electrons, and using a constrained pairing function g n with n = N d /2.

F. Implicit multiconfigurational character of the AGP
In the previous sections, it has been mentioned several times that the Pfaffian and the AGP ansätze have a multiconfigurational character, despite they can be evaluated at the cost of a single determinant. In this section, we expand the AGP ansatz in terms of Slater determinants, to show this aspect explicitly. In order to simplify the derivation, we will consider here a simplified case, while the most general case could be studied with a similar approach but involving more cumbersome expressions.
We consider the AGPs ansatz (Eq. 53) for an unpolarized system of N = 2N p electrons (i.e., N u = N d = N/2). The symmetry implies that the twin molecular orbitals φ k (i),φ k (i) appearing in the RHS of Eq. 61 have the same spatial part, that we denoteφ k (r), modulus a sign (because they are orthonormal), and they have opposite spin part. Without loss of generality we can assume that φ k is spin up andφ k is spin down. Given this convention, the scalar product l k between the spatial parts of φ k andφ k with either be +1 or -1. We defineλ k = l k λ k , where λ k is the same as the one in the RHS of Eq. 61. Notice that λ k = λ k , so λ are ranked in decreasing order of their absolute value (i.e., λ k ≥ λ k+1 ). The pairing function then can be written as This expression is useful for comparing with the standard CI expansion.
It is convenient now to use the second quantization notations in order to simplify the derivation. In particular, we indicate withâ † k,↑ (â † k,↓ ) the operator that creates an electron of spin up (down) in the orbitalsφ k , and satisfies the canonical anticommutation relations. We can write the pairing functions g(i, j) ≡ i, j|ĝ|0 , where |0 is the vacuum, |i, j is the WF of a system with one electron with coordinates (r i , σ i ) and another with coordinates (r j , σ j ), andĝ is the operator: whereλ k is the one defined above. Using this notation, it is easy to show (see Appendix of Ref. 66) that the pairing function g is equivalent to the complete active space of 2 electrons on the L molecular orbitals φ k .
Within this notation, the AGPs WF is where p ≡ N/2 is the number of electron pairs. If we substitute inĝ p the expansion forĝ in Eq. 69, after having conveniently defined the operatorb k ≡â † k,↑â † k,↓ which created an electron pair on the orbitalφ k , and having noticed that b kbl =b lbk andb 2 k = 0 (as following from the anticommutation relations of theâ † k operators) we obtain that: The ranking of theλ k coefficients implies that the leading term inĝ p is given by the term with (i 1 , . . . , i p ) = (1, . . . , p). 72 Thus, we can rewriteĝ p in terms of this leading term and excitations with respect to it: In the above expansion, we shall recognize that the first element is a single Slater determinant of the molecular orbitals φ i i=1,..., N/2 ; the second element is the summation of all the possible double excitations both going from an orbital j ≤ N/2 to an orbital q > N/2; the third element is the summation of a subset of the possible quadruple excitations, and so on. In other terms, AGPs can be written as the zero seniority 73 subset of the CI expansion, having some constraints on the coefficients of the expansion.
It is worth mentioning that AGPs is also reliable in cases where a single Slater determinant is not a good reference, as for instance, in the case of a broken covalent bond, where the HOMO and LUMO orbitals are degenerate (see Fig. 5). In fact, we can have that λ p = λ p+1 , so the two determinants p i=1b i |0 and p−1 i=1b i b p+1 |0 have the same weight and they both are the leading terms of the constrained zero-seniority expansion in Eq. 72.
Notice that in the AGPn ansatz we can perform an expansion similar to Eq. 72, but the excitations stops at the orbital n rather than L. Thus, it is straightforward to see that for n = N/2 there are no excitations and the only term is a single Slater determinant.

G. Atomic basis set for the pairing function and the Jastrow factor
TurboRVB employs localized atomic orbitals such as the Gaussian type: or the Slater type: where the real and the imaginary part (m > 0) of the spherical harmonic function Y l,m,I (θ, ϕ) centered at R I are (is) taken and rewritten in Cartesian coordinates in order to work with real defined and easy to compute orbitals, l is the corresponding angular momentum and m ≥ 0 is its projection number along the z−quantization axis. The localized atomic orbitals are also present in the onebody and three/four body Jastrow parts (i.e., denoted as χ (r) in Eq. 30 and Eq. 34). One can use standard basis sets with exponents ζ (also coefficients for a contracted basis set) taken from the available standard database such as the Basis Set Exchange 74 or from other more specific references when using pseudo potential [75][76][77][78][79][80] . A python wrapper, named Turbo-Genius, makes the above procedure much easier, as shown later.

H. Pseudopotential
TurboRVB supports pseudopotential calculations both in VMC and LRDMC calculations. Many ECPs have been generated and successfully used in quantum chemistry codes, but they are usually tuned to match Density Functional Theory (DFT) or Hartree-Fock (HF) all-electron (AE) calculations, which are not expected to be optimal for state of the art manybody techniques. Recently some progress has been made in this direction and pseudopotentials determined by correlated many-body techniques are also available [75][76][77][78][79][80][81][82] . All the pseudopotentials used in QMC employ the standard semi-local form: where r i,I = |r i − R I | is the distance between the i-th electron and the I-th ions, l max is the maximum angular momentum of the ion I, and l max l=0 l m=−l Y l,m Y l,m is a projection operator on the spherical harmonics centered at the ion I. In TurboRVB, the angular momentum projector is calculated by using standard polyhedral quadrature formulas for the angular integrations 83 . As it is now becoming a common practice not only in QMC, both the local V I loc r i,I and the non-local V I l r i,I functions, are expanded over a simple Gaussian basis parametrized by coefficients (e.g., effective charge Z eff and other simple constants), multiplying simple powers of r, and a corresponding gaussian term: where α k,l , β k,l (usually small positive integers), and γ k,l are the parameters obtained by appropriate fitting. Several published pseudopotentials have already been tabulated in Tur-boRVB. Of course, one can also use any pseudopotential employing the semi-local form in the mentioned Gaussian basis with a straightforward little extra work for the input preparation.

I. Contraction of the primitive atomic basis
The contraction of atomic orbitals has been widely used in quantum chemistry and DFT codes, and was originally introduced to define a pseudo-Slater orbital by combining several primitive Gaussian orbitals. This is also important in the QMC context because it decreases the number of variational parameters by a large factor, as it will be shown in this section. Although one can obtain a contracted basis set directly from a database, TurboRVB allows us to prepare a high-quality hybrid basis set (contraction) starting from a given primitive basis, by the so-called "geminal embedding scheme". 84 To this purpose, let us decompose the genimal function in terms of atomic contribution, as far as the dependence over i (the left argument) is concerned: where I represents an atom in a system, U I proj (i, j) is the pairing function projected on the atom I, A I µ,ν is A µ,ν if µ ∈ I, otherwise A I µ,ν = 0, where A µ,ν is assumed to be given for the system under consideration, e.g., obtained by a standard DFT calculation, where in this case, by Eq. 63: where c k µ andc k µ are the coefficients of the DFT molecular or- Quite generally, the projected pairing function can be expanded in a truncated space spanned by q terms: In other words, the Schmidt decomposition is applied to the matrixŨ I proj (i, j) describing the coupling between a given atom I and the enviroment, within the geminal ansatz. This procedure defines the so-called geminal embedded orbitals (GEOs), that are determined in terms of an expansion over all the atomic orbitals used for the atom I: where ψ GEO k (i) are orthonormal. Following the Schmidt decomposition, it is possible to determine the best GEOs by minimizing the Euclidian distance between the original and the truncated geminal functions: Considering all possible unconstrained functions ψ Env k (j) and employing the steady condition δd 2 δψ Env k (j) = 0, d 2 reads: where U I is the density matrix that reads: Since the GEOs ψ GEO k (r) are orthonormal and the density matrix is hermitian, Eq. 82 becomes minimum when the GEOs coincides with the p eigenvectors of the density matrix with the maximum eigenvalues (denoted as w i ). The original atomic basis ψ µ (i) is usually non-orthogonal, so the problem turns into the generalized eigenvalue equation: The truncation error is readily estimated by the summation of the eigenvalues d 2 = U I proj 2 − p i=1 w p . Since the eigenvalues are sorted in ascending order, a suitably chosen value of p allows the user to neglect the most irrelevant vector components with small eigenvalues w i and work with enough accuracy even with a few GEOs per atom, thus minimizing the number of variational parameters necessary to describe well the system, as it will be shown below. With this construction a new geminal is defined in the GEO basis, namely: where the matrix coefficientÃ µ,ν are given by maximizing the normalized overlap (Q) between the original and the new geminals: where g|g = didjg (i, j) * g (i, j). It turns out that the overlap remains large even for small GEO basis set size p; implying that, by using this scheme, one can decrease the number of variational parameters corresponding to the matrix A, i.e., from 4L 2 to 4p 2 ≪ 4L 2 .

J. Conversion of the WF
As described in Sec. III, TurboRVB implements different types of ansatz: ( i ) the Pfaffian (Pf), (ii) the Pfaffian with constrained number of molecular orbitals (Pfn) (iii) the Antisymmetrized Geminal Power (AGP), (iv) the Antisymmetrized Geminal Power with constrained number of molecular orbitals (AGPn), and (v) the single Slater determinant (SD). One can choose a proper ansatz depending on a target system, considering the computational cost of a chosen ansatz and the relevant physical and chemical properties of a target material. During the simulation, a user can go back and forth between the ansatz using modules implemented in TurboRVB, with/without losing the information of an optimized ansatz (Fig. 4): The first case is to add molecular orbitals to an ansatz, i.e., JAGP ⇒ JSD, JAGP ⇒ JAGPn, or JPf ⇒ JPfn. In TurboRVB, this is obtained by rewriting the expansion of the geminal in terms of molecular orbitals (see sec. III D and sec. III E). The corresponding tool is convert-fort10mol.x. The second important case is to convert an ansatz among the available ones, i.e., JSD, JAGP, or JAGPn ⇒ JAGP. In TurboRVB, this is the purpose of the convertfort10.x tool and is achieved by maximizing the overlap between the two WFs, one of them being the input (fort.10 in) and the other being the type (fort.10 out) to be filled by new geminal matrix coefficients (result written in fort.10 new). In more details, in TurboRVB, the following overlap between two geminals is maximized: in order to obtain new geminal matrix coefficients A new µ,ν , defining the new pairing function as: while the original geminal was given in terms of the parameter matrix A ori µ,ν : Notice that 0 Q 1; therefore the larger is Q, the better is the conversion, and Q approaches the unit value if the conversion is perfect. For this type of conversion, one can also apply the geminal embedding scheme to construct a hybrid basis set, as described in the sec.III I. The final case is to convert a JAGP ansatz to JPf. Since the JAGP ansatz is a special case of the JPf one, where only G ud and G du terms are defined as described in the section III C, the conversion can be realized just by direct substitution. Therefore, the main challenge is to find a reasonable initialization for the two spin-triplet sectors G uu and G dd that are not described in the JAGP and that otherwise have to be set to 0. There are two possible approaches 61,68 : ( i ) for polarized systems, we can build the G uu block of the matrix by using an even number of unpaired orbitals {φ i } and build an antisymmetric g uu by means of Eq. 63, where the eigenvalues λ k are chosen to be large enough to occupy certainly these unpaired states, as in the standard Slater determinant used for our initialization. Again, we emphasize that this works only for polarized systems. (ii) The second approach that also works in a spin-unpolarized case is to determine a standard broken symmetry single determinant ansatz (e.g., by the TurboRVB built-in DFT within the LSDA) and modify it with a global spin rotation. Indeed, in the presence of finite local magnetic moments, it is often convenient to rotate the spin moments of the WF in a direction perpendicular to the spin quantization axis chosen for our spin-dependent Jastrow factor, i.e., the z quantization axis. In this way one can obtain reasonable initializations for G uu and G dd . TurboRVB allows every possible rotation, including an arbitrary small one close to the identity. A particularly important case is when a rotation of π/2 is applied around the y direction. This operation maps One can convert from a AGP the pairing function that is obtained from a VMC optimization: to a Pf one: This transformation provides a meaningful initialization to the Pfaffian WF that can be then optimized for reaching the best possible description of the ground state within this ansatz.

IV. BULK SYSTEMS
The application of TurboRVB is not limited to open systems such as atoms and molecules. TurboRVB can also simulate bulk systems in large supercells with arbitrary twisted boundary conditions. These are used to minimize finite-size effects, and represent quite an important approach 85,86 in order to reach a meaningful and accurate thermodynamic limit.

A. CRYSTAL basis set
For periodic system calculations, the many-body WF should satisfy the many-body Bloch condition 87,88 : which follows from the property that the many-body Hamiltonian is invariant under the translation of any electron coordinate by a simulation-cell vector T s , where T s = la + mb + nc is determined by arbitrary integers l, n, m and the three vectors a, b and c define the supercell.
In TurboRVB, a single-particle basis set satisfies the following condition: where k s is a twist vector (k s = k x s , k y s , k z s ), and T s represents an arbitrary simulation cell vector. Notice that the use of a non vanishing twist vector generally makes a many-body WF complex. TurboRVB implements the CRYSTAL periodic basis 15,89,90 : where ψ l,m,I is a non-periodic real atomic orbital such as Gaussian (Eq. 73) and Slater (Eq. 74). The use of Gaussian or Slater orbitals that rapidly decay far from nuclei guarantees that the above summation converges fast with a finite small number of T s . Notice that, in TurboRVB, T s are not limited to orthorhombic ones (e.g., rhombohedral, hexagonal) but has been recently extended to include all possible crystal translation groups (e.g., triclinic). The same procedure is applied to the basis set for the Jastrow part, although using simple periodic boundary conditions 15 , because the twists do not affect the Jastrow part of the WF, namely: which satisfies χ PBC l,m,I (r + T s ; ζ) = χ PBC l,m,I (r; ζ). Moreover, the homogeneous one-body and two-body Jastrow factors have to be appropriately periodized because they are not defined in terms of the above periodic basis. Namely, the homogeneous one-body Jastrow part (Eq. 30) should satisfy:J and the two-body Jastrow part (Eq. 32) should fulfill: In order to satisfy the above constraints, we consider the relative electron-nuclei or electron-electron coordinate differences r d , necessary to evaluate J 1 and J 2 , respectively, and expand them as: where r a , r b and r c are appropriate transformed coordinates, that are conveniently defined within a cube of unit length, because of the assumed periodicity of the supercell, namely |r a |, |r b |, |r c | < 1/2. As a consequence, this mapping makes the physical electron-electron and electron-ion distance periodic by definition (i.e., they refer to the minimum distance image of the supercell). However, there may be divergences or singularity at the boundaries of this unit cube. Therefore, before computing the distance corresponding to r d , these coordinates are transformed (r a , r b , r c ) → (r a ,r b ,r c ) = (p(r a ), p(r b ), p(r c )) by use of an appropriate function p(x), with at least continuous first derivative for |x| < 1/2. This function is chosen to preserve the physical meaning at short distances, i.e., p(x) = x in these cases, and being nonlinear elsewhere, in order to satisfy not only the periodicity but also the requirement of continuous first derivatives of the many-body WF Ψ k s . We have, therefore, defined p(x) as follows: Indeed, though the modified relative distance diverges (i.e., |r d | → ∞) at the edges of the Wigner-Seitz cell (e.g., r = ± 1 2 a, ± 1 2 b, ± 1 2 c), the exponential (Eq. 31) and the Padé (Eq. 33) functions remain finite, u (r) → 1/2b ea and v σ i ,σ j r i, j → 1/4b para ee or 1/2b anti ee , respectively, thus preserving with continuity the periodicity of the one-/two-body homogeneous parts of the Jastrow factor. However, for the Padé form, one has to change the expression for p(x) in order to satisfy also the continuity in the WF derivatives, even when the modified relative distances diverge, i.e.,: In this case the region where p(x) = x further shrinks. Thus, it is often more convenient to use the exponential form to obtain a more accurate variational WF, because the long-range part is implicitly corrected by the inhomogeneous terms in Eq. 34. Finally, we remark that the many-body WF also obeys the second Bloch condition 87,88 , namely: where T p represents a unit-cell (not supercell) vector, and k p is the crystal momentum. This comes from the property that the many-body Hamiltonian is invariant under the simultaneous translation of all-electron coordinates by a unit-cell vector T p . Within TurboRVB, this condition can be employed by imposing the intra-unit cell translational symmetries on the Jastrow and the pairing function, as simple linear constraints in the variational parameters. However, this option is restricted to the case k p = 0. On the other hand, it is possible to use different twists on each spin component, that has proven very effective for implementing the mentioned translation symmetries within pairing WFs. 91

B. Finite-size effects
The systematic error induced by a finite simulation cell is a long-standing issue in the ab initio QMC calculation. There are two types of finite-size errors in QMC calculations; one is the so-called one-body effect that arises from the kinetic energy term of the Hamiltonian, and the other one is the so-called two-body effect that arises from the periodic Ewald contribution resulting from the electron-electron interaction. Notice that, in the independent-particle calculations (i.e., DFT), only the former is present, which can be readily evaluated by k summation in the first Brillouin zone, but the two-body finite-size effects are not present because the exchange-correlation energy used in DFT usually derives from DMC results extrapolated to the infinite simulation cell size. To correct the one-body finite-size error, one can use the twisted averaged boundary conditions 85 , special k-points methods 87,88 , or the exact special twist (EST) method 86 , all of which have been implemented in TurboRVB. Within the Tur-boRVB implementation, the Jastrow part is independent of twists (i.e., TurboRVB uses a common Jastrow for all twists), As emphasized in Ref. 85, at variance with DFT, the QMC computational effort is independent of the number of k-points used. For the two-body finite-size effects, which cannot be alleviated by the above remedy, one can employ the model periodic Coulomb interaction (MPC) [92][93][94] . This method has not been implemented in TurboRVB yet. Nevertheless, simpler alternatives exist. For instance one can alleviate the twobody finite size effects by directly increasing the supercell size, or one can estimate these effects by employing the KZK exchange-correlation function 95 at the DFT level. Moreover, it is also possible to employ systematic finite-size corrections based on the knowledge of the density structure factor in momentum space 96 , that can be readily computed within Tur-boRVB, with a short postprocessing computation.

V. BUILT-IN DENSITY FUNCTIONAL THEORY (DFT) CODE
Although most QMC codes load their trial WFs from available DFT/quantum chemistry codes such as Gaussian 97 , CRYSTAL 90 , and Quantum Espresso 98 , TurboRVB is a selfconsistent many-body package and does not require input from other codes. Indeed, in TurboRVB, an original DFT code is implemented that has the important feature to work in exactly the same basis and pseudopotentials (if any) used for the QMC many-body WF, and for instance can work with mixed Slater-Gaussian basis or with arbitrary contraction, i.e., mixing angular momenta. On the other hand, the drawback is that only very few functionals are available so far because, within the TurboRVB approach, the DFT has to be used only to have a reasonable initialization of the antisymmetric part of the many-body WF. This is then optimized at best by direct minimization of the energy in the presence of the Jastrow factor and should very weakly depend on the initialization, or for instance, on the DFT functional used. At present, TurboRVB supports Perdew-Zunger (PZ81) Local-Spin Density Approximation (LSDA) 99 , and those combined with KZK exchangecorrelation energy 95 , but, for instance, does not allow LDA + U calculations. This will be the goal of a possible future work.

A. Stabilization algorithm for large localized basis sets
An important issue in any DFT algorithm based on localized basis sets, such as Gaussian, is how to reach the so-called complete basis set limit. In the plane-wave method, this is systematically achieved by increasing the number of plane waves (e.g., kinetic cutoff energy). However, with a localized basis set, by increasing its dimension, numerical instabilities may occur as a consequence of the redundancy and non-orthonormality of the basis. On the other hand, the use of a small basis set implies too much biased results. In Tur-boRVB, an original and systematic algorithm to remove this redundancy 100 has been implemented. When the basis set is non-orthogonal, the key self consistent step for solving the DFT Kohn-Sham equations, turns to a generalized eigenvalue problem: where H is the single-particle Hamiltonian matrix, c is the vector of the coefficients of the orbitals, E is the single-particle energy, and the S is the N b × N b overlap matrix, where N b is the dimension of the basis set: The overlap matrix is strictly positive definite, namely all its eigenvalues are positive. The problem mentioned above shows up when the overlap matrix S becomes ill-conditioned. If the condition number, namely the ratio between the largest (denoted as s N b ) and the smallest eigenvalues (denoted as s 1 ), becomes very large, the CBS limit cannot be achieved by the standard procedure due to numerical instabilities. Therefore, in TurboRVB, the small eigenvalues and the corresponding eigenvectors are disregarded. Indeed, the original basis is modified as follows: where s i and υ i j are the i-th eigenvalues and eigenvectors of the overlap matrix, ε mach is usually set to the machine precision, j runs from 1 to N b , i runs from 1 to M b , and M b represents the number of eigenvalues satisfying the above inequality. If infinite numerical accuracy were available, the matrix S in the new basis would be just the identity and the generalized eigenvalue problem in Eq. 104 would turn to a standard one in this basis. Within finite precision arithmetic it is better to iterate further the stabilization procedure and define a new basis: wheres i andυ i j are the i-th eigenvalues and eigenvectors of the recomputed overlap matrixs i, j = e i | S |e j (now well conditioned since very close to the identity), respectively. For using the latter basis set, the N b × M b global transformation matrix is stored, as it takes into account the projection from the original basis to the final modified one: On this basis set, the generalized eigenvalue problem becomes a well-conditioned one with a corresponding M b × M b Hamiltonian matrixH = UHŨ. In this way, one can avoid the numerical instabilities induced by a too redundant large basis set. This stabilization introduces an error at most √ ε mach · N b , which is typically negligible compared with the target chemical accuracy.

B. Electron-ion cusp condition
Another feature of the DFT code implemented in Tur-boRVB, is that electron-nuclei cusp conditions are exactly fulfilled for any basis (i.e., Gaussian orbital) even within the DFT framework. This is achieved by an appropriate modification of the standard basis sets commonly used (e.g., ccpVTZ) 101 for WF based calculations: the new basis is obtained by multiplying each element of the original basis by a suitably chosen one-body Jastrow factor introducing the correct cusps, whereJ 1 (r) is the same as in Eq. 29 and the parameter b in Eq. 31 is optimized by direct minimization of the chosen DFT energy functional. In this way each element of the reshaped basis set satisfies the so-called electron-ion cusp conditions, namely that, when r is close to any atomic position R b , the following relation holds: for all a, b. This formulation allows us to reach the CBS limit extremely fast in DFT calculations, especially within all electrons. Indeed, we have experienced that a given target accuracy can be obtained with a smaller basis, e.g., our modified ccpVDZ basis is typically equivalent in accuracy to the much larger ccpVTZ (and the ccpVTZ equivalent to the ccpVQZ, and so on and so forth). For QMC application, this tool is particularly useful because the Slater determinant obtained with this new basis does not have divergences in the local energy, that are instead present in the original basis, when electron positions are close to the nuclear ones. This feature emphasizes further the clear advantage of TurboRVB to allow the use of a special purpose DFT code, just devoted to the optimal initialization of a QMC many-body electronic WF.

C. Double-grid remedy for all-electron calculations
As described in the previous section, TurboRVB employs standard atomic orbitals, modified by means of an appropri-ately chosen one-body factor. Thus the electron-electron integrations cannot be evaluated analytically even when the Gaussian atomic basis set is employed. Therefore, TurboRVB calculates the electron-electron Hartree potential by solving the Poisson's equation with the fast Fourier transform (FFT) on a Cartesian grid both in open and periodic systems. The numerical solution becomes problematic in all-electron calculations when the atomic number Z becomes large. Indeed, in this case, the Kohn-Sham Hamiltonian matrix elements have to be computed in the presence of a rapidly varying electron density in the vicinity of nuclei, and therefore require an exceedingly large FFT mesh, with an almost prohibitive computational cost for its evaluation. TurboRVB alleviates this drawback by a double-grid scheme 102 , in which the Hartree potential is calculated with standard FFT convolution on a coarse mesh. Then, these values are interpolated on a much finer mesh in the vicinity of the nuclei. Although the DFT energy obtained with the above approximation is not exactly consistent with the one corresponding to a very dense uniform mesh, the corresponding VMC energies and the variances are almost indistinguishable each other, at a much cheaper computational cost.

VI. DERIVATIVES OF ENERGY
Derivatives of total energies with respect to variational parameters represent an essential ingredient for optimizing a many-body WF. Forces (derivative with respect to atomic positions) are also essential for performing structural optimization or molecular dynamics. However, in a complex code, and especially in QMC, the evaluation of the functional derivatives, necessary for the WF optimization, are very complicated, mainly for the complexity of the algorithm that, in turn, may lead to a very inefficient implementation, though recent progress has been done 103 . For instance , a simple approach is to compute them with finite difference expressions, leading to a too large computational time, because obviously proportional to the number of targeted derivatives.
Algorithmic differentiation (AD) is a method capable of solving all the above problems, essentially by a smart application of the differentiation chain rule. There are two types of algorithms, the forward algorithmic differentiation (FAD) that implements the chain rule straightforwardly from the beginning to the end of the algorithm and, the adjoint algorithmic differentiation (AAD), that uses the chain rule starting from the end of the algorithm (also known as backward propagation). When the number of input parameters is much larger than the corresponding output ones, AAD is much more efficient than FAD and indeed allows the calculation of all possible derivatives in a computational time proportional (with a small prefactor, see e.g., Fig. 6) to the one for computing the target function (i.e., the energy or the WF value). The former is therefore the ideal method for Quantum Monte Carlo when the variational WF contains several variational parameters. The AAD was applied for the calculation of atomic forces in Ref. 39. To our knowledge, this was the first time that AAD was used within QMC. Now, TurboRVB implements all derivatives such as those of the local energy ( ∂ ∂α k e L (x i )) or those corresponding to the WF logarithm ( ∂ ∂α k ln Ψ (x i )) using the AAD, which drastically improves the efficiency (Fig. 6) 39 and reliability of the calculation.

A. Derivatives with respect to variational parameters
The derivatives of the energy with respect to a given real variational parameter α k (one complex parameter can be thought to be composed by two real ones, its real and imaginary part, respectively) is represented as a generalized force: In variational Monte Carlo, the derivative can be evaluated using M configurations of electron coordinates 15 : where e L (x) is the local energy, O k (x) is the logarithmic derivative of the WF (i.e., O k (x) = ∂ ln Ψ α (x) ∂α k ), andŌ k is its aver- . In TurboRVB, the logarithmic derivative (O k (x)) is computed very efficiently by using the AAD algoritm. Notice that the derivatives of the local energy are not needed here because the Hamiltonian does not depend on any variational parameter. Instead, these terms are necessary in order to calculate ionic forces (i.e., derivatives of the total energies with respect to atomic postions). If the WF is an exact eigenstate of the hamiltonian, the generalized forces f k exactly vanish without statistical errors because the local energy is no longer dependent on x. In other words, the derivatives have the zero-variance property and represent therefore the fundamental ingredients for an efficient WF optimization, as described sec. VII.
In practice, during an optimization, the code monitors the variational energy (E (α)) and the maximum value of the signal to noise ratio among all the force components, which is denoted as devmax in the code: where σ f k represents the estimated error bar of the force f k . This value is used in TurboRVB as one of the convergence criteria of the optimization. To our experience, after a reasonable number of iterations devmax stabilizes to a value ≤ 4.0, when a correct (possibly local) minumum is being approached.

B. Atomic forces
In TurboRVB, a very useful space warp coordinate transformation (SWCT) 39 is employed for the calculation of atomic forces. This scheme was introduced several years ago in Ref. 104, in order to decrease the statistical errors of forces. The finite-difference expression for ionic forces calculations is: The SWCT is used to mimic the displacement of charges around the nucleus,r where κ (r) is a function that decays sufficiently fast for large r because the charges far from the nuclei should not be affected by the SWCT, and ω → 0, as a consequence of this requirement. It turns out that the original 104 simple choice κ (r) = 1/r 4 works very well and is indeed adopted in Tur-boRVB.
Starting from the finite-difference expression, one can straightforwardly derive the corresponding differential expression 39 : where J is the Jacobian of the above transformation, and the brackets indicates a Monte Carlo average over the trial WF. All the terms above can be written by the partial derivatives of the local energy and those of the logarithm of the WF 39 : (119) In order to evaluate these differential expressions, 6N + 6N at components have to be evaluated, namely, ∂ ∂r i e L , ∂ ∂r i ln Ψ , ∂ ∂R a ln e L , ∂ ∂R a ln Ψ 39 . These values are very efficiently computed in TurboRVB, by using the aforementioned AAD, that works even in presence of pseudopotentials.
The SWCT significantly decreases the statistical errors, but the forces still have infinite variance properties because ∂e L and ∂ ln (Ψ) diverge in the vicinity of the nodal surfaces. Attaccalite and Sorella 105 developed a reweighting method to address this issue in an unbiased way. Within this scheme, the QMC sampling is not driven by the chosen trial WF Π T (x) = Ψ 2 T (x), but by a slightly different guiding function Ψ G (x) defined by where R (x) vanishes more weakly than the trial function does, when the configuration x approaches the nodal surface, whereas R ε is an appropriately chosen regularization of R, depending on ε, that never vanishes. In this way, the guiding function is larger when approaching the nodal surface, and this singular region can be more accurately sampled in order to avoid infinite variance problems in the calculation of forces, as it will be shown later. For this purpose, R (x) is defined in a way to vanish as the antisymmetric part of the WF to some power 2θ R ≤ 1. However, in order to avoid too large fluctuations, a simple relation is used that Pf(G) ≃ 1 the pairing function (Eq.45), here referred to the general AGP case (the other cases straightforwardly follows upon some restriction of the matrix G). Indeed, if the Pfaffian goes to zero, most of the matrix elements of the inverse of G diverge inversely proportional to it. Thus the quantity satisfies the requirement without depending explicitly on the full Pfaffian that has fluctuations exponentially large in the number of particles and would lead to an inefficient regularization 105 . Here, we have introduced also the scaling factor S , that takes into account that the WF may vanish also when a single particle is going to infinity, and not because it is approaching the nodal surface. In this case, the scaling factor has to vanish in a proper way so as to allow a non-vanishing R. In this way, the guiding function can decay sufficiently fast in this limit and can be normalized in an open system. This represents a necessary condition to have a stable simulation, otherwise all electrons will be kicked out at infinity, after a long Markov chain. Therefore, the factor S is chosen to vanish as long as one raw of the pairing function is vanishing, because it is corresponding to an electron going to infinity (all matrix elements containing such electron coordinate have to vanish): This definition is also useful, because the regularization proposed is scale-invariant, namely R (x) remains unchanged if G is scaled by an arbitrary constant and therefore is adopted as it is also for bulk systems, and it is particular important when, in these cases, there exist large regions of almost negligible electronic density. 106 After that R ε (x) is defined as: that is much simpler than the original proposal reported in Ref. 105. By using the new probability, forces (the HF and the Pulay) can be evaluated with finite variance as: where the W (x) is the new weight: (125) and F (x), according to Eq. 117, indicates the average of an appropriate random variable containing either the derivative of the local energy (the HF contribution) or the product of the local energy and the log derivative of the WF (the Pulay one), both contributions diverging as ≃ 1 δ 2 , where δ is the distance from the nodal surface. The point is that, since in the vicinity of the nodal surface W (x) is proportional to δ 4θ R , whereas the probability vanishes much slower as Π G (x) ≃ δ 2−4θ R , the reweighting scheme solves the infinite variance problem for 1/4 < θ R ≤ 1/2 because the random variable W (x)F (x) remains with finite variance. This follows after simple inspection of its integrability in the small δ region, i.e., dδδ 4θ R −2 < ∞. Thus, both the HF and the Pulay terms can be evaluated with finite variances in the TurboRVB evaluation of atomic forces. In the most recent versions of TurboRVB, the default value of θ R is taken to be at the middle of the interval of stability, i.e., θ R = 3/8, because the value θ = 1/2, previously chosen 105 , leads to some instability, namely, during the QMC sampling, the nodal surface is approached too much closely, and the determinants or Pfaffians become very ill-conditioned with obvious problems of numerical accuracy. Moreover, the value of ε is automatically selected during the first steps of the simulation in a way that the average reweighting factor W (x) ≃ 0.8, that represents empirically an almost optimal setting.
Finally, we remark that this method to compute the energy derivatives with finite variance is applied by default in Tur-boRVB even for the optimization of the WF variational parameters, and that a similar approach can be used when computing atomic forces within the LRDMC algorithm.

A. Stochastic reconfiguration
Once the energy derivatives can be computed, the most straightforward strategy to optimize a WF is to employ the steepest descent method, where the WF parameters are iteratively updates as follows: where ∆ is a small constant and f k ≡ − ∂E ∂α k is the generalized force already defined in Eq. 111. However, it does not work well when optimizing highly non-linear WF parameters because a small change of a given variational parameter may produce a very different WF, whereas another parameter change may weakly affect the WF. Of course, one can use more sophisticated methods such as the Newton-Raphson method, the conjugate gradient, the quasi-Newton method, but the straightforward implementation of these optimizations do not work efficiently within a stochastic approach like QMC. In order to overcome this difficulty, a more efficient change in the variational parameters has been defined by means of a positive-definite preconditioning matrix S and the generalized force vector f: where the matrix S is stochastically evaluated by means of M configuration samples x = {x 1 , x 2 , . . . x M }: where The resulting approach is the so-called stochastic reconfiguration (SR) method 40 . Mazzola et al. 107 pointed out that the matrix S is essentially a metric for the parameter space, measuring the distance of the underlying normalized WF. Therefore, Eq. 128 is simply the steepest descent in this curved manifold. This observation connects the SR method with the so-called natural gradient method, widely used in the context of deep learning 108 . In this context, for each parameter α, Ô k −Ō k and S k,k ′ can be interpreted as the score function (i.e., the gradient of the log-likelihood function) and the Fisher information matrix (FIM), respectively, while the WF square |Ψ (x)| 2 plays the role of the likelihood function. In this sense, the stochastic reconfiguration method is essentially identical to the natural gradient optimization with the FIM that has been intensively used in the machine-learning community.
The straightforward implementation of the SR method is not stable mainly because the statistical noise sometimes makes the matrix S ill-conditioned, which deteriorates the efficiency of the optimization method 41 . Therefore, in practice, the diagonal elements of the preconditioning matrix S are shifted by a small positive parameter (ε) as: This modification improves the efficiency of the optimization by several orders of magnitude, as shown in Fig. 7. Finally, the variational parameters are updated as: In practice, the user should provide two input parameters, namely, ∆ (denoted as tpar in the code input variables) and ε (denoted as parr).

B. Linear method
In TurboRVB, the state-of-the-art QMC optimization method is also implemented, namely the so-called linear method [42][43][44] . In this scheme, a many-body WF is expanded up the linear order (i.e., considering only the first derivatives for the variational parameters): where α k is the k-th variational parameter, and Ψ k (x) is the first derivative with respect to the k-th variational parameter. Within this expansion, the expectation value of the energy reads: In order to obtain the vector z that minimizes this expectation value, one should solve the generalized eigenvalue problem: where This algorithm is more conveniently implemented in the socalled semi-orthogonal basis 15 : where which can be readily evaluated by a Monte Carlo sampling, using not only the random variables O k (x i ) necessary for the simpler SR technique, but also the parameter derivatives of the local energy ∂ α k e L (x i ), that can be directly computed by AAD, and thus allow the calculation of the above Hamiltonian matrix elements, by straightforward algebra: Notice that the code does not always take the eigenvector corresponding to the lowest eigenvalue (that in turn may be also complex, as in the linear method the matrix H k,k ′ is not Hermitian), but takes the one that maximizes |z 0 | 2 (i.e., the coefficient of the zeroth-order term) in order to have the most stable parameter change, namely the updated WF being as close as possible to the old one. Finally, the code updates the variational parameters by using the obtained z and an input parameter ∆ (denoted as tpar) as: where ∆ ≃ 1/3 is the default TurboRVB choice that is much more stable than the original algorithm (∆ = 1), but it is approximately 1/3 slower. The linear method is usually rather unstable for large systems and several parameters, and also, in such cases, each iteration requires a diagonalization of a huge matrix (a task that can be often prohibitive). In Tur-boRVB, a practical compromise has been devised by using the linear method for a restricted variational space containing up to a maximum number (npbra) of variational parameters with largest signal to noise ratio (optionally larger than parcutpar, see Eq. 113) and/or a number (ncg) of global line parameter directions spanned by the natural gradient ones (e.g., ∆ is one variational parameter of this form in Eq. 131) calculated at the given optimization iteration or the previous closest ones 41 , optionally neglecting the ones with smaller signal-noise ratio (< parcutmin).

C. Practical Rule
In most optimizations, the number of samplings in VMC is much larger than the number of optimized variational parameters and optimization is stable and efficient. In the opposite case, the preconditioning matrix S becomes rank-deficient singular; therefore the optimization does not work properly (c.f. the Cramér-Rao inequality). A practical rule is to set the number of samples M such that 15 : where p is the number of optimized variational parameters. However, one can also employ very small number of samplings with sufficiently large ε in the stochastic reconfiguration method because the regularization of the matrix S in Eq. 130 can remove all instabilities as long as ε is sufficiently large, as in the infinite ε limit one recovers the very powerful stochastic gradient technique well-known within machine learning. This method seems promising for optimizing even a huge amount of variational parameters and could certainly represent an important development in QMC. At present, the inversion in Eq. 131 can be done in Tur-boRVB without storing explicitly the matrix 109 as the associated linear problem is solved iteratively using conjugate gradients and implicit small matrix-vector operations, optimally distributed within the MPI protocol. Thus the cost of this inversion becomes negligible in the large ε limit, and linear with the number of variational parameters.

VIII. MOLECULAR DYNAMICS
In TurboRVB, several types of ab initio molecular dynamics (MD) have been implemented for both classical and quantum nuclei. The MD is driven by ionic forces . . , R N }, and the potential energy landscape V(R) is evaluated by VMC, namely In Eq. 141, |Ψ R is the QMC WF that, according to the Born-Oppenheimer (BO) approximation, minimizes the expectation value of H(R) at each ionic position R. The expression for the ionic forces is quite complex. Computing these forces with finite variance and in a fast way, as done in Tur-boRVB, (Sec. VI B), is of paramount importance to make a QMC-based MD possible. Moreover, in the QMC framework the evaluation of V(R) and its corresponding force estimators are intrinsically noisy. The statistical noise must be kept under control, if one wants to have an unbiased sampling of the phase space during the propagation of the trajectory. This issue has been solved by resorting to a Langevin type of molecular dynamics in both classical and quantum formulations, where the QMC noise becomes a controlled source of thermalization at the target temperature T . Two types of Langevin dynamics (LD) have been implemented in TurboRVB: the second-order LD, in its classical (Sec. VIII A) and quantum (Sec. VIII B) variants, and the firstorder LD accelerated by the covariance matrix of QMC forces (Sec. VIII C).

A. Second-order Langevin dynamics
The equations of motion of the second-order LD read 105 : where R, υ, f, η are the 3N at -dimensional vectors representing atomic positions, velocities, deterministic and stochastic forces of N at atoms, written in mass-scaled coordinates: for i = 1, . . . , N at . The stochastic forces are related to the friction matrix γ through the fluctuation-dissipation theorem, namely: with the temperature T expressed in atomic units. TurboRVB exploits the freedom in Eq. 146, by assuming: where s 0 and ∆ 0 are constants, I is the identity, and S QMC (R) the covariance matrix of QMC forces: In the above Equation, · · · refers to the average over the Monte Carlo sampling at the given MD step. The friction, a quantity that controls the sampling efficiency, is, therefore, position-dependent now. Luo et al. has shown that the force covariance matrix S QMC i, j (R) is, within a good approximation, proportional to the dynamical matrix (Fig. 8). Therefore, the choice in Eq. 148 realizes the nearly optimal damping of the high-frequency modes. After discretizing Eqs. 142 and 143 in the time interval t n − τ/2 < t < t n + τ/2, where the index n denotes the time slice t n = nτ, with time step τ, and by integrating the above Equations one obtains: υ n+1 (t) = e −γ n τ · υ n + Γ n (F n +η) with F n ≡ F (R n ), γ n ≡ γ (R n ), andη = γ n 2 sinh (γ n τ/2) t n +τ/2 By using Eq. 152, it is easy to show that the correlator defining the discrete noise is given by the following 3N at × 3N at matrix: In this way,η fulfill the fluctuation-dissipation theorem required by the Langevin thermostat. However, part of this noise is already present in the QMC force evaluation (η QMC ), as quantified by Eq. 148. Hence, TurboRVB adds the external random forceη ext , such thatη =η ext +η QMC . This is generated according to the following correlator: in order to fulfill the original fluctuation-dissipation theorem.
Since the correlation matrix in Eq. 154 is positive definite as long as τ < ∆ 0 , the correlated external force can be readily generated by standard algorithms (e.g., the Box-Muller's method 111 combined with the matrix diagonalization). This noise correction plays an important role in realizing a stable target temperature even in the presence of statistical errors in the QMC forces. The original numerical integration detailed here has been further improved by G. Mazzola with τ β = β/L denotes the imaginary time step. We have replaced here the true quantum particles by fictitious classical ring polymers, consisting of L replicas (beads) of the system. The beads in these necklaces are connected to each other by harmonic springs, evolving in the imaginary time. Without loss of generality, we can extend the definition of classical vectors to the quantum case, by including all L replicas, such that the resulting vector q ≡ R (1) , . . . , R (i) , . . . , R (L) T is 3N at L-dimensional. If interpreted classically, the partition function Z in Eq. 155 describes a system at the effective temperature LT . The Hamiltonian H L corresponding to this quantum-to-classical isomorphism reads as written in mass-scaled coordinates Eq. 145, withω L = L/β , and subjected to the ring boundary condition R (0) i ≡ R (L) i . The system is thermalized by a Langevin thermostat, such that the related Liouville operator L can be written as: in mass-scaled coordinates, with ∂ q i ≡ ∂ . L is built upon the Hamiltonian propagation, driven by the first two terms, and the Langevin thermostat, represented by the last two, deriving from the Fokker-Planck equation. Hereafter in this Section, i and j run over all particles and beads indexed together. In Eq. 157, F i ≡ F BO i + F harm i is the total force acting on each replica, comprising the BO (intra-replica) and harmonic contributions (inter-replicas), where F BO i ≡ −∂ q iṼ and 1 , .., R ( j) N at ). The quantum-to-classical isomorphism Hamiltonian in Eq. 156 includes very different energy scales. To cope with this, we split the Liouvillian in Eq. 157 into two operators, one containing only the physical (BO) modes, the other depending exclusively on the harmonic modes of the rings. To do so, we first separate the friction matrix into two contributions: We can then rewrite the total Liouvillian as the sum of two terms, L = L BO + L harm , where in such a way that we can break up the evolution operator via a Trotter factorization 114 to get the following symmetric propagator 115,116 : The equations of motion corresponding to the propagator in Eq. 161 have been implemented in TurboRVB, where now F BO i can be computed in a QMC framework, and therefore, it can be affected by an intrisic noise, as it was in the classical scheme (Sec. VIII A).
We note that the equations of motion generated by L harm are linear in both p and q, thus exactly solvable in an analytic closed form. This is because the non-linear BO components of the total force have been put in the L BO factor. The dynamics generated by L harm is a quantum Ornstein-Uhlenbeck process, which describes the motion of Brownian quantum particles under the influence of friction.
The exact integration of L harm without further splitting the Langevin thermostat from the harmonic modes of the rings is a major achievement, which gives the algorithm an enhanced stability and better scaling with respect to the time step size τ and the number of replica L. Moreover, the time step error is significantly smaller than the one of other PILD methods 117 . Thus, this is a key feature of the PI algorithm implemented in TurboRVB. This is also why the implemented method is called path integral Ornstein-Uhlenbeck dynamics (PIOUD).
If we now look at the L BO factor in Eq. 161, the corresponding equations of motion are the ones in Eq. (12) of Ref. 110.
Indeed, only the momenta are evolved in L BO , and the resulting equations are equal to those of the simple classical Langevin algorithm introduced in Ref. 105, restricted to p, which in the present case is a 3N at L-dimensional vector. As in the classical second-order QMC-driven LD, the QMC noise affecting the BO forces is controlled thanks to the friction γ BO and the fluctuation-dissipation theorem. A noise reduction scheme similar to the one described in Eqs. 153 and 154 is also applied in this case, yielding a thermalization towards a steady temperature ensemble even for quantum particles.
We refer the interested reader to Ref. 113, where we report a detailed description of the PIOUD algorithm and the integrated equations of motion.

C. First-order Langevin dynamics
TurboRVB also features an accelerated first-order dynamics to sample the (classical) ionic canonical distribution. This molecular dynamics scheme has been introduced recently by G. Mazzola and S. Sorella 118 . Here, the ions formally undergo the following dynamicṡ where T is the temperature, F (R) and η are the deterministic and stochastic forces, S is the covariance matrix of the QMC forces as introduced in Eq. 148. This dynamics generalizes the Newton method, which is valid for structural relaxation, to finite-temperature. Indeed, the use of the matrix S, which is empirically proportional to the Hessian matrix, can drastically decrease the autocorrelation time of the Markov chain. In other words, this metric plays an important role in treating the different time scales of a given system. In this way, one can still use a quite large integration time step τ, which is beneficial in reducing the autocorrelation time of quantities that depend on slowly varying degrees of freedom (e.g., molecular diffusion), without introducing any bias in the sampling.
Discretization of Eq. 162 for a finite time step τ at finite temperature is rather involved and has been described in the Ref. 118, and reads with As the effective temperature depends on the finite integration time, the convergence to the target temperature T target for τ → 0 can be improved by a renormalization of the temperature T used in the simulation, namely by using the above equations with an appropriate choice of T : where a = 1/2 in the case of the harmonic potential leads to the exact result for any discretization time τ, as long as S is the exact elastic matrix. In the general case, the parameter a should be tuned to speed up convergence to the τ → 0 limit. Conversely, the target temperature of the simulation can be measured a posteriori, as it is often done in the case of secondorder Langevin dynamics, by measuring the average squared velocities of the particles. The accelerated first-order Langevin dynamics has been successfully employed in the most recent studies of dense liquid hydrogen, providing equilibrated simulations even close to phase transitions 118,119 . Moreover, nuclear quantum effects can also be computed within this first-order LD, in exactly the same way as it was implemented for the second-order case discussed in the previous subsection. However, so far, it is not clear which method is more efficient and further studies are also necessary to clarify what is the optimal choice of the acceleration matrix S (R) within the QMC framework. Table II shows the main modules implemented in Tur-boRVB and a brief description of their functionalities. All programs reported in Table II are coded in Fortran90. The flow-chart of a typical QMC calculation (see Fig. 1) is detailed as follows: starting from the geometry and the chemical composition, the user first chooses the basis set for both Jastrow and antisymmetric factors, and generates a JAGP-type ansatz using makefort10.x (fort.10 is the wave function filename in TurboRVB).

IX. IMPLEMENTATION AND OPERATIONS
The antisymmetric part of the WF is initially determined at the DFT mean-field level. Since the DFT calculation is based on a single Slater determinant, the generated JAGP WF should be converted to a JSD one. This is performed by the convertfort10mol.x module, which adds M molecular orbitals to the WF, where M = N ↑ . 120 The DFT molecular orbitals are then used to construct the initial trial WF for subsequent QMC calculations. In QMC, the user can choose to work with five different ansätze, namely, JSD, JAGP, JAGPn, JPfn, and JPf. The initial antisymmetric part of the JSD and JAGPn wave functions can be directly obtained from the DFT orbitals. In the JSD case, only the occupied orbitals are imported, while in the JAGPn WF also the empty DFT orbitals will be used in the geminal expansion, up to the n-th orbital. Another possibility is to employ a JAGP ansatz. In that case, one has to convert the initial JSD trial WF obtained by DFT to the JAGP one using convertfort10.x.
JAGPn can also be obtained by applying convert-fort10mol.x to a previously determined JAGP ansatz. Anal-ogously, a JPf WF can be converted from a JAGP WF using convertpfaff.x. JPfn can be obtained by applying convert-fort10mol.x to the JPf ansatz. These conversions can be done either before or after a QMC optimization. The possibilities for ansatz conversion are schematically drawn in Fig. 4.
After the initial determination of the antisymmetric part of the WF, one performs the WF optimization at the VMC level. This is an important step used to determine the best variational parameters of the Jastrow factor at fixed determinant, or to fully relax the WF by optimizing both Jastrow and determinant parameters.
Finally, the user can proceed to VMC and LRDMC calculations using the optimized WF. In LRDMC, in order to remove the lattice discretization error, it is better to repeat the calculations with several lattice spaces a, and to extrapolate the results with a quartic function E (a) = E 0 + k 1 a 2 + k 2 a 4 , where E 0 is the extrapolated (a → 0) energy.
The python tool Turbo-Genius makes all the above steps more straightforward, as discussed later.

X. PARALLELIZATION AND BENCHMARKS
TurboRVB has been parallelized using the MPI and OpenMP libraries, also supporting the hybrid MPI/OpenMP protocol. As well known, the QMC algorithm is readily parallelized by employing several independent random walks (called walkers), initialized with different seeds for random number generators. Fig. 9 shows the strong and the weak scaling performances of TurboRVB, which are compared with the ideal scaling. The benchmarks have been measured using the conventional 2 × 2 × 2 diamond supercell containing 64 carbon atoms for a simulation with 256 electrons on Marconi KNL nodes (Cineca supercomputer/SCAI, Intel Xeon Phi 7250 CPU, 1.40 GHz, 68 cores per a node without hyperthreading). For the weak scaling case, we set the number of walkers M to the minimum possible one (1 walker per a core) and increased the number of cores from 68 (1 node) to 27200 (400 nodes). On the other hand, for the strong scaling case, we fixed the number of walkers M = 2176 and increased the number of cores from 544 (8 nodes) to 8704 (128 nodes) with the hybrid MPI/OpenMP protocol (the number of openMP threads was set to four). For both cases, the performances are excellent, because all techniques are off from the ideal scaling by ∼ 5% at worst in all cases performed here. At present, the code does not support the use of GPU accelerators. This feature is currently under development.

A. Correlated Samplings
A very efficient and easy-to-use correlated sampling technique has been implemented in TurboRVB. This is very useful in evaluating total energy differences between two very similar WFs. One of its practical sides is to investigate whether

Module
Description makefort10.x Generates a JAGP WF from a given basis set and structure. assembling pseudo.x Prepares pseudo potentials. convertfort10mol.x Adds molecular orbitals. If the number of molecular orbitals is equal to (larger than) the half of the number of electrons in a system, the resultant WF becomes JSD (JAGPn). convertfort10.x Converts a JSD/JAGP/JAGPn WF to a JAGP one. It also converts an uncontracted atomic basis set to a hybrid (contracted) one using the geminal embedding scheme. convertfortpfaff.x Converts a JAGP WF to JPf one. prep.
x Performs a DFT calculation. turborvb.x Performs VMC optimization, VMC evaluation, LRDMC, structural optimization and molecular dynamics. readforward.x Performs correlated samplings, and calculates various physical properties. a conversion of a WF (e.g., JSD → JAGP) has been successful. A user should run a short VMC calculation using one of the two mentioned WFs, during which the walkers' history is stored in an appropriate scratch area (turborvb.scratch). After that, a second run should be executed (by using the readforward.x module) to determine at the end the correlated energy difference and the overlap between the two WFs. Correlated samplings for LRDMC or for the difference in other quantities, e.g., forces, have not been implemented yet. The latter could be very important for the frozen phonon calculation based on the ionic force difference method, which represents a possible subject of research.

B. Physical properties
TurboRVB enables us to calculate various properties using the optimized many-body WF at both the VMC and LRDMC levels by averaging local observables and applying the forward walking technique 15 . The observables presently available in the readforward.x module can be categorized as follows: ( i ) Charge density (ρ (r) = ρ ↑ (r) + ρ ↓ (r)) and spin den-sity (ρ σ (r)), fundamental properties obtained by a stochastic evaluation of multi-dimensional integrals, involving the many-body WF and improved estimators, such as those proposed by Assaraf, Caffarel, and Scemama 121 . Fig. 10 shows the charge density of the square four hydrogen (H 4 ) depicted by C. Genovese et al. (ii) Electronic correlation functions, such as the charge-charge and spin-spin structure factors, useful to study the critical behavior close to a phase transition, or the physical properties of a given phase. (iii) The expectation values of the spin component along the quantization axis (S z ) and the spin square (S 2 ) operators inside a sphere centered on each atom 61,68 , whose radius can be specified by the user. (iv) The Berry phase: where the complex polarization is computed for the α = {x, y, z} component r α i of the position vector r i , and L α is the supercell length in the same direction α 122 . This property can be used to characterize the metal-insulator transition based on the many-body WF, as it has been done in the one-dimensional hydrogen chain by L. Stella et al. (Fig. 11).

XII. APPLICATIONS
Since the start of the TurboRVB project in 2003 by M. Casula and S. Sorella 64 , the code has been applied to study atomic species, molecular systems, and various materials, including challenging systems such as large complexes, surfaces, liquids, and so on. Here, we provide a brief review of the applications done so far.

A. Molecular systems
TurboRVB has been employed to study the properties of several molecular systems, including small diatomic systems (such as Li 2 , 126 Fig. 13 report the binding energies for the molecular dimers of the first-row atoms, evaluated using LRDMC with the JSD, JAGP, or JPf WF ansätze as guiding functions. Some of them are close or reach the chemical accuracy (i.e., 1 kcal/mol or ∼0.04 eV). A systematic improvement going from the JSD → JAGP → JPf can be appreciated. It fol-lows from the employment of a more flexible parametrization of the WF, which improves both the variational energy and the quality of the nodes. However, the improvement given by the JPf ansatz over JAGP and JSD is even more fundamental than what Table III and Fig. 13 indicate, where the reported binding energies are evaluated as the difference between the energy of the molecule and twice the energy of the atom. This evaluation does not highlight that JPf is a size-consistent ansatz, while both the JSD and the JAGP have a size-consistency issue with most of the molecules considered. 136 Let us consider the oxygen molecule as an example: the O 2 ground state is a spin triplet ( 3 Σ − g ) which, in either JSD or JAGP, can be described using two unpaired electrons with the same spin. However, as we take the oxygen atoms far apart (say, the O-O system), they should both be in the ground state for the oxygen atom, which is also a spin triplet ( 3 P). Neither JSD nor JAGP have the flexibility to describe the O-O system by the same WF parametrization used for the O 2 , 137 thus E O-O > 2E O . In contrast, JPf has the flexibility to describe the dissociation limit Fig. 14 shows that JPf is indeed able to provide a reliable dissociation curve of the oxygen molecule, whereas JSD and JAGP are reliable only in proximity of the minimum. Notice that the FN projection scheme alleviates in part the limitations of the underlying ansatz, as both JSD and JAGP have a smaller size-inconsistency at the LRDMC than at the VMC level, but the correction is incomplete and the inconsistency is sizable (> 1 eV) for all but the JPf ansatz. Therefore, in general, JPf will yield a better description of the overall potential energy surface (PES), including binding energies, vibrational properties, zero-point motion, transition states, and so on.  The applications of TurboRVB are not limited to single point evaluations. As discussed in Sec. VI, a distintive feature of TurboRVB among most QMC codes is the efficent evaluation of atomic forces. This capability has been extensively exploited to investigate several properties of the PES, such as vibrational properties, and to perform geometry optimizations and molecular dynamics simulations at finite temperature. It is important to remember that QMC atomic forces are, as the QMC energy, affected by a stochastic error. In a first attempt 143 to evaluate vibrational frequencies of a triatomic molecule (namely, water), a grid of points close to the guessed minimum geometry was used to evaluate energy and forces. The results were used to fit the PES around the minimum with a truncated Taylor expansion, yielding evaluations of the structural minimum and of both the harmonic and anharmonic vibrational frequencies (using the second-order perturbation theory). It was observed that the inclusion of the estimates of forces reduces the stochastic error on the vibrational frequencies by around an order of magnitude with respect to the fit using only the energy evaluations (on the same grid and using the same sampling size on each single point evaluation). This corresponds to about a 25 times larger speedup. 144 Moreover, it was observed that the optimal grid is a tradeoff between maximizing the precision on the vibrational properties and minimizing the errors on the fit. A systematic study of the quality of the variational WF ansätze in relation to the vibrational (and other) properties of the water molecule is reported in Zen et al. 132 . Afterward, Y. Luo et al. 110,145 automated the process of finding the optimal grid by performing fits on samples generated by a finite temperature molecular dynamics, and applied the approach to H 2 O, H 2 S, SO 2 , NH 3 , and PH 3 . The obtained vibrational frequencies are in good agreement with the experimental values.
Zen et al. 134 applied the second-order LD to study the structural properties of liquid water. They obtained radial distribution functions in good agreement with recent neutron TABLE III. Binding energy of the first-row dimers calculated by TurboRVB with LRDMC and different types of ansätze as a guiding function (GF), at the experimental bond length (see. Ref. 138). Js denotes the symmetric (i.e., spin-independent) Jastrow factor. LRDMC calculations are performed with several lattice spaces a and extrapolated for a → 0. The spin-orbit coupling and the zero point energies 70  scattering and X-ray experiments, especially for the position of the oxygen-oxygen peak, which is very sensitive to the quality of the water description, in particular to the inclusion of dispersive vdW forces. This is a stringent test for theories as the shape of the peak can be directly compared with neutron scattering data.
F. Mouhat et al. developed a path integral LD formalism particularly suited for QMC and applied it to the Zundel ion (H 5 O + 2 ), with the aim of studying nuclear quantum effects 113 . The outcome of their study shows how essential is the inclusion of nuclear quantum effects to properly deal with proton transfer in water and aqueous systems, even at room temperature.

B. Weakly bonded systems and non-covalent interactions
An important application for QMC is the description of systems interacting through weak van der Waals (vdW) dispersive interactions, which are fundamental in supramolecular chemistry, layered and porous materials, molecular crystals, adsorption of molecules on surfaces, and so on. vdW interactions play a major role also in non-covalent intra-molecular bonds, which result from the interplay of forces of different nature, including electrostatic contributions. An accurate evaluation of the overall strength of these non-covalent bonds requires each component of the interaction to be modeled correctly. This is a big challenge for all ab initio approaches, and only a few methods prove reliable. QMC and CCSD(T) are considered to yield very accurate estimates, and there is a general agreement among the predictions from these two theories. Thus, they are typically used to produce benchmark values. Traditionally, mean-field schemes have difficulties in describing dispersive interactions, because they arise from nonlocal electron correlations, completely missed by local density functionals. However, in the context of DFT, there has been an intensive effort in developing new exchange correlation functionals 146 to include vdW interactions. This is usually done in a semi-empirical way, often building on the benchmark provided by CCSD(T) or QMC.
The non-covalent interactions are typically a tiny fraction of the total energy of a compound. Similarly to the cases discussed in Sec. XII A, also in the evaluation of weak interactions it is important to use a size-consistent approach. However, very often, the molecules form a closed shell, and most of the times the intermolecular interactions that we need to evaluate are among fragments having zero spin. In this case, a single Slater determinant is a size-consistent ansatz, and so is the JSD. The AGP (without Jastrow) is not sizeconsistent, due to the presence of unphysical charge fluctuations 41,147 , but a very flexible (i.e., large) Jastrow factor allows damping the charge fluctuations and recovering the sizeconsistency in the JAGP. This is of course true also in a subsequent FN-DMC projection, which is indeed equivalent to applying an infinitely flexible Jastrow factor to the trial WF. At the variational level, the JSD ansatz leads sometimes to better results than the JAGP. Even though JAGP leads to a lower total energy, it is indeed difficult to use a Jastrow factor large enough to damp any charge fluctuation. So, there is a worse error cancellation when energy differences are considered in JAGP. This was observed, for instance, in the water dimer, 134 although at the FN level JSD and JAGP lead to the same binding energy. It should also be considered that in the widely used DMC scheme 148 a size-consistency issue was present in any simulation with finite time step τ, until the modification of the "standard" algorithm has been introduced by Zen et al. 149 , whose solution removes almost completely the size-consistency error. At variance with DMC, the LRDMC algorithm does not suffer from this bias, and it is always size consistent.
The DMC scheme is commonly used to project a single Slater determinant WF where only the Jastrow factor has been optimized, while the antisymmetric part is filled with DFT orbitals 150 (JDFT). The resulting FN energy is usually biased by the nodal surface, kept frozen from previous DFT or lower level mean-field calculations. Optimizing the Jastrow factor and determinant together, which is routinely done in Tur-boRVB, improves significantly the total energy. However, in the evaluation of non-covalent interactions, there is not yet a clearcut indication of the advantage of the full optimization over JDFT. This is due to the fact that the error cancellation plays a major role. Indeed, it is expected that, for non-covalent interactions, JDFT has almost the same bias in the interacting and non-interacting systems. When pseudopotentials are used, DLA 53 can be used to remove any bias given by the Jastrow optimization.
TurboRVB with JDFT, JSD, and JAGP guiding functions has been applied to several interesting systems, such as molecular hydrogen adsorbed on benzene 151 , benzene dimers 41 , graphite layers 152 , water or methanol molecules adsorbed on a clay surface (see Fig. 15) 153

C. Strongly correlated and superconducting materials
Another important QMC application is the study of strongly correlated materials, that represents a challenge for any theoretical method, because in several cases a fully consistent description remains, so far, elusive. In this context, QMC methods have been widely used for solving, as accurately as possible, strongly correlated lattice Hamiltonians, like, for instance, the Hubbard model. Real-space ab initio QMC methods could fill the gap between realistic Hamiltonians and lattice models, by providing an accuracy similar to the one reached in lattice Hamiltonians, while keeping the complexity of real materials. TurboRVB has been originally conceived in this perspective, with the idea that one of the most flexible WF ansätze tested on correlated lattice models [156][157][158][159][160] , namely the generalized version of the resonating valence bond (RVB) WF, could be translated into a successful variational WF for ab initio systems. One of the most appealing features of the RVB WF is that it allows several broken-symmetry phases, such as antiferromagnetic insulators and superconductors, which are present -and sometimes coexisting -in the extremely rich phase diagrams of strongly correlated materials. As introduced in Sec. III, the RVB WF employed in TurboRVB takes the static and dynamic correlation effects into consideration beyond the commonly used Jastrow-Slater determinant (JSD) WF. Therefore, the code is expected to describe strongly correlated systems more accurately. In particular, among the various families of correlated WFs which go beyond the simplest JSD form, the JAGP/JPfaffian families are very suitable to describe extended systems, as they keep a compact and still strongly correlated form even in the thermodynamic limit. Indeed, as detailed in Sec. III, the AGP part is the particle number conserving version of the BCS WF. In the same way, the JAGP variational form can be seen as efficient formulation of a RVB WF. Thus, superconducting materials including cuprates 161 and iron pnictides and selenides 162 are prominent applications of the JAGP/RVB WFs. By running extensive VMC calculations for CaCuO 2 , a parent compound of cuprate high-temperature superconductors, M. Marchi et al. validated the JAGP description by correctly finding the expected d-wave symmetry of the pairing function 163 . M. Casula and S. Sorella applied the RVB WF to FeSe, one of the iron-based high-temperature superconductors 164 , whose electronic structure and pairing mechanism have been unclear and intensively debated. They determined the symmetry of the superconducting order parameter and the size of its gap entirely from first principles, by analyzing the AGP pairing function (Fig. 16). B. Busemeyer et al. applied DMC and LRDMC to study the structural and magnetic properties of the normal state of FeSe under pressure, and reported that collinear spin configurations are energetically more favorable than other spin patterns, such as ferromagnetic, checkerboard, and staggered dimer, over a large range of pressures 165 .

D. Aromatic and honeycomb lattice compounds
Aromatic compounds are also successful examples of the RVB theory, because they are characterized by resonating C-C bonds that can be efficiently represented by antisymmetrized singlet pairing functions (i.e., antisymmetrized product of geminals, or AGP). M. Casula et al. applied the AGP WF combined with a Jastrow factor to the simplest aromatic compound, that is, the benzene molecule (C 6 H 6 ) for the first time 126,130 . They reported that the inclusion of the resonance between the two possible Kekulé states significantly lowers the VMC energy compared with the one obtained by a single Kekulé WF. They also showed the importance of adding a three-body Jastrow factor to improve the AGP description, just in the spirit of the RVB framework. By selectively switching on the intersite resonances, they proved that the Dewar contributions (i.e., including the third nearest neighbor carbons) improves the description of the resonating valence bond. C. Genovese et al. have applied a more general AGP WF, i.e., the Pfaffian, to the benzene molecule 68 . They reported that the obtained atomization energies are reasonably consistent with the experimental value at the LRDMC level. On the other hand, when starting from JDFT, i.e., from frozen DFT nodes, the obtained atomization energy is severely underestimated in most cases, implying that the optimization of the nodal surface at the VMC level is essential to obtain a correct value. As mentioned in Sec. XII B, S. Sorella et al. applied the RVB WF to calculate the binding energy of the face-to-face and displaced parallel benzene dimers 41 , and obtained the binding energy very close to the experimental value in the displaced parallel case. N. Dupuy and M. Casula applied the same RVB WF to study the ground-state properties of the oligoacene series, from anthracene up to the nonacene. 131 They found that the ground state obtained by the RVB WF has a weak diradical or polyradical instabilities until the nonacene, which is in contrast with the results previously found by lower-level theories, such as Hartree-Fock, hybrid functionals, and CASSCF in a restricted active space. TurboRVB has been applied not only to simple aromatic molecules but also to honeycomb lattice compounds, in order to study the role of RVB correlations, in other words, to investigate whether the resonance energy is sizable also in extended systems. M. Marchi et al. studied the nature of chemical bonds in graphene using the RVB WF 163 , and found that the RVB energy gain becomes extremely small in the thermodynamic limit i.e., in the infinite lattice. Conversely, S. Sorella et al. applied the RVB WF to isotropically strained graphene 166 (Fig. 17) and found that the RVB effect is crucial in the Kekulé-like dimerized (DIM) phase that has been proposed to become stable under the isotropic strain (Fig. 18). They also reported that the JAGP energy gain with respect to JSD is negligible in the perfect honeycomb structure -as previously found by M. Marchi et al. -but is extremely important in the DIM phase.

E. Organic compounds: diradicals and conjugated systems
Diradical and conjugated organic compounds are other interesting applications of the RVB WF, because both static and dynamic correlation effects play important roles in these sys-  65 . They also investigated the chargetransfer and diradical nature of the electronic states of the retinal minimal model (penta-2,4-dieniminium cation) at the VMC and LRDMC levels 66 . They proved that the dynamical electronic correlation is key to get a reliable ground state energy surface in proximity of the conical intersection between the electronic ground-and first excited-state of the molecule. Their obtained energy landscapes significantly different from the CASSCF ones, inverting the relative stability of several torsional paths, similar to what was obtained by other correlated approaches. M. Barborini et al. applied VMC and LRDMC to calculate the torsion barriers of 1,3-butadiene and the ring-opening barrier of cyclobutene 167 , which are in good agreement with the experimental results. M. Barborini and E. Coccia investigated the potential energy curves of the two spin states of tetramethyleneethane (TME) molecule and the two anionic states of TME as a function of the torsion of the central dihedral angle 168 . Through ab initio geometrical optimizations at the VMC level, they proposed possible structural interconversions between the states, which are in agreement with the ion photoelectron spectroscopy experiments.
The bond length alternation (BLA), namely, the impact on the geometry of the difference between the single and double carbon bonds, is one of the key structural descriptors in conjugated organic compounds. M. Barborini et al.  consistency between BLAs obtained by CCSD(T)-CBS (the most accurate) and VMC-JAGP calculations, but the reason for the discrepancy is still unclear. They also investigated the structural properties of polyacetylene chains H-(C 2 H 2 ) n -H up to N = 12 acetylene units 170 , in which they revealed that the BLA obtained by the extrapolation to n → ∞ is 0.0910(7) Å, which is compatible with the experimental data. E. Coccia et al. applied VMC to study geometries of the retinal model, the chromophore of the rhodopsin involved in the mechanism of vision, and discussed the effects of the electronic correlation on the BLA 171,172 by comparing with DFT and quantum chemistry methods. Wave-function and geometry optimizations of the retinal model were also performed in the presence of the protein field, classically described via electrostatic and mechanical interactions: the complex environment dramatically affects the BLA of the chromophore, which consistently increases with respect to the gas-phase case. They also performed VMC geometry optimization of the peridinin chromophore in order to verify the interplay between the BLA and the optical properties of a large conjugated moiety, involved in the light-harvesting step of photosynthesis in the peridinin-chlorophyll-protein complex 173 . Comparison with DFT and wave-function approaches shows that the combined application of RVB-based VMC and Bethe-Salpeter methods for geometry and excitations, respectively, gives an accurate estimation of the electronic transition energies of peridinin. They recently applied VMC optimization to obtain a ground state structure of keto-1, enol-1, and enol-2 forms of oxyluciferin 174 , which is used to calculate its S 1 ← S 0 absorption energy based on the Bethe-Salpeter formalism. E. Coccia et al. carried out a systematic basis-set analysis on polarizability, quadrupole moment and electronic density of C 2 H 2 molecule, with the aim to test VMC and LRDMC calculations in the presence of an external electric field 175 . In particular, a relatively small hybrid (Gaussian + Slater functions) basis set was seen to be accurate enough to properly reproduce reference values.

F. Metal-organic complexes
Recently, TurboRVB has also been used to study metalorganic complexes. S. Chu et al. studied free-energy reaction barriers of water splitting reactions, using a simplified computational model based on the cobalt ion 176 (Fig. 19). They found that the total free-energy differences of the water oxidation, computed at the QMC level, fairly agree with the experimental reference values, in a surprisingly better way than the corresponding CCSD(T) ones. M. Barborini and L. Guidoni applied VMC to tackle the geometry optimization of a Fe 2 S 2 (SH) 2 − model complex (High-Spin and Broken Symmetry states) 177 based on the extended broken-symmetry (EBS) approach. The number of applications is still limited in this category, but the recent development of sophisticated transition-metal pseudopotentials will make it more feasible in the near future.

G. Crystals polymorphism
Ab initio electronic calculations allow one to study the phase stability between crystals polymorphs, by comparing the enthalpies or free energies, and by computing their equation of states (EOSs). As mentioned in the introduction, DFT's predictive power strongly depends on the choice of the exchange and correlation functionals, sometimes causing qualitatively different phase stability predictions. On the other hand, VMC and LRDMC calculations do not suffer from this drawback in principle. S. Sorella et al. studied the pressureinduced metal-insulator transition in bulk silicon through EOS calculations 178 . They obtained a reasonable transition pressure at room temperature by VMC, however still a little far from the experimental value. N. Devaux et al. studied the α-γ phase transition in elemental cerium at the VMC and LRDMC levels 179 . They revealed that the volume-collapse transition could be understood as a conventional first-order transition of electronic origin, induced by the p-f hybridization, and by the strong local repulsion affecting the f states. H. Hay et al. investigated the effect of long-range van der Waals forces in the SiO 2 polymorphs using DFT and QMC 180 . Their LRDMC calculations showed that the obtained energy differences in quartz-cristobalite and quartz-stishovite agree well with the experimental values within sub-chemical accuracy, implying that QMC is a promising method in describing both conventional and high-pressure SiO 2 polymorphs. Notice that, at present, only the zero-temperature energy (without the entropic term) is available at the VMC and LRDMC levels. Therefore, the entropic effects are estimated by DFT-PBE calculations 178 . In the near future, it will be possible to estimate the entropic corrections directly within the VMC framework, via phonon calculations under the quasi-harmonic approximation. As described in Sec. IV, considering twist boundary conditions is essential to correct the one-body finite-size errors for periodic systems, especially in a metallic case. Although the twisted average approach is commonly used to correct the error, one can also use the exact special twist (EST) method 86 developed by M. Dagrada et al. to save some computational cost. EST works even in realistic systems, namely, the bcchydrogen, the bcc-lithium, and the silicon in the high-pressure β-tin phase 86 .

H. One-dimensional chains
One-dimensional chains of atoms have been intensively studied based on both effective and ab initio Hamiltonians, because they embody many important open questions in modern condensed matter physics, despite their simplicity. The finite or infinite hydrogen chains are ideal systems for this purpose, because it is not necessary to use a pseudopotential or to consider the relativistic effects. They are therefore benchmark systems. L. Stella et al. investigated the metal-insulator transition in the hydrogen chain by explicitly calculating the complex polarization 122 , as shown in Sec. XI. They found that the model has not a metal-insulator transition and always behaves as an insulating 1D Hubbard model at half-filling, at least for not too small interatomic distance (R > 1). M. Motta and his Simons collaborators on the many-electron problem 181 undertook a comprehensive benchmark study in the finite and infinite linear hydrogen chains using state-of-the-art manybody methods, including VMC and LRDMC with extrapolations to the thermodynamic and the complete-basis-set limits.
They provide accurate potential energies curves online, which will be very useful for benchmarking other methods. They have recently investigated the ground states of the hydrogenchain more in details, namely, the insulator-to-metal transition, dimerization, and magnetic phases 182 . They revealed a fascinating phase diagram, with several emergent quantum phases, depending on the interproton distance.

I. Liquid and solid hydrogen/helium
Hydrogen behavior at high pressures is still not fully understood and subject of intense research. In 1935, Wigner and Huntington predicted its metallization upon compression 183 , while Ashcroft also proposed its high-T c superconductivity 184 . Besides its fundamental interest as the simplest realistic condensed matter system in Nature, calculating its equilibrium properties remains of fundamental importance for planetary science applications 185 . This system is particularly challenging for DFT-based simulations, due to the presence of strong correlation effects and a subtle interplay between structural and electronic phase transitions. A straightforward way to study a phase diagram under pressure and at finite temperature in QMC is to apply the Langevin molecular dynamics described in Sec. VIII. First-order phase transitions can be identified from a discontinuous behavior in the calculated equation of state, and of pair correlation functions. C. Attaccalite and S. Sorella have formulated a second-order LD suitable for QMC simulation for the first time, as described in Sec. VIII. They applied it to assess the relative stability of solid and liquid hydrogen at high pressure and room temperature 105 . Years later, in a series of works, G. Mazzola et al. applied this technique to study the proposed liquid-liquid transition (LLT) from a molecular insulating to an atomic conducting phase of hydrogen 112,186 . G. Mazzola and S. Sorella have also formulated the accelerated first-order LD 118 , as described in Sec. VIII, to study such an elusive phase transition with properly equilibrated simulations 118 . Finally, under the same framework the first QMC molecular dynamics simulations of an hydrogenhelium mixture has been performed at Jupiter's interior conditions 119 . The study revealed that mixing He has a significant influence on the H metallization pressure as shown in Fig. 20, and provided useful benchmark for the equation of states currently adopted by the planetary science community. M. Dagrada et al. revealed that their developed EST technique also works well in the LD simulation for liquid hydrogen, but a perfect agreement with the twisted-average technique is achieved only when relatively large supercells are used 86 . Therefore, a careful finite size scaling analysis is needed when the EST technique is employed, though the thermodynamic limit extrapolation is much smoother at the special twist, saving some computational cost. Moreover, EST allows for the determinant/AGP optimization.

J. Excited states
TurboRVB allows studying not only the ground-state electronic properties but also excited states. It is trivial to obtain excitation properties when they are characterized by different particles or different spin projection, because these quantities are clearly conserved. For example, M. Barborini et al. calculated the 3 B 1u ← 1 A g vertical triplet excitation of the ethylene molecule and its adiabatic excitation 3 A 1 ← 1 A g 187 . One can also deal with nontrivial cases based on a flexible symmetryadapted RVB WF within its MO representation, which was proposed by N. Dupuy et al. in 2015 188 . Instead of dealing with a linear combination of symmetry-adapted configuration state functions, they construct a symmetry-adapted RVB WF that accurately describes the targeted symmetry. They reported that their developed constrained (fixed-rank) minimization of the geminal expansion is stable without symmetry contamination of the starting excited state symmetry; as a result, the vertical ionization energies and the electron affinities of anthracene obtained by VMC and LRDMC calculations agree with the experimental values.  190 . The STM image of [CuCl 4 ] 2− obtained by this approach is shown in Fig. 21. It reveals that the electronic correlation has a significant effect on the QPWF in this system. Indeed, the resulting STM image clearly differs from that obtained by uncorrelated methods (i.e., the Hartree-Fock counterpart). In general, the electronic calculation should be useful to interpret the experimental results theoretically. We believe therefore that Tur-boRVB can become useful also for state-of-the-art research in this field.
FIG. 21. The STM images of in-plane [CuCl 4 ] 2− simulated using UHF and QMC calculations. The left column shows the isosurfaces of the square modulus of the eQPWF for a fixed value of probability density of 9.7 × 10 −5 Å −3 , while in the right column the isosurfaces take the value of 6.7 × 10 −6 Å −3 . The top panels (a and a ′ ) display the LUMO+1 UHF molecular orbitals, while the middle panels (b and b ′ ) display the QMC images. Reprinted with permission from Barborini et al. 189 , J. Chem. Theory Comput., published by ACS in 2016.

XIII. PYTHON MODULE:TURBO-GENIUS
In general, a computational study involves many complicated operations, such as preparing input files, searching for optimal parameters, performing benchmark studies, transferring information coming from previous calculations to a subsequent simulation's input file, and analyzing output files. Automatizing these operations can significantly save researcher's burden, avoid trivial human errors, and is beneficial for accelerating the distribution of a computational package to a much wider community. We are now developing a python-based workflow system suitable for TurboRVB named Turbo-Genius, like Nexus 76 and QMC-SW 191 . The workflow system has been implemented in Python 3 in an object-oriented fashion, so the modules are highly extensible. Notice that, since TurboRVB is an all-in-one package, any other commercial code is not necessary to run the code. There are two main modules, turbo-genius-serial.py and turbo-geniussequential.py, which manage other modules/classes. The former allows one to perform a single job specified by an option (e.g., turbo-genius-serial.py -j makefort10), and the latter does a sequential job (e.g., makefort10.x → convertfort10mol.x → prep.x → turborvb.x). By using Turbo-Genius, one can also analyze a simulation output. For example, one can plot an optimization history, average the variational parameters after optimization, and search for an optimal number of equilibration steps and a reblocking length. 38 Turbo-Genius is still at an early stage of development. However, its main features are already there. In the near future, the Turbo-Genius modules will be capable of sophisticated TurboRVB jobs automatizations, as well as setting optimal parameters depending on target systems.

XIV. CONCLUSIONS
In this paper, we have reviewed TurboRVB, an ab initio quantum Monte Carlo package featuring two well established QMC algorithms for electronic structure: variational Monte Carlo, and diffusion Monte Carlo. The beginning of the Tur-boRVB development dates back to 2003, with the Ph.D thesis project of M. Casula supervised by S. Sorella. Since then, TurboRVB has rapidly grown up, and it has been applied to the study of several molecular and condensed matter systems, from benzene to high-temperature superconducting materials. Central to this package are variational parameterizations of correlated electronic WFs in the product form Ψ = Φ AS * exp J, made of an antisymmetric part (Φ AS ), and a Jastrow factor (J).
Concerning the antisymmetric part, TurboRVB implements five ansätze (in order of decreasing variational flexibility): ( i ) the Pfaffian (Pf), (ii) the Pfaffian with constrained number of molecular orbitals (Pfn), (iii) the Antisymmetrized Geminal Power (AGP), (iv) the Antisymmetrized Geminal Power with constrained number of molecular orbitals (AGPn), and (v) the single Slater determinant (SD). Notably, the user can freely navigate between these ansätze. Indeed, the package includes conversion modules, which allow the user to choose the most suitable WF form for the target system.
All the above ansätze can be optimized by state-of-the-art algorithms, namely, the stochastic reconfiguration and the linear method, thus achieving highly accurate variational and fixed-node energies, with a computational cost remaining at the single-determinant level, thanks to efficient algorithmic developments.
The stochastic optimization of several variational parameters (so far up to the order of 10 5 ) is feasible in TurboRVB thanks to an efficient evaluation of energy derivatives using the AAD. This algorithmic scheme, which gives the code an original architecture, drastically decreases the computational cost of ionic forces calculations, keeping them of order N 3 . Thus, it paves the way for efficient structural optimizations and molecular dynamics simulations at the VMC level of theory. Indeed, by means of the TurboRVB package, it is at present possible to perform QMC molecular dynamics simulations of dense hydrogen (with up to 256 ions) and liquid water (with up to 64 water molecules), with a remarkably large number of atoms, not far from the DFT state of the art. Thanks to TurboRVB, nuclear quantum effects are also accessible in a correlated electronic environment provided by the VMC WF.
Large scale calculations are possible with TurboRVB not only because of cutting-edge algorithmic developments but also by virtue of an efficient parallelization based on a hybrid MPI and OpenMP protocol, that is ideal for employing GPU accelerators as well, recently available in the most advanced HPC infrastructures.
A python wrapper, named Turbo-Genius, has recently been developed to make the code user-friendly, and to allow both beginners and professional users to handle it more efficiently. We will continue developing TurboRVB and Turbo-Genius to implement new QMC algorithms and to extend their range of applications for the current developers and expected new upcoming contributors and users.