Julien Toulouse
email: julien.toulouse@upmc.fr

C J Umrigar
email: cyrusumrigar@cornell.edu

Full optimization of Jastrow-Slater wave functions with application to the first-row atoms and homonuclear diatomic molecules

We pursue the development and application of the recently-introduced linear optimization method for determining the optimal linear and nonlinear parameters of Jastrow-Slater wave functions in a variational Monte Carlo framework. In this approach, the optimal parameters are found iteratively by diagonalizing the Hamiltonian matrix in the space spanned by the wave function and its first-order derivatives, making use of a strong zero-variance principle. We extend the method to optimize the exponents of the basis functions, simultaneously with all the other parameters, namely the Jastrow, configuration state function and orbital parameters. We show that the linear optimization method can be thought of as a so-called augmented Hessian approach, which helps explain the robustness of the method and permits us to extend it to minimize a linear combination of the energy and the energy variance. We apply the linear optimization method to obtain the complete ground-state potential energy curve of the C2 molecule up to the dissociation limit, and discuss size consistency and broken spin-symmetry issues in quantum Monte Carlo calculations. We perform calculations of the first-row atoms and homonuclear diatomic molecules with fully optimized Jastrow-Slater wave functions, and we demonstrate that molecular well depths can be obtained with near chemical accuracy quite systematically at the diffusion Monte Carlo level for these systems.

I. INTRODUCTION

Quantum Monte Carlo (QMC) methods (see e.g. Refs. [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF][START_REF]Quantum Monte Carlo Methods in Physics and Chemistry[END_REF][3] constitute an alternative to standard quantum chemistry approaches for accurate calculations of the electronic structure of atoms, molecules and solids. The two most commonly used variants, variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC), use a flexible trial wave function, generally consisting for atoms and molecules of a Jastrow factor multiplied by a short expansion in configuration state functions (CSFs), each consisting of a linear combination of Slater determinants of orbitals expanded in a localized one-particle basis. To fully benefit from the considerable flexibility in the form of the wave function, it is crucial to be able to efficiently optimize all the parameters in these wave functions.

In recent years, a lot of effort has been devoted to developing efficient methods for optimizing a large number of parameters in QMC wave functions. One the most effective approaches is the linear optimization method of Refs. 4 and 5. This is an extension of the zero-variance generalized eigenvalue equation approach of Nightingale and Melik-Alaverdian [6] to arbitrary nonlinear parameters, and it permits a very efficient and robust energy minimization in a VMC framework. This method was applied successfully to the simultaneous optimization of Jastrow, CSF and orbital parameters of Jastrow-Slater wave functions for some all-electron atoms and molecules in Refs. 4 and 7, and for the C 2 and Si 2 molecules with pseudopotentials in Ref. 5. It has also been applied in Ref. 8 for optimizing Jastrow, CSF and backflow parameters to obtain very accurate wave functions for the firstrow atoms.

In this paper, we pursue the development and application of the linear optimization method. We extend the method to optimize the exponents of the basis functions, simultaneously with all the other parameters, a capacity rarely available in standard quantum chemistry methods (see, however, Refs. [9][10][11]. This uses a slight generalization of the parametrization employed for Jastrow-Slater wave functions to allow nonorthogonal orbitals. Also, we show that the linear optimization method can be thought of as a so-called augmented Hessian approach. This allows us to clearly establish the connection between the linear optimization method and a stabilized Newton optimization method, helping to explain the robustness of the approach. Moreover, this formulation permits us to extend the method to minimize a linear combination of the energy and the energy variance, as done in Ref. 12 with the Newton method. We then apply the linear optimization method to obtain the complete ground-state potential energy curve of the C 2 molecule up to the dissociation limit, and discuss size consistency and broken spin-symmetry issues in QMC calculations. Finally, we perform calculations of the first-row atoms and homonuclear diatomic molecules with fully optimized Jastrow-Slater wave functions, and we demonstrate that molecular well depths can be obtained with near chemical accuracy quite systematically at the DMC level for these systems.

When not explicitly indicated, Hartree atomic units are assumed throughout this work.

II. WAVE FUNCTION OPTIMIZATION A. Wave function parametrization

To optimize a large number of parameters of a wave function, it is important to use a convenient and efficient parametrization that eliminates redundancies. We use an N -electron Jastrow-Slater wave function parametrized as (see Ref. 4)

|Ψ(p) = Ĵ(α)e κ(κ, ζ) NCSF I=1 c I |C I (ζ) , (1)
where Ĵ(α) is a Jastrow operator, e κ(κ, ζ) is an orbital transformation operator, and |C I (ζ) are configuration state functions (CSFs). The parameters to be optimized, collectively designated by p = (α, c, κ, ζ), are the Jastrow parameters α, the CSF parameters c, the orbital rotation parameters κ and the basis exponent parameters ζ.

We use a flexible Jastrow factor consisting of the exponential of the sum of electron-nucleus, electronelectron and electron-electron-nucleus terms, written as systematic polynomial and Padé expansions [13] (see also Refs. 14,15). Each CSF is a short linear combination of products of spin-up and spin-down Slater determinants

|D ↑ k (ζ) and |D ↓ k (ζ) |C I (ζ) = k d I,k |D ↑ k (ζ) |D ↓ k (ζ) , (2)
where the coefficients d I,k are fully determined by the spatial and spin symmetries of the state considered.

The N ↑ -electron and N ↓ -electron Slater determinants are generated from the set of nonnecessarily orthonormal orbitals at the current optimization step,

|D ↑ k (ζ) = â0 † k1↑ (ζ)â 0 † k2↑ (ζ) • • • â0 † kN ↑ ↑ (ζ)|vac and |D ↓ k (ζ) = â0 † kN ↑ +1 ↓ (ζ)â 0 † kN ↑ +2↓ (ζ) • • • â0 † kN ↓ (ζ)|vac
, where â0 † k σ (ζ) (with σ =↑, ↓) is the fermionic creation operator for the orbital |φ 0 k (ζ) (the superscript 0 referring to the current optimization step) in the spin-σ determinant, and |vac is the vacuum state of second quantization. The (occupied and virtual) orbitals are written as linear combinations of N bas basis functions

|χ µ (ζ µ ) with current coefficients λ 0 k,µ |φ 0 k (ζ) = N bas µ=1 λ 0 k,µ |χ µ (ζ µ ) . (3)
Specifically, in this work, we use Slater basis functions whose expression in position representation, using spherical coordinates r = (r, θ, φ) around an atom position r a , is

r|χ µ (ζ µ ) = N nµ (ζ µ ) r nµ-1 e -ζµr S lµ,mµ (θ, φ), (4)
where N n (ζ) = (2ζ) 2n+1 /(2n!) is the radial normalization constant and S l,m (θ, φ) are the normalized real spherical harmonics. Some parameters in the Jastrow factor are fixed by imposing cusp conditions [16] on the wave function; the other Jastrow parameters are varied freely. Due to the arbitrariness of the overall normalization of the wave function, only N CSF -1 CSF coefficients need be varied, e.g. the coefficient of the first CSF is kept fixed. The only restrictions that we impose on the exponents are the equality of the exponents of the basis functions composing the 2l+1 components of spherical harmonics having the same l (e.g., p x , p y and p z ) and, of course, the equality of the exponents of symmetry-equivalent basis functions on equivalent atoms. Finally, the parametrization of the orbital coefficients through the orbital rotation parameters κ, to which we shall come next, allows one to conveniently eliminate the redundancies due to the invariance properties of determinants under elementary row operations.

Because straightforward variation of the exponents ζ of the basis functions results in orbitals being nonorthogonal, in this work we use an orbital optimization formalism that applies to nonorthogonal orbitals. The general idea of this formalism appears also in valence bond theory (see Refs. [17][18][START_REF]Valence Bond Theory[END_REF] and is a direct generalization of the formalism for optimizing orthonormal orbitals used in standard multiconfiguration self-consistent field (MCSCF) theory (see, e.g., Ref. [START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF], and, in a QMC context in Ref. 4. At each optimization step, the orbitals are transformed using the operator e κ(κ, ζ) where κ(κ, ζ) is the total real singlet orbital excitation operator

κ(κ, ζ) = kl κ kl Êkl (ζ), (5)
where the parameters κ kl are nonzero only for nonredundant orbital pairs (see below) and Êkl (ζ) is the singlet excitation operator from orbital l to orbital k

Êkl (ζ) = â0 † k↑ (ζ) b0 l↑ (ζ) + â0 † k↓ (ζ) b0 l↓ (ζ), (6)
where b0 lσ = q (O -1 ) lq â0 qσ is the dual orbital annihilation operator at the current optimization step written in terms of the usual annihilation operators â0 qσ and the inverse of the overlap matrix O of the orbitals with elements O lq = φ 0 l |φ 0 q . The operators â † kσ and blσ satisfy the canonical anti-commutation relations

{â 0 † kσ , â0 † lσ ′ } = 0, { b0 kσ , b0 lσ ′ } = 0 and {â 0 † kσ , b0 lσ ′ } = δ kl δ σσ ′ .
The action of the operator Êkl (ζ) on a spin-σ Slater determinant is simply to replace the l spin-σ orbital by the k spin-σ orbital. Thus, in practice, the calculation of the orbital overlap matrix O is not needed. In contrast to the case of orthonormal orbitals, the operator κ(κ, ζ) is not anti-Hermitian and thus the operator e κ(κ, ζ) is not unitary. The action of this operator on a Slater determinant is seen by inserting e -κ(κ, ζ) e κ(κ, ζ) = 1 after each orbital creation operator making up the determinant and using e κ(κ, ζ) |vac = |vac ; this leads to a new Slater determinant made with the transformed orbital creation operators

â † kσ (κ, ζ) = e κ(κ, ζ) â0 † kσ (ζ) e -κ(κ, ζ) , (7)
and, accordingly, the corresponding transformed orbitals are

|φ k (κ, ζ) = e κ(κ, ζ) |φ 0 k (ζ) = l (e κ ) lk |φ 0 l (ζ) , (8)
where the sum is over all (occupied and virtual) orbitals, and (e κ ) lk are the elements of the transformation matrix e κ constructed as the exponential of the matrix κ with elements κ kl . The nonredundant orbital excitations l → k have already been discussed in Ref. 4. For a single-determinant wave function, the nonredundant excitations are: closed → open, closed → virtual and open → virtual. For a multi-configuration complete active space (CAS) [START_REF] Roos | [END_REF] wave function, the nonredundant excitations are: inactive → active, inactive → secondary and active → secondary. For both single-determinant and multi-determinant CAS wave functions, if l → k is an allowed excitation in Eq. ( 5) then the action of the reverse excitation k → l is zero, and one can choose to impose the condition κ lk = -κ kl . In this case, κ is a real anti-symmetric matrix and thus e κ is an orthogonal matrix that simply rotates the orbitals. For a general multiconfiguration wave function (not CAS), some active → active excitations must also be included and the action of the corresponding reverse excitation is generally not zero. If the reverse excitation is independent, one does not have to enforce the orthogonality condition κ lk = -κ kl , in which case the transformation of Eq. ( 8) is no longer a rotation. Of course, in addition to these restrictions, only excitations between orbitals of the same spatial symmetry have to be considered. We note that this orbital optimization formalism is well suited for the use of localized orbitals, although we do not explore this possibility in this work. We note that it is also possible to keep the orbitals exactly orthonormal when varying the basis exponents (see Refs. 9, 10). If one starts from orthonormal orbitals and uses basis functions that are for instance symmetrically orthonormalized [22,23]

| χµ = N bas ν=1 (B -1/2 ) νµ |χ µ , ( 9
)
where B is the overlap matrix of the basis functions with elements B νµ = χ ν |χ µ , then the orthonormality of the orbitals is preserved during the optimization. However, this complicates the calculation of the derivatives of the wave function with respect to the exponent parameters (in particular, one needs to calculate the derivatives of the matrix B -1/2 as done in Ref. 24), and the computational effort of the optimization is significantly increased. In practice, as discussed in Section IV, we have found that using orthogonalized basis functions does not significantly reduce the number of iterations needed to reach convergence. We denote by N p the total number of parameters to be optimized. The parameters at the current optimization step are denoted by p 0 = (α 0 , c 0 , κ 0 = 0, ζ 0 ) and the corresponding current wave function by

|Ψ 0 = |Ψ(p 0 ) = Ĵ(α 0 ) NCSF I=1 c 0 I |C I (ζ 0 ) . (10)
B. First-order derivatives of the wave function

We now give the expressions for the first-order derivatives of the wave function |Ψ(p) of Eq. ( 1) with respect to the parameters p i at p = p 0

|Ψ i = ∂|Ψ(p) ∂p i p=p 0 , (11)
which collectively designate the derivatives with respect to the Jastrow parameters

|Ψ αi = ∂ Ĵ(α 0 ) ∂α i NCSF I=1 c 0 I |C I (ζ 0 ) , (12)
with respect to the CSF parameters

|Ψ cI = Ĵ(α 0 )|C I (ζ 0 ) , (13)
with respect to the orbital parameters

|Ψ κ kl = Ĵ(α 0 ) NCSF I=1 c 0 I Êkl (ζ 0 ) |C I (ζ 0 ) , (14)
and with respect to the exponent parameters

|Ψ ζµ = Ĵ(α 0 ) NCSF I=1 c 0 I ∂|C I (ζ 0 ) ∂ζ µ . (15)
The derivatives with respect to the orbital parameters in Eq. ( 14) are thus simply generated by single excitations of orbitals out of the CSFs. In the derivatives with respect to the exponent parameters in Eq. ( 15), the orbital transformation operator e κ(κ, ζ) in Eq. 1 does not contribute since the orbitals are transformed at each step so that we always have κ 0 = 0.

C. Linear optimization method

We use the linear optimization method of Refs. 4 and 5 to optimize all the parameters in our wave functions. This is an extension of the zero-variance generalized eigenvalue equation approach of Nightingale and Melik-Alaverdian [6] to arbitrary nonlinear parameters, and it permits a very robust and efficient energy minimization in a VMC context. We review here this approach from a somewhat different perspective and show how the method can be extended to minimize a linear combination of the energy and the energy variance.

At each step of the optimization, the quantummechanical averages are computed by sampling the probability density of the current wave function Ψ 0 (R) = R|Ψ 0 in the N -electron position representation |R = |r 1 , r 2 , • • • , r N . We will denote the statistical average of a local quantity, f (R), by

f (R) = (1/M ) M k=1 f (R k ) with M electron configurations R k .
∆p j |Ψ j , (16)
where |Ψ 0 = |Ψ(p 0 ) is the normalized current wave function, and

|Ψ j = 1 Ψ 0 |Ψ 0 |Ψ j - Ψ 0 |Ψ j Ψ 0 |Ψ 0 |Ψ 0 , (17)
(iii) update the current parameters as p 0 → p 0 + ∆p. The energy minimization step (ii) can be written in matrix notation as

E lin = min ∆p 1 ∆p T E 0 g T /2 g/2 H 1 ∆p 1 ∆p T 1 0 T 0 S 1 ∆p , (19)
where E 0 = Ψ 0 | Ĥ|Ψ 0 is the current energy, g is the gradient of the energy with respect to the N p parameters with components g i = 2 Ψ i | Ĥ|Ψ 0 , H is the Hamiltonian matrix in the basis consisting of the N p wave function derivatives with elements H ij = Ψ i | Ĥ|Ψ j , and S is the overlap matrix in this basis with elements

S ij = Ψ i |Ψ j .
Clearly, the minimization of Eq. ( 19) is equivalent to solving the following (N p + 1)-dimensional generalized eigenvalue equation

E 0 g T /2 g/2 H 1 ∆p = E lin 1 0 T 0 S 1 ∆p , (20)
for its lowest eigenvector. Each optimization step thus consists of a standard Rayleigh-Ritz approach. When applied to the specific case of the orbital rotation parameters, this linear optimization method is known in the quantum chemistry literature as the super configuration interaction or generalized Brillouin theorem approach [25][26][27][28]. It can also be seen as an instance of the class of augmented Hessian optimization methods with a particular choice of the matrix H, sometimes also called rational function optimization methods, which are based on a rational quadratic model of the function to optimize at each step rather than a quadratic one, and which are known to be very powerful for optimizing MCSCF wave functions and molecular geometries (see, e.g., Refs. [29][30][31][32][33][34][35][36][37][38].

We note that after solving Eq. ( 20), and before updating the current parameters, the parameter variations ∆p can be advantageously transformed according to Eq. ( 32 Implementation in variational Monte Carlo. Nightingale and Melik-Alaverdian [6] showed how to realize efficiently this energy minimization approach on a finite Monte Carlo (MC) sample, which is not as obvious as it may seem. The procedure makes use of a strong zerovariance principle and can be described as follows. If the current wave function and its first-order derivatives with respect to the parameters |Ψ 0 , |Ψ 1 , • • • , |Ψ Np form a complete basis of the Hilbert space considered (or, less stringently, span an invariant subspace of the Hamiltonian Ĥ), then there exist optimal parameters variations ∆p so that the linear wave function of Eq. ( 16) is an exact eigenstate of Ĥ, and therefore satisfies the local Schrödinger equation for any electron configuration

R k R k | Ĥ   |Ψ 0 + Np j=1 ∆p j |Ψ j   = R k | E lin   |Ψ 0 + Np j=1 ∆p j |Ψ j   , (21)
where the ∆p j 's are independent of R k . Multiplying this equation by Ψ i (R k )/Ψ 0 (R k ) 2 (where i = 0, 1, ..., N p ) and averaging over the M points R k sampled from |Ψ 0 (R)| 2 leads to the following stochastic version of the generalized eigenvalue equation of Eq. ( 20)

M E 0 M g T R /2 M g L /2 M H 1 ∆p = E lin 1 0 T 0 M S 1 ∆p , (22)
whose lowest eigenvector solution gives the desired optimal parameter variations ∆p independently of the MC sample, i.e. with zero variance. In Eq. ( 22), M E 0 = E L (R) is the average of the local energy E L (R) = R| Ĥ|Ψ 0 /Ψ 0 (R) (the superscript M denoting the dependence on the MC sample), M g L and M g R are two estimates of the energy gradient with components

M g L,i = 2 Ψ i (R) Ψ 0 (R) R| Ĥ|Ψ 0 Ψ 0 (R) = 2 Ψ i (R) Ψ 0 (R) E L (R) - Ψ i (R) Ψ 0 (R) E L (R) , (23)
where

Ψ i (R)/Ψ 0 (R) = Ψ i (R)/Ψ 0 (R) -Ψ i (R)/Ψ 0 (R) has been used, M g R,j = 2 R| Ĥ|Ψ j Ψ 0 (R) = 2 Ψ j (R) Ψ 0 (R) E L (R) - Ψ j (R) Ψ 0 (R) E L (R) + E L,j (R) , (24)
where

E L,j (R) = R| Ĥ|Ψ j /Ψ 0 (R) - [Ψ j (R)/Ψ 0 (R)] E L (R)
is the the derivative of the local energy with respect to the parameter p j (which is zero in the limit of an infinite sample), M H is the following nonsymmetric estimate of the Hamiltonian matrix

M H ij = Ψ i (R) Ψ 0 (R) R| Ĥ|Ψ j Ψ 0 (R) = Ψ i (R) Ψ 0 (R) Ψ j (R) Ψ 0 (R) E L (R) - Ψ i (R) Ψ 0 (R) Ψ j (R) Ψ 0 (R) E L (R) - Ψ j (R) Ψ 0 (R) Ψ i (R) Ψ 0 (R) E L (R) + Ψ i (R) Ψ 0 (R) Ψ j (R) Ψ 0 (R) E L (R) + Ψ i (R) Ψ 0 (R) E L,j (R) - Ψ i (R) Ψ 0 (R) E L,j (R) , (25)
and M S is the estimated overlap matrix

M S ij = Ψ i (R) Ψ 0 (R) Ψ j (R) Ψ 0 (R) = Ψ i (R) Ψ 0 (R) Ψ j (R) Ψ 0 (R) - Ψ i (R) Ψ 0 (R) Ψ j (R) Ψ 0 (R) . (26)
Now, in practice for nontrivial problems, the basis

|Ψ 0 , |Ψ 1 , • • • , |Ψ Np
is never complete, and consequently solving Eq. ( 22) actually gives an eigenvector solution M ∆p and associated eigenvalue M E lin that do depend on the MC sample. But this solution is not the solution that would be obtained by naively minimizing the energy of the MC sample

M E lin = min ∆p 1 ∆p T M E 0 M g T R /2 M g L /2 M H 1 ∆p 1 ∆p T 1 0 T 0 M S 1 ∆p , (27)
which would yield instead a generalized eigenvalue equation similar to Eq. ( 22) but with symmetrized analogues of the energy gradients M g L , M g R and the Hamiltonian matrix M H. In fact, solving the generalized eigenvalue equation of Eq. ( 22) leads to parameter variations with statistical fluctuations about one or two orders of magnitude smaller than the parameter fluctuations obtained using the symmetrized eigenvalue equation resulting from the minimization of Eq. ( 27). Of course, in the limit of an infinite sample M → ∞, the generalized eigenvalue equation of Eq. ( 22) and the minimization of Eq. ( 27) become equivalent.

Minimization of a linear combination of the energy and the energy variance

In some cases, it is desirable to minimize a linear combination of the energy E and the energy variance V :

(1 -q)E + qV where 0 ≤ q ≤ 1. For instance, it has been shown that mixing in a small fraction of the energy variance (e.g., q = 0.05) in the optimization can significantly decrease the variance while sacrificing almost nothing of the energy [12]. Also, there is a theoretical argument suggesting that if one wants to minimize the number of MC samples needed to obtain a fixed statistical uncertainty on the average energy in DMC, then one should minimize in VMC a linear combination of the energy and the energy variance (usually, the energy dominates by far in this linear combination) [39,40]. Indeed, in practice, the number of MC samples needed in DMC is often reduced by a few percent by using q = 0.05 rather then q = 0 [12]. Finally, the energy variance may be more sensitive to some parameters than the energy.

The formulation of the linear optimization method as an augmented Hessian approach shows clearly how to introduce minimization of the energy variance in the method. Suppose that, at each optimization step, we have some quadratic model of the energy variance to minimize

V min = min ∆p V 0 + g T V • ∆p + 1 2 ∆p T • h V • ∆p , (28)
where

V 0 = Ψ 0 |( Ĥ -E 0 ) 2 |Ψ 0
is the energy variance of the current wave function |Ψ 0 , g V is the gradient of the energy variance with components g V,i = 2 Ψ i |( Ĥ -E 0 ) 2 |Ψ 0 and h V is some approximation to the Hessian matrix of the energy variance. Then, alternatively, one could minimize the following rational quadratic model

V min = min ∆p 1 ∆p T V 0 g T V /2 g V /2 h V /2 + V 0 S 1 ∆p 1 ∆p T 1 0 T 0 S 1 ∆p , (29)
which agrees with the quadratic model in Eq. ( 28) up to second order in ∆p, and which leads to the following generalized eigenvalue equation

V 0 g T V /2 g V /2 h V /2 + V 0 S 1 ∆p = V min 1 0 T 0 S 1 ∆p . ( 30
)
On a finite MC sample, the overlap matrix S is estimated as before by the expression given in Eq. ( 26), the energy variance is estimated as

M V 0 = E L (R) 2 -E L (R) 2 , its gradient is calculated as [12] M g V,i = 2 E L,i (R)E L (R) -E L,i (R) E L (R) + Ψ i (R) Ψ 0 (R) E L (R) 2 - Ψ i (R) Ψ 0 (R) E L (R) 2 -M g L,i E L (R) , (31)
and its Hessian can be approximated by the (positivedefinite) Levenberg-Marquardt approximation [12] M h

V,ij = 2 E L,i (R)E L,j (R) -M g L,i E L,j (R) -M g L,j E L,i (R) + M g L,i M g L,j . (32)
We have also tested the use of the exact Hessian of the energy variance of the linear wave function of Eq. ( 16), but this Hessian containing a quadricovariance term tends to be more noisy than the simpler Hessian of Eq. ( 32), is not guaranteed to be positive definite, and leads to a less efficient optimization.

It is clear that an arbitrary linear combination of the energy and the energy variance can be minimized by combining the augmented Hessian matrix of the energy on the left-hand-side of Eq. ( 22) with the augmented Hessian matrix of the energy variance on the left-hand-side of Eq. ( 30) with the estimators of Eqs. (31) and (32). We note however that this procedure destroys the zerovariance principle described in the previous section which holds if only the energy is minimized. In practice, introducing only a small fraction (∼ 5%) of the energy variance does not adversely affect the benefit gained from the zero-variance principle, and in fact in most cases makes the optimization converge more rapidly. We note that if we were to undertake the additional computational effort of computing R| Ĥ2 |Ψ i then it is possible to formulate an optimization method that obeys the strong zero-variance principle even when optimizing the energy variance or a linear combination of the energy and the energy variance.

In fact, following this procedure, any penalty function imposing some constraint, for which we have estimates of the gradient and of some approximation to the Hessian, can be added to the energy and optimized with the linear method, hopefully without spoiling very much the benefit of the zero-variance principle for the energy.

Robustness of the optimization

The linear optimization method has been found to be somewhat more robust than the Newton optimization method of Ref. 41 using an approximate Hessian, and even of that of Ref. 12 using the exact Hessian, for optimizing QMC wave functions. On a finite MC sample, when minimizing the energy, the zero-variance principle of Nightingale and Melik-Alaverdian [6] is certainly a major ingredient in the robustness of the method. But even in the limit of an infinite MC sample, augmented Hessian approaches are known in the quantum chemistry literature to be more robust than the simple Newton method. In fact, it has been shown that augmented Hessian methods can have a greater radius of convergence than the Newton method [32]. This can be understood by rewriting the generalized eigenvalue equation ( 20) as (see Refs. 32, 34)

h -2∆E S • ∆p = -g, (33a)
2∆E = g T • ∆p, (33b)
where h = 2(H -E 0 S) is a N p -dimensional matrix and ∆E = E lin -E 0 < 0 is the energy stabilization obtained from going from the current wave function |Ψ 0 to the linear wave function |Ψ lin (which is necessarily negative within statistical noise). Equations (33) show that the linear optimization method is equivalent to a Newton method with an approximate Hessian matrix h which is level-shifted by the positive-definite matrix -2∆E S, which acts as a stabilizer. This method is thus closely related to the stabilized, approximate Newton method of Ref. 41. The advantage of the present approach, beside the previously discussed zero-variance principle, is that the stabilization constant -2∆E is automatically determined from the solution of the generalized eigenvalue equation. Clearly, from Eq. (33b), -2∆E tends tends to be large far from the minimum, and tends to zero at convergence. In some cases, when the initial parameters are very bad or the MC sample is not large enough, it is necessary to further stabilize the optimization. A variety of different stabilization schemes are conceivable. In practice, we have found that adding a positive constant a diag to the diagonal of the Hamiltonian matrix, i.e. H → H + a diag I where I is the identity matrix, works well. The value of a diag is adjusted at each optimization step by performing three very short MC calculations using correlated sampling with wave function parameters obtained with three trial values of a diag and predicting by parabolic interpolation the value of a diag that minimizes the energy, with some bounds imposed. In addition, a diag is forced to increase if the norm of linear wave function variation || Np j=1 ∆p j |Ψ j || or the norm of the parameter variations ||∆p|| exceed some chosen thresholds, or if some parameters exit their allowed domain of variation (e.g., if a basis exponent becomes negative).

III. COMPUTATIONAL DETAILS

We now give some details of the calculations that we have performed on the first-row atoms and homonuclear diatomic molecules in their ground states.

We start by generating a standard ab initio wave function using the quantum chemistry program GAMESS [42], typically a restricted Hartree-Fock (RHF) wave function or a MCSCF wave function with a complete active space generated by distributing n valence electrons in m valence orbitals [CAS(n,m)]. We use the CVB1 Slater basis of Ema et al. [43], each Slater function being actually approximated by a fit to 14 Gaussian functions [44][45][46] in GAMESS.

This standard ab initio wave function is then multiplied by a Jastrow factor, imposing the electron-electron cusp condition, but with essentially all other free parameters chosen to be initially zero to form our starting trial Jastrow-Slater wave function, and QMC calculations are performed with the program CHAMP [START_REF] Champ | a quantum Monte Carlo program written[END_REF] using the true Slater basis set rather than its Gaussian expansion. The Jastrow, CSF, orbital and exponents parameters are simultaneously optimized with the linear energy minimization method in VMC, using an accelerated Metropolis algorithm [START_REF] Umrigar | [END_REF][START_REF] Umrigar | Quantum Monte Carlo Methods in Physics and Chemistry[END_REF]. We usually start with 10000 MC iterations for the first optimization step and then this number is progressively increased at each step (typically by a factor 1.5 to 4) during the optimization until the energy is converged to 10 -4 (for the lighter systems) or 10 -3 (for the heavier systems) Hartree within a statistical accuracy of 5.10 -5 or 5.10 -4 , respectively. The optimization typically converges in less than 10 steps. Once the trial wave function has been optimized, we perform a DMC calculation within the short-time and fixed-node (FN) approximations (see, e.g., Refs. 50-54). We use an imaginary time step of usually τ = 0.01 Hartree -1 in an Slater determinant multiplied by a Jastrow factor. Crude initial exponents have been intentionally chosen as the integers nearest to the exponent values of the CVB1 basis [43]. In this calculation, the number of MC sample was initially 10000 and this number was progressively increased at each iteration until a final statistical uncertainty of 0.5 mHartree was reached.

efficient DMC algorithm featuring very small time-step errors [55]. For Ne and Ne 2 we computed the DMC energies at four time steps, 0.020, 0.015, 0.01, and 0.005 Hartree -1 , and performed an extrapolation to zero time step.

We found it convenient to start from the CVB1 basis exponents and from CSF and orbitals coefficients generated by GAMESS, but in fact the ability to optimize all these parameters in QMC allows us to also start from cruder starting parameters without relying on a external quantum chemistry program. In this case, a larger number of optimization iterations are needed to achieve convergence.

IV. OPTIMIZATION OF THE BASIS FUNCTION EXPONENTS

In Ref. 4, the optimization of the Jastrow, CSF and orbital parameters has been discussed in detail. We complete here the discussion with the optimization of the exponent parameters.

Figure 1 shows the convergence of the VMC total energy of the all-electron C 2 molecule during the optimization of the 12 exponent parameters in a wave function composed of a Jastrow factor multiplied by a single Slater determinant where the Jastrow and orbital parameters have been previously optimized (with the exponents fixed at the CVB1 values). Crude initial exponents have been intentionally chosen as the integers nearest to the exponent values of the CVB1 basis. One sees that the linear energy minimization method yields a fast convergence of the energy in about three iterations, typically as fast as when optimizing the other parameters. The simultaneous optimization of the Jastrow, CSF, orbital and exponent parameters generally converges as fast as the simultaneous optimization of only the Jastrow, CSF and orbital parameters reported in Ref. 4.

When optimizing the exponents without simultaneous optimization of the orbitals, we have found the optimization very stable, the introduction of the stabilization constant a diag often being unnecessary. However, when optimizing simultaneously the orbital and exponents parameters, the optimization tends to be less stable because of near redundancies between these two sets of parameters, and a diag typically increases up to 10 -3 -10 -4 to retain stability.

Tests on a few atoms have shown that, for the optimization of exponents only, the use of the orthonormalized basis functions of Eq. ( 9) tends to be a bit more stable. Typically, the overlap matrix S of the wave function derivatives have eigenvalues that span about 7 orders of magnitude from 1 to 10 -7 for (unnormalized or normalized) nonorthogonalized functions whereas, for orthonormalized basis functions, the eigenvalues span only about 4 orders of magnitude from 1 to 10 -4 . Thus, orthogonalization of the basis functions reduces the range of the eigenvalues, attenuating near redundancies among the exponent parameters. However, the wave function derivatives with respect to the exponents take significantly longer to compute when using orthonormalized basis functions, and in addition, when also optimizing the orbitals, the near redundancies between some orbital and exponents parameters make a diag increase anyway. Thus, we have not found it worthwhile for our purpose to use orthonormalized basis functions.

For the first-row atoms, optimizing the exponents (simultaneously with the other parameters) rather than using the exponents of the CVB1 basis typically yields improvements of the total VMC energies between 0.1 and 1 mHartree, which is at the edge of statistical significance and accuracy of the optimization. Thus, at the accuracy that we are concerned with, the exponents of the CVB1 basis for these atoms are nearly optimal for Jastrow-Slater wave functions. As expected, larger improvements are obtained for the first-row diatomic molecules. The largest improvement of the total energy is observed for the O 2 molecule with a gain of about 3.4 mHartree in VMC and 1.2 mHartree in DMC with a Jastrow-Slater single-determinant wave function. The largest improvement of the standard deviation of the energy is obtained for the Be 2 dimer using a Jastrow-Slater singledeterminant wave function, with a gain of 0.25 Hartree. (The standard deviation of the energy is σ = 0.57 Hartree with the CVB1 exponents and σ = 0.32 Hartree with the reoptimized exponents using pure energy minimization. A mixed minimization with an energy weighting of 0.95 and a variance weighting of 0.05 and reoptimized exponents results in the same energy but a σ of 0.32, whereas a pure variance minimization yields a variational energy that is higher by 1 mHartree and a σ of only 0.24.) Properties other than the energy might be more sensitive to the basis exponents. We note that the optimization method that we are using is designed to find local minima, and we cannot be sure that we have found the global minimum for the form of the trial wave function considered. In particular, optimization of the exponent parameters typically leads to multiple local minima. We have found that by optimizing the exponents it is possible to reduce the size of the basis, without sacrificing the energy or the energy variance, but the results in this paper were all obtained using a basis size corresponding to the CVB1 basis. A smaller basis has also the advantage of having fewer local minima.

As noted in Ref. 10, because the virial theorem within the Born-Oppenheimer approximation at the equilibrium nuclear geometry holds if the energy is stationary with respect to scaling of the electron coordinates, optimization of the basis exponents along with optimization of the scaling factors of the interelectron coordinates in the Jastrow factor permits one to satisfy exactly the virial theorem in VMC in the limit of infinite sample size.

V. POTENTIAL ENERGY CURVE OF THE C2 MOLECULE

In Fig. 2, we show the total energy curve of the all-electron C 2 molecule as a function of the interatomic distance R calculated in (plot a) RHF, VMC and DMC with a fully optimized Jastrow × singledeterminant wave function [J×SD], and (plot b) MC-SCF CAS (8,8), VMC and DMC with a fully optimized Jastrow × multi-determinant CAS (8,8) wave function [J×CAS (8,8)]. In each case, the horizontal line represents twice the energy of an isolated atom calculated with the same method, and provides a check of the size consistency of the method. We stress that the wave functions have been optimized by energy minimization rather than variance minimization, and in appendix A we present an argument suggesting that, as regards size consistency, energy-optimized wave functions are to be preferred over variance-optimized wave functions. For comparison, we also plot a Morse potential [56], E Morse (R) = E exact (R e )+D e (1-e -x 2 ) 2 where x = 2β(R-R e )/R e and β = ω e /( 4√ B e D e ), using an estimate of the exact energy at equilibrium E exact (R e ) = -75.9265 Hartree [57] and accurate spectroscopic constants: equilibrium distance R e = 2.3481 [14], well-depth D e = 6.44 eV [57], first vibrational frequency ω e = 1855 cm -1 [START_REF] Linstrom | NIST Chemistry WebBook, NIST Standard Reference Database Number 69[END_REF] and rotational constant B e = 1.81984 cm -1 [START_REF] Huber | Molecular Spectra Molecular Structure -IV. Constants of Diatomic Molecules[END_REF]. For analysis, we report in Table I the distribution of the four π electrons among the two carbon atoms A and B in the dissociation limit, the remaining eight σ electrons being unimportant for the study of the dissociation.

We note that Sorella et al. [START_REF] Sorella | [END_REF] have also reported recently QMC calculations of the potential energy curve of the C 2 molecule, using a pseudopotential and Jastrow (8,8), VMC and DMC with a fully optimized Jastrow × multi-determinant CAS (8,8) wave function [J × CAS(8,8)], using the CVB1 Slater basis form [43]. In each case, the horizontal line represents twice the energy of an isolated atom calculated with the same method, and provides a check of the size consistency of the method. For comparison, a Morse potential [56] using accurate spectroscopic constants is also shown (see text). antisymmetrized geminal power wave functions. At very large interatomic distances, lack of ergodicity in the QMC calculations may be an issue, as electrons tend to remain stuck around an atom, and nonequilibrated results can be obtained. In VMC calculations it is always possible to make large moves of the electrons between the two atoms, as done for example in Ref. 60, but this is not possible in DMC calculations where dynamics of the moves is specified and becomes exact only in the small time-step limit. One can nevertheless avoid being deceived by using a large population of walkers (thereby improving the sampling of the configuration space), looking at the evolution of the results as the bond is stretched, and performing several runs with different starting locations of the walkers.

We first discuss the single-determinant case. RHF is of course not size consistent, and leads to a large percentage of incorrect ionic dissociations (62.5%). Our Jastrow factor has a multiplicatively separable form (at dissociation, it reduces to the product of the Jastrow factors employed for the isolated atoms), so that fulfillment of size consistency in VMC calculations is only dependent on the determinantal part of the wave function. The VMC calculation using a single determinant is not size consistent, but the Jastrow factor reduces the size-consistency error and decreases the percentage of ionic dissociations to about 31%. Interestingly, the DMC calculation using the nodes of a non-size-consistent single-determinant trial wave function appears to be size consistent within the accuracy of the calculation, ionic dissociations being absent. Moreover, examination of the distribution of electrons in the DMC calculation shows that the distribution A(↑↓) + B(↑↓) has a vanishing probability in the dissociation limit (see Table I). Only the distributions A(↑↑) + B(↓↓) and A(↓↓) + B(↑↑) remain at dissociation. In appendix B, we show that this implies that the singlet-spin symmetry of the ground-state is broken with an expectation value of the total spin operator Ŝ2 over the FN wave function of 2. In quantum chemistry, it is well known that spin (and/or spatial) symmetry breaking frequently occurs in unrestricted Hartree-Fock or unrestricted Kohn-Sham calculations at large interatomic II).

distances where electron correlation gets stronger. Our results show that spin-symmetry breaking can also occur in DMC calculations even using an unbroken-symmetry trial wave function, meaning that its nodal surface does not impose the spin symmetry. One may wonder how the repeated application of a spin-independent Green function to an initial trial wave function that is a spin eigenstate (neglecting the small spin contamination that can be introduced by imposing spin-dependent electronelectron cusp conditions in the Jastrow factor [61]) can result in a wave function that is not a spin eigenstate. In fact, breaking of spin symmetry is possible in DMC calculations because the Green function is not applied exactly but only by finite sampling. After all, the fermion-sign problem can also been seen as resulting from the breaking of the antisymmetry of the trial wave function due to finite sampling.

We now discuss the multi-determinant case. The MC-SCF calculation in a full valence CAS, i.e. CAS (8,8) for the C 2 molecule, is size consistent, the corresponding atomic calculation being taken as a MCSCF CAS(4,4) calculation as in Refs. 4 and 5. Not surprisingly, the corresponding VMC and DMC calculations are also size consistent. The DMC energy curve agrees closely with the reference Morse potential. For these three calculations, at the dissociation limit, the distributions A(↑↓)+ B(↑↓), A(↑↑) + B(↓↓) and A(↓↓) + B(↑↑) are obtained with equal weights, which is expected for a proper spin-singlet wave function describing two dissociated carbon atoms, as noted in Ref. [START_REF] Sorella | [END_REF].

VI. RESULTS ON FIRST-ROW ATOMS AND HOMONUCLEAR DIATOMIC MOLECULES

In Table II, we report total energies of the firstrow atoms and homonuclear diatomic molecules at their experimental bond length using several computational methods: RHF, MCSCF in a full valence complete active space (FVCAS), VMC with Jastrow × single determinant [J×SD] and Jastrow × multi-determinant FV-CAS [J×FVCAS] wave functions (where the Jastrow, CSF, orbital and exponent parameters have been simultaneously optimized), and DMC with the same J×SD and J×FVCAS wave functions. For atoms, the active space of FVCAS wave functions consists of the 2s and 2p orbitals. For the molecules, it consists of all the orbitals coming from the n = 2 atomic shells, i.e. 2σ g 2σ u 3σ g 1π u,x 1π u,y 1π g,x 1π g,y 3σ u (this is the energy ordering of the Hartree-Fock orbitals for 5 molecules out of the 8 molecules). For the atoms Li, N, O, F and Ne, and for the dimer Ne 2 , orbital occupations and symmetry constraints imply that the FVCAS wave functions contain only a single determinant. Thus, for these systems, the FVCAS MCSCF wave functions are identical to the RHF wave functions, and the J×FVCAS wave functions are identical to J×SD wave functions. The well depths (dissociation energy + zero-point energy) have been calculated consistently by using single-determinant wave functions for both the molecule and the atom, or multi-determinant FVCAS wave functions for both the molecule and the atom. The errors of the computed well depths are plotted in Fig 3.

The largest errors of the DMC total energy using J×FVCAS wave functions are obtained for the heaviest systems and are of the order of 15 mHartree for the atoms and 30 mHartree for the molecules. Of course, one can always improve the total energy by increasing the number of CSFs as done for example by Brown et al. [8], but good well depths are already obtained with J×FVCAS wave functions due to a compensation of errors between the atoms and the molecule. DMC calculations using J×FVCAS wave functions give well depths with near chemical accuracy (1 kcal/mol ≈ 0.04 eV), the largest absolute error being of about 0.1 eV for the N 2 molecule. In particular, we note that, although the Be 2 dimer is unbound at the RHF, MCSCF and VMC level, the weak bond is well reproduced at the DMC level. Because of the extremely weak van der Waals bond of the Ne 2 dimer we computed the DMC energies of Ne and Ne 2 at four time steps τ = 0.020, 0.015, 0.010, and 0.005 Hartree -1 and extrapolated to zero time step. The time-step error at τ = 0.01 for Ne 2 was -0.0023 Hartree whereas that for Ne was -0.00068 Hartree.

VII. CONCLUSIONS

To summarize, we have extended our earlier published linear optimization method to allow for nonorthogonal orbitals. This then makes it possible to optimize all the parameters in the wave function, including the basis exponents. Moreover, by noting that the linear optimization method can be seen as an augmented Hessian method, we have shown that it is possible to minimize a TABLE II: Total energies (in Hartree) and well depths (in eV) of first-row atoms and homonuclear diatomic molecules at their experimental bond lengths R0 (in Bohr) using several computational methods (see text). The RHF and MCSCF calculations have been performed with the Slater CVB1 basis set [43], expanding each Slater function into 14 Cartesian Gaussian functions. The QMC calculations have been performed with the true Slater basis set rather than its Gaussian expansion. The Jastrow, CSF, orbital and exponent parameters of the Jastrow-Slater wave functions have been optimized in VMC, and the resulting wave functions have been used in DMC. The DMC energies are for time step τ = 0.01 Hartree -1 , with the exception of Ne and Ne2 for which the energies extrapolated to τ = 0 are given. The well depths have been calculated consistently using single-determinant (SD) wave functions for both the molecule and the atom, or, full valence complete active space (FVCAS) multi-determinant wave functions for both the molecule and the atom. linear combination of the energy and the energy variance with the linear optimization method. We have applied the method to the calculation of the full ground-state potential energy curve of the C 2 molecule, and we have shown that although a VMC calculation using a spinrestricted single-determinant Jastrow-Slater wave function is not size consistent, the corresponding DMC calculation using the same trial wave function is size consistent within statistical uncertainty. The price to pay for this size consistency is the breaking of the spin-singlet symmetry at dissociation: the fixed-node DMC wave function has an expectation value of 2 for the total spin operator Ŝ2 , although the spin-singlet trial wave function is an eigenstate of Ŝ2 with eigenvalue 0. Of course, using multi-determinant FVCAS Jastrow-Slater wave functions, both the VMC and the DMC calculations are size consistent without breaking of spin symmetry. Finally, we have performed calculations on the first-row atoms and homonuclear diatomic molecules and showed that well depths can be computed with near chemical accuracy using just fully optimized multi-determinant FV-CAS Jastrow-Slater wave functions.

Atoms

) Be2 ( 1 Σ + g ) B2 ( 3 Σ - g ) C2 ( 1 Σ + g ) N2 ( 1 Σ + g ) O2 ( 3 Σ - g ) F2 ( 1 Σ + g ) Ne2 ( 1 Σ + g ) Interatomic distances (
where the coefficients c ij are also consistently determined by imposing the stationarity of the energy

∂E AB ∂c ij = 2 i A j B | ĤAB -E AB |AB AB|AB = 0. (A8)
It is then easy to see that the product wave function |AB P = ψA ψB |vac makes the corresponding energy E P AB stationary

∂E P AB ∂c ij = 2 i A | ĤA -E A |A A|A j B |B B|B + 2 j B | ĤB -E B |B B|B i A |A A|A = 0, (A9)
since both terms vanish according to Eqs. (A6). The product wave function is thus a possible solution, and, if this is the actual solution given by the method, then the method is size consistent.

As an aside, we note that the energy of the product wave function will converge exponentially to the sum of the constituent energies as a function of the interfragment distance because the overlap of the fragment wave functions decays exponentially. On the other hand, the true wave function has an energy that converges only as an inverse power to the sum of the constituent energies (Van der Waals interaction).

Lack of multiplicative separability of variance-optimized linear wave functions

In the case of variance-optimized linear wave functions, the coefficients c iA and c jB of the fragment wave functions of Eq. (A5) are determined by requiring the stationarity of the energy variance

∂V A ∂c iA = 2 i A |( ĤA -E A ) 2 -V A |A A|A = 0, (A10a)
∂V B ∂c jB = 2 j B |( ĤB -E B ) 2 -V B |B B|B = 0. (A10b)
Correspondingly, the coefficients c ij of the composite wave function of Eq. (A7) are also determined by imposing stationarity of the energy variance

∂V AB ∂c ij = 2 i A j B |( ĤAB -E AB ) 2 -V AB |AB AB|AB = 0.
(A11) In contrast to the case of energy-optimized wave functions, the product wave function |AB P = ψA ψB |vac now does not make the corresponding energy variance

V P AB stationary ∂V P AB ∂c ij = 2 i A |( ĤA -E A ) 2 -V A |A A|A j B |B B|B + 2 j B |( ĤB -E B ) 2 -V B |B B|B i A |A A|A + 4 i A | ĤA -E A |A A|A j B | ĤB -E B |B B|B = 0, ( A12
)
since the last term in Eq. (A12) is the product of the energy gradients of the fragments which do not now generally vanish. Thus, the wave function minimizing the energy variance of the composite system is not a product wave function. This suggests that variance minimization does not generally yield additively separable energies, E AB = E A + E B . Further, since the variance is additive for the product wave function, if the method minimizes over a space that includes the product wave function, then V AB ≤ V A + V B . This happens by having anticorrelated energy fluctuations on the two fragments. Of course, in the limit of exact wave functions, the gradients of the energy and of the energy variance simultaneously vanish, and size consistency is ensured. Consequently, the magnitude of the violation of size consistency of variance-optimized linear wave functions is expected to become smaller as the wave function becomes more accurate.

where ψ † σ (r) is the fermionic field creation operator at point r and spin σ. In Eq. (B1), Ψ(r 1 , r 2 , • • • , r N ) is taken as antisymmetric under the exchange of two samespin electron space coordinates and normalized to unity, i.e. dr 1 dr 2 • • • dr N |Ψ(r 1 , r 2 , • • • , r N )| 2 = 1, implying that |Γ is also normalized to unity, Γ|Γ = 1. Even if the wave function |Ψ is not antisymmetric under the exchange of two opposite-spin electrons (product of spinup and spin-down determinants), the total wave function |Γ is always fully antisymmetric.

To study the dissociation of the C 2 molecule, it is sufficient to consider only the four π electrons. The total FN wave function |Γ FN corresponding to the spin-assigned real-space FN wave function Ψ FN (r 1 , r 2 , r 3 , r 4 ) is thus written as Examination of the real-space location of the electrons during DMC calculations using a spin-restricted single-determinant trial wave function shows that, in the dissociation limit, the mixed distribution Ψ(r 1 , r 2 , r 3 , r 4 )Ψ FN (r 1 , r 2 , r 3 , r 4 ) vanishes whenever two electrons of opposite spins are in the neighborhood of the same C nucleus. In contrast, we know from VMC calculations that the trial wave function Ψ(r 1 , r 2 , r 3 , r 4 ) does not forbid dissociation with two opposite-spin electrons around the same nucleus, thus we conclude that it is the FN wave function Ψ FN (r 1 , r 2 , r 3 , r 4 ) that vanishes for this type of dissociation. In other words, it means that at dissociation Ψ FN (r 1 , r 2 , r 4 , r 3 ) can be written as (assuming that inversion symmetry is preserved) Ψ FN (r 1 , r 2 , r 4 , r 3 ) = f A (r 1 , r 2 )g B (r 3 , r 4 ) +g B (r 1 , r 2 )f A (r 3 , r 4 ), (B4)

|Γ
where f A and g B are antisymmetric two-electron functions localized around nuclei A and B, respectively. We now investigate the spin symmetry of this FN wave function. First, using the anticommutation rules of the field operators, it is easy to show that the wave function |Γ FN is an eigenstate of the spin-projection operator In conclusion, we have shown that the singlet-spin symmetry of the ground-state of the C 2 molecule is broken in the FN wave function at dissociation using the nodes of a spin-restricted single-determinant wave function. The expectation value of Ŝ2 over the FN wave function is 2, which is identical to the value found for the lowest broken symmetry solution of the unrestricted Hartree-Fock [67] or unrestricted Kohn-Sham equations with usual approximate density functionals [68].

1 .

 1 Minimization of the energy Deterministic optimization method. The idea of the method is to iteratively: (i) expand the normalized wave function |Ψ(p) = |Ψ(p) / Ψ(p)|Ψ(p) to first order in the parameter variations ∆p = pp 0 around the current parameters p 0 |Ψ lin (p) = |Ψ 0 + Np j=1

  ) of Ref. 4 which corresponds to changing the normalization of the wave function |Ψ(p) [69].

FIG. 1 :

 1 FIG.1: Convergence of the VMC total energy of the allelectron C2 molecule during the optimization of the 12 exponent parameters in a wave function composed of a single Slater determinant multiplied by a Jastrow factor. Crude initial exponents have been intentionally chosen as the integers nearest to the exponent values of the CVB1 basis[43]. In this calculation, the number of MC sample was initially 10000 and this number was progressively increased at each iteration until a final statistical uncertainty of 0.5 mHartree was reached.

FIG. 2 :

 2 FIG.2: Total energy of the all-electron C2 molecule as a function of the interatomic distance R calculated in (plot a) RHF, VMC and DMC with a fully optimized Jastrow × single determinant wave function [J × SD], and (plot b) MCSCF CAS(8,8), VMC and DMC with a fully optimized Jastrow × multi-determinant CAS(8,8) wave function [J × CAS(8,8)], using the CVB1 Slater basis form[43]. In each case, the horizontal line represents twice the energy of an isolated atom calculated with the same method, and provides a check of the size consistency of the method. For comparison, a Morse potential[56] using accurate spectroscopic constants is also shown (see text).

FIG. 3 :

 3 FIG.3: Error in well depths of the first-row homonuclear diatomic molecules using several computational methods (see text and TableII).

FN = 1 2 dr 1 dr 2 dr 3 dr 4 Ψ

 4 FN (r 1 , r 2 , r 3 , r 4 ) ψ † ↑ (r 1 ) ψ † ↑ (r 2 ) ψ † ↓ (r 3 ) ψ † ↓ (r 4 )|vac , (B2)with the antisymmetry constraintsΨ FN (r 2 , r 1 , r 3 , r 4 ) = -Ψ FN (r 1 , r 2 , r 3 , r 4 ),(B3a)Ψ FN (r 1 , r 2 , r 4 , r 3 ) = -Ψ FN (r 1 , r 2 , r 3 , r 4 ). (B3b)

  Ŝz = (1/2) dr[ ψ † ↑ (r) ψ↑ (r)-ψ † ↓ (r) ψ↓ (r)] with eigenvalue zero Ŝz |Γ FN = 0, (B5) for any function Ψ FN (r 1 , r 2 , r 4 , r 3 ). The action of the total spin operator Ŝ2 = Ŝ+ Ŝ-+ Ŝz ( Ŝz -1) with Ŝ+ = dr ψ † ↑ (r) ψ↓ (r) and Ŝ-= dr ψ † ↓ (r) ψ↑ (r) on |Γ FN givesŜ2 |Γ FN = 1 2 dr 1 dr 2 dr 3 dr 4 [2Ψ FN (r 1 , r 2 , r 3 , r 4 ) -4Ψ FN (r 1 , r 3 , r 2 , r 4 )] ψ † ↑ (r 1 ) ψ † ↑ (r 2 ) ψ † ↓ (r 3 ) ψ † ↓ (r 4 )|vac ,(B6)where the anticommutation rules of the field operators and permutations of electron space coordinates have been used. Equation (B6) shows that |Γ FN is generally not an eigenstate of Ŝ2 . At dissociation, it is nevertheless possible to calculate the expectation value of Ŝ2 over|Γ FN Γ FN | Ŝ2 |Γ FN = 2 -dr 1 dr 2 dr 3 dr 4 Ψ FN (r 1 , r 2 , r 3 , r 4 ) ×Ψ FN (r 1 , r 3 , r 2 , r 4 ) = 2,(B7)since the last integral vanishes due to the localized form of Ψ FN (r 1 , r 3 , r 2 , r 4 ) given in Eq. (B4).

TABLE I :

 I Distribution of the electrons among the two carbon atoms A and B of the C2 molecule in the dissociation limit. For the neutral dissociations, only the distribution of the four π electrons are considered since the remaining eight σ electrons are unimportant for the study of the dissociation. The same methods used in Fig.2are compared. For the RHF and MCSCF wave functions, the percentages of the distributions can be calculated analytically.

			neutral dissociation		ionic dissociation
		A(↑↓) + B(↑↓) A(↑↑) + B(↓↓) A(↓↓) + B(↑↑)
	RHF	25 %	6.25 %	6.25 %	62.5 %
	VMC J×SD	≈ 43 %	≈ 13 %	≈ 13 %	≈ 31 %
	DMC J×SD	≈ 0 %	≈ 50 %	≈ 50 %	≈ 0 %
	MCSCF CAS(8,8)	33.33 %	33.33 %	33.33 %	0 %
	VMC J×CAS(8,8)	≈ 33 %	≈ 33 %	≈ 33 %	≈ 0 %
	DMC J×CAS(8,8)	≈ 33 %	≈ 33 %	≈ 33 %	≈ 0 %

  Ref. 62, b Ref. 59, c Ref. 58, d Ref. 14, e Ref. 63, f Ref. 57, g Ref. 64.

						Bohr)
		5.051 b	4.65 c	3.005 d	2.3481 d	2.075 b	2.283 b	2.668 b	5.84 g
				Numbers of CSFs in FVCAS wave functions
		8	38	137	165	107	30	8	1
					Total energies (Hartree)
	RHF	-14.87127	-29.13148	-49.08961	-75.40154	-108.98650 -149.65881	-198.76323 -257.09105
	MCSCF FVCAS -14.89758	-29.22111	-49.22009	-75.63991	-109.13585 -149.76453	-198.84307
	VMC J×SD	-14.98255(5) -29.29768(4) -49.3457(5) -75.8088(5) -109.4520(5) -150.2248(5) -199.4209(5) -257.80956(2)
	VMC J×FVCAS -14.99229(5) -29.33180(5) -49.3916(2) -75.8862(2) -109.4851(3) -150.2436(2) -199.4443(3)
	DMC J×SD	-14.99167(2) -29.31895(5) -49.38264(9) -75.8672(1) -109.5039(1) -150.2872(2) -199.4861(2) -257.84707(5)
	DMC J×FVCAS -14.99456(1) -29.33736(2) -49.4067(2) -75.9106(1) -109.5206(1) -150.29437(9) -199.4970(1)
	Estimated exact -14.995(1)	-29.3380(4) -49.415(2)	-75.9265 f	-109.5427 f -150.3274 f	-199.5304 f -257.8753
					Well depths (eV)
	RHF	0.159	-0.395	0.858	0.668	5.042	1.021	-1.510	-0.00189
	MCSCF FVCAS 0.875	-0.331	2.521	6.105	9.106	3.897	0.662
	VMC J×SD	0.726(3)	-0.048(3)	2.367(3)	4.75(1)	8.88(1)	4.20(1)	0.55(1)	-0.050(5)
	VMC J×FVCAS 0.991(3)	-0.042(3)	2.814(6)	6.369(6)	9.78(1)	4.713(8)	1.19(1)
	DMC J×SD	0.9679(8)	0.125(1)	2.798(3)	5.656(3)	9.583(3)	4.992(7)	1.349(6)	0.004(2)
	DMC J×FVCAS 1.0465(6)	0.0767(8)	2.906(3)	6.482(3)	10.037(3)	5.187(5)	1.645(4)
	Estimated exact 1.06(4) c	0.09(1) c	2.92(6) e	6.44(2) f	9.908(<3) f 5.241(<3) f	1.693(5) f	0.00365 g
	a

Acknowledgments

We thank Peter Nightingale, Sandro Sorella, Andreas Savin, Frank Petruzielo, Roland Assaraf, Benoît Braïda and Paola Gori-Giorgi for valuable discussions. We also thank Alexander Kollias, Peter Reinhardt and Roland Assaraf for providing us with the Gaussian fits of Slater basis functions of Ref. 46. This work was supported in part by a European Marie Curie Outgoing International Fellowship (039750-QMC-DFT), the National Science Foundation (EAR-0530301) and the DOE (DE-FG02-07ER46365). Most of the calculations were performed on the Intel cluster at the Cornell Nanoscale Facility (a member of the National Nanotechnology Infrastructure Network supported by the National Science Foundation) and at the Cornell Theory Center.

Appendix A: A remark on size consistency and variance minimization

In this appendix, we briefly review the concept of size consistency of an electronic-structure method (see, e.g., Refs. 20, 65, 66 for more details), and we give an argument for preferring energy-optimized wave functions over variance-optimized wave functions as regards size consistency.

Definition of size consistency

Consider an electronic system AB made of two noninteracting fragments A and B (e.g., a diatomic molecule at dissociation). This system has a Hamiltonian ĤAB = ĤA + ĤB , (A1)

where ĤA and ĤB are the Hamiltonians of the fragments, commuting with each other. If E A and E B are the (approximate) energies of the fragments given by some method, and E AB is the (approximate) energy of the composite system given by the same method, then this method is said to be size consistent if and only if

i.e., the energy is additive.

In particular, if |A = ψA |vac and |B = ψB |vac are the (approximate) wave functions given by the method under consideration where ψA and ψB are secondquantized wave operators (commuting or anticommuting with each other), then a sufficient condition for size consistency of the method is that it leads to an (approximate) wave function for the composite system of the product form

i.e., the wave function is multiplicatively separable. However, this is not a necessary condition, as exemplified by perturbation theory (see, e.g., the discussion in Ref. 66). Also, in general, due to the nonlocality of the total spin operator Ŝ2 , one has in fact to consider a sum of products of degenerate spin-multiplet component wave functions of the fragments to accommodate non-singlet spin symmetry, but it is sufficient to take the simple product form of Eq. (A3) for our purpose. If the wave function of the system AB has this product form, then it is easy to show that the energy variance is also additively separable

where

is the energy variance of the system AB, and

|B / B|B are the energy variances of the fragments.

Multiplicative separability of energy-optimized linear wave functions

Before discussing variance-optimized wave functions, it is useful to repeat briefly the proof of the multiplicative separability of energy-optimized linear wave functions, given for instance in Ref. [START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF].

Consider that the fragments are described by the following linearly parametrized (approximate) wave functions

and

where |i A = ψiA |vac and |j B = ψjB |vac are some many-body basis states, and the coefficients c iA and c jB are determined by requiring the stationarity of the energy of each fragment

Correspondingly, consider that the composite system is described by a linearly parametrized wave function in the product basis

Appendix B: Spin-symmetry breaking in DMC for the C2 molecule at dissociation

In this appendix, we show how the information on the real-space location of the electrons in DMC calculations of the C 2 molecule at dissociation using a spin-restricted single-determinant trial wave function reveals that the spin-singlet symmetry of the exact ground state is broken in the FN wave function.

For that, we need to determine the expectation value of the total spin operator Ŝ2 over the FN wave function. As we use in the QMC calculation a real-space spin-assigned wave function, we first need to reconstitute the corresponding total wave function. Reference 61 shows how to do so using straightforward first quantization. Here, we use the alternative formalism of real-space second quantization.

In this formalism, the total wave function |Γ in abstract Hilbert space corresponding to a N -electron realspace spin-assigned wave function Ψ(r 1 , r 2 , • • • , r N ) = r 1 ↑, r 2 ↑, • • • , r N ↓ |Ψ with N ↑ spin-up followed by N ↓ spin-down electrons is written as


Keywords: deforestaton, development, credit, public spendings. Codes JEL / JEL classifcaton: C21, O13, Q33

Improving access to capital through credit and public spendings is an important step toward development and poverty alleviaton. At the same tme, deforestaton-related actvites, like agricultural expansion, can be seen as relying on natural capital, through the depleton of forest resources and the use of land in an extensive way. It is then important to beter understand how a beter access to capital infuence the use of land as a natural capital. This paper assesses the relatonship between fnancial development, public spendings and deforestaton. Are they substtute or complement? Our econometric analysis shows that deforestaton is positvely correlated to access to credit and public spendings, which gives some evidence that natural capital is a complement to credit and public spendings.

Avertssement :

Les commentaires et analyses développés n'engagent que leurs auteurs qui restent seuls responsables des erreurs et insufsances.

Introduction

While deforestation is a crucial environmental concern worldwide, it also represents a source of income for many agents and communities, since it results in agricultural expansion and timber harvest [START_REF] Celentano | Welfare outcomes and the advance of the deforestation frontier in the brazilian amazon[END_REF]). In this sense, deforestation can be considered as the use of natural capital (land and forest resources) on the development path by agents lacking of other capital opportunities [START_REF] Azqueta | Valuing nature: From environmental impacts to natural capital[END_REF]). Indeed, in developing countries natural capital represents a larger part of their wealth than in developed ones [START_REF] Ruta | Handbook of sustainable development., chapter The capital approach to sustainability[END_REF]World Bank (The) 2006).

Using natural capital for development purposes may be related to the issue of complementarity or substitutability between deforestation and man-made capital that has been investigated and debated in the literature (see among others, [START_REF] Barbier | The role of natural resources in economic development[END_REF]). Optimal growth models have been extended by some authors (e.g. [START_REF] Dasgupta | The optimal depletion of exhaustible resources[END_REF][START_REF] Solow | Intergenerational equity and exhaustible resources[END_REF][START_REF] Stiglitz | Growth with exhaustible natural resources: Efficient and optimal growth paths[END_REF] where the conditions of technical progress and substitution between natural capital and material capital are analyzed. Their approach implicitly suggests that substitution between natural and material capital can create the conditions of unbounded economic development, even in a world with exhaustible resources i.e. finite natural capital. This view has been challenged by proponents of strong sustainability who after Daly and Georgescu-Roegen argue in favor of maintaining each type of capital (Daly and Townsend 1993;[START_REF] Daly | Georgescu-Roegen versus Solow/Stiglitz[END_REF]. For instance, [START_REF] Jouvet | La croissance verte : de l'intention à la mise en oeuvre[END_REF] argue that natural capital depletion can affect the marginal productivity of the whole economy. The Meadows Reports or more recently Diamond gave examples of unsustainable and collapsing societies which have relied too heavily on natural capital depletion (Meadows et al. 1972;[START_REF] Meadows | The Limits to Growth: The 30-year Update[END_REF][START_REF] Diamond | Collapse: How Societies Choose to Fail or Survive[END_REF]).

In this paper, we intend to go further the weak versus strong sustainability debate by considering how man-made capital that encompasses credit and public spending has an effect on deforestation.

Improving access to credit may be seen as a key issue to help rural communities to develop, as credit constraints are an important barrier to development [START_REF] Anjini | An empirical investigation of rationing constraints in rural credit markets in india[END_REF]) and the development of credit markets an important factor of growth [START_REF] Hung | Credit market development and economic growth[END_REF]). Several papers have been considering empirically the impact of the development of credit banking on the efficiency of investments (Odedokun (1996)). [START_REF] Ge | Financial development, bank discrimination and trade credit[END_REF] suggest that informal capital can be a substitute to formal capital when credit markets are poorly developed. Following [START_REF] Rudel | The paths to rain forest destruction: Crossnational patterns of tropical deforestation, 1975-1990[END_REF], [START_REF] Geist | What Drives Tropical Deforestation? A Meta Analisys of Proximate and Undeliyng Causes of Deforestation Based on Subnational Case Study Evidence[END_REF] argue that deforestation may be capital driven. Some studies in Latin America do find evidence that access to credit favors more deforestation-related activities than the others [START_REF] Barbier | Economic analysis of deforestation in mexico[END_REF]Pfaff, 1997). In contrast, [START_REF] Godoy | Household determinants of deforestation by amerindians in honduras[END_REF] consider in their studies on Bolivia and Honduras that families with a better access to credit are less forest-dependent than others. Finally, [START_REF] Monela | Tropical Rainforest Deforestation, Biodiversity Benefits,and Sustainable Land Use : Analysis of Economic and Ecological Aspects Related to the Nguru Mountains,Tanzania[END_REF] finds a positive relationship between credit availability and forest clearing in Tanzania.

Government expenditures have also often been thought to foster deforestation (e.g. [START_REF] Angelsen | Rethinking the Causes of Deforestation: Lessons from Economic Models[END_REF][START_REF] Geist | What Drives Tropical Deforestation? A Meta Analisys of Proximate and Undeliyng Causes of Deforestation Based on Subnational Case Study Evidence[END_REF]. Several authors have evidenced the role of transport infrastructures (e.g. [START_REF] Pfaff | What Drives Deforestation in the Brazilian Amazon?: Evidence from Satellite and Socioeconomic Data[END_REF]. In a similar vein, [START_REF] Bulte | On the gains of committing to inefficiency: Corruption, deforestation and low land productivity in Latin America[END_REF] show that rural subsidies towards large farmers triggered deforestation in Latin America. At last, [START_REF] Galinato | The role of government spending on deforestation and carbon dioxide emissions from land use change[END_REF] evidence a positive impact of government spending on deforestation. Interestingly they do not confirm the effect of a composition effect: the provision of public goods has no effect on deforestation. Perhaps, one possible explanation of the results is whether there is an omitted variable.

A crucial question is whether public spending and access to credit jointly affect deforestation process. In this paper, we try to fill a gap in the literature by exploring the simultaneous effect of credit and public spending, which both contribute to man-made capital, on deforestation. Does this capital substitute or complement to natural capital?

We argue that assessing the impact of public spending (access to credit) on deforestation cannot be made without considering access to credit (public spending). Indeed, the link between credit and public spending is well established in the macro-economic literature. On the one hand, there could exist a crowding in effect à la Barro (Barro 1990). Under that hypothesis, private production depends on both credit and productive government expenditures which considered as an input to private production. Marginal productivity of private capital is therefore positively affected by public spending. On the other hand, a crowding out effect can be put forward according to which an increase in public spending dries out the credit available to private agents. The rest of this paper is organized as follows. We define first the kinds of mechanism that relate deforestation to credit and public spending. Then we assess the nature of this relationship, using panel-data analysis relating credit development and public spending to deforestation. Section 2 presents a simple model of substitutability and complementarity between deforestation and man-made capital. Section 3 presents our econometric analysis, and section 4 concludes.

A simple model of capital input and net income maximization

In this section, we investigate the channels by which a better access to private capital or public spendings may be related to deforestation. Consider a representative agent, that maximizes its net income:

max D,K,δ I(D, K, δ) = A(D, δK) + O((1 -δ)K) -c(D) -c(K) (1)
Net income is derived from the agents' economic activities, that depend on its access to two types of assets: natural capital, through deforestation D, and made-man capital (henceforth called capital) K. By this term, we mean here capital requiring both credit and public expenses that benefit to agents1 . Income A(D, δK) comes from activities related to deforestation (agricultural expansion, timber harvesting, fuelwood collection...), while income O((1 -δ)K) comes from activities not related to deforestation (subsidies to direct consumption, agricultural intensification, trade and business, work to the city...). δ is the share of capital allocated to deforestation-related activities. We assume here standard properties2 :

A i ≥ 0, A ii ≤ 0, ∀i = D, K; O K ≥ 0, O KK ≤ 0.
The costs function that are defined in equation 1 are increasing and convex: :

c i ≥ 0, c ii ≥ 0, ∀i = D, K.
The cost of access to natural capital is related to property rights, as well as to the availability of the resource: Are forests open access? Are property rights safe and well defined?

Are forests a scarce resource? The cost of credit is related to the level of financial development, credit supply and public expenses in the country: Is credit easily accessible to rural communities?

Is credit supply abundant? Which citizens are targeted by redistributive policies?

First-order conditions implicitly give the level of deforestation D * , capital K * and share of capital allocated to deforestation-related activities δ * :

X = A D -c D = 0 (2) Y = δA K + (1 -δ)O K -c K = 0 (3) Z = KA K -KO K = 0 (4)
Equation 2 indicates that deforestation D * increases with its marginal productivity and decreases with its marginal cost. Equation 3 shows that capital K * is chosen in order to equal marginal productivity and marginal cost of capital. Equation 4 defines the allocation of capital δ * as the equalization of the marginal productivity of capital between deforestation-related and deforestation-unrelated activities.

In this context, if there is some development of the financial system, an increase in credit supply or an increase in public expenses that result in a decrease of the marginal cost of capital, it is likely to have an indirect impact on deforestation. First, looking at equation 3 and applying the implicit function theorem, the decrease in marginal cost of capital tends to increase the access to capital for the agent:

K * c K = - Y c K Y K * = 1 δ * A KK + (1 -δ * )O KK -c KK < 0 (5)
Second, from equation 4, this increased access to capital changes its allocation between deforestation-related and deforestation-unrelated activities:

δ * K * = - Z K * Z δ * = - δ * A KK -(1 -δ * )O KK KA KK + KO KK (6)
The denominator being unambiguously negative, the sign of equation 6 is the one of the nominator. Therefore, we can see that the share of capital allocated to deforestation-related activities increases if the marginal productivity of those activities decreases more slowly than the one of deforestation-unrelated activities. In other words, the agents allocate a larger part of capital to deforestation-related activities if the related additional amount of capital is more productive in those activities than in deforestation-unrelated activities.

Third, using equation 2, it is possible to see what could be the end-of-pipe effect of an increase in access to capital on deforestation:

D * (δ * K * ) = - X (δ * K * ) X D * = - δ * A DK -c DK A DD -c DD (7)
Here again, the denominator is unambiguously negative. Therefore, an increase in capital allocated to deforestation-related activities tends to increase deforestation if the net marginal productivity of deforestation is increasing in capital allocated.

Overall, this two-steps relationship defines the impact of an increased access to capital on deforestation. Deforestation may be positively related to better access to capital if :

• an increased access to capital increases the profitability of deforestation-related activities more than the one of deforestation-unrelated activities and capital and deforestation are complement:

δ * A KK -(1 -δ * )O KK > 0, δ * A DK -c DK > 0;
• an increased access to capital increases the profitability of deforestation-unrelated activities more than the one of deforestation-related ones and capital and deforestation are substitute:

δ * A KK -(1 -δ * )O KK < 0, δ * A DK -c DK < 0. Although this case is theoretically possible,
it is quite unlikely in real life, as an increased in access to capital is likely to increase the marginal profitability of deforestation, mainly through machinery.

Since the first case is more likely than the second one, we consider that observing a positive relationship between deforestation and access to capital describes a situation of complementarity between natural capital and made-man capital, while a situation of substitutability will be considered in the case of a negative relationship.

Within this very simple model, the term "capital" encompasses both credit and public expenses.

However, those two types of asset may have diverse relationships with deforestation. Moreover, as better access to credit and better management of public spendings is a crucial interest in the implementation of development policies, it is important to investigate those relationships, whether deforestation tends to be positively or negatively related to access to those two other forms of capital.

3 The empirical analysis of deforestation, credit, and public spendings Data and econometric specifications are described first, then the results are presented and discussed. Regressions are run on a sample of forested developing countries in Asia, Latin America as well as Sub-Saharan Africa excepted small islands of which forested areas are negligible.3

Data set

The rate of deforestation from FAOstat databases is the dependent variable, measured over the 1990 2010 period. Other variables are either from the World Bank Indicators or the PRS group for variables measuring institutional quality. 4 All variables in the model are computed over five years periods (1990-1194,1995-1999,2000-2004,2005-2010) to catch delays in the deforestation adjustment processes.

The two interest variables which approach the theoretical notion of man-made capital derived from our model are domestic credit provided by the banking sector over GDP (Credit) and government consumption over GDP (GovCons).

Three groups of control variables are considered, following the usual literature on deforestation [START_REF] Angelsen | Rethinking the Causes of Deforestation: Lessons from Economic Models[END_REF] and [START_REF] Barbier | The Economics of Tropical Deforestation and Land Use: An Introduction to the Special Issue[END_REF], among others).

Deforestation is influenced by structural factors. Rates of deforestation may be characterized by inertia thus inducing the effect of lagged on current deforestation rates (Lagged Def or). According to the Environmental Kuznets Curve (EKC) assumption, GDP per capita and squared GDP per capita in logs (Lgdppc and Lgdppcsq) should respectively influence deforestation positively and negatively. Nevertheless, as shown in Choumert et al. ( 2013) existing studies provide ambiguous results: [START_REF] Bhattarai | Institutions and the environmental kuznets curve for deforestation: A crosscountry analysis for latin america, africa and asia[END_REF] or [START_REF] Culas | Deforestation and the environmental kuznets curve: An institutional perspective[END_REF] find an EKC contrary to [START_REF] Meyer | Institutional, social and economic roots of deforestation: a cross-country comparison[END_REF]. Concerning population in logs (Lpop), the expected impact is also ambiguous. Indeed, if population fuels the demand for arable lands, fuelwood or charcoal [START_REF] Cropper | The interaction of population growth and environmental quality[END_REF], it could also promote the demand for forest products [START_REF] Foster | Economic growth and the rise of forests[END_REF] and thereby having a negative impact on deforestation rates. Whether urbanization (U rbpop)

drives deforestation is also an empirical question. On the one hand, pressure on forestland would decline with the rural-urban migration but on the other hand, urbanization could raise the level the demand for agricultural products.5

The macroeconomic environment matters in the deforestation process. We assume that the real exchange rate (Reer) and its instability (ReerInst) could lead to deforestation. The first variable catches the competitiveness of the export sector: an increase in the real exchange rate means a real appreciation. Real appreciation is expected to dampen deforestation since exported forest products are negatively affected [START_REF] Arcand | Deforestation and the real exchange rate[END_REF]. The second variable is a proxy for the macroeconomic instability which can boost deforestation. Indeed risk adverse agents could mitigate the influence of macroeconomic instability relying more on the depletion of natural resources, including forests. We also assume that the higher the debt service (DebtBurden), the higher the deforestation. Indeed, [START_REF] Kahn | Third-world debt and tropical deforestation[END_REF] outline a significant positive association between public debt and deforestation in LDCs.

Institutions could also affect the deforestation process. Following [START_REF] Bhattarai | Institutions and the environmental kuznets curve for deforestation: A crosscountry analysis for latin america, africa and asia[END_REF] and [START_REF] Ferreira | Deforestation, Property Rights, and International Trade[END_REF], better institutions preserve forests. In this paper, we approach the institutional quality with an index extracted from the International Country Risk Guide (ICRG) database. More precisely, we consider two aspects of institutional quality: political stability and public integrity as put forward by [START_REF] Ferreira | Governance and timber harvests[END_REF] with three measures that are alternatively introduced:

government stability (Gov Stab),6 internal conflicts (Int Conf ) and corruption (Cor), of which an increase means an improvement. We therefore assume that these variables have a negative impact on deforestation.

Econometric specification

The estimations are made using country-specific and period-specific fixed-effects to control respectively omitted variables constant over time (e.g. country characteristics) and omitted variables common to the countries (e.g. international prices). The Arellano and Bond's GMM estimator for dynamic panel is implemented for two reasons [START_REF] Arellano | Some tests of specification for panel data: Monte carlo evidence and an application to employment equations[END_REF]. Firstly, the OLS estimator is inconsistent since the lagged dependent variable is introduced besides country fixed-effects [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF]. Secondly the GMM estimator controls for the potential endogeneity of the two interest variables due to measurement error, reverse causality or omission of relevant variables.The equations are estimated in first differences and the lagged levels of the explanatory variables are used as instruments. 7 The validity of the instruments is tested by the Sargan/Hansen overidentification test and by the serial correlation test AR(2), where the null hypothesis is that the error term does not exhibit second order auto-correlation.

Results

The results related to our control variable (table 1) are consistent with the literature: ( 1 Burden has a positive impact on deforestation (equations 6 to 9). A real appreciation protects the forest while its instability is detrimental to it (equation 10). ( 5) Institutional quality variables (equations 7 to 10) have the expected negative signs.

As regards our interest variables, Credit and Gov Cons have weakly significant impacts on deforestation when introduced separately (equations 1 and 2). This result is interpreted as a symptom of omitted variable bias: other things held equal, Gov Cons seems to crowd out Credit.

Interestingly, marginal effects are lower than those reported in equation 3 when Credit and Gov Cons are introduced simultaneously. Equation 3 therefore illustrates a complementarity effect between man-made capital and natural capital. These positive signs are robust to the introduction of the controls (equations 5 to 10). When interacted (equation 4), Gov Cons x Credit have a positive though weakly significant effect on deforestation which may interpreted as a crowding in impact of Gov Cons (resp. Credit) on Credit (resp. Gov Cons) à la Barro. Moreover, according to equation 3, the impact of an increase in our interest variables on deforestation rates can be simulated. A one standard deviation increase in Gov Cons (Credit) leads to an increase in the deforestation rate, which represents 21% (28%) of its average value. This is a rather important impact highlights the intensity in the complementarity effect between man-made capital and natural capital. This paper draws attention on conflicting environmental and development objectives: capital as a mean of development may be built at the expenses of natural assets like tropical forests.

Indeed we find evidence that an increased access to credit and public spendings seem to boost more deforestation-related activities than others. In this sense, man-made capital and natural capital may be considered as complement more than substitute in this context.

An important research recommendation here is then to deepen the analysis of this complementarity and to find potential sources of substitutability. Indeed it is important to understand in which context a better access to credit -an crucial development requisite-is not achieved at the expense of forests depletion. Disentangling the access to man-made capital as a way to boost depletion of natural capital is then essential.

In terms of policy recommendation, our results can be interpreted as evidencing the existence of a Tinbergen rule which states that several objectives can be fulfilled only with a sufficient number of instruments: reducing the threat on forests is not likely to indirectly come from better access to credit and local development, unless further effort are made to orientate the investment in deforestation-unrelated activities. Environmental instruments need to be added to traditional economic instruments in the achievement for instance of the Millennium Development Goals.

Those suggestions have to be considered in the light of development policies. According to [START_REF] Humphreys | Logjam: Deforestation and the crisis of global governance[END_REF], the World Bank has long considered forests as an abundant and under-harvested asset, and focused on the need of developing countries to build access to capital in order to develop forest-related activities. In contrast, since the 1990, the World Bank seemed to took into consideration this complementarity between natural assets and human-made assets. In this sense, the WB strategy around forests has been to focus on the development of markets for goods and services (NTFP, carbon sequestration, eco-tourism), that could be described in our model as deforestationunrelated activities. At the same time, the focus has been put in order to decreases the expected profitability of deforestation-related activities (timber certification for instance). A next step of

  ) The dynamic panel estimator is not rejected by the aforementioned tests. The coefficient of the lagged dependent variable is significantly positive which involves an AR(1) process of deforestation; (2) The Environmental Kuznets Curve is not validated by our data; (3) Population and urban population (equations 5 to 10) respectively affect positively and negatively the rate of deforestation; (4) The macroeconomic environment variables have the expected signs. Debt service as measured by Debt

  Building capital is an important goal of economic development. Financial development provides means to achieve this objective. Increasing public spending may also be eligible for the provision of capital in developing countries. At the same time, deforestation and land-use change can be considered as the use of natural capital to achieve development. It follows that assessing the relationship between those different forms of capital is important, in a world in which poverty alleviation through access to credit and deforestation are both important development and environment issues.

	4 Concluding Remarks
		1	2	3	4	5	6	7	8	9	10
	Lagged Defor	0.0614	0.0437	0.0630	0.0697	0.0758	0.0665	0.0869	0.0713	0.0645	0.0640
		(2.5695)**	(1.8398)*	(3.5838)***	(7.6934)***	(4.4244)***	(3.6970)***	(5.7899)***	(4.1028)***	(2.8216)***	(3.8277)***
	Lgdppc	-0.0015	0.0013	0.0014	0.0030	0.0008	0.0006	0.0012	0.0007	0.0003	-0.0011
		(-0.5220)	(0.4054)	(0.6832)	(2.1805)**	(0.3938)	(0.2516)	(0.5432)	(0.2658)	(0.1073)	(-0.4628)
	Lgdppc sq	3.92E-0.5	-0.0001	-0.0002	-0.0003	-0.0002	-7.58E-05	-9.78E-05	-6.58E-05	-3.05E-05	4.94E-05
		(0.1898)	(-0.6181)	(-1.2975)	(-2.7286)***	(-1.0552)	(-0.3983)	(-0.5479)	(-0..3403)	(-0.1507)	(0.2937)
	L Pop	0.0050	0.0044	0.0025	0.0024	0.0037	0.0012	0.0016	0.0020	0.0008	0.0038
		(2.2303)**	(1.8569)**	(1.5279)	(1.9036)*	(1.9462)*	(0.7957)	(0.9664)	(1.1809)	(0.5357)	(1.9976)**
	Urb Pop					-6.55E-05	-5.94E-05	-7.89E-05	-6.30E-05	-7.11E-05	-0.0001
						(-1.3352)	(-1.2585)	(-1.9693)*	(-1.3401)	(-1.6085)*	(-2.7497)***
	GovCons		4.87E-05	8.70E-05	9.43E-05	8.70 E-05	6.71E-05	6.24 E-05	5.65E-05	7.87E-05	3.88E-05
			(1.6150)	(3.5642)***	(4.2523)***	(3.4939)***	(2.2840)**	(2.0197)**	(1.8413)*	(2.4337)**	(1.9472)*
	Credit	8.91E-06		1.38E-05	-2.90E-06	1.70E-05	1.65E-05	1.71 E-05	1.72E-05	1.43E-05	6.90E-06
		(1.2175)		(3.9166)***	(-0.6344)	(4.4122)***	(3.7301)***	(3.7009)***	(3.3848)***	(3.0944)***	(2.1663)**
	GovCons x Credit				4.35E-07
					(1.6790)*
	Debt Burden						3.87E-05	2.72E-05	3.67E-05	3.58E-05
							(2.0682)**	(2.1016)**	(1.8743)*	(2.7503)***
	Gov Stab							-0.0016
								(-2.3565)**
	Int Conf								-0.0001
									(-2.2379)**
	Cor									-0.0003	-0.0002
										(-2.0472)**	(-2.6252)***
	Reer										-3.29E-06
											(-1.6605)*
	Reer Inst										9.81E-06
											(3.4689)***
	Obs	230	230	230	230	230	216	216	216	216	223
	Periods included	4	4	4	4	4	4	4	4	4	4
	Countries included	61	61	61	61	61	57	57	57	57	59
	Number of instruments	27	27	39	59	40	41	42	42	42	51
	OID test p value	0.70	0.69	0.52	0.47	0.51	0.74	0.67	0.75	0.79	0.81
	AR(2) p value	0.25	0.22	0.25	0.29	0.24	0.27	0.25	0.28	0.28	0.23

Table 3 :

 3 VariablesReerReal effective exchange rate, base 100 in 1990. Weights determined by the country's first ten partners (imports and exports), oil countries excluded or not. An increase means a real appreciation. from the International Country Risk Guide database. It assesses political violence in the country and its actual or potential impact on governance. Ranks from 0 to 12. An increase means an improvement i.e. a decrease in the intensity of political violence.

	Variables	Description	Source	Average (Standard deviation)	Median
	Dependent
	variable:
	Defor	Average annual deforestation rate	FAOstat. Four periods: 1991-95; 1996-00; 2001-05; 2006-	0.0021 (0.0039)	0.0018
			10. FAO production yearbook provided instruments for
			lagged values of deforestation rates
	Explanatory	All variables are period averages	Four periods: 1991-95; 1996-00; 2001-05; 2006-10.
	variables:
	Interest variables
	Credit	Domestic Credit provided by the Banking Sector, % of	World Bank Data: http://data.worldbank.org	42.7307 (40.0493)	30.5907
		GDP.
	Gov cons	General government final consumption expenditure (% of	idem	13.3456 (5.1048)	12.1443
		GDP)
	Control variables
	GDPPC	GDP per capita, constant 2000 USD.	idem	2395.261 (3436.476)	1024.6370
	Pop	Population, total.		69.07E+06 (2.12E+08)	13.17E+06
	UrbPop	Urban population, % of total.	idem	47.1430 (21.5917)	44.9071
	DebtBurden	Total debt service (% of exports of goods, services and	idem	5.4792 (5.1698)	4.3196
		income).
			Authors' calculations, World Bank Data	111.9092 (35.5919)	102.8959
	Reer Inst	Instability of the Real effective exchange rate. It is	Authors' calculations	4.8109 (14.4654)	2.7615
		calculated as the squared residuals of the country-
		individual regressions of Reer on one-year lagged Reer and
		a deterministic trend
	Corruption	Corruption index from the International Country Risk	http://www.prsgroup.com/ICRG.aspx	2.5161 (0.8849)	2.4792
		Guide database. Ranks from 0 to 6. An increase means an
		improvement i.e. a decrease in corruption.
	Int Conf	Internal conflicts index idem	8.6084 (1.6650)	8.8333
	Gov Stab	Government stability index from the International Country	idem	8.0008 (1.7687)	8.4458
		database. It both assesses the government's ability to carry
		out its declared program(s), and its ability to stay in office.
		Ranks from 0 to 12. An increase means an improvement

Those two kinds of capital will be distinguished in the empirical section.

The subscripts refer to first and second derivatives.

The list of countries included in regression analysis is given in the appendix

Variables are described in the appendix.

[START_REF] Defries | Deforestation driven by urban population growth and agricultural trade in the twenty-first century[END_REF] stress the positive influence of urbanization on deforestation with recent data which contrast with results obtained in previous studies for the 1980s and the 1990s.

The variable is introduced in logs. In order to deal with zero values, we add 1 to the variable before computing logs.

The lagged deforestation rate and the interest variables are instrumented by their lagged values. The controls variables are considered as exogenous.To deal with the problem of instruments proliferation, the matrix of instruments is collapsed to ensure that the number of instruments does not exceed the number of countries[START_REF] Roodman | A note on the theme of too many instruments[END_REF].

Acknowledgements / Remerciements

This paper has been presented at the 2012 CERDI conference on ''Environment and Natural Resources Management in Developing and Transition Economies'' (Clermont Ferrand, France, October 2012) and at the 2013 EAERE conference (Toulouse, France, June 2013). We thank participants for their valuable comments.

Appendix