TurboRVB: A many-body toolkit for ab initio electronic simulations by quantum Monte Carlo

Nakano, Kousuke; Attaccalite, Claudio; Barborini, Matteo; Capriotti, Luca; Casula, Michele; Coccia, Emanuele; Dagrada, Mario; Genovese, Claudio; Luo, Ye; Mazzola, Guglielmo; Zen, Andrea; Sorella, Sandro

doi:10.1063/5.0005037

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 (WFs), capable of describing several materials with very high accuracy, even when standard mean-field approaches [e.g., density functional theory (DFT)] fail. The electronic 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 power with spin-singlet pairing or as a Pfaffian, including both singlet and triplet correlations. This WF can be viewed as an efficient implementation of the so-called resonating valence bond (RVB) Ansatz, first proposed by Pauling and Anderson in quantum chemistry [L. Pauling, The Nature of the Chemical Bond (Cornell University Press, 1960)] and condensed matter physics [P.W. Anderson, Mat. Res. Bull 8, 153 (1973)], 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-of-the-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, which is also an ideal environment for exploiting the computational power of modern Graphics Processing Unit accelerators.

I. INTRODUCTION

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 ago, 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¹ and materials informatics.² 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³ 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 non-interacting electron system, defined by an effective mean-field potential, to be determined self-consistently.⁴ While DFT is an exact theory in principle, the exact form for the exchange-correlation 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,⁵ partly because there is no established strategy for systematic improvement,⁶ 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 wave function (WF)-based approaches popular in quantum chemistry applications.⁹ In these methods, electronic correlations are captured either variationally or perturbatively by post-Hartree–Fock theories, such as the Møller–Plesset perturbation theory (MP),¹⁰ configuration interaction (CI) and Full-CI (FCI),¹¹ multi-configurational self-consistent field (MCSCF),¹² and Coupled-Cluster (CC) theory,¹³ to name a few. 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 a reasonable balance between accuracy and computational affordability (despite its cost growing 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, which are 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,¹⁶ stochastic approaches to numerical algorithms had a pervasive influence in a wide range of fields, from physics and engineering to finance. In this respect, the QMC framework is qualitatively different from the deterministic one 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).¹⁴ 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 correlated 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,¹⁵ which, 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, can be determined by an accurate variational optimization.

So far, several groups have implemented the above algorithms and established excellent QMC codes such as QMCPACK,¹⁷ CASINO,¹⁸ QWALK,¹⁹ CHAMP,²⁰ and HANDE-QMC.²¹ We also remark that other QMC algorithms have also been developed recently, such as the auxiliary field quantum Monte Carlo (AF-QMC),^22,23 full-configuration interaction quantum Monte Carlo (FCI-QMC),^24,25 coupled cluster Monte Carlo (CCMC),^26,27 density matrix quantum Monte Carlo (DMQMC),^28,29 model space quantum Monte Carlo (MSQMC),^30–32 clock quantum Monte Carlo (CQMC),³³ or driven-dissipative quantum Monte Carlo (DDQMC).³⁴

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 appropriate variational wave functions, 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, first proposed by Pauling and Anderson in quantum chemistry³⁵ and condensed matter physics,³⁶ respectively. It includes statical and dynamical correlation effects beyond the commonly used Slater determinant, while keeping the computational cost at the single-determinant 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.³⁷

FIG. 1.

View large Download slide

Schematic picture of some of the possible tools available in TurboRVB. The arrows indicate typical workflow calculations.

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 dynamics 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; and 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 review of QMC¹⁴ for details. TurboRVB also implements an original DFT engine to generate trial WFs, which is 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 as

⟨E⟩ = \frac{\int d x Ψ^{2} (x) \cdot \hat{H} Ψ (x) / Ψ (x)}{\int d x Ψ^{2} (x)} = \int d x e_{L} (x) π (x),

(1)

where $x = (r_{1} σ_{1}, r_{2} σ_{2}, \dots, r_{N} σ_{N})$ here and henceforth is a shorthand notation for the N electron coordinates and their spins, whereas

e_{L} (x) \equiv \frac{\hat{H} Ψ (x)}{Ψ (x)} and π (x) \equiv \frac{Ψ^{2} (x)}{\int d x^{'} Ψ^{2} (x^{'})}

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^38,39) Metropolis method and by averaging the obtained local energies $e_{L} (x_{i})$ ⁠,

E_{VMC} = {⟨e_{L} (x)⟩}_{π (x)} \approx \frac{1}{M} \sum_{i = 1}^{M} e_{L} (x_{i}),

(2)

which has an associated statistical error of $\sqrt{Var [e_{L} (x_{i})] / \tilde{M}}$ ⁠, where Var[e_L(x_i)] is the variance of the sampled local energies and $\tilde{M}$ 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 is worth to note that if Ψ(x) is an eigenfunction of $\hat{H}$ ⁠, say with eigenvalue E₀, then e_L(x) = E₀ for each x, implying that the variance of the local energy is zero and E_VMC = E₀ 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) = Ψ_{G}^{2} (x) / \int d x Ψ_{G}^{2} (x)$ and estimate a generic local observable $O (x)$ either by using $Ō_{VMC} = {⟨O (x)⟩}_{π (x)}$ or

Ō_{VMC} = \frac{{⟨O (x^{'}) W (x^{'})⟩}_{π^{'} (x^{'})}}{{⟨W (x^{'})⟩}_{π^{'} (x)}} \approx \frac{\sum_{i = 1}^{M^{'}} O ({x^{'}}_{i}) W ({x^{'}}_{i})}{\sum_{i = 1}^{M^{'}} W ({x^{'}}_{i})},

(3)

where $W (x) = | Ψ (x) / Ψ_{G} (x) |^{2}$ ⁠, and the points x′_i are distributed according to π′. This reweighting scheme is very important when evaluating atomic forces, as discussed in Sec. VI B. Since the evaluations of the standard deviations are nontrivial in this case, TurboRVB employs the bootstrap and jackknife methods in order to estimate the mean value and the statistical error,¹⁵ 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.⁴⁰

One can optimize the WF based on the variational theorem by introducing a set of parameters $(α_{1}, α_{2}, \dots, α_{p})$ to the WF $Ψ (x, α)$ ⁠,

E_{VMC} (α) = \int d x e_{L} (x, α) π (x, α) ⩾ E_{exact} .

(4)

However, the optimization of a many-body WF remains a difficult challenge not only because optimizing a cost function containing many 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} = - \frac{\partial E}{\partial α_{i}}$ ⁠)⁴¹ and very efficient optimization methods, the so-called “stochastic reconfiguration”^42,43 and “the linear method.”^44–46 These methods are discussed in Secs. VII A and VII B.

B. Lattice regularized Monte Carlo

Lattice regularized diffusion Monte Carlo (LRDMC), which was initially proposed by Casula et al.,⁴⁷ is a projection technique that allows us to improve a variational Ansatz systematically. This method is based on Green’s function Monte Carlo (GFMC),^48–50 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

|Ψ_{T}⟩ = \sum_{n} a_{n} |ϒ_{n}⟩,

(5)

where a_n is the coefficient for the nth eigenvectors (Υ_n). Therefore, by applying $G^{M} = {(Λ - \hat{H})}^{M}$ ⁠, one can obtain

\begin{matrix} |ϒ_{0}⟩ & \propto & lim_{M \to \infty} {(Λ - \hat{H})}^{M} |Ψ_{T}⟩ \\ = & lim_{M \to \infty} {(λ - E_{0})}^{M} [a_{0} |ϒ_{0}⟩ + \sum_{n \neq 0} {(\frac{λ - E_{n}}{λ - E_{0}})}^{M} a_{n} |ϒ_{n}⟩], \end{matrix}

(6)

where Λ is a diagonal matrix with Λ_x′,x = λδ_x′,x (λ should be sufficiently large to obtain the ground state) and E_n is the nth eigenvalue of $\hat{H}$ ⁠. Since $\frac{λ - E_{n}}{λ - E_{0}} < 1$ ⁠, the projection filters out the ground state WF Υ₀ from a given trial WF $|Ψ_{T}⟩$ ⁠, as long as the trial WF is not orthogonal to the true ground state (i.e., a₀ ≡ ⟨Υ₀|Ψ_T⟩ ≠ 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 ${\hat{H}}^{a}$ is then defined such that ${\hat{H}}^{a} \to \hat{H}$ for a → 0. Namely, the kinetic part is approximated by a finite difference form

\begin{matrix} Δ_{i} f (x_{i}, y_{i}, z_{i}) & \approx & Δ_{i}^{a} f (x_{i}, y_{i}, z_{i}) \\ \equiv & \frac{1}{a^{2}} \{[f (x_{i} + a) - f (x_{i})] + [f (x_{i} - a) - f (x_{i})]\} \\ \leftrightarrow & y_{i} \leftrightarrow z_{i}, \end{matrix}

(7)

and the potential term is modified as

V^{a} (x) = V (x) + \frac{1}{2} [\frac{\sum_{i} (Δ_{i}^{a} - Δ_{i}) Ψ_{G} (x)}{Ψ_{G} (x)}] .

(8)

The corresponding Green’s function matrix elements are

G_{x^{'}, x} = ⟨x^{'} |(Λ - {\hat{H}}^{a})| x⟩ \equiv {(Λ - H^{a})}_{x^{'}, x}

(9)

and the single LRDMC iteration step is given by the following equation:

Ψ_{n + 1} (x^{'}) = \sum_{x} G_{x^{'}, x} Ψ_{n} (x) .

(10)

The sketch of the LRDMC algorithm, a Markov chain that evolves the many-body WF according to Eq. (10), is as follows:⁴⁷ (STEP 1) Prepare a walker with configuration x and weight w (w₀ = 1). (STEP 2) A new configuration x′ is generated by the transition probability,

p_{x^{'}, x} = G_{x^{'}, x} / b_{x},

(11)

where

b_{x} = \sum_{x^{'}} G_{x^{'}, x}

(12)

is a normalization factor. By applying the discretized Hamiltonian to a given configuration (⁠ ${\hat{H}}^{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.⁵¹ This allows the evaluation of the normalization factor b_x in Eq. (12) even in a continuous model. Note that 6N comes from the diffusion of each electron in two directions (±a) and the remaining 1 stands for the starting configuration x before the possible hoppings (all N electrons) (i.e., x′ = x). (STEP 3) Finally, update the weight with w_n+1 = w_nb_x and return to STEP I. After a sufficiently large number of iterations (the Markov process is equilibrated), one can calculate the ground state energy E₀,

E_{0} \approx \frac{⟨w_{n} e_{L} (x_{n})⟩}{⟨w_{n}⟩},

(13)

where $⟨\dots⟩$ denotes the statistical average over many independent samples generated by the Markov chain and $e_{L} (x)$ is called the (bare) local energy that reads

e_{L} (x) = \sum_{x^{'}} H_{x^{'}, x} = λ - b_{x} .

(14)

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

G_{n}^{p} = \prod_{j = 1}^{p} b_{n - j} .

(15)

Then, the ground state energy E₀ can be estimated by

E_{0} \approx \frac{\sum_{n} G_{n}^{p} e_{L} (x_{n})}{\sum_{n} G_{n}^{p}} .

(16)

This straightforward implementation of the above simple method is not suitable for realistic simulations due to fluctuations of weights, large correlation times, 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,¹⁵ in which the original Green’s function is modified using the so-called guiding function Ψ_G as

{\tilde{G}}_{x^{'}, x} = G_{x^{'}, x} \frac{Ψ_{G} (x^{'})}{Ψ_{G} (x)},

(17)

and the projection is modified as

Ψ_{G} (x^{'}) Ψ_{n + 1} (x^{'}) = \sum_{x} {\tilde{G}}_{x^{'}, x} Ψ_{G} (x) Ψ_{n} (x) .

(18)

In practice, the guiding function is prepared by a VMC calculation. The modified Green’s function for importance sampling ${\tilde{G}}_{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

{\tilde{b}}_{x} = \sum_{x^{'}} {\tilde{G}}_{x^{'}, x},

(19)

and the local energy with importance sampling is

{\tilde{e}}_{L} (x) = \frac{⟨ Ψ_{G} | \hat{H} | x ⟩}{⟨ Ψ_{G} | x ⟩} = \sum_{x^{'}} H_{x^{'}, x} \frac{Ψ_{G} (x^{'})}{Ψ_{G} (x)} .

(20)

Equation (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 in 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)].⁵¹ 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¹⁵ in order to avoid the sign problem. Indeed, the Hamiltonian is modified using the spin-flip term $V_{s f} (x) = \sum_{x^{'} : s_{x^{'}, x} > 0} H_{x^{'}, x} Ψ_{G} (x^{'}) / Ψ_{G} (x)$ ⁠,

H_{x^{'}, x}^{FN, γ} = \{\begin{matrix} H_{x, x} + (1 + γ) V_{S F} (x) f o r x^{'} = x \\ H_{x^{'}, x} f o r x^{'} \neq x, s_{x^{'}, x} < 0 \\ - γ H_{x^{'}, x} f o r x^{'} \neq x, s_{x^{'}, x} > 0, \end{matrix})

(21)

where $s_{x^{'}, x} = Ψ_{G} (x^{'}) H_{x^{'}, x} / Ψ_{G} (x)$ and γ ≥ 0 is a real parameter. The use of the fixed-node Green’s function

{\tilde{G}}_{x^{'}, x}^{FN} = {(Λ - H^{FN})}_{x^{'}, x} \frac{Ψ_{G} (x^{'})}{Ψ_{G} (x)}

(22)

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⁴⁹ in TurboRVB) scheme for a more efficient computation.¹⁵ The code performs the branching as follows: (1) Set the new weights equal to the average of the old ones,

{w^{'}}_{α, n} = \bar{w} \equiv \frac{1}{N_{w}} \sum_{β} w_{β, n} .

(23)

(2) Select the new walkers from the old ones with a probability that is proportional to the old walkers’ weights,

p_{α, n} = \frac{w_{α, n}}{\sum_{β} w_{β, n}},

(24)

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 τ, which 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,

G_{n}^{p} = \prod_{j = 1}^{p} {\bar{w}}_{n - j},

(25)

and calculate the ground state energy,

E_{0} \approx \frac{\sum_{n} G_{n}^{p} e_{L} (x_{n})}{\sum_{n} G_{n}^{p}},

(26)

where $e_{L} (x_{n})$ is the mean local energy averaged over the walkers, which reads

e_{L} (x_{n}) = \frac{\sum_{α} w_{α, n} e_{L} (x_{α, n})}{\sum_{α} w_{α, n}}

(27)

and is evaluated just before each reconfiguration. Note that p is also an input parameter, which 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 $\frac{M}{Λ} = τ$ fixed. In the M → ∞ limit, the projection ${(Λ - \hat{H})}^{M}$ is equal to the imaginary time evolution $\exp (- τ \hat{H})$ ⁠, apart from an irrelevant constant Λ^M. Thus, the user should specify only τ as the input parameter. Indeed, the walker weight is updated by $w \to w \exp (- τ_{ξ} e_{L} (x))$ and the imaginary time is updated by τ_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 after each projection time of length τ 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.¹⁵ 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.⁴⁷ 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. Second, 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,^37,47 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,³⁷ 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 recently been reported.³⁷ Finally, the inclusion of non-local pseudopotentials in this framework is straightforward by means of an additional spherical grid defined in an appropriate mesh.⁵¹ As described in Refs. 47 and 52, LRDMC provides an upper bound for the true ground-state energy and allows the estimation of E_FN, even in the presence of non-local pseudopotentials. Note that this variational property has also been extended to the standard DMC framework,⁵³ with the introduction of the so-called T-moves. Moreover, the recently introduced determinant locality approximation (DLA)⁵⁴ to deal with non-local pseudopotentials is also implemented in TurboRVB and can be optionally used in LRDMC.

FIG. 2.

View large Download slide

Fixed-node energies for the all-electron carbon atom computed within DMC, single-grid LRDMC (one lattice), and double-grid LRDMC. The lattice discretization parameter a is mapped to the time step τ as $a = \sqrt{τ}$ ⁠. Reprinted with permission from M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett. 95, 100201 (2005). Copyright 2005 APS.

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 trade-off 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 Sec. II, we have seen that QMC algorithms, both at the variational and fixed-node levels, 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 system size.

TurboRVB employs a many-body WF Ansatz Ψ, which can be written as the product of two terms,

Ψ = Φ_{AS} \times \exp J,

(28)

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 inter-electronic distances, is applied to Φ_SD, the resulting Ansatz Ψ_JSD = Φ_SD * exp J often provides over 70% of the correlation energy⁵⁵ 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.⁵⁶ A Jastrow correlated SD (JSD) function yields a computational cost for QMC simulations—both VMC and FN—about ∝N³, 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,⁵⁷ which is a remapping of the electronic configurations that enter 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 levels. However, it can be used at a price to increase significantly an already large computational cost.⁵⁸ Another possibility is to improve Ψ_SD similar 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 pre-defined 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 a constrained number of molecular orbitals (Pfn), (iii) the Antisymmetrized Geminal Power (AGP), (iv) the Antisymmetrized Geminal Power with a constrained number of molecular orbitals (AGPn), and (v) the single Slater determinant. It is interesting to observe that the latter four WFs are all 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 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 Sec. VII. In this section, we will describe 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).

FIG. 3.

View large Download slide

Ansatz hierarchy. The output of Hartree–Fock (HF) or DFT simulations with different exchange-correlation functionals are special instances of the SD Ansatz.

FIG. 4.

View large Download slide

Ansatz conversion.

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 conditions.⁵⁹ TurboRVB implements the Jastrow term composed of one-body, two-body, and three/four-body factors (J = J₁ + J₂ + 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 electron pair contribution. The one-body Jastrow factor J₁ is the sum of two parts, the homogeneous part (enforcing the electron–ion cusp condition),

J_{1}^{h} (r_{1}, \dots, r_{N}) = \sum_{i = 1}^{N} \sum_{a = 1}^{N_{at}} (- {(2 Z_{a})}^{3 / 4} u_{a} ({(2 Z_{a})}^{1 / 4} |r_{i} - R_{a}|)),

(29)

and the corresponding inhomogeneous part,

J_{1}^{i n h} (r_{1} σ_{1}, \dots, r_{N} σ_{N}) = \sum_{i = 1}^{N} \sum_{a = 1}^{N_{at}} (\sum_{l} M_{a, l}^{σ_{i}} χ_{a, l} (r_{i})),

(30)

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_{a, l}^{σ_{i}}}$ are variational parameters, and $u_{a} (r)$ is a simple bounded function. In TurboRVB, the most common choice for u_a is

u_{a} (r) = \frac{1}{2 b_{e a}} (1 - e^{- r b_{e a}}),

(31)

depending on a single variational parameter b_ea, which may be optimized independently for each atomic species.

The two-body Jastrow factor is defined as

J_{2} (r_{1} σ_{1}, \dots, r_{N} σ_{N}) = \sum_{i < j} v_{σ_{i}, σ_{j}} (|r_{i} - r_{j}|),

(32)

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:

v_{σ_{i}, σ_{j}} (r_{i, j}) = \{\begin{matrix} \frac{r_{i, j}}{4} \cdot {(1 + b_{e e}^{p a r a} \cdot r_{i, j})}^{- 1} & (σ_{i} = σ_{j}) \\ \frac{r_{i, j}}{2} \cdot {(1 + b_{e e}^{a n t i} \cdot r_{i, j})}^{- 1} & (σ_{i} \neq σ_{j}), \end{matrix})

(33)

where $r_{i, j} = |r_{i} - r_{j}|$ ⁠, and $b_{e e}^{p a r a}$ and $b_{e e}^{a n t i}$ are variational parameters.

The three/four-body Jastrow factor reads

J_{3 / 4} (r_{1} σ_{1}, \dots, r_{N} σ_{N}) = \sum_{i < j} (\sum_{a, l} \sum_{b, m} M_{a, l, b, m}^{σ_{i}, σ_{j}} χ_{a, l} (r_{i}) χ_{b, m} (r_{j})),

(34)

where the indices l and m again indicate different orbitals centered on corresponding atoms a and b, and ${M_{a, l, b, m}^{σ_{i}, σ_{j}}}$ 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.^60–62 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,

i \leftrightarrow (r_{i}, σ_{i}),

(35)

corresponding to the ith 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

Φ_{AS} (1, \dots, N) = A \{g (1, 2) g (3, 4) \dots g (N - 1, N)\},

(36)

where $A$ is the antisymmetrization operator,

A \equiv \frac{1}{N!} \sum_{P \in S_{N}} ϵ_{P} \hat{P},

(37)

S_N is the permutation group of N elements, $\hat{P}$ is the operator corresponding to the generic permutation P, and ϵ_P is its sign.

Let us define G the N × N matrix with elements G_i,j = g(i, j). Note that

g (i, j) = - g (j, i) (and G_{i, j} = - G_{j, i}),

(38)

as a consequence of the statistics of fermionic particles; thus, G is skew-symmetric (i.e., G^T = −G, ^T being the transpose operator), so the diagonal is zero and the number of independent entries is $\sum_{n = 1}^{N - 1} n = N (N - 1) / 2$ ⁠.

The Pfaffian⁶³ of a N × N skew-symmetric matrix G is defined as

Pf (G) \equiv \frac{1}{2^{N / 2} (N / 2)!} \sum_{P \in S_{N}} ϵ_{P} G_{P (1), P (2)} \dots G_{P (N - 1), P (N)}

(39)

if N is even, and it is zero if N is odd.⁶⁴ Therefore, the WF Φ_AS defined in the right-hand side (RHS) of Eq. (36) equals $\frac{1}{N!!} Pf (G)$ ⁠, where the semifactorial $N!! \equiv \frac{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

Φ_{Pf} = Pf (G) .

(40)

Note 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

(41)

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_{d u} = - {G_{u d}}^{T}$ ⁠, i.e., g_du(i, j) = −g_ud(j, i). g_uu describes the pairing between a pair of electrons with spin-up,

g_{u u} (i, j) = f_{u u} (r_{i}, r_{j}) |↑ ↑⟩,

(42)

where the function f_uu describes the spatial dependence on the coordinates r_i, r_j for i, j ≤ N_u. Note 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,

g_{d d} (i, j) = f_{d d} (r_{i}, r_{j}) |↓ ↓⟩

(43)

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⟩ = \frac{|↑ ↓⟩ - |↓ ↑⟩}{\sqrt{2}}$ or a triplet $|1,0⟩ = \frac{|↑ ↓⟩ + |↓ ↑⟩}{\sqrt{2}}$ ⁠, in the general case, the pairing function g_ud will be a linear combination of the two components,

g_{u d} (i, j) = f_{S} (r_{i}, r_{j}) \frac{|↑ ↓⟩ - |↓ ↑⟩}{\sqrt{2}} + f_{T} (r_{i}, r_{j}) \frac{|↑ ↓⟩ + |↓ ↑⟩}{\sqrt{2}},

(44)

where f_S(r_i, r_j) = f_S(r_j, r_i) describes the spatial dependence of the singlet part of g_ud and f_T(r_i, r_j) = −f_T(r_j, r_i) describes the spatial dependence of the triplet part. Therefore, the generic pairing function g(i, j) is the sum of all four components mentioned above, namely,

\begin{matrix} g (i, j) & = & f_{S} (r_{i}, r_{j}) |0,0⟩ + f_{T} (r_{i}, r_{j}) |1,0⟩ \\ + f_{u u} (r_{i}, r_{j}) |1, + 1⟩ + f_{d d} (r_{i}, r_{j}) |1, - 1⟩ . \end{matrix}

(45)

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

f (r_{i}, r_{j}) = \sum_{l, m, a, b} A_{\{a, l\}, \{b, m\}} ψ_{a, l} (r_{i}) ψ_{b, m} (r_{j}),

(46)

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 $A_{\{a, l\}, \{b, m\}} = A_{\{b, m\}, \{a, l\}}$ because f_S(r_i, r_j) = f_S(r_j, r_i), whereas $A_{\{a, l\}, \{b, m\}} = - A_{\{b, m\}, \{a, l\}}$ for f_T, f_uu, and f_dd.

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 N physical electrons. The extra matrix elements are easily computed,

\begin{matrix} g (N + 1, N + 1) = 0, \\ g (i, N + 1) = - g (N + 1, i) = ϕ (i) f o r i \leq N, \end{matrix}

(47)

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

Φ_{AS} (1, \dots, N) = A \{g (1, 2) g (3, 4) \dots g (N, N + 1)\},

(48)

where the antisymmetrization operator acts on the (now even) N + 1 electrons. 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 and (ii) it does not depend on the extra N + 1th 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

Φ_{Pf} = Pf (\tilde{G}),

(49)

where

(50)

is a (N + 1) × (N + 1) skew-symmetric matrix and Φ is a N-dimensional vector whose elements are ϕ(i), i = 1, …, N. Its Pfaffian is quite generally nonvanishing for at least some configuration, provided that 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 electrons 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 electrons coupled by the matrix G 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 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

(51)

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-mentioned derivation, it is possible to generalize further the Pfaffian WF Ansatz by introducing fictitious particles and relaxing the constraints on the $\tilde{G}$ matrix, even in the lower block diagonal part, which 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 AGP).

Let us consider first an unpolarized system, having an even number N of 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 Sec. III B are both zero matrices of size N/2 × N/2, and the matrix G_ud has only the contribution coming from the singlet, which we dub G_S. The antisymmetrization operator implies the computation of

(52)

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,

Φ_{AGPs} = det (G_{S}) .

(53)

It should be noted that it is not necessary that the matrix G_ud be 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

Φ_{AGP} = det (G_{u d})

(54)

is not an eigenstate of the spin. In other terms, there is a spin contamination, similar 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 with 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

A \{g (1, 2) \dots g (2 p - 1, 2 p) ϕ_{1} (2 p + 1) \dots ϕ_{k} (N)\} .

(55)

As mentioned above, we assume that the pairing function is zero for like-spin 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 was done at the end of Subsection III B, 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.,

(56)

with the N_u × N_u matrix $\tilde{G} = [G_{u d} | Φ]$ ⁠, 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

Φ_{AGP} = det (\tilde{G}) .

(57)

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^65–67 (see Sec. 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 Casula and Sorella⁶⁵ for the first time in 2003; then, it has also been implemented in other QMC codes.

D. AGP and Pf with a constrained number of molecular orbitals (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 of the 2L atomic orbitals used, L for each spin component. Therefore,

g (i, j) = \sum_{μ, ν}^{2 L} A_{μ, ν} ψ_{μ} (i) ψ_{ν} (j),

(58)

and clearly, the symmetry of g implies that A_μ,ν = −A_ν,μ. The coefficients A_μ,ν define a 2L × 2L skew-symmetric matrix A. If we define the 2L dimensional vector $Ψ_{i} = {(ψ_{1} (i), \dots, ψ_{2 L} (i))}^{T}$ ⁠, Eq. (58) rewrites as $g (i, j) = Ψ_{i}^{T} A Ψ_{j}$ ⁠. Moreover, the overlaps S_μ,ν ≡ ⟨ψ_μ|ψ_ν⟩ between atomic orbitals define the overlap matrix S, which in the case of a spin-dependent basis is block diagonal,

(59)

with S_uu and S_dd being 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,

{\tilde{ψ}}_{μ} (i) = \sum_{ν} S_{μ, ν}^{- 1 / 2} ψ_{ν} (i),

(60)

where S^−1/2 is well defined since S is strictly positive definite.⁶⁸ The matrix elements of g can be recasted in the new orthonormal basis $\{{\tilde{ψ}}_{μ} (i)\}$ ⁠, yielding $g (i, j) = \sum_{μ, ν}^{2 L} Ã_{μ, ν} {\tilde{ψ}}_{μ} (i) {\tilde{ψ}}_{ν} (j) = {\tilde{Ψ}}_{i}^{T} Ã {\tilde{Ψ}}_{j}$ ⁠, with the matrix Ã ≡ S^1/2AS^1/2 that is skew symmetric.

At this point, from the spectral theory of skew symmetric matrices, it is possible to perform the Youla decomposition of Ã (see Ref. 69), 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.⁷⁰ So, the pairing function g(i, j) can be written as $Φ_{i}^{T} Σ Φ_{j}$ ⁠, where $Φ_{i} = Q^{T} {\tilde{Ψ}}_{i}$ for each i. This defines a basis of L MOs $\{ϕ_{k} (i)\}$ and corresponding twinned ones $\{{\bar{ϕ}}_{k} (i)\}$ ⁠, forming together a basis of 2L mutually orthonormal elements for which the original geminal function reads

g (i, j) = \sum_{k = 1}^{L} λ_{k} [ϕ_{k} (i) {\bar{ϕ}}_{k} (j) - {\bar{ϕ}}_{k} (i) ϕ_{k} (j)]

(61)

with λ_k ≥ 0. After these transformations, these MOs can be finally written in the chosen (hybrid) atomic basis,

\begin{matrix} ϕ_{k} (i) & = & \sum_{ν = 1}^{2 L} P_{μ, l} ψ_{μ} (i), \\ {\bar{ϕ}}_{k} (i) & = & \sum_{ν = 1}^{2 L} {\bar{P}}_{μ, l} ψ_{μ} (i) \end{matrix}

(62)

by appropriate p × 2L rectangular matrices P and $\bar{P}$ ⁠. Then, with no loss of generality, we can assume that the molecular orbitals $\{ϕ_{k} (i), {\bar{ϕ}}_{k} (i)\}$ are ranked such that λ₁ ≥ λ₂ ≥…≥ λ_L ≥ 0. 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 orbitals with k > n, yielding the pairing function

g_{n} (i, j) = \sum_{k = 1}^{n} λ_{k} [ϕ_{k} (i) {\bar{ϕ}}_{k} (j) - {\bar{ϕ}}_{k} (i) ϕ_{k} (j)],

(63)

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⁷¹ 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 increasing the variational energy expectation value, energy differences may actually improve by an appropriate choice of n. In the mentioned work,⁷¹ this promising approach was applied with an AGP containing only singlet correlations, but the binding energies were defined without using a rigorous size consistent criterion. This drawback can now be 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 corresponding matrix $A_{μ, ν}^{n} = \sum_{k = 1}^{n} λ_{k} [P_{μ, k} {\bar{P}}_{ν, k} - P_{ν, k} {\bar{P}}_{μ, k}]$ ⁠. According to Eq. (63), an arbitrary small variation δg_n of the constrained pairing function g_n reads

\begin{matrix} δ g_{n} (i, j) & = & \sum_{k = 1}^{n} δ λ_{k} [ϕ_{k} (i) {\bar{ϕ}}_{k} (j) - {\bar{ϕ}}_{k} (i) ϕ_{k} (j)] \\ + \sum_{k = 1}^{n} λ_{k} [δ ϕ_{k} (i) {\bar{ϕ}}_{k} (j) - δ {\bar{ϕ}}_{k} (i) ϕ_{k} (j)] \\ + \sum_{k = 1}^{n} λ_{k} [ϕ_{k} (i) δ {\bar{ϕ}}_{k} (j) - {\bar{ϕ}}_{k} (i) δ ϕ_{k} (j)] \end{matrix}

(64)

and therefore satisfies the following property, as will be shown later,

(Î - L) δ g_{n} (Î - R) = 0,

(65)

where Î is the identity operator, L and R are projection operators over the occupied MOs, i.e., L²(i, j) = ∫ dkL(i, k)L(k, j) = L(i, j), and similarly, R² = R, where here and henceforth the shorthand integration symbol $\int d k = \sum_{σ_{k}} \int d r_{k}$ contains implicitly also the spin summation. These operators are then defined as follows:

\begin{matrix} L (i, j) & = & \sum_{k = 1}^{n} [ϕ_{k} (i) ϕ_{k}^{*} (j) + {\bar{ϕ}}_{k} (i) {\bar{ϕ}}_{k}^{*} (j)], \\ R (i, j) & = & \sum_{k = 1}^{n} [ϕ_{k}^{*} (i) ϕ_{k} (j) + {\bar{ϕ}}_{k}^{*} (i) {\bar{ϕ}}_{k} (j)] . \end{matrix}

(66)

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. Note 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:

δ g_{n} = δ g - (Î - L) δ g (Î - R) .

(67)

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)² = (Î − R) and (Î − L)² = (Î − L), from which Eq. (67) fulfills Eq. (65). Equation (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, $\bar{P}$ ⁠, and the overlap matrix S that can be easily and efficiently implemented.¹⁵

This linear relation between A and Aⁿ can be therefore easily implemented together with the corresponding derivatives necessary for the optimization of the energy⁷² 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.¹⁵

E. Single Slater determinant (SD)

An important special case of the AGPn and Pfn Ansätze discussed in Sec. III D 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,71 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 Secs. III B and III C, it has been mentioned that the Pfaffian and the AGP Ansätze have a multiconfigurational character despite the fact that 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 real 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), {\bar{ϕ}}_{k} (i)\}$ appearing in the RHS of Eq. (61) have the same spatial part, which we denote ${\tilde{ϕ}}_{k} (r)$ ⁠, modulus a sign (because they are orthonormal), and they have an opposite spin part. Without the loss of generality, we can assume that ϕ_k is spin up and ${\bar{ϕ}}_{k}$ is spin down. Given this convention, the scalar product l_k between the spatial parts of ϕ_k and ${\bar{ϕ}}_{k}$ will either be +1 or −1. We define ${\tilde{λ}}_{k} = l_{k} λ_{k}$ ⁠, where λ_k is the same as the one in the RHS of Eq. (61). Note that $|{\tilde{λ}}_{k}| = λ_{k}$ ⁠, so $\{\tilde{λ}\}$ are ranked in decreasing order of their absolute value (i.e., $|{\tilde{λ}}_{k}| \geq |{\tilde{λ}}_{k + 1}|$ ⁠). The pairing function can then be written as

g (i, j) = \sum_{k = 1}^{L} {\tilde{λ}}_{k} {\tilde{ϕ}}_{k} (r_{i}) {\tilde{ϕ}}_{k} (r_{j}) (|↑ ↓⟩ - |↓ ↑⟩) .

(68)

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 ${\tilde{ϕ}}_{k}$ and satisfies the canonical anticommutation relations. We can write the pairing functions $g (i, j) \equiv ⟨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

ĝ = \sum_{k = 1}^{L} {\tilde{λ}}_{k} â_{k, ↑}^{†} â_{k, ↓}^{†},

(69)

where ${\tilde{λ}}_{k}$ is the one defined above. Using this notation, it is easy to show (see Appendix of Ref. 67) that the pairing function g is equivalent to the complete active space of two electrons on the L molecular orbitals $\{{\tilde{ϕ}}_{k}\}$ ⁠.

Within this notation, the AGP WF is

| Φ_{AGPs} (⟩ = ĝ^{p} | 0 (⟩,

(70)

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 operator ${\hat{b}}_{k} \equiv â_{k, ↑}^{†} â_{k, ↓}^{†}$ ⁠, which created an electron pair on the orbital ${\tilde{ϕ}}_{k}$ ⁠, and having noted that ${\hat{b}}_{k} {\hat{b}}_{l} = {\hat{b}}_{l} {\hat{b}}_{k}$ and ${\hat{b}}_{k}^{2} = 0$ (as following from the anticommutation relations of the $â_{k}^{†}$ operators), we obtain that

ĝ^{p} = p! \sum_{1 \leq i_{1} < i_{2} < \dots < i_{p} \leq L} {\tilde{λ}}_{i_{1}} {\tilde{λ}}_{i_{2}} \dots {\tilde{λ}}_{i_{p}} {\hat{b}}_{i_{1}} {\hat{b}}_{i_{2}} \dots {\hat{b}}_{i_{p}} .

(71)

The chosen order of the ${\tilde{λ}}_{k}$ coefficients implies that the leading term in ĝ^p is given by the term with (i₁, …, i_p) = (1, …, p).⁷³ Thus, we can rewrite ĝ^p in terms of this leading term and excitations with respect to it,

\begin{matrix} ĝ^{p} \frac{|0⟩}{p! \prod_{i = 1}^{p} {\tilde{λ}}_{i}} & = & \prod_{i = 1}^{p} {\hat{b}}_{i} | 0 (⟩ + \sum_{j = 1}^{p} \sum_{q = p + 1}^{L} \frac{{\tilde{λ}}_{q}}{{\tilde{λ}}_{j}} (\underset{i \neq j}{\prod_{i = 1}^{p}} {\hat{b}}_{i}) {\hat{b}}_{q} | 0 (⟩ \\ + \sum_{1 \leq j < k \leq p} \sum_{p < q < r \leq L} \frac{{\tilde{λ}}_{q} {\tilde{λ}}_{r}}{{\tilde{λ}}_{j} {\tilde{λ}}_{k}} (\underset{i \neq k \land i \neq j}{\prod_{i = 1}^{p}} {\hat{b}}_{i}) {\hat{b}}_{q} {\hat{b}}_{r} | 0 (⟩ + \dots . \end{matrix}

(72)

In the above expansion, we shall recognize that the first element is a single Slater determinant of the molecular orbitals ${\{{\tilde{ϕ}}_{i}\}}_{i = 1, \dots, N / 2}$ ⁠; the second element is the summation of all possible double excitations 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⁷⁴ subset of the CI expansion, having some constraints on the coefficients of the expansion.

It is worth mentioning that AGPs are 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 highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) are degenerate (see Fig. 5). In fact, we can have that $|{\tilde{λ}}_{p}| = |{\tilde{λ}}_{p + 1}|$ ⁠, so the two determinants $\prod_{i = 1}^{p} {\hat{b}}_{i} |0⟩$ and $(\prod_{i = 1}^{p - 1} {\hat{b}}_{i}) {\hat{b}}_{p + 1} |0⟩$ have the same weight, and they both are the leading terms of the constrained zero-seniority expansion in Eq. (72).

FIG. 5.

View large Download slide

The twisted ethylene molecule (C₂H₄) is a prototypical system requiring a multireference method. The ground state geometry is planar (torsional angle equal to zero); C atoms have a double bond, and the HOMO and LUMO orbitals are separated. As the molecule is twisted, one bond is broken and the HOMO–LUMO orbitals become degenerate for a torsional angle of 90°. The upper panel shows evaluations of the torsional barrier as obtained using JSD and JAGP at the VMC and LRDMC level and using a complete active space (CAS) of 12 orbitals in 12 electrons. JAGP provides reliable results, whereas JSD (which neglects the LUMO orbital) overestimates the barrier both at the variational and diffusion levels. The middle panel shows the ratio ${\tilde{λ}}_{LUMO} / {\tilde{λ}}_{HOMO}$ in the JAGP Ansatz and the corresponding value in the CAS calculations. The bottom panel shows the HOMO and LUMO orbitals for the planar and twisted geometries. Readapted with permission from Zen et al., J. Chem. Theory Comput. 10, 1048 (2014). Copyright 2014 ACS.

Note that in the AGPn Ansatz, we can perform an expansion similar to Eq. (72), but the excitations stop 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,

ψ_{l, \pm m, I}^{Gaussian} (r; ζ) = {|r - R_{I}|}^{l} e^{- ζ {|r - R_{I}|}^{2}} \cdot R [{(- i)}^{\frac{1 \pm 1}{2}} Y_{l, m, I} (θ, φ)],

(73)

or the Slater type,

ψ_{l, \pm m, I}^{Slater} (r; ζ) = {|r - R_{I}|}^{l} e^{- ζ |r - R_{I}|} \cdot R [{(- i)}^{\frac{1 \pm 1}{2}} Y_{l, m, I} (θ, φ)],

(74)

where the real and the imaginary part (m > 0) of the spherical harmonic function $Y_{l, m, I} (θ, φ)$ centered at R_I are 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 one-body and three/four-body Jastrow parts [i.e., denoted as $χ (r)$ in Eqs. (30) and (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⁷⁵ or from other more specific references when using pseudopotential.^76–81 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 many-body techniques. Recently, some progress has been made in this direction, and pseudopotentials determined by correlated many-body techniques are also available.^76–83 All the pseudopotentials used in QMC employ the standard semi-local form

{\hat{V}}_{pp}^{I} (r_{i}) = V_{loc}^{I} (r_{i, I}) + \sum_{l = 0}^{l_{max}} V_{l}^{I} (r_{i, I}) \sum_{m = - l}^{l} |Y_{l, m}⟩ ⟨Y_{l, m}|,

(75)

where r_i,I = |r_i − R_I| is the distance between the ith electron and the Ith ions, l_max is the maximum angular momentum of the ion I, and $\sum_{l = 0}^{l_{max}} \sum_{m = - l}^{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.⁸⁴ As it is now becoming a common practice not only in QMC, both the local $V_{loc}^{I} (r_{i, I})$ and the non-local $V_{l}^{I} (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

r^{2} V_{l} (r) = \sum_{k} α_{k, l} r^{β_{k, l}} \exp (- γ_{k, l} r^{2}),

(76)

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 TurboRVB. 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 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.”⁸⁵ To this purpose, let us decompose the geminal function in terms of atomic contribution, as far as the dependence over i (the left argument) is concerned,

g (i, j) = \sum_{I} U_{proj}^{I} (i, j) = \sum_{I} \sum_{μ, ν} A_{μ, ν}^{I} ϕ_{μ} (i) ϕ_{ν} (j),

(77)

where I represents an atom in a system; $U_{proj}^{I} (i, j)$ is the pairing function projected on the atom I; $A_{μ, ν}^{I}$ is A_μ,ν if μ ∈ I and $A_{μ, ν}^{I} = 0$ otherwise, 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),

A_{μ, ν} = \sum_{k = 1}^{n} (c_{μ}^{k} {\bar{c}}_{ν}^{k} - {\bar{c}}_{μ}^{k} c_{ν}^{k}),

(78)

where $c_{μ}^{k}$ and ${\bar{c}}_{μ}^{k}$ are the coefficients of the DFT molecular orbitals $ϕ_{k} (i) = \sum_{μ} c_{μ}^{k} ϕ_{μ} (i)$ and ${\bar{ϕ}}_{k} (i) = \sum_{ν} {\bar{c}}_{ν}^{k} ϕ_{ν} (i)$ in the atomic basis expansion, respectively.

Quite generally, the projected pairing function can be expanded in a truncated space spanned by q terms,

Ũ_{proj}^{I} (i, j) = \sum_{k = 1}^{q} σ_{k} ψ_{k}^{GEO} (i) ψ_{k}^{Env} (j) .

(79)

In other words, the Schmidt decomposition is applied to the matrix $Ũ_{proj}^{I} (i, j)$ describing the coupling between a given atom I and the environment, within the geminal Ansatz. This procedure defines the so-called geminal embedded orbitals (GEOs), which are determined in terms of an expansion over all atomic orbitals used for the atom I,

ψ_{k}^{GEO} (i) = \sum_{μ \in I} μ_{k, μ}^{GEO} ψ_{μ} (i),

(80)

where $ψ_{k}^{GEO} (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,

d^{2} = {|U_{proj}^{I} - Ũ_{proj}^{I}|}^{2} .

(81)

Considering all possible unconstrained functions $ψ_{k}^{Env} (j)$ and employing the steady condition $\frac{δ d^{2}}{δ ψ_{k}^{Env} (j)} = 0$ ⁠, d² reads

d^{2} = {|U_{proj}^{I}|}^{2} - \sum_{k = 1}^{p} \int d i d j D_{proj}^{I} (i, j) ψ_{k}^{GEO} (i) ψ_{k}^{GEO} (j),

(82)

where ${|U_{proj}^{I}|}^{2} = \int d i d j {|U_{proj}^{I} (i, j)|}^{2}$ and $D_{proj}^{I} (i, j)$ is the density matrix that reads

D_{proj}^{I} (i, j) = \int d k U_{proj}^{I} {(i, k)}^{*} U_{proj}^{I} (j, k) .

(83)

Since the GEOs $ψ_{k}^{GEO} (r)$ are orthonormal and the density matrix is Hermitian, Eq. (82) becomes minimum when the GEOs coincide 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,

\sum_{j \in J} {[(A^{I *} S A^{I}) S]}_{i, j} μ_{k, j}^{GEO} = ω_{k} μ_{k, i}^{GEO} .

(84)

The truncation error is readily estimated by the summation of the eigenvalues $d^{2} = {|U_{proj}^{I}|}^{2} - \sum_{i = 1}^{p} 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 will be shown below. With this construction, a new geminal is defined in the GEO basis, namely,

\tilde{g} (i, j) = \sum_{μ, ν} Ã_{μ, ν} ψ_{μ}^{GEO} (i) ψ_{ν}^{GEO} (j),

(85)

where the matrix coefficients Ã_μ,ν are given by maximizing the normalized overlap (Q) between the original and the new geminals,

Q = \frac{{⟨ \tilde{g} | g ⟩}^{2}}{⟨ \tilde{g} | \tilde{g} ⟩ ⟨ g | g ⟩},

(86)

where $⟨ \tilde{g} | g ⟩ = \int d i d j \tilde{g} {(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² to 4p² ≪ 4L².

J. Conversion of the WF

As described in Sec. III, TurboRVB implements different types of Ansätze: (i) the Pfaffian (Pf), (ii) the Pfaffian with a constrained number of molecular orbitals (Pfn), (iii) the Antisymmetrized Geminal Power (AGP), (iv) the Antisymmetrized Geminal Power with a 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 (see 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 Secs. III D and III E). The corresponding tool is convertfort10mol.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:

max Q = \frac{{⟨ g^{new} | g^{ori} ⟩}^{2}}{⟨ g^{new} | g^{new} ⟩ ⟨ g^{ori} | g^{ori} ⟩},

(87)

in order to obtain new geminal matrix coefficients $A_{μ, ν}^{new}$ ⁠, defining the new pairing function as

g^{new} (i, j) = \sum_{μ \cdot ν} A_{μ \cdot ν}^{new} ψ_{μ}^{new} (i) ψ_{ν}^{new} (j),

(88)

while the original geminal was given in terms of the parameter matrix $A_{μ, ν}^{ori}$ ⁠,

g^{ori} (i, j) = \sum_{μ, ν} A_{μ, ν}^{ori} ψ_{μ}^{ori} (i) ψ_{ν}^{ori} (j) .

(89)

Note that 0 ⩽ Q ⩽ 1; therefore, the larger the Q, the better 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 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 Sec. 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:^62,69: (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 for 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

| ↑ (⟩ \to \frac{1}{\sqrt{2}} (| ↑ (⟩ + | ↓ (⟩) and | ↓ (⟩ \to \frac{1}{\sqrt{2}} (| ↑ (⟩ - | ↓ (⟩) .

(90)

One can convert from an AGP the pairing function that is obtained from a VMC optimization

g_{u d} (i, j) = f_{S} (r_{i}, r_{j}) \frac{|↑ ↓⟩ - |↓ ↑⟩}{\sqrt{2}} + f_{T} (r_{i}, r_{j}) \frac{|↑ ↓⟩ + |↓ ↑⟩}{\sqrt{2}}

(91)

to a Pf one

\begin{matrix} g_{u d} (i, j) \to g (i, j) = f_{S} (r_{i}, r_{j}) \frac{|↑ ↓⟩ - |↓ ↑⟩}{\sqrt{2}} + f_{T} (r_{i}, r_{j}) (|↑ ↑⟩ - |↓ ↓⟩) . \end{matrix}

(92)

This transformation provides a meaningful initialization to the Pfaffian WF that can then be 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^86,87 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,^88,89

Ψ_{k_{s}} (r_{1}, \dots, r_{i} + T_{s}, \dots, r_{N}) = e^{i k_{s} \cdot T_{s}} Ψ_{k_{s}} (r_{1}, \dots, r_{i}, \dots, r_{N}),

(93)

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:

ψ_{l, m, I}^{PBC} (r + T_{s}; ζ) = e^{i k_{s} \cdot T_{s}} ψ_{l, m, I}^{PBC} (r; ζ),

(94)

where k_s is a twist vector [ $k_{s} = (k_{s}^{x}, k_{s}^{y}, k_{s}^{z})$ ] and T_s represents an arbitrary simulation cell vector. Note that the use of a non-vanishing twist vector generally makes a many-body WF complex. TurboRVB implements the CRYSTAL periodic basis,^15,90,91

ψ_{l, m, I}^{PBC} (r; ζ) = \sum_{T_{s}} ψ_{l, m, I} (r + T_{s}; ζ) e^{- i k_{s} \cdot T_{s}},

(95)

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. Note that, in TurboRVB, T_s are not limited to orthorhombic ones, but has recently been extended to include all possible crystal translation groups (e.g., rhombohedral, hexagonal triclinic).

The same procedure is applied to the basis set for the Jastrow part, although using simple periodic boundary conditions,¹⁵ because the twists do not affect the Jastrow part of the WF, namely,

χ_{l, m, I}^{PBC} (r; ζ) = \sum_{T_{s}} χ_{l, m, I} (r + T_{s}; ζ),

(96)

which satisfies $χ_{l, m, I}^{PBC} (r + T_{s}; ζ) = χ_{l, m, I}^{PBC} (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

{\tilde{J}}_{1} (r + T_{s}) = {\tilde{J}}_{1} (r),

(97)

and the two-body Jastrow part [Eq. (32)] should fulfill

\begin{matrix} J_{2} (r_{1} σ_{1}, \dots, (r_{i} + T_{s}) σ_{i}, \dots, r_{N} σ_{N}) = J_{2} (r_{1} σ_{1}, \dots, r_{i} σ_{i}, \dots, r_{N} σ_{N}) . \end{matrix}

(98)

In order to satisfy the above constraints, we consider the relative electron–nuclei or electron–electron coordinate differences r_d, necessary to evaluate J₁ and J₂, respectively, and expand them as

r_{d} = r_{a} a + r_{b} b + r_{c} c,

(99)

where r_a, r_b, and r_c are appropriate transformed coordinates, which 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}) \to ({\bar{r}}_{a}, {\bar{r}}_{b}, {\bar{r}}_{c}) = (p (r_{a}), p (r_{b}), p (r_{c}))$ by use of an appropriate function p(x), with at least a 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 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:

p (x) = \{\begin{matrix} x & (- \frac{1}{4} < x < \frac{1}{4}) \\ - \frac{1}{8 (1 + 2 x)} & (- \frac{1}{2} ⩽ x ⩽ - \frac{1}{4}) \\ \frac{1}{8 (1 - 2 x)} & (\frac{1}{4} ⩽ x ⩽ \frac{1}{2}) . \end{matrix})

(100)

Indeed, although the modified relative distance diverges (i.e., $|r_{d}| \to \infty$ ⁠) at the edges of the Wigner–Seitz cell (e.g., $r = \pm \frac{1}{2} a, \pm \frac{1}{2} b, \pm \frac{1}{2} c$ ⁠), the exponential [Eq. (31)] and the Padé [Eq. (33)] functions remain finite, $u (r) \to 1 / 2 b_{e a}$ and $v_{σ_{i}, σ_{j}} (r_{i, j}) \to 1 / 4 b_{e e}^{p a r a}$ or $1 / 2 b_{e e}^{a n t i}$ ⁠, 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.,

p (x) = \{\begin{matrix} x & (- \frac{1}{6} < x < \frac{1}{6}) \\ - \frac{1}{54 {(1 / 2 + x)}^{2}} & (- \frac{1}{2} ⩽ x ⩽ - \frac{1}{6}) \\ \frac{1}{54 {(1 / 2 - x)}^{2}} & (\frac{1}{6} ⩽ x ⩽ \frac{1}{2}) . \end{matrix})

(101)

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,^88,89 namely,

Ψ_{k_{p}} (\{r_{i} + T_{p}\}) = e^{i k_{p} \cdot T_{p}} Ψ_{k_{p}} (\{r_{i}\}),

(102)

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, which has proven very effective for implementing the mentioned translation symmetries within pairing WFs.⁹²

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 is the so-called two-body effect that arises from the periodic Ewald contribution resulting from the electron–electron interaction. Note 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,⁸⁶ special k-points methods,^88,89 or the exact special twist (EST) method,⁸⁷ all of which have been implemented in TurboRVB. Within the TurboRVB implementation, the Jastrow part is independent of twists (i.e., TurboRVB uses a common Jastrow for all twists),

Ψ_{k_{s}} (r_{1}, \dots, r_{N}) = J (r_{1}, \dots, r_{N}) Φ_{A S}^{k_{s}} (r_{1}, \dots, r_{N}) .

(103)

As emphasized in Ref. 86, 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 (MPC) interaction.^93–95 This method has not been implemented in TurboRVB yet. Nevertheless, simpler alternatives exist. For instance, one can alleviate the two-body finite size effects by directly increasing the supercell size, or one can estimate these effects by employing the KZK exchange-correlation function⁹⁶ 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,⁹⁷ which can be readily computed within TurboRVB, 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,⁹⁸ CRYSTAL,⁹¹ and Quantum Espresso,⁹⁹ TurboRVB is a self-consistent 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)¹⁰⁰ and those combined with KZK exchange-correlation energy⁹⁶ 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 TurboRVB, an original and systematic algorithm to remove this redundancy¹⁰¹ 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,

H c = E S c,

(104)

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,

S_{i, j} = ⟨ ϕ_{i} | ϕ_{j} ⟩ = \int d r ϕ_{i}^{*} (r) ϕ_{j} (r) .

(105)

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₁), 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:

e_{j}^{i} = \frac{1}{\sqrt{s_{i}}} υ_{j}^{i} for s_{i} / s_{N} ⩾ ε_{mach},

(106)

where s_i and $υ_{j}^{i}$ are the ith 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,

{\tilde{e}}_{j}^{i} = \frac{1}{\sqrt{{\tilde{s}}_{i}}} {\tilde{υ}}_{j}^{i},

(107)

where ${\tilde{s}}_{i}$ and ${\tilde{υ}}_{j}^{i}$ are the ith eigenvalues and eigenvectors of the recomputed overlap matrix ${\tilde{s}}_{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,

U_{i, j} = \sum_{k} e_{j}^{k} {\tilde{e}}_{k}^{i} .

(108)

On this basis set, the generalized eigenvalue problem becomes a well-conditioned one with a corresponding M_b × M_b Hamiltonian matrix $\tilde{H} = U H \tilde{U}$ ⁠. In this way, one can avoid the numerical instabilities induced by a too redundant large basis set. This stabilization introduces an error at most $\sqrt{ε_{mach} \cdot 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 TurboRVB 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)¹⁰² 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, namely,

{\tilde{ϕ}}_{j}^{b} (r - R_{b}) = ϕ_{j}^{b} (r - R_{b}) {\tilde{J}}_{1} (r),

(109)

where ${\tilde{J}}_{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:

lim_{r \to R_{b}} \frac{\nabla {\tilde{ϕ}}_{j}^{a}}{{\tilde{ϕ}}_{j}^{a}} = - Z_{b} \frac{r - R_{b}}{| r - R_{b} |}

(110)

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, which 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 Sec. V B, TurboRVB employs standard atomic orbitals, modified by means of an appropriately 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 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,¹⁰³ 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 to 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, is very complicated, mainly for the complexity of the algorithm that, in turn, may lead to a very inefficient implementation, although recent progress has been done.¹⁰⁴ For instance, a simple approach is to compute them with finite difference expressions, leading to a too large computational time, as 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 many variational parameters. The AAD was applied for the calculation of atomic forces in Ref. 41. 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 [ $\frac{\partial}{\partial α_{k}} e_{L} (x_{i})$ ] or those corresponding to the WF logarithm [ $\frac{\partial}{\partial α_{k}} \ln Ψ (x_{i})$ ] using the AAD, which drastically improves the efficiency (Fig. 6)⁴¹ and reliability of the calculation.

FIG. 6.

View large Download slide

Ratio of CPU time required to compute energies and all force components referenced to the one required for the simple energy calculation within VMC. The calculations refer to 1, 2, 4, and 32 water molecules. The inset is an expansion of the lower part of the plot. Reprinted with permission from S. Sorella and L. Capriotti, J. Chem. Phys. 133, 234111 (2010). Copyright 2010 AIP Publishing LLC.

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,

f_{k} = - \frac{\partial E (α)}{\partial α_{k}} = - \frac{\partial}{\partial α_{k}} \frac{⟨ Ψ_{α} | \hat{H} | Ψ_{α} ⟩}{⟨ Ψ_{α} | Ψ_{α} ⟩} .

(111)

In variational Monte Carlo, the derivative can be evaluated using M configurations of electron coordinates,¹⁵

\begin{matrix} f_{k} & = & - 2 R [\sum_{x} \frac{e_{L}^{*} (x) (O_{k} (x) - Ō_{k}) {|Ψ_{α} (x)|}^{2}}{\sum_{x} {|Ψ_{α} (x)|}^{2}}] \\ \approx & - 2 R [\frac{1}{M} \sum_{i = 1}^{M} e_{L}^{*} (x_{i}) (O_{k} (x_{i}) - Ō_{k})], \end{matrix}

(112)

where $e_{L} (x)$ is the local energy, $O_{k} (x)$ is the logarithmic derivative of the WF [i.e., $O_{k} (x) = \frac{\partial \ln Ψ_{α} (x)}{\partial α_{k}}$ ], and $Ō_{k}$ is its average over M samples [i.e., $Ō_{k} = \frac{1}{M} \sum_{i = 1}^{M} O_{k} (x_{i})$ ]. In TurboRVB, the logarithmic derivative [ $O_{k} (x)$ ] is computed very efficiently by using the AAD algorithm. Note 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 positions). 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 therefore represent the fundamental ingredients for an efficient WF optimization, as described in 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,

devmax \equiv max_{k} (|\frac{f_{k}}{σ_{f_{k}}}|),

(113)

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) minimum is being approached.

B. Atomic forces

In TurboRVB, a very useful space warp coordinate transformation (SWCT)⁴¹ is employed for the calculation of atomic forces. This scheme was introduced several decades ago in Ref. 105, in order to decrease the statistical errors of forces. The finite-difference expression for ionic force calculations is

F_{a} = - \frac{Δ}{Δ R_{a}} E = - \frac{E (R_{a} + Δ R_{a}) - E (R_{a})}{Δ R_{a}} + O (Δ R_{a}) .

(114)

The SWCT is used to mimic the displacement of charges around the nucleus,

{\bar{r}}_{i} = r_{i} + Δ R_{a} ω_{a} (r_{i}),

(115)

ω_{a} (r) = \frac{κ (|r - R_{a}|)}{\sum_{b = 1}^{M} κ (|r - R_{b}|)},

(116)

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¹⁰⁵ simple choice $κ (r) = 1 / r^{4}$ works very well and is indeed adopted in TurboRVB.

Starting from the finite-difference expression, one can straightforwardly derive the corresponding differential expression,⁴¹

\begin{matrix} F_{a} = - ⟨\frac{d}{d R_{a}} e_{L}⟩ + 2 (⟨e_{L}⟩ ⟨\frac{d}{d R_{a}} \ln (J^{\frac{1}{2}} Ψ)⟩ - ⟨e_{L} \frac{d}{d R_{a}} \ln (J^{\frac{1}{2}} Ψ)⟩), \end{matrix}

(117)

where $J$ is the Jacobian of the above transformation and the brackets indicate 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,⁴¹

\frac{d}{d R_{a}} e_{L} = \frac{\partial}{\partial R_{a}} e_{L} + \sum_{i} ω_{a} (r_{i}) \frac{\partial}{\partial r_{i}} e_{L},

(118)

\begin{matrix} \frac{d}{d R_{a}} \ln (J^{\frac{1}{2}} Ψ) = \frac{\partial}{\partial R_{a}} \ln (Ψ) + \sum_{i} [ω_{a} (r_{i}) \frac{\partial}{\partial r_{i}} \ln (Ψ) + \frac{1}{2} \frac{\partial}{\partial r_{i}} ω_{a} (r_{i})] . \end{matrix}

(119)

In order to evaluate these differential expressions, 6N + 6N_at components have to be evaluated, namely, $⟨\frac{\partial}{\partial r_{i}} e_{L}⟩$ ⁠, $⟨\frac{\partial}{\partial r_{i}} \ln Ψ⟩$ ⁠, $⟨\frac{\partial}{\partial R_{a}} \ln e_{L}⟩$ ⁠, and $⟨\frac{\partial}{\partial R_{a}} \ln Ψ⟩$ ⁠.⁴¹ These values are very efficiently computed in TurboRVB, by using the aforementioned AAD, which works even in the presence of pseudopotentials.

The SWCT significantly decreases the statistical errors, but the forces still have infinite variance properties because ∂e_L and $\partial \ln (Ψ)$ diverge in the vicinity of the nodal surfaces. Attaccalite and Sorella¹⁰⁶ 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) = Ψ_{T}^{2} (x)$ but by a slightly different guiding function $Ψ_{G} (x)$ defined by

Ψ_{G} (x) = \frac{R^{ε} (x)}{R (x)} Ψ_{T} (x),

(120)

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 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 $P f (G) ≃ \frac{1}{G_{i, j}^{- 1}}$ ⁠, where G_i,j is the pairing function [Eq. (45)], here referred to as 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

R (x) = {(S \sum_{i, j} {|{G_{i, j}}^{- 1}|}^{2})}^{- θ_{R}}

(121)

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.¹⁰⁶ Here, we have also introduced 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 vanishes because it is corresponding to an electron going to infinity (all matrix elements containing such an electron coordinate have to vanish),

S = {m i n}_{i} \sum_{j} | G_{i j} |^{2} .

(122)

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 is also for bulk systems, and it is particularly important when, in these cases, there exist large regions of almost negligible electronic density.¹⁰⁷

After that, $R^{ε} (x)$ is defined as

R^{ε} (x) = m a x [R (x), ε],

(123)

which is much simpler than the original proposal reported in Ref. 106. By using the new probability, forces (the HF and the Pulay) can be evaluated with finite variance as

{⟨W (x) F (x)⟩}_{Π_{G} (x)} / {⟨W (x)⟩}_{Π_{G} (x)},

(124)

where $W (x)$ is the new weight,

W (x) = {(Ψ_{T} (x) / Ψ_{G} (x))}^{2} \equiv {(\frac{R (x)}{m a x [R (x), ε]})}^{2}

(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 $≃ \frac{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., $\int d δ δ^{4 θ_{R} - 2} < \infty$ ⁠. 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,¹⁰⁶ 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$ ⁠, which 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 TurboRVB 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.

VII. OPTIMIZATION OF WFs

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 updated as follows:

α_{k} \to α_{k}^{'} = α_{k} + δ α_{k},

(126)

δ α_{k} = Δ f_{k},

(127)

where Δ is a small constant and $f_{k} \equiv - \frac{\partial E}{\partial α_{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 in 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, and the quasi-Newton method, but the straightforward implementation of these optimizations does not work efficiently within a stochastic approach such as 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,

α_{k} \to α_{k} + Δ \cdot {(S^{- 1} f)}_{k},

(128)

where the matrix $S$ is stochastically evaluated by means of M configuration samples $x = \{x_{1}, x_{2}, \dots, x_{M}\}$ ⁠,

S_{k, k^{'}} = [\frac{1}{M} \sum_{i = 1}^{M} {(O_{k} (x_{i}) - Ō_{k})}^{*} (O_{k^{'}} (x_{i}) - Ō_{k^{'}})],

(129)

where $O_{k} (x_{i}) = \frac{\partial \ln Ψ (x_{i})}{\partial α_{k}}$ and $Ō_{k} = \frac{1}{M} \sum_{i = 1}^{M} O_{k} (x_{i})$ ⁠. The resulting approach is the so-called stochastic reconfiguration (SR) method.⁴² Mazzola et al.¹⁰⁸ 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.¹⁰⁹ 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.⁴³ Therefore, in practice, the diagonal elements of the preconditioning matrix $S$ are shifted by a small positive parameter (ε) as

s_{i, i}^{'} = s_{i \cdot i} (1 + ε) .

(130)

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

α_{k} \to α_{k} + Δ \cdot {({S^{'}}^{- 1} f)}_{k} .

(131)

In practice, the user should provide two input parameters, namely, Δ (denoted as tpar in the code input variables) and ε (denoted as parr).

FIG. 7.

View large Download slide

Optimization of the variational WF in the simple one-dimensional Heisenberg model $H = J \sum_{i} {\vec{S}}_{i} \cdot {\vec{S}}_{i + 1}$ with the standard SR (ε = 0, open circles) and with the present regularization (ε = 0.001, open triangles). The evolution of the nearest neighbor spin–spin (S_z) Jastrow parameter is plotted. The figure clearly shows that the SR method with regularization is several orders of magnitude more efficient than the standard SR for determining the variational parameter with a given statistical accuracy. The inset shows the first few iterations. Reprinted with permission from S. Sorella, M. Casula, and D. Rocca, J. Chem. Phys. 127, 014105 (2007). Copyright 2007 AIP Publishing LLC.

B. Linear method

In TurboRVB, the state-of-the-art QMC optimization method is also implemented, namely, the so-called linear method.^44–46 In this scheme, a many-body WF is expanded up the linear order (i.e., considering only the first derivatives for the variational parameters),

| Ψ_{α + δ α} (⟩ = | Ψ_{α} (⟩ + \sum_{k = 1} δ α_{k} \frac{\partial}{\partial α_{k}} | Ψ_{α} (⟩ \equiv | Ψ_{0} (⟩ + \sum_{k = 1} δ α_{k} | Ψ_{k} (⟩,

(132)

where α_k is the kth variational parameter and $Ψ_{k} (x)$ is the first derivative with respect to the kth variational parameter. Within this expansion, the expectation value of the energy reads

E_{α + z} = \frac{\sum_{k, k^{'} = 0}^{p} {z_{k}}^{*} z_{k^{'}} ⟨ Ψ_{k} | \hat{H} | Ψ_{k^{'}} ⟩}{\sum_{k, k^{'} = 0}^{p} {z_{k}}^{*} z_{k^{'}} ⟨ Ψ_{k} | Ψ_{k^{'}} ⟩} .

(133)

In order to obtain the vector z that minimizes this expectation value, one should solve the generalized eigenvalue problem,

H z = E S z,

(134)

where $H_{k, k^{'}} = ⟨ Ψ_{k} | \hat{H} | Ψ_{k^{'}} ⟩$ and $S_{k, k^{'}} = ⟨ Ψ_{k} | Ψ_{k^{'}} ⟩$ ⁠.

This algorithm is more conveniently implemented in the so-called semi-orthogonal basis,¹⁵

| Ψ_{α + δ α} (⟩ = z_{0} | {\tilde{Ψ}}_{α} (⟩ + \sum_{k = 1} z_{k} (Ô_{k} - Ō_{k}) | {\tilde{Ψ}}_{α} (⟩,

(135)

where $| {\tilde{Ψ}}_{α} (⟩ = | Ψ_{α} (⟩ / ∥ Ψ_{α} ∥$ ⁠, $⟨ x | Ô_{k} | x^{'} ⟩ = δ_{x, x^{'}} O_{k} (x)$ ⁠, $O_{k} (x) = \frac{\partial}{\partial α_{k}} \ln |Ψ_{α} (x)|$ ⁠, and $Ō$ _k = ⟨Ψ_α| $Ô$ _k|Ψ_α⟩/⟨Ψ_α|Ψ_α⟩. TurboRVB calculates z by solving the generalized equation [Eq. (134)] with the matrices,

H_{k, k^{'}} \approx \frac{1}{M} \sum_{i = 1}^{M} {(O_{k} (x_{i}) - Ō_{k})}^{*} \frac{⟨ x_{i} | \hat{H} (\hat{O_{k^{'}}} - Ō_{k^{'}}) | Ψ_{α} ⟩}{⟨ x_{i} | Ψ_{α} ⟩},

(136)

S_{k, k^{'}} \approx \frac{1}{M} \sum_{i = 1}^{M} {(O_{k} (x_{i}) - Ō_{k})}^{*} (O_{k^{'}} (x_{i}) - Ō_{k^{'}}),

(137)

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 $\partial_{α_{k}} e_{L} (x_{i})$ ⁠, which can be directly computed by AAD, and thus allow the calculation of the above Hamiltonian matrix elements, by straightforward algebra,

\partial_{α_{k}} e_{L} (x_{i}) = \partial_{α_{k}} \frac{⟨ x_{i} | \hat{H} | Ψ_{α} ⟩}{⟨ x_{i} | Ψ_{α} ⟩} = \frac{⟨ x_{i} | \hat{H} (\hat{O_{k^{'}}} - Ō_{k^{'}}) | Ψ_{α} ⟩}{⟨ x_{i} | Ψ_{α} ⟩} .

(138)

Note 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

α_{k} \to α_{k} + Δ \cdot z_{k} / z_{0},

(139)

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 many parameters, and also, in such cases, each iteration requires a diagonalization of a huge matrix (a task that can be often prohibitive). In TurboRVB, 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 the 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,⁴³ optionally neglecting the ones with a smaller signal-noise ratio (<parcutmin).

C. Practical rule

In most optimizations, the number of samplings in VMC is much larger than the number of variational parameters and the optimization is stable and efficient. In the opposite case, the preconditioning matrix S becomes rank-deficient singular; therefore, the optimization does not work properly (cf. the Cramér–Rao inequality). A practical rule is to set the number of samples M such that¹⁵

M ⩾ 5 \sim 10 \times p,

(140)

where p is the number of variational parameters. However, one can also employ a 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 TurboRVB without storing explicitly the matrix¹¹⁰ 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 F = −∇_RV(R), where F ≡ {F₁, …, F_N}, R ≡ {R₁, …, R_N}, and the potential energy landscape V(R) is evaluated by VMC, namely,

V (R) = \frac{⟨ Ψ_{R} | \hat{H} (R) | Ψ_{R} ⟩}{⟨ Ψ_{R} | Ψ_{R} ⟩} .

(141)

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 TurboRVB (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 first-order 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¹⁰⁶

\dot{υ} = - γ (R) \cdot υ + F (R) + η (t),

(142)

\dot{R} = υ,

(143)

⟨η_{i} (t) η_{j} (t^{'})⟩ = S_{i, j} (R) δ (t - t^{'}),

(144)

where R, υ, f, and η are the 3N_at-dimensional vectors representing atomic positions, velocities, and deterministic and stochastic forces of N_at atoms, written in mass-scaled coordinates,

\begin{aligned} R_{i} & = R_{i}^{0} \sqrt{m_{i}}, \\ F_{i} & = F_{i}^{0} / \sqrt{m_{i}}, \\ η_{i} & = η_{i}^{0} \sqrt{m_{i}} \end{aligned}

(145)

for i = 1, …, N_at. The stochastic forces are related to the friction matrix γ through the fluctuation–dissipation theorem, namely,

S (R) = 2 T γ (R),

(146)

with the temperature T expressed in atomic units. TurboRVB exploits the freedom in Eq. (146), by assuming

S (R) = s_{0} I + Δ_{0} S^{QMC} (R),

(147)

where s₀ and Δ₀ are constants, $I$ is the identity, and $S^{QMC} (R)$ is the covariance matrix of QMC forces,

S_{i, j}^{QMC} (R) = ⟨(F_{i} (R) - ⟨F_{i} (R)⟩)⟩ ⟨(F_{j} (R) - ⟨F_{j} (R)⟩)⟩ .

(148)

In the above equation, $⟨\dots⟩$ refers to the average over the QMC sampling at the given MD step. The friction, a quantity that controls the nuclear sampling efficiency, is, therefore, position-dependent now. Luo et al. have shown that the force covariance matrix $S_{i, j}^{QMC} (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} τ} \cdot υ_{n} + Γ_{n} (F_{n} + \tilde{η}),

(149)

R_{n + 1} = R_{n} + τ υ_{n + 1},

(150)

with $F_{n} \equiv F (R_{n})$ ⁠, $γ_{n} \equiv γ (R_{n})$ ⁠,

Γ_{n} = γ_{n}^{- 1} (I - e^{- γ_{n} τ}),

(151)

and

\tilde{η} = \frac{γ_{n}}{2 \sinh (γ_{n} τ / 2)} \int_{t_{n} - τ / 2}^{t_{n} + τ / 2} e^{- γ_{n} (t - t_{n})} \cdot η (t) d t .

(152)

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:

⟨{\tilde{η}}_{i} {\tilde{η}}_{j}⟩ \equiv S^{*} = T γ_{n}^{2} \coth (γ_{n} τ / 2) .

(153)

In this way, $\tilde{η}$ fulfill the fluctuation–dissipation theorem required by the Langevin thermostat. However, part of this noise is already present in the QMC force evaluation (⁠ ${\tilde{η}}^{QMC}$ ⁠), as quantified by Eq. (148). Hence, TurboRVB adds the external random force ${\tilde{η}}^{ext}$ such that $\tilde{η} = {\tilde{η}}^{ext} + {\tilde{η}}^{QMC}$ ⁠. This is generated according to the following correlator:

⟨{\tilde{η}}_{i}^{ext} {\tilde{η}}_{j}^{ext}⟩ = S^{*} - S_{QMC},

(154)

in order to fulfill the original fluctuation–dissipation theorem. Since the correlation matrix in Eq. (154) is positive definite as long as τ < Δ₀, the correlated external force can be readily generated by standard algorithms (e.g., the Box–Muller method¹¹² combined with the matrix diagonalization). This noise correction plays an important role in targeting a stable temperature even in the presence of statistical errors in the QMC forces. The original numerical integration detailed here has been further improved by Mazzola and Sorella (see the supplementary notes of Ref. 113).

FIG. 8.

View large Download slide

Eigenvectors of the 3 × 3 correlation matrix of QMC forces $S^{QMC}$ in the water monomer. The obtained eigenvectors correspond to the well-known three vibrational modes: bending (red arrows), symmetrical (blue arrows), and asymmetrical (green arrows) stretching. The smaller eigenvalue corresponds to the lowest-frequency vibrational mode. The eigenvalues in the plot are all rescaled by the lowest one. Reprinted with permission from Y. Luo, A. Zen, and S. Sorella, J. Chem. Phys. 141, 194112 (2014). Copyright 2014 AIP Publishing LLC.

B. Path integral Ornstein–Uhlenbeck dynamics

Mouhat et al. have recently introduced the nuclear quantum effects into the second-order LD based on a path integral molecular dynamics (PIMD) approach,¹¹⁴ driven by QMC nuclear forces. As a starting point of the Feynman path integral (PI) theory, the quantum partition function $Z = Tr [e^{- β H}]$ (with β = 1/k_BT) can be written as¹¹⁵

Z ≃ \frac{1}{{(2 π ℏ)}^{3 N_{a t} L}} \int d^{N_{a t} L} p d^{N_{a t} L} q e^{- τ_{β} H_{L} (p, q)},

(155)

with τ_β = β/L denoting 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 the 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 \equiv {(R^{(1)}, \dots, R^{(i)}, \dots, R^{(L)})}^{T}$ is 3N_atL-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

\begin{matrix} H_{L} (p, q) & = & \sum_{i = 1}^{N_{a t}} \sum_{j = 1}^{L} (\frac{1}{2} {[p_{i}^{(j)}]}^{2} + \frac{1}{2} {\tilde{ω}}_{L}^{2} {(R_{i}^{(j)} - R_{i}^{(j - 1)})}^{2}) \\ + \sum_{j = 1}^{L} V (R_{1}^{(j)}, \dots, R_{N_{a t}}^{(j)}), \end{matrix}

(156)

written in mass-scaled coordinates [Eq. (145)], with ${\tilde{ω}}_{L} = L / β ℏ$ ⁠, and subjected to the ring boundary condition $R_{i}^{(0)} \equiv R_{i}^{(L)}$ ⁠. The system is thermalized with a Langevin thermostat such that the related Liouville operator $L$ can be written as

L = \sum_{i = 1}^{3 N_{a t} L} (F_{i} \partial_{p_{i}} + p_{i} \partial_{q_{i}} - \sum_{j = 1}^{3 N_{a t} L} γ_{i j} (\partial_{p_{i}} p_{j} + k_{B} T M \partial_{p_{i}} \partial_{p_{j}}))

(157)

in mass-scaled coordinates, with $\partial_{q_{i}} \equiv \frac{\partial}{\partial_{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} \equiv F_{i}^{B O} + F_{i}^{h a r m}$ is the total force acting on each replica, comprising the BO (intra-replica) and harmonic contributions (inter-replicas), where $F_{i}^{B O} \equiv - \partial_{q_{i}} \tilde{V}$ and $\tilde{V} = \sum_{j = 1}^{L} V (R_{1}^{(j)}, \dots, R_{N_{a t}}^{(j)})$ ⁠.

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 and the other depending exclusively on the harmonic modes of the rings. To do so, we first separate the friction matrix into two contributions,

γ = γ^{B O} + γ^{h a r m} .

(158)

We can then rewrite the total Liouvillian as the sum of two terms, $L = L^{B O} + L^{h a r m}$ ⁠, where

L^{h a r m} = \sum_{i} (F_{i}^{h a r m} \partial_{p_{i}} + p_{i} \partial_{q_{i}} - \sum_{j} γ_{i j}^{h a r m} (\partial_{p_{i}} p_{j} + k_{B} T M \partial_{p_{i}} \partial_{p_{j}})),

(159)

L^{B O} = \sum_{i} (F_{i}^{B O} \partial_{p_{i}} - \sum_{j} γ_{i j}^{B O} (\partial_{p_{i}} p_{j} + k_{B} T M \partial_{p_{i}} \partial_{p_{j}})),

(160)

in such a way that we can break up the evolution operator via a Trotter factorization¹¹⁵ to get the following symmetric propagator:^116,117

e^{- L δ t} ≃ e^{- L^{B O} δ t / 2} e^{- L^{h a r m} δ t} e^{- L^{B O} δ t / 2} .

(161)

The equations of motion corresponding to the propagator in Eq. (161) have been implemented in TurboRVB, where now $F_{i}^{B O}$ can be computed in a QMC framework, and therefore, it can be affected by an intrinsic noise, as was in the classical scheme (Sec. VIII A).

We note that the equations of motion generated by $L^{h a r m}$ 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^{B O}$ factor. The dynamics generated by $L^{h a r m}$ is a quantum Ornstein–Uhlenbeck process, which describes the motion of Brownian quantum particles under the influence of friction.

The exact integration of $L^{h a r m}$ 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.¹¹⁸ 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^{B O}$ factor in Eq. (161), the corresponding equations of motion are the ones in Eq. (12) of Ref. 111. Indeed, only the momenta are evolved in $L^{B O}$ ⁠, and the resulting equations are equal to those of the simple classical Langevin algorithm introduced in Ref. 106, restricted to p, which in the present case is a 3N_atL-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 toward a steady temperature ensemble even for quantum particles.

We refer the interested reader to Ref. 114, 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 Mazzola and Sorella.¹¹⁹ Here, the ions formally undergo the following dynamics:

\dot{R} = S^{- 1} (R) \cdot F (R) + η,

(162)

⟨η_{i} η_{j}⟩ = 2 T δ (t - t^{'}) S_{i, j}^{- 1} (R),

(163)

where T is the temperature, $F (R)$ and η are the deterministic and stochastic forces, and $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 Ref. 119 and reads

\begin{matrix} R (t + τ) & = & R (t) + \sqrt{2 T τ} \cdot z (t) + τ S^{- 1} (R) F (R) \\ - S^{- 1} (R) (\frac{S (R (t - τ)) - S (R)}{2}) (R (t - τ) - R (t)), \end{matrix}

(164)

with

⟨z_{i} (t) z_{j} (t)⟩ = S_{i, j}^{- 1} (R (t)) .

(165)

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,

T = T_{target} (1 - a τ),

(166)

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 is often done in the case of second-order 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.^119,120 Moreover, nuclear quantum effects can also be computed within this first-order LD, in exactly the same way as was implemented for the second-order case discussed in Subsection VIII B. 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.

IX. IMPLEMENTATION AND OPERATIONS

Table I shows the main modules implemented in TurboRVB and a brief description of their functionalities. All programs reported in Table I 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).

TABLE I.

Main modules in TurboRVB.

Module	Description
makefort10.x	Generates a JAGP WF from a given basis set and structure
assembling_pseudo.x	Prepares pseudopotentials
convertfort10mol.x	Adds molecular orbitals. If the number of molecular orbitals is equal to (larger than) half the number of electrons in a system, the resultant WF is the JSD (JAGPn)
convertfort10.x	Converts a wave function type (fort.10_in), i.e., JSD/JAGP/JAGPn/JPf/JPfn, to another one with a different type (fort.10_out). The output (fort.10_new) is as close as possible to the input and may also include an optimal contracted (hybrid) basis set.
convertfortpfaff.x	Converts a JAGP WF to a 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

Module	Description
makefort10.x	Generates a JAGP WF from a given basis set and structure
assembling_pseudo.x	Prepares pseudopotentials
convertfort10mol.x	Adds molecular orbitals. If the number of molecular orbitals is equal to (larger than) half the number of electrons in a system, the resultant WF is the JSD (JAGPn)
convertfort10.x	Converts a wave function type (fort.10_in), i.e., JSD/JAGP/JAGPn/JPf/JPfn, to another one with a different type (fort.10_out). The output (fort.10_new) is as close as possible to the input and may also include an optimal contracted (hybrid) basis set.
convertfortpfaff.x	Converts a JAGP WF to a 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

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^↑.¹²¹ 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 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 unoccupied orbitals will be used in the geminal expansion, up to the nth 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 convertfort10mol.x to a previously determined JAGP Ansatz. Analogously, a JPf WF can be converted from a JAGP WF using convertpfaff.x. JPfn can be obtained by applying convertfort10mol.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₀ 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 many independent random walks (called walkers), initialized with different seeds for random number generators. Figure 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 node without hyperthreading). For the weak scaling case, we set the number of walkers M to the minimum possible one (one walker per core) and increased the number of cores from 68 (one node) to 27 200 (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 Graphics Processing Unit (GPU) accelerators. This feature is currently under development.

FIG. 9.

View large Download slide

(a) The weak and (b) strong scalings of TurboRVB. The benchmarks were measured using the conventional 2 × 2 × 2 diamond with the ccECP pseudopotential (256 electrons in the simulation cell) on Marconi KNL nodes (Cineca supercomputer/SCAI).

XI. POSTPROCESSING

A. Correlated samplings

A very efficient and easy-to-use correlated sampling technique has been implemented in TurboRVB. This is very useful in evaluatingtotal energy differences between two very similar WFs. One of its practical sides is to investigate whether 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.¹⁵ The observables presently available in the readforward.x module can be categorized as follows: (i) Charge density [ $ρ (r) = ρ^{↑} (r) + ρ^{↓} (r)$ ] and spin density [ $ρ^{σ} (r) = (ρ^{↑} (r) - ρ^{↓} (r)) / 2$ ], 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.¹²² Figure 10 shows the charge density of the square four hydrogen (H₄) depicted by 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²) operators inside a sphere centered on each atom,^62,69 whose radius can be specified by the user. (iv) The Berry phase: $z_{α} = ⟨ Ψ | e^{i (2 π / L^{α}) \sum_{i = 1}^{N} r_{i}^{α}} | Ψ ⟩$ ⁠, 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 α.¹²³ This property can be used to characterize the metal–insulator transition based on the many-body WF, as has been done in the one-dimensional hydrogen chain by Stella et al. (Fig. 11).

FIG. 10.

View large Download slide

The xy-plane charge density of the square four hydrogen (H₄), calculated using a JAGP WF. Reprinted with permission from C. Genovese, A. Meninno, and S. Sorella, J. Chem. Phys. 150, 084102 (2019).¹²⁴ Copyright 2019 AIP Publishing LLC.

FIG. 11.

View large Download slide

(a) Electronic localization length λ_N divided by the interatomic distance a as a function of a, for the one-dimensional hydrogen chain, where λ_N is defined as $λ_{N} = (L / 2 π) \sqrt{- \ln {|z_{N}|}^{2} / N}$ ⁠, N is the number of electrons (hydrogen atoms), and L is the length of the simulation cell along the one dimensional chain. (b) Modulus of the complex polarization z_N as a function of the interatomic distance of the hydrogen atoms. According to the previous works,^125,126 in the thermodynamic limit (N → 0), a metal can be characterized by a vanishing modulus of the complex polarization (⁠ $|z_{N}| \to 0$ ⁠, λ_N → +∞), while it becomes unity in the insulating case (⁠ $|z_{N}| \to 1$ ⁠, λ_N → 0). Thus, one can discuss the metal–insulator transition. Reprinted with permission from Stella et al., Phys. Rev. B 84, 245117 (2011). Copyright 2011 APS.

XII. APPLICATIONS

Since the start of the TurboRVB project in 2003 by Casula and Sorella,⁶⁵ the code has been applied to study atomic species, molecular systems, and various materials, including challenging systems such as large complexes, surfaces, and liquids. 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₂,¹²⁷ Be₂,^71,128 B₂,⁷¹ C₂,^62,71 N₂,⁷¹ O₂,^{69,71,127,129} F₂,⁷¹ Na₂,¹⁰³ LiF,⁷¹ CN,⁷¹ and Fe₂¹³⁰), reactive oxygen species (singlet O₂, $O_{2}^{-}$ ⁠, OH^•, OH⁻, NO^•, NO⁻, HOO^•, HOO⁻, and cis and trans HOOO^•),¹²⁹ aromatic molecules (benzene^127,131 and oligoacene series,¹³² see Sec. XII D), the water molecule,^107,127 dimer^133,134 and hexamer clusters,¹³⁴ the Zundel ion (⁠ $H_{5} O_{2}^{+}$ ⁠),^114,135 H₂S, SO₂, NH₃, PH₃,¹¹¹ and others. Most of the studies report both variational and FN diffusion Monte Carlo results, and several of the WF Ansätze available in TurboRVB are tested. Often, it is shown that quite good results are already obtained at the variational level, and the FN DMC further improves the accuracy of the outcomes. For instance, see the case of the sodium dimer in Fig. 12.

FIG. 12.

View large Download slide

PES of the sodium dimer. The dashed lines show the obtained Murrell–Sorbie function. The error bars are within the markers. The obtained LRDMC PES agrees with the CCSD(T) and experimental ones. Reprinted with permission from K. Nakano, R. Maezono, and S. Sorella, J. Chem. Theory Comput. 15, 4044 (2019). Copyright 2019 ACS.

Table II and 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 follows 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 II and Fig. 13 indicate, where the reportedbinding 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.¹³⁶ Let us consider the oxygen molecule as an example: the O₂ ground state is a spin triplet (⁠ $^{3} Σ_{g}^{-}$ ⁠) which, in both JSD and 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 (³P). Neither JSD nor JAGP has the flexibility to describe the O–O system by the same WF parametrization used for O₂;¹³⁷ thus, E_O-O > 2E_O. In contrast, JPf has the flexibility to describe the dissociation limit correctly, so $E_{O - O}^{JPf} = 2 E_{O}^{JPf}$ ⁠. Figure 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. Note 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, and transition states.

TABLE II.

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⁷¹ are not included in these evaluations.

Binding energy (eV)	Li₂	Be₂	B₂	C₂	N₂	O₂	F₂
LRDMC (GF = JSD)	0.976(6)	0.144(7)	2.83(1)	5.74(2)	9.67(1)	4.94(3)	1.27(1)
LRDMC (GF = JsAGPs)	0.9812(12)	−0.0270(12)	2.66(1)	6.01(1)^a	9.91(1)^a	5.06(2)^a	1.56(1)
LRDMC (GF = JPf)	1.0580(12)	0.0304(22)	2.75(1)	6.31(1)^a	9.97(1)^a	5.127(6)^a	1.64(1)
Est. exact	1.06(4)^b	0.1153(3)^c	2.91(6)^d	6.43(2)^e	9.902(3)^e	5.233(3)^e	1.693(5)^e

Binding energy (eV)	Li₂	Be₂	B₂	C₂	N₂	O₂	F₂
LRDMC (GF = JSD)	0.976(6)	0.144(7)	2.83(1)	5.74(2)	9.67(1)	4.94(3)	1.27(1)
LRDMC (GF = JsAGPs)	0.9812(12)	−0.0270(12)	2.66(1)	6.01(1)^a	9.91(1)^a	5.06(2)^a	1.56(1)
LRDMC (GF = JPf)	1.0580(12)	0.0304(22)	2.75(1)	6.31(1)^a	9.97(1)^a	5.127(6)^a	1.64(1)
Est. exact	1.06(4)^b	0.1153(3)^c	2.91(6)^d	6.43(2)^e	9.902(3)^e	5.233(3)^e	1.693(5)^e

^a

See Ref. 69.

^b

See Ref. 139.

^c

See Ref. 140.

^d

See Ref. 141.

^e

See Ref. 142.

FIG. 13.

View large Download slide

Binding energies of the first-row dimers at their experimental bond length (see Table II). The vertical axis shows the difference between the LRDMC and the estimated exact binding energies. The results obtained from the DMC using a Jastrow correlated full-valence-CAS as trial WF [DMC(JFVCAS)] are taken from Ref. 138.

FIG. 14.

View large Download slide

O₂ dissociation curve as obtained with TurboRVB at the variational (VMC) and fixed-node projected (LRDMC) levels of theory and using JSD with DFT orbitals (JDFT), JAGP, and JPf WF Ansätze.

The applications of TurboRVB are not limited to single point evaluations. As discussed in Sec. VI, a distinctive feature of TurboRVB among most QMC codes is the efficient evaluation of atomic forces. This capability has been extensively exploited to investigate several properties of the PES, such as vibrational frequencies, 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¹⁴³ 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.¹⁴⁴ Moreover, it was observed that the optimal grid is a trade-off 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 the work of Zen et al.¹⁰⁷ Afterward, Luo et al.^111,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₂O, H₂S, SO₂, NH₃, and PH₃. The obtained vibrational frequencies are in good agreement with the experimental values.

Zen et al.¹³⁴ applied the second-order LD to study the structural properties of liquid water. They obtained radial distribution functions in good agreement with recent neutron 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 and, in particular, to the inclusion of dispersive vdW forces. This is a stringent test for ab initio approaches as the shape of the peak can be directly compared with neutron scattering data.

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.¹¹⁴ 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 non-local 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¹⁴⁶ 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. Similar 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 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 size-consistent due to the presence of unphysical charge fluctuations,^43,147 but a very flexible (i.e., large) Jastrow factor allows damping the charge fluctuations and recovering the size-consistency 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 occasionally leads 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,¹³⁴ 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,¹⁴⁸ 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.,¹⁴⁹ 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¹⁵⁰ (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 TurboRVB, 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⁵⁴ 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,¹⁵¹ benzene dimers,⁴³ graphite layers,¹⁵² water or methanol molecules adsorbed on a clay surface (see Fig. 15),¹⁵³ and metallic clusters,¹⁵⁴ and to evaluate the cohesive energies of B₂O₃ polymorphs.¹⁵⁵

FIG. 15.

View large Download slide

Adsorption of water and methanol on the hydroxyl-terminated and the silicate-terminated faces of kaolinite (side view in the first row and top view in the second row). The adsorbed molecule on kaolinite is depicted in cyan and gray, and the H-bonds are represented by the blue dashed lines. Reprinted with permission from Zen et al., J. Phys. Chem. C 120, 26402 (2016). Copyright 2016 ACS.

C. Strongly correlated and superconducting materials

Another important QMC application is the study of strongly correlated materials, which 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, such as, 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–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, thecode is expected to describe strongly correlated systems more accurately. In particular, among the various families of correlated WFs that 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 an efficient formulation of a RVB WF. Thus, superconducting materials including cuprates¹⁶¹ and iron pnictides and selenides¹⁶² are prominent applications of the JAGP/RVB WFs. By running extensive VMC calculations for CaCuO₂, a parent compound of cuprate high-temperature superconductors, Marchi et al. validated the JAGP description by correctly finding the expected d-wave symmetry of the pairing function.¹⁶³ Casula and Sorella applied the RVB WF to FeSe, one of the iron-based high-temperature superconductors,¹⁶⁴ 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). 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.¹⁶⁵

FIG. 16.

View large Download slide

(a) Even–even and (b) odd–even components of the AGP pairing function of FeSe at 0 GPa obtained by VMC calculations, where even and odd refer to the parities of the orbitals to the reflection through the Fe plane. The contour plots show $Φ (R_{center}, r)$ with R_center set to be the iron lattice site at the center of the supercell, while r spans the plane defined by the 4 × 4 lattice. Red (yellow) balls are iron (selenium) sites. Arbitrary units of blue (red) intensity indicate negative (positive) regions with corresponding magnitude. Reprinted with permission from M. Casula and S. Sorella, Phys. Rev. B 88, 155125 (2013). Copyright 2013 APS.

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). Casula et al. applied the AGP WF combined with a Jastrow factor to the simplest aromatic compound, that is, the benzene molecule (C₆H₆) for the first time.^127,131 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. Genovese et al. have applied a more general AGP WF, i.e., the Pfaffian, to the benzene molecule.⁶⁹ 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 for the atomization energy. As mentioned in Sec. XII B, Sorella et al. applied the RVB WF to calculate the binding energy of the face-to-face and displaced parallel benzene dimers⁴³ and obtained the binding energy very close to the experimental value in the displaced parallel case. Dupuy and Casula applied the same RVB WF to study the ground-state properties of the oligoacene series, from anthracene up to the nonacene.¹³² 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. Marchi et al. studied the nature of chemical bonds in graphene using the RVB WF¹⁶³ and found that the RVB energy gain becomes extremely small in the thermodynamic limit, i.e., in the infinite lattice. Conversely, Sorella et al. applied the RVB WF to isotropically strained graphene¹⁶⁶ (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 Marchi et al.—but is extremely important in the DIM phase.

FIG. 17.

View large Download slide

The graphene structures studied in Ref. 166. (a) SEM honeycomb, semimetallic; (b) AFI honeycomb antiferromagnetic insulator; (c) DIM dimerized Kekulé-like insulator; and (d) HEX distorted hexagonal insulator. There are two carbons per unit cell in (a) and (b) and six in (c) and (d). t_A, t_B, and t_C schematically denote different hopping integrals magnitudes. Reprinted with permission from Sorella et al., Phys. Rev. Lett. 121, 066402 (2018). Copyright 2018 APS.

FIG. 18.

View large Download slide

Energy gains by the resonance effect, measured by the energy per atom difference between the single determinant Ansatz (Jastrow–Slater determinant WF) and the corresponding multideterminant JAGP WF. The largest energy gain occurs in the DIM state, underlining its resonating valence-bond nature, actually increasing for large strain ε. Small negative values at small strain are finite-size effects. Inset: finite-size scaling of this correlation energy gain in the DIM state at ε = 15%. Reprinted with permission from Sorella et al., Phys. Rev. Lett. 121, 066402 (2018). Copyright 2018 APS.

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 systems. Zen et al. found that the RVB WF gives a correct torsion barrier of the ethylene (C₂H₄) and the triple-singlet gap in the methylene (CH₂).⁶⁶ They also investigated the charge-transfer and diradical nature of the electronic states of the retinal minimal model (penta-2,4-dieniminium cation) at the VMC and LRDMC levels.⁶⁷ 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 are significantly different from the CASSCF ones, inverting the relative stability of several torsional paths, and similar to what was obtained by other correlated approaches. Barborini and Guidoni applied VMC and LRDMC to calculate the torsion barriers of 1,3-butadiene and the ring-opening barrier of cyclobutene,¹⁶⁷ which are in good agreement with the experimental results. Barborini and Coccia investigated the potential energy curves of the two spin states of the tetramethyleneethane (TME) molecule and the two anionic states of TME as a function of the torsion of the central dihedral angle.¹⁶⁸ 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. Barborini et al. investigated the effects of the static and dynamical electronic correlations on the BLA of the 1,3-butadiene¹⁶⁹ using a VMC-JAGP Ansatz and other quantum chemistry calculations. They found inconsistency 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₂H₂)_n-H up to N = 12 acetylene units,¹⁷⁰ in which they revealed that the BLA obtained by the extrapolation to n → ∞ is 0.0910(7) Å, which is compatible with the experimental data. 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.¹⁷³ 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,¹⁷⁴ which is used to calculate the S¹ ← S⁰ absorption energy based on the Bethe–Salpeter formalism. Coccia et al. carried out a systematic basis-set analysis on polarizability, quadrupole moment, and electronic density of the C₂H₂ molecule, with the aim to test VMC and LRDMC calculations in the presence of an external electric field.¹⁷⁵ In particular, a relatively small hybrid (Gaussian + Slater functions) basis set was seen to be able to properly reproduce reference values.

F. Metal-organic complexes

Recently, TurboRVB has also been used to study metal-organic complexes. Chu et al. studied free-energy reaction barriers of water splitting reactions, using a simplified computational model based on the cobalt ion¹⁷⁶ (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₂S₂(SH) $_{2}^{-}$ model complex (high-spin and broken symmetry states)¹⁷⁷ 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.

FIG. 19.

View large Download slide

Reaction cycle involving the cobalt ion model, studied in Ref. 176 by VMC and LRDMC. The pink, red, and white balls represent the cobalt, oxygen, and hydrogen atoms, respectively. Reprinted with permission from Chu et al., J. Chem. Theory Comput. 12, 5803 (2016). Copyright 2016 ACS.

G. Crystal polymorphism

Ab initio electronic calculations allow one to study the phase stability between crystal 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. Sorella et al. studied the pressure-induced metal–insulator transition in bulk silicon through EOS calculations.¹⁷⁸ They obtained a reasonable transition pressure at room temperature by VMC, however still a little far from the experimental value. Devaux et al. studied the α-γ phase transition in elemental cerium at the VMC and LRDMC levels.¹⁷⁹ 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. Hay et al. investigated the effect of long-range van der Waals forces in the SiO₂ polymorphs using DFT and QMC.¹⁸⁰ 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₂ polymorphs. Note 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.¹⁷⁸ 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⁸⁷ developed by Dagrada et al. to save some computational cost. EST works even in realistic systems, namely, the bcc-hydrogen, the bcc-lithium, and the silicon in the high-pressure β-tin phase.⁸⁷

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. Stella et al. investigated the metal–insulator transition in the hydrogen chain by explicitly calculating the complex polarization,¹²³ as shown in Sec. XI. They found that the model has no metal–insulator transition and always behaves as an insulating 1D Hubbard model at half-filling, at least for a large enough interatomic distance (R > 1). Motta et al. (Simons Collaboration on the Many-Electron Problem)¹⁸¹ undertook a comprehensive benchmark study in the finite and infinite linear hydrogen chains using state-of-the-art many-body methods, including VMC and LRDMC with extrapolations to the thermodynamic and the complete-basis-set limits. They provided accurate potential energy curves online, which will be very useful for benchmarking other methods. They have recently investigated the ground states of the hydrogen-chain more in detail, namely, the insulator-to-metal transition, dimerization, and magnetic phases.¹⁸² 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 is subject of intense research. In 1935, Wigner and Huntington predicted its metallization upon compression,¹⁸³ while Ashcroft also proposed its high-T_c superconductivity.¹⁸⁴ 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.¹⁸⁵ 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. Attaccalite and 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.¹⁰⁶ Years later, in a series of works, 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.^113,186 Mazzola and Sorella have also formulated the accelerated first-order LD,¹¹⁹ as described in Sec. VIII, to study such an elusive phase transition with properly equilibrated simulations.¹¹⁹ Finally, under the same framework, the first QMC molecular dynamics simulations of a hydrogen–helium mixture has been performed at Jupiter’s interior conditions.¹²⁰ 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. 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.⁸⁷ Therefore, a careful finite size scaling analysis is needed when the EST technique is employed, although the thermodynamic limit extrapolation is much smoother at the special twist, saving some computational cost. Moreover, EST allows for the determinant/AGP optimization.

FIG. 20.

View large Download slide

Phase diagram of dense hydrogen and a hydrogen–helium mixture. “This work” refers to the results obtained by VMC-Langevin dynamics simulations. The shaded areas represent the liquid–liquid transitions (LLTs) between the insulating-molecular and the metallic-atomic fluids. The solid symbols refer to their QMC LLTs for the pure hydrogen (blue circles) and for the hydrogen–helium mixture (red triangles), and the blue and red solid lines indicate first order LLTs. The empty (solid) left (right) triangles correspond to simulations displaying a clear atomic (molecular) behavior; on the other hand, the red diamonds represent an intermediate behavior. A calculated Jupiter’s adiabat is shown as a gray line. Pure hydrogen first-order LLT predicted by QMC (cyan) and DFT (solid and dashed lines) are also shown. The light green and dark green triangles and the brown stars represent experimental results. Reprinted with permission from G. Mazzola, R. Helled, and S. Sorella, Phys. Rev. Lett. 120, 025701 (2018). Copyright 2018 APS.

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, Barborini et al. calculated the ³B_1u ← ¹A_g vertical triplet excitation of the ethylene molecule and its adiabatic excitation ³A₁ ←¹A_g.¹⁸⁷ One can also deal with nontrivial cases based on a flexible symmetry-adapted RVB WF within its MO representation, which was proposed by Dupuy et al. in 2015.¹⁸⁸ 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.

K. Connection to experimental spectroscopies

Barborini et al. proposed a new ab initio approach to calculate quasiparticle WFs (QPWFs) of isolated systems and resolved in momentum space, based on the correlated sampling technique. They applied the new approach to simulate Scanning Tunneling Microscopy (STM) images¹⁸⁹ and Angle-Resolved Photoemission Spectroscopy (ARPES).¹⁹⁰ The STM image of ${[{C u C l}_{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 TurboRVB can become useful also for state-of-the-art research in this field.

FIG. 21.

View large Download slide

The STM images of in-plane ${[{C u C l}_{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⁻⁵ Å⁻³, while in the right column the isosurfaces take the value of 6.7 × 10⁻⁶ Å⁻³. Panels (a) and (a′) display the LUMO+1 UHF molecular orbitals, while panels (b) and (b′) display the QMC images. Reprinted with permission from Barborini et al., J. Chem. Theory Comput. 12, 5339 (2016). Copyright 2016 ACS.

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, such as Nexus¹⁹¹ and QMC-SW.¹⁹² The workflow system has been implemented in Python 3 in an object-oriented fashion, so the modules are highly extensible. Note 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-genius-sequential.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.⁴⁰ 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 job 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 TurboRVB development dates back to 2003, with the Ph.D. thesis project of Casula supervised by 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 a constrained number of molecular orbitals (Pfn), (iii) the Antisymmetrized Geminal Power (AGP), (iv) the Antisymmetrized Geminal Power with a 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 many variational parameters (so far up to the order of 10⁵) 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 force calculations, keeping them of order N³. 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, which 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.

DATA AVAILABILITY

The data that support the findings of this study are available within this article and the cited references. TurboRVB is available from S.S.’s web site (https://people.sissa.it/∼sorella/) upon request.

ACKNOWLEDGMENTS

The authors acknowledge all contributors for the development of the code TurboRVB over the last 20 years. We would especially like to acknowledge very fruitful collaborations with F. Becca, G. Carleo, C. Cavazzoni, N. Devaux, N. Dupuy, L. Guidoni, H. Hay, R. Hlubina, Y. Iqbal, R. Maezono, F. Mouhat, A. Parola, K. Seki, T. Shirakawa, L. Spanu, L. Stella, L. F. Tocchio, and S. Yunoki.

K.N. is grateful for computational resources from PRACE, Project No. 2019204934, and those from the facilities of the Research Center for Advanced Computing Infrastructure at Japan Advanced Institute of Science and Technology (JAIST). K.N. also acknowledges financial support from the Simons Foundation and that from Grant-in-Aid for Scientific Research on Innovative Areas (No. 16H06439). C.A. acknowledges funding from the European Union Seventh Framework Program under Grant Agreement No. 785219 Graphene Core 2 and COST Action TUMIEE CA17126. M.C. is grateful to the French grand équipement national de calcul intensif (GENCI) for the computational time provided through these years under Project No. 0906493; for the funded PRACE Project Nos. 2012061116, 2015133179, and 2016163936; and for the access to the Hokusai and K-computer granted by the Institute of Physical and Chemical Research (RIKEN). Y.L. was supported by the Argonne Leadership Computing Facility, which is a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-06CH11357. A.Z.’s work was sponsored by the Air Force Office of Scientific Research, Air Force Material Command, U.S. Air Force, under Grant No. FA9550-19-1-7007. S.S. acknowledges financial support from PRIN 2017BZPKSZ and computational resources from PRACE, Project No. 2019204934. S.S. is also grateful to his wife L. Urgias for bearing the too many weekends spent on the development of TurboRVB.

We dedicate this paper to the memory of P. W. Anderson, remembering him as one of the most influential contributors to physics and chemistry of the past century, and in particular for deeply inspiring this work with his RVB theory of matter.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-access-plan.

NOMENCLATURE

$A$: antisymmetrization operator
a, b, I, J: indices of atoms
F_a: atomic force for an atom a, $F_{a} \equiv - \frac{d E}{d R_{a}}$
f_i: general force for the ith variational parameter, $f_{i} \equiv - \frac{\partial E}{\partial α_{i}}$
$f (r_{i}, r_{j})$: spatial part of the pairing function
G: N × N matrix with elements G_i,j = g(i, j)
$g (i, j)$: pairing function between electrons i ↔ (r_i, σ_i) and j ↔ (r_j, σ_j)
i, j: compact notation for (r_i, σ_i) and (r_j, σ_j), respectively
i, j: indices of electrons/orbitals/Monte Carlo sampling
J: Jastrow factor, composed of J₁, J₂, and J_3/4
l, m: indices of orbitals
M: the number of configurations of electrons (the number of Monte Carlo samplings)
N: the number of electrons
N_at: the number of atoms
N_d: the number of electrons with spin down
N_u: the number of electrons with spin up
$O_{k} (x)$: the logarithmic derivative of a many-body WF $O_{k} (x) \equiv \frac{\partial \ln Ψ (x)}{\partial α_{k}}$
R_a: ath ion coordination
r_i: ith electron coordination
S: overlap matrix with elements S_μ,ν ≡⟨ψ_μ|ψ_ν⟩
$S$: variance–covariance matrix of the logarithmic derivative of a many-body WF
$S^{QMC} (R)$: variance–covariance matrix of QMC forces
x: N set of electron coordinations including spins, the same as x_i
x_i: N set of electron coordinations including spins $\equiv {(r_{1} σ_{1}, r_{2} σ_{2}, \dots, r_{N} σ_{N})}_{i}$ ⁠. The index i refers to the Monte Carlo sampling index
Φ: N × k matrix with elements $Φ_{i, k} = ϕ_{k} (i)$
Φ_AGP: antisymmetrized geminal power (AGP)
Φ_AS: antisymmetric (AS) part
Φ_Pf: Pfaffian (Pf)
Φ_SD: single Slater determinant (SD)
Ψ: many-body WF Ansatz, Ψ = Φ_AS × exp J
Ψ_G: many-body guiding function
Ψ_T: many-body trial WF
α₁, …, α_p: variational parameters of the trial WF
α_k: kth variational parameter
μ, ν: indices of orbitals
σ_i: spin of ith electron
$χ_{i} (r)$: primitive atomic orbital for the Jastrow part
$ψ_{i} (r)$: primitive atomic orbital for the antisymmetric part
$ϕ_{k} (r)$: kth molecular orbital for the antisymmetric part
∫di: compact notation for $\sum_{σ_{i}} \int d r_{i}$

REFERENCES

1.

R.

Ramprasad

,

R.

Batra

,

G.

Pilania

,

A.

Mannodi-Kanakkithodi

, and

C.

Kim

,

npj Comput. Mater.

3

,

54

(

2017

).

https://doi.org/10.1038/s41524-017-0056-5

Crossref

2.

K.

Rajan

,

Mater. Today

8

,

38

(

2005

).

https://doi.org/10.1016/s1369-7021(05)71123-8

Crossref

3.

W.

Kohn

and

L. J.

Sham

,

Phys. Rev.

140

,

A1133

(

1965

).

https://doi.org/10.1103/physrev.140.a1133

Crossref

4.

R. M.

Martin

,

Electronic Structure: Basic Theory and Practical Methods

(

Cambridge University Press

,

2004

).

Google Scholar

Crossref

5.

M. G.

Medvedev

,

I. S.

Bushmarinov

,

J.

Sun

,

J. P.

Perdew

, and

K. A.

Lyssenko

,

Science

355

,

49

(

2017

).

https://doi.org/10.1126/science.aah5975

Crossref

PubMed

6.

J. P.

Perdew

and

K.

Schmidt

,

AIP Conf. Proc.

577

,

1

–

20

(

2001

).

https://doi.org/10.1063/1.1390175

Crossref

7.

K.

Burke

,

J. Chem. Phys.

136

,

150901

(

2012

).

https://doi.org/10.1063/1.4704546

Crossref

PubMed

8.

A. J.

Cohen

,

P.

Mori-Sánchez

, and

W.

Yang

,

Chem. Rev.

112

,

289

(

2012

).

https://doi.org/10.1021/cr200107z

Crossref

PubMed

9.

A.

Szabo

and

N.

Ostlund

,

Modern Quantum Chemistry

(

Macmillan Publishing Co., Inc.

,

1982

).

Google Scholar

10.

C.

Møller

and

M. S.

Plesset

,

Phys. Rev.

46

,

618

(

1934

).

https://doi.org/10.1103/physrev.46.618

Crossref

11.

P. J.

Knowles

and

N. C.

Handy

,

Chem. Phys. Lett.

111

,

315

(

1984

).

https://doi.org/10.1016/0009-2614(84)85513-x

Crossref

12.

B. O.

Roos

,

P. R.

Taylor

, and

P. E. M.

Sigbahn

,

Chem. Phys.

48

,

157

(

1980

).

https://doi.org/10.1016/0301-0104(80)80045-0

Crossref

13.

J.

Čížek

,

J. Chem. Phys.

45

,

4256

(

1966

).

https://doi.org/10.1063/1.1727484

Crossref

14.

W.

Foulkes

,

L.

Mitas

,

R.

Needs

, and

G.

Rajagopal

,

Rev. Mod. Phys.

73

,

33

(

2001

).

https://doi.org/10.1103/revmodphys.73.33

Crossref

15.

F.

Becca

and

S.

Sorella

,

Quantum Monte Carlo Approaches for Correlated Systems

(

Cambridge University Press

,

2017

).

Google Scholar

Crossref

16.

N.

Metropolis

and

S.

Ulam

,

J. Am. Stat. Assoc.

44

,

335

(

1949

).

https://doi.org/10.1080/01621459.1949.10483310

Crossref

PubMed

17.

J.

Kim

,

A. D.

Baczewski

,

T. D.

Beaudet

,

A.

Benali

,

M. C.

Bennett

,

M. A.

Berrill

,

N. S.

Blunt

,

E. J. L.

Borda

,

M.

Casula

,

D. M.

Ceperley

,

S.

Chiesa

,

B. K.

Clark

,

R. C.

Clay

,

K. T.

Delaney

,

M.

Dewing

,

K. P.

Esler

,

H.

Hao

,

O.

Heinonen

,

P. R. C.

Kent

,

J. T.

Krogel

,

I.

Kylänpää

,

Y. W.

Li

,

M. G.

Lopez

,

Y.

Luo

,

F. D.

Malone

,

R. M.

Martin

,

A.

Mathuriya

,

J.

McMinis

,

C. A.

Melton

,

L.

Mitas

,

M. A.

Morales

,

E.

Neuscamman

,

W. D.

Parker

,

S. D. P.

Flores

,

N. A.

Romero

,

B. M.

Rubenstein

,

J. A. R.

Shea

,

H.

Shin

,

L.

Shulenburger

,

A. F.

Tillack

,

J. P.

Townsend

,

N. M.

Tubman

,

B.

Van Der Goetz

,

J. E.

Vincent

,

D. C.

Yang

,

Y.

Yang

,

S.

Zhang

, and

L.

Zhao

,

J. Phys.: Condens. Matter

30

,

195901

(

2018

).

https://doi.org/10.1088/1361-648x/aab9c3

Crossref

PubMed

18.

R. J.

Needs

,

M. D.

Towler

,

N. D.

Drummond

, and

P. L.

Ríos

,

J. Phys.: Condens. Matter

22

,

023201

(

2009

).

https://doi.org/10.1088/0953-8984/22/2/023201

Crossref

PubMed

19.

L. K.

Wagner

,

M.

Bajdich

, and

L.

Mitas

,

J. Comput. Phys.

228

,

3390

(

2009

).

https://doi.org/10.1016/j.jcp.2009.01.017

Crossref

20.

C. J.

Umrigar

and

C.

Filippi

, Cornell-Holland Ab-initio Materials Package (CHAMP), http://pages.physics.cornell.edu/∼cyrus/champ.html; accessed

27 October 2019

.

21.

J. S.

Spencer

,

N. S.

Blunt

,

S.

Choi

,

J.

Etrych

,

M.-A.

Filip

,

W. M. C.

Foulkes

,

R. S.

Franklin

,

W. J.

Handley

,

F. D.

Malone

,

V. A.

Neufeld

et al,

J. Chem. Theory Comput.

15

,

1728

(

2019

).

https://doi.org/10.1021/acs.jctc.8b01217

Crossref

PubMed

22.

S.

Zhang

,

J.

Carlson

, and

J. E.

Gubernatis

,

Phys. Rev. B

55

,

7464

(

1997

).

https://doi.org/10.1103/physrevb.55.7464

Crossref

23.

M.

Motta

and

S.

Zhang

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1364

(

2018

).

https://doi.org/10.1002/wcms.1364

Crossref

24.

G. H.

Booth

,

A. J.

Thom

, and

A.

Alavi

,

J. Chem. Phys.

131

,

054106

(

2009

).

https://doi.org/10.1063/1.3193710

Crossref

PubMed

25.

G. H.

Booth

,

A.

Grüneis

,

G.

Kresse

, and

A.

Alavi

,

Nature

493

,

365

(

2013

).

https://doi.org/10.1038/nature11770

Crossref

PubMed

26.

A. J. W.

Thom

,

Phys. Rev. Lett.

105

,

263004

(

2010

).

https://doi.org/10.1103/physrevlett.105.263004

Crossref

PubMed

27.

J. S.

Spencer

and

A. J.

Thom

,

J. Chem. Phys.

144

,

084108

(

2016

).

https://doi.org/10.1063/1.4942173

Crossref

PubMed

28.

N. S.

Blunt

,

T. W.

Rogers

,

J. S.

Spencer

, and

W. M. C.

Foulkes

,

Phys. Rev. B

89

,

245124

(

2014

).

https://doi.org/10.1103/physrevb.89.245124

Crossref

29.

F. D.

Malone

,

N.

Blunt

,

J. J.

Shepherd

,

D.

Lee

,

J.

Spencer

, and

W.

Foulkes

,

J. Chem. Phys.

143

,

044116

(

2015

).

https://doi.org/10.1063/1.4927434

Crossref

PubMed

30.

S.

Ten-no

,

J. Chem. Phys.

138

,

164126

(

2013

).

https://doi.org/10.1063/1.4802766

Crossref

PubMed

31.

Y.

Ohtsuka

and

S.

Ten-no

,

J. Chem. Phys.

143

,

214107

(

2015

).

https://doi.org/10.1063/1.4936412

Crossref

PubMed

32.

S.

Ten-no

,

J. Chem. Phys.

147

,

244107

(

2017

).

https://doi.org/10.1063/1.5003222

Crossref

PubMed

33.

J. R.

McClean

and

A.

Aspuru-Guzik

,

Phys. Rev. A

91

,

012311

(

2015

).

https://doi.org/10.1103/physreva.91.012311

Crossref

34.

A.

Nagy

and

V.

Savona

,

Phys. Rev. A

97

,

052129

(

2018

).

https://doi.org/10.1103/physreva.97.052129

Crossref