Kohn-Sham orbitals and potentials from quantum Monte Carlo molecular densities.

Kohn-Sham orbitals and potentials from quantum Monte Carlo molecular densities Daniele Varsano, Matteo Barborini, and Leonardo Guidoni Citation: The Journal of Chemical Physics 140, 054102 (2014); doi: 10.1063/1.4863213 View online: http://dx.doi.org/10.1063/1.4863213 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/5?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Quantum Monte Carlo study of the first-row atoms and ions J. Chem. Phys. 134, 084105 (2011); 10.1063/1.3554625 Quantum Monte Carlo study of first-row atoms using transcorrelated variational Monte Carlo trial functions J. Chem. Phys. 126, 164109 (2007); 10.1063/1.2715581 On the universality of the long-/short-range separation in multiconfigurational density-functional theory J. Chem. Phys. 126, 074111 (2007); 10.1063/1.2566459 Quantum Monte Carlo calculations of the dissociation energies of three-electron hemibonded radical cationic dimers J. Chem. Phys. 124, 024318 (2006); 10.1063/1.2150818 Correlated Monte Carlo electron-pair density for the atoms helium to neon J. Chem. Phys. 109, 7075 (1998); 10.1063/1.477390

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 146.189.194.69 On: Tue, 07 Apr 2015 07:57:32

THE JOURNAL OF CHEMICAL PHYSICS 140, 054102 (2014)

Kohn-Sham orbitals and potentials from quantum Monte Carlo molecular densities Daniele Varsano,1,a) Matteo Barborini,2,3 and Leonardo Guidoni3,b) 1

Dipartimento di Fisica, Sapienza-Università di Roma, P.le Aldo Moro 5, 00185 Roma, Italy Dipartimento di ingegneria e scienze dell’informazione e matematica, Università degli studi dell’Aquila, Via Vetoio, Coppito, 67100 L’Aquila, Italy 3 Dipartimento di Scienze Fisiche e Chimiche, Università degli studi dell’Aquila, Via Vetoio, Coppito, 67100 L’Aquila, Italy 2

(Received 22 October 2013; accepted 13 January 2014; published online 3 February 2014) In this work we show the possibility to extract Kohn-Sham orbitals, orbital energies, and exchange correlation potentials from accurate Quantum Monte Carlo (QMC) densities for atoms (He, Be, Ne) and molecules (H2 , Be2 , H2 O, and C2 H4 ). The Variational Monte Carlo (VMC) densities based on accurate Jastrow Antisymmetrised Geminal Power wave functions are calculated through different estimators. Using these reference densities, we extract the Kohn-Sham quantities with the method developed by Zhao, Morrison, and Parr (ZMP) [Phys. Rev. A 50, 2138 (1994)]. We compare these extracted quantities with those obtained form CISD densities and with other data reported in the literature, finding a good agreement between VMC and other high-level quantum chemistry methods. Our results demonstrate the applicability of the ZMP procedure to QMC molecular densities, that can be used for the testing and development of improved functionals and for the implementation of embedding schemes based on QMC and Density Functional Theory. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4863213] I. INTRODUCTION

In the last years the application of Quantum Monte Carlo (QMC) methods1 to the study of molecular systems has become more and more popular due to the great accuracy, the good scalability with the number of electrons and the excellent parallelization of the algorithms. Beside ground state total energy calculations, we have assisted to the progressive development of QMC methods to obtain electronic properties,2, 3 electronic excitations,4, 5 vibrational frequencies,6–10 accurate relaxed molecular geometries,5, 6, 11–16 and chemical reaction barriers:17, 18 recently also the possibility to use QMC methods for non-adiabatic molecular dynamics has been explored.19 Moreover, thanks to the great precision that it can achieve, this method has been a very important and essential tool for the development of Density Functional Theory,20 both for benchmark issues and for the development of exchange and correlation functionals, like in the emblematic case of the Local Density Approximation.21, 22 In this paper, we show the possibility to use very accurate electronic densities calculated via QMC in order to obtain Kohn-Sham (KS) orbitals and exchange-correlation potentials that reproduce those target densities. Since the birth of Density Functional Theory, several methods have been developed to extract Kohn-Shan quantities from accurate ab initio densities,23–31 that are usually obtained through post

a) Electronic mail: [email protected]. Current address: Centro S3,

CNR Istituto di Nanoscienze Via Campi 213/A, 41125 Modena, Italy.

b) Electronic mail: [email protected]

0021-9606/2014/140(5)/054102/15/$30.00

Hartee-Fock methods like Configuration Interaction (CI) or Coupled Cluster (CC). Due to the computational limitations of the aforementioned methods, accurate KS orbitals and exchange correlation potentials have been calculated only for atoms and small molecules.29, 32–37 These limitations arise mainly because the CI and CC methods scale as the 5th–7th power of the number of electrons and very large basis sets are required to obtain accurate descriptions of the electronic densities. On the other hand the QMC techniques are correlated wave function methods which scale at most as the 4th power of the number of electrons in the system, becoming computationally favoured for large molecules. For this reason it would be highly desirable to define a procedure to obtain KS quantities from QMC densities. To the best of our knowledge the only attempt in this direction has been made by Filippi, Gonze, and Umrigar38 for atomic systems (Be and Ne atoms) exploiting the radial symmetry of the density. In this work we aim to develop a scheme to obtain reliable potentials and Kohn-Sham quantities from Quantum Monte Carlo densities for molecular systems on cubic grids, where it is not possible to take advantage of the spherical symmetry. One of the main problems in using QMC for this purpose is that electronic densities, as all the other observables in QMC, are affected by stochastic errors that might be problematic for the numerical instability of the inversion algorithm. Here we extend the procedure to reconstruct the exchange-correlation potential from QMC densities also to molecular systems using very accurate Jastrow Antisymmetrised Geminal Power (JAGP) wave functions39 and an improved estimator for the evaluation of the electronic density.3, 40 The extraction of the

140, 054102-1

© 2014 AIP Publishing LLC


054102-2

Varsano, Barborini, and Guidoni

exchange-correlation potentials and Kohn-Sham orbitals is performed using the method developed by Zhao, Morrison, and Parr27 (ZMP) implemented in a mixed Gaussian and plane-wave approach.41 The calculated vxc [ρ QMC ] potentials are compared with those obtained from CISD densities built with the largest computationally affordable basis sets. The definition of a solid and stable procedure to obtain exchange-correlation potentials by many-body QMC densities will be verified for several important test cases including small molecules. Due to its favourable scaling with the number of electrons, QMC calculations will therefore be able to provide accurate reference densities and potentials for larger molecular systems, where post Hartree-Fock methods become unaffordable. The extracted KS orbitals and eigenvalues can be used as fundamental ingredients to benchmark and develop exchange-correlation potentials of real molecular test cases. Another important application could be found in the calculation of the kinetic potentials42 used in embedding approaches43 with the aim of coupling accurate wave function based methods, such as QMC, with Density-Functional Theory (DFT). This article is organized as follows: in Sec. II we outline the computational methods used for the QMC procedure, the adopted variational wave functions, the calculation of the electronic density and we shortly recall the method used to extract the related Kohn-Sham quantities. In Sec. III we provide the computational details and implementations for both the QMC calculations and the ZMP procedure, and in Sec. IV we report our results for He, Be and Ne atoms, calculated using both logarithmic radial grids and cubic regular grids, and for H2 , Be2 , H2 O, and C2 H4 molecules, discussing the open problems and future improvements. Finally, our conclusions are summarized in Sec. V.

II. COMPUTATIONAL METHODS

In this section, we briefly describe the schemes used for the calculation of the QMC densities, and the methodology used to extract the KS orbitals and the exchange correlation potentials that reproduce these densities. In Sec. II A we summarize the QMC methods used, in Sec. II B we give details on the variational wave functions we have used, in Sec. II C we illustrate the algorithm used to compute the densities and finally in Sec. II D we shortly recall the method developed in Ref. 27 to obtain the KS orbitals and potentials reproducing the target densities.

A. Quantum Monte Carlo methods

The Quantum Monte Carlo methods are stochastic approaches to evaluate the mean values of quantum physics observables over a defined wave function. Here we limit ourselves to a short description of the principal QMC method used in this work, which is the Variational Monte Carlo (VMC) described in Ref. 1 together with other QMC algorithms. In the VMC approach the energy

J. Chem. Phys. 140, 054102 (2014)

functional ¯ ¯ R})] = E[T ({α,

¯ Hˆ T (¯x; {α, ¯ x¯ T (¯x; {α, ¯ R}) ¯ R})d , (1) ¯ 2 d x¯ |T (¯x; {α, ¯ R})|

which is the integral over the set of x¯ = {xi , . . . , x3N } electronic spin and Cartesian coordinates, is evaluated through a stochastic procedure, with respect to a trial wave function ¯ depending on a set of parameters α. ¯ R}) ¯ To stochasT (¯x; {α, tically solve the integral in Eq. (1), the integrand is rewritten as the product of two functions ¯ ¯ R})] = Eloc (¯x)(¯x)d x¯ , (2) E[T ({α, which are, respectively, the local energy, Eloc (¯x) = Hˆ T (¯x; ¯ ¯ x; {α, ¯ R}), that is the energy of an elec{α, ¯ R})/ T (¯ tronic configuration x¯ associated with the trial wave func¯ ¯ R}), and the probability density to find tion T (¯x; {α, ¯ 2/ ¯ R})| the electrons in that configuration (¯x) = |T (¯x; {α, 2 ¯ ¯ R})| d x¯ . |T (¯x; {α, The principle is to sample a certain number N of independent electronic configurations x¯ distributed according to (¯x), evaluating for each of them the local energy: in this way the EV MC energy defined as the mean value of Eloc over these configurations, ¯ E[T ({α, ¯ R})] ≈ EV MC =

N 1 Eloc (x¯i ), N i=1

(3)

gives an estimation √ of the energy functional, with an error that decreases like 1/ N . The great advantage of these methods lies in the fact that, no matter how sophisticated the parametrization α¯ of the trial wave function, it is possible to obtain an evaluation of the 3Ndimensional integral (Eq. (3)) within a certain stochastic accuracy. Next, the energy functional is minimized with respect to the set of parameters α, ¯ giving a variational estimation of the ground state energy, and thus finding an optimized wave function to describe the ground state: for this purpose in this work we have used the linear method described in Refs. 44 and 45. It is thus clear that the accuracy of the VMC calculations is strictly related to the functional form of the adopted trial wave function. B. Variational wave function

In the last years in the QMC community many efforts have been made to define compact and highly correlated wave functions. Some of them are based on the Resonating Valence Bond (RVB)46–48 theory introduced by Pauling49 and already applied to closed shell molecules in the early 1950s by Pople and Hurley.50 In this work we have adopted one of the most successful wave function based on this approach: the Jastrow Antisymmetrised Geminal Power (JAGP).39 This wave function has already been applied to the study of various molecular systems,8, 39 transition metals,9 graphene sheets,51 and weak van deer Waals,44 and hydrogen bonded systems,52 showing its reliability in predicting correct energies and structural parameters even at the VMC level.


054102-3


J. Chem. Phys. 140, 054102 (2014)

The JAGP39 wave function, is built as the product between an Antisymmetric Geminal Power (AGP)50 and a Jastrow factor J (¯r). In the case of closed shell atomic systems of N electrons in a spin singlet state, i.e., N/2 = N↑ = N↓ , the determinantal AGP part of the wave function is written as the antisymmetrised product:

By substituting this equation in the traditional QMC estimator, and integrating by parts, we obtain the simple ACS formulation40 N 1 ∇r2i T2 (r1 , . . . , rN ) 1 ρACS (r) = − . 4π i |ri − r| T2 (r1 , . . . , rN ) (¯r)

N/2

↑ ↓ G xi ; xi ,

(4)

λμa νb ψμa (ri ) ψνb (rj ) |0, 0 .

(5)

AGP (¯x) = Aˆ

i=1

of geminal functions G (xi ; xj ) =

M a,b=1 μ,ν

These G (xi ; xj ) are defined as a linear combination of products of two atomic orbitals, of quantum numbers μ, ν = (n, l, lz ) and centred on the ath and bth atoms, in a spin singlet state, i.e., |0, 0 = 1/2(|↑1 |↓2 − |↓1 |↑2 ). The Jastrow factor we have used, described in Ref. 8, satisfies both the nuclear and electronic cusp conditions, and also treats the dynamical correlation between electrons. C. VMC density estimators

The calculation of the electronic density through the VMC method is essentially the mean value of the three dimensional delta function: N (3) ρ(r) = δ (ri − r) , (6) i=1

(¯r)

(¯r)

over the sampled probability (¯r), and so it is based on the stochastic evaluation of the integral:

(3) ¯ |T (¯r)|2 N i=1 δ (ri − r)d r ρ(r) = . (7) 2 |T (¯r)| d r¯ In practice, it is obtained by dividing the three dimensional real space r in a discrete grid of finite volume intervals, in which the stochastic sampling is accumulated and afterwards renormalized to the total electronic charge of the system. Unfortunately, this straightforward procedure has two disadvantages related to the space discretization: (i) in those regions in which the charge varies rapidly, i.e., the nuclear cusps, it is usually required the use of a very fine grid to describe the correct distribution; (ii) the use of fine grids, however, may lead to big statistical fluctuations due to the low sampling, especially in those regions, like the density tails, in which the number of samplings is particularly low. These drawbacks make the QMC electronic densities very difficult to compare with those obtained with other post-HF methods that are calculated by the analytic integration of the wave functions, even though these densities usually converge very slowly with the size of the basis set. To overcome this drawback Assaraf, Caffarel, and Scemama40 (ACS) have recently defined an improved charge density estimator. The ACS estimator is based on the definition of the Green’s function for the Poisson equation: δ (3) (ri − r) = −

1 1 2 ∇ri . 4π |ri − r|

(9) Thanks to this definition, the density in each point can be calculated using all the electronic sampling reweighed with respect to the coulomb interaction and the second derivative of the wave function, reducing in this way the stochastic fluctuations of the traditional estimator and providing a continuous density. Unfortunately, two main drawbacks remain. The first is that the variance of the density calculated on the nuclei is unbounded as an electron ri approaches the nucleus;40 second, the density can assume non-physical negative values in those regions which are poorly sampled (usually the tails of the charge density distribution), because of the particular reweighing factor that appears in the estimator. To partially cure both aspects, Assaraf et al. have proposed the introduction in Eq. (9) of two auxiliary functions f(ri , r) and g(r),40 so that the final formulation for the ACS estimator is N 1 1 − g(r) ρACS (r) = − 4π i |ri − r| ∇i2 [f (ri ; r)T2 (r1 , . . . , rN )] . (10) × T2 (r1 , . . . , rN )

(8)

The function f(ri , r), which is only required to satisfy the condition that f(ri = r; r) = 1,40 is split in a short- and a longregime functions, respectively, fSR and fLR , defined as fSR = 2Za (|ri − Ra | − |r − Ra |)

(11)

fLR = (1 + l|ri − r|)e−l|ri −r| .

(12)

and

While fSR is used to cure the divergence of the variance near the nuclear cusps, fLR is needed to reintroduce the correct exponential decay of the density through the choice of the parameter l. Assaraf et al. suggested to choose the switching between the two regimes according to the variance of the density in each point of the space. Recently, some of us have proposed in Ref. 3 a more compact expression for the f(ri , r) function, coupling equations (11) and (12) together through Hill functions depending on the distance |r − Ra | between the nucleus and the position in which the density is evaluated. This expression has been generalized also for molecular applications and is based on the choice of two parameters regulating the switching point, and the steepness of the Hill functions.3 The choice of all the parameters is discussed in Sec. III A. Finally, the function g(r): g(r) =

1 |r − Ra |

(13)


054102-4


J. Chem. Phys. 140, 054102 (2014)

is introduced in Eq. (10) to reduce density fluctuations in the |r| → +∞ regime and to correct the negative density values in the poorly sampled regions.3, 40 D. Kohn-Sham orbitals and potentials from QMC density

i

In order to extract Kohn-Sham orbitals, orbital energies, and exchange correlation potentials from the QMC densities obtained as described above, we have used the method developed by ZMP.27 The ZMP procedure is based on the Levy constrained-search method53 for minimizing the noninteracting Kinetic energy23, 54 and given a target density ρ 0 (r) permits to calculate the Kohn-Sham orbitals such that |φi (r)|2 = ρ0 (r), (14) i∈occ

by solving self-consistently the differential equations: 1

− ∇r2 + vext (r) + vconstr (r) φi (r) = i φi (r), 2 where vext (r) is the external nuclear potential, and ρ(r ) − ρ0 (r )

dr vconstr (r) = |r − r |

(15)

(16)

is an effective potential arising from having considered a constraint of the form: 1 [ρ(r) − ρ0 (r)][ρ(r ) − ρ0 (r )] drdr = 0 C[ρ, ρ0 ] = 2 |r − r | (17) in the minimization of the kinetic energy, and is a Lagrange multiplier for the constraint C. As suggested in Ref. 27 we add to the left side of Eq. (15) the Fermi-Amaldi potential55, 56 defined as (1 − N1 )vH (r), where N is the number of electrons and vH (r) is the classical Hartree potential: ρ(r ) vH (r) = dr . (18) |r − r | In this way we arrive to the ZMP equation:

1 1

− ∇r2 + vext (r) + 1 − vH (r) + vconstr (r) φi (r) 2 N = i φi (r).

the Kohn-Sham kinetic energy. This can be achieved by solving the Kohn-Sham equations with the exchange correlation potential of Eq. (20) and inverting the density-functional expression for the total energy: i + EH [ρ] + ρ(r)vxc (r)dr, (21) Exc [ρ] = E[ρ] −

(19)

As demonstrated in Ref. 27 the solution of Eq. (19) in the limit

→ ∞ provides the Kohn-Sham energies and Kohn-Sham orbitals that satisfy Eq. (14). Moreover, the exchange correlation potential vxc (r) that allows to reproduce the target density in a Kohn-Sham equation, can be numerically calculated as 1

(20) vxc (r) = lim vconstr (r) − vH (r) .

→∞ N The insertion of the Fermi-Amaldi term in Eq. (19) permits to speed up convergence, and to ensure the correct long range behaviour −1/r in the exchange correlation potential (Eq. (20)). The calculation of the vxc (r) potential from ab initio densities allows us also to provide an estimation of the exchange-correlation energy Exc [ρ], defined as Exc [ρ] = (Vee [ρ] − EH [ρ]) + (T [ρ] − Ts [ρ]) where Vee is the total electron-electron repulsion, EH is the Hartree energy and Ts is

where i are the eigenvalues of the KS equations and E[ρ] is the total ab initio energy associated with the target density ρ QMC . III. COMPUTATIONAL DETAILS A. Quantum Monte Carlo density calculations

The QMC study of the electronic densities for atoms and molecules is carried out using all-electron JAGP wave functions, fully optimized at the variational level using the linear method with Hessian acceleration described in Ref. 44, and implemented in the TurboRVB quantum Monte Carlo package57 by Sorella. For each of the three atoms considered, He, Be, and Ne, basis set convergence studies have been performed for both the fermionic AGP part of the wave function, and for the Jastrow three-body term. For the fermionic AGP part the reduced aug-cc-pVDZ (aDZ) and the reduced aug-cc-pVTZ (aTZ) basis sets are used as starting points for the wave function optimizations. The biggest exponents (Z > 100) of the contracted s orbitals of the two basis sets have been discarded as the electron-nucleus cusp condition is already described through the one body Jastrow factor. Also all the orbitals with angular momentum l > 2 were discarded. For the Ne atom we have also used as starting point the primitive Gaussians, with Z < 300, of the fully uncontracted UGBS1V basis set without d orbitals. To optimize this basis set made of 12 Gaussian primitive for s and p orbitals we have chosen to follow the constrained procedure described in Ref. 8. In this procedure the AGP wave function is projected in a subspace built through a fixed number of molecular orbitals. For the case of the neon atom we have used the smallest set of doubly occupied molecular orbitals, i.e., 5, corresponding to a Jastrow Single Determinant wave function. The basis sets used for the three-body Jastrow factor of the three atoms were built from contracted orbitals of primitive Gaussians, respectively (3s2p)/[2s1p], and (3s2p)/[2s2p], and for the more complex Neon atom also the fully uncontracted Gaussian basis set of (3s2p1d) orbitals was adopted. The basis sets for the two diatomic molecules H2 , Be2 have been built following the same scheme used for the three atoms, although considering only fully uncontracted Gaussian basis sets for the Jastrow factor. Energy convergence studies for the fermionic part of the H2 molecule were done using the cc-pVDZ (DZ), cc-pVTZ (TZ), and aTZ basis sets reduced as described above, while for the three body Jastrow factor only the (2s1p) uncontracted basis set was used. As for the H2 molecule, also for the Be2 dimer the AGP basis set convergence has been tested using reduced aDZ and aTZ basis sets, while two Jastrow basis sets were used with (3s2p) and (3s2p1d) orbitals. For the H2 O molecule, we have


054102-5


used two different wave functions that have been optimized in the recent work by Zen et al.58 where the structural, vibrational and electronic properties of this molecule have been calculated. The first wave function is the traditional JAGP described above. In this case the AGP part of the wave function is built starting form the aug-cc-pVDZ (aDZ) basis set for both the hydrogens and the oxygen atoms, reduced as described above. The basis set used for the Jastrow factor is built using (4s3p1d) uncontracted Gaussian orbitals for the oxygen atom, and of (2s1p) Gaussian orbitals for the hydrogens. Also the second wave function is of the AGP type although the basis set functions that appear in the geminal expression (Eq. (5)) are not atomic orbitals, but hybrid orbitals,58 built as linear combinations of the basis set centred on the atoms. We used 13 hybrid orbitals for the oxygen and 4 for the hydrogens. These hybrids are built from a linear combination of (12s9p4d3f) uncontracted Gaussian orbitals for the oxygen atom and (9s4p3d) for hydrogen. The uncontracted Gaussian basis set used for the Jastrow factor is built of (4s2p1d) primitives for the oxygen atom and of (2s1p) for the hydrogens.58 For the C2 H4 molecule we have used the wave function optimized by Barborini et al.5 with the constrained procedure described in Ref. 8. In this case we have projected the AGP ansatz on 12 doubly occupied molecular orbitals5 in order to recover the property of size-consistency. For the carbon and hydrogen atoms we have used a basis set of (9s5p2d)/[4s3p2d] and (5s2p)/[3s2p] contracted Gaussian orbitals, respectively. The basis set used for the Jastrow factor was composed of (4s4p)/[2s2p] contracted Gaussian orbitals for the carbon atoms, and (3s3p)/[1s1p] for the hydrogens. The wave function optimization procedure is the same for all the studied systems and was carried out via three consequent steps. First only the one and two body Jastrow factors together with the λij of Eq. (5) and the contracted orbitals’ coefficients of the atomic basis sets were optimized. The second step consisted in optimizing only the parameters of the complete Jastrow factor, keeping fixed all the parameters of the fermionic AGP part, and finally all the parameters of the wave function, including the coefficients and exponents of the basis set, were optimized together. Using the most accurate basis sets obtained through the stochastic optimization procedures we have calculated accurate electronic densities through the traditional (Eq. (6)) and the ACS (Eq. (10)) estimators. To build the QMC densities on uniform cubic grids for all the atoms and molecules, we have used the traditional estimator (Eq. (6)) accumulating 2.048 × 109 MC steps per electron, while to evaluate the linear atomic densities of He, Be, and Ne over a radial grid of 400 points logarithmically distributed up to 35 bohrs from the nucleus, we have used the ACS estimator. Due to the computational cost of the ACS procedure, only 1.024 × 109 MC steps per electron were stored. To calculate these densities we have defined the Hill functions that regulate the switching between the long and short range f functions defined in Eqs. (11) and (12) following Ref. 3. In this work we have used the switching distance K = 0.1 au fixing the steepness of the Hill functions to n = 3 (see Eq. (31) of Ref. 3) in order to guarantee a sufficiently smooth transition between the two descriptions.

J. Chem. Phys. 140, 054102 (2014)

As discussed in Sec. II C the ACS densities have a low dependency on the parameter l that appears in Eq. (12), anyway this dependency is attenuated when accurate wave functions are used and an extensive/sufficiently large VMC sampling is performed. This parameter affects the exponential decay of the ACS density, and in principle should affect the value of the HOMO eigenvalue, which is related to the ionization energy. We have observed that this dependency is more evident for the lighter He atom, but not dramatic when also the tails of the electronic distribution are sufficiently sampled, as in our case. For this reason, we have chosen the l parameter’s values that yield the smaller fluctuations in the density tail: we have found that the best values were l = 2.0, l = 1.4, and l = 2.5, for the He, Be, and Ne atoms, respectively. Since the computational cost of the ACS estimator depends both on the number of the sampling points and on the total number of grid points, it becomes quite expensive in the calculation of densities on 3D grids. In the case of the molecular systems, the ACS estimator has been used to evaluate densities only in the regions around the nuclei, overcoming the difficulty of the traditional estimator in describing the nuclear cusp region.40, 59 Restricting the use of such estimator in the region close to the nuclei has also the additional advantage of avoiding the dependence of the calculated densities at large distance from the nuclei on the parameter l. As the nuclear region has the highest probability of finding an electronic configuration, only 2.048 × 108 VMC samplings per electron were used to obtain these core densities. The full CISD densities were calculated with the Gaussian 0960 computer code. Very large basis sets were used in order to have well converged densities. In particular we used the fully uncontracted aug-cc-pV5Z basis set for the He atom, while for Be and Ne we used the fully uncontracted UGBS1V basis. The two correlated consistent basis set were reduced neglecting the g and h shells. Concerning the molecular systems, we used the aug-cc-pV5Z for the H2 molecule, the augcc-pVQZ for the Be2 dimer and the aug-cc-pVTZ for the H2 O and C2 H4 molecules. B. Details on the calculations for the KS orbitals and νxc potentials

The method illustrated in Sec. II D has been implemented in the development version 2.5 of the CP2K code (starting from the revision number 13384).41, 61 Equation (19) is solved using a dual basis of atom centred Gaussian orbitals and plane waves (GPW)41, 61, 62 where the wave functions and the Hamiltonian matrix are expanded in Gaussian-type orbitals while the density is also represented in a plane wave basis. The in is calculated on a real-space put target density used in vconstr cubic grid defined by the plane wave cut-off, in this way, the

terms of Eq. (19) are calHartree, Fermi-Amaldi, and vconstr culated in reciprocal space via Fast Fourier Transform. For the calculations of the atoms we have also used the atomic CP2K code (only Gaussian-type basis set) adopting a radial grid built of 400 points distributed logarithmically from 0 to 35 bohrs from the nucleus. In order to evaluate the Kohn-Sham orbitals, energies and the corresponding exchange-correlation potential, Eq. (19) is


054102-6


solved for different values of . All the quantities of interests can be then extrapolated to → ∞ as a power series of 1/ or alternatively, as the vxc is the only quantity that depends on , it is possible to extrapolate Eq. (20) and solve the KS equation with the obtained potential. To reproduce the CISD atomic densities within the ZMP procedure implemented in the CP2K atomic code we have used the same basis sets used to create those densities (Sec. III A). On the other hand, for the ACS densities of He and Be we have used a totally uncontracted basis set of 40 primitive Gaussians for the s and p shells. For the Ne atom we have used a basis set of 80 primitive Gaussians for s and p shells. The extrapolation in of Eq. (20) using the CISD densities, was done considering values of in the interval between 100 and 1000, and at = 1000, we observed that the vxc (r) was converged. For the ACS densities, due to the previous considerations on the incompleteness of the basis set, we observed that non-physical fluctuations appeared in vxc (r) for

> 300 and for this reason we performed the extrapolation considering the interval between = 50 and = 300. The ZMP calculations on cubic grids were done using the CP2K code with the GPW description. The cell size of the cubic grids has been chosen requiring that the charge density at the edge of the cell was of the order of 10−10 a.u.: we used cubic cells of sides 5, 10, and 7 bohrs for He, Be, and Ne, respectively, and we found that grid spacings of 0.026, 0.031, and 0.019 bohrs corresponding to plane wave cut-offs of 1800, 1300, 3600 hartree were enough to get converged results. For the H2 , Be2 , H2 O, and C2 H4 molecular systems we have used the uncontracted correlation consistent basis sets: the H, O, and C atoms where described through the uncontraction of the aug-cc-pV5Z basis set, while for Be we have used the uncontracted aug-cc-pVQZ. To avoid non-physical fluctuations in the exchange correlation potential33 the ZMP quantities for the molecules were extrapolated through = {50, 100, 150, 200, 250, 300}.

J. Chem. Phys. 140, 054102 (2014) TABLE I. Ground state and ionization potentials of He, Be, and Ne. The VMC ionization energies are calculated as total energy differences EN − EN −1 . J3, 4

VMC aDZ VMC aDZ VMC aTZ VMC aTZ VMCa CISD SEb VMC aDZ VMC aDZ VMC aTZ VMC aTZ VMCc VMCd VMC(SD)e VMC(FVCAS)e CISD SEb VMC aDZ VMC aDZ VMC aDZ VMC aTZ VMC aTZ VMC aTZ VMC UGBS1V CISD (CBS)f VMCa VMCc VMCd VMC(SD)e CISD SEb

E (hartree)

Helium atom (3s2p)/[2s1p] − 2.903667(2) (3s2p)/[2s2p] − 2.903669(2) (3s2p)/[2s1p] − 2.903700(1) (3s2p)/[2s2p] − 2.903711(1) − 2.903527(9) − 2.9031007044 − 2.903724 Beryllium atom (3s2p)/[2s1p] − 14.66583(1) (3s2p)/[2s2p] − 14.66609(1) (3s2p)/[2s1p] − 14.66607(1) (3s2p)/[2s2p] − 14.66622(1) − 14.66681(1)

IP (hartree)

0.903667(2) 0.903669(2) 0.903700(1) 0.903711(1)

0.9037 0.34328(1) 0.34295(1) 0.34208(1) 0.34182(1) 0.3424(2)

− 14.64972(5) − 14.66668(5) − 14.661478521 − 14.66736 Neon atom (3s2p)/[2s1p] − 128.8703(1) (3s2p)/[2s2p] − 128.8919(1) (3s2p1d) − 128.8930(1) (3s2p)/[2s1p] − 128.8764(1) (3s2p)/[2s2p] − 128.8931(1) (3s2p1d) − 128.8942(1) (3s2p1d) − 128.89437(6)

0.3424 0.7862(1) 0.7937(1) 0.7949(1) 0.7952(1) 0.7958(1) 0.7942(1)

− 128.91982 − 128.891(5) − 128.9029(3) 0.8008(3) − 128.9057(1) − 128.86173911 − 128.9376

0.7923

a

Calculations from Ref. 63 from numerical orbitals. Estimated exact semi-empirical results from Refs. 64 and 65, based on optical spectroscopy data from Ref. 66. c Calculations from Ref. 67, with configuration state function (CSF) with the addition of a 4-electron correlation factor in the Jastrow part. d QMC calculations from Ref. 68 with a multiconfigurational wave function for Be and a single determinant wave function for Ne. e Variational results with Single determinant (SD) of full valence complete active space (FVCAS) wave functions with an additional Jastrow factor, from Ref. 69. f Calculations from Ref. 70. b

IV. RESULTS AND DISCUSSION

As anticipated above, to study the feasibility of using QMC densities for extracting KS quantities of molecular systems, we have first applied the ZMP procedure to the He, Be, and Ne atomic densities calculated on a radial logarithmic grid. Afterwards, the obtained results have been compared with those obtained on a cubic uniform three dimensional grid, to test the procedure that will be used for the molecular systems. Finally, KS orbitals, energies and potentials based on QMC densities are calculated for selected molecules, and in order to test their quality, these have been compared with the ones obtained using full CISD target densities and with other results present in the literature.

A. Assessment of the QMC wave functions

The first important step to achieve very accurate QMC densities to be used as targets in the ZMP procedure, is the optimization of the variational wave functions through which

those densities are obtained. In order to verify the convergence of the total energies and electronic properties of the atoms and molecules considered, all the variational parameters of the JAGP wave functions, including the coefficients and exponents of the atomic orbitals, were fully optimized, as described in Sec. III A. The total energies of the all-electron atoms for different fermionic and Jastrow basis sets are reported in Table I. As expected for the two lighter He and Be atoms the JAGP gives total energies in perfect agreement with experimental results within errors of 10−5 and 10−3 hartree, respectively. Moreover, the calculated total energies and ionization potentials are in excellent agreement with previous VMC


054102-7


J. Chem. Phys. 140, 054102 (2014)

TABLE II. Total energies and electronic properties of C2 H4 , H2 O, Be2 , and H2 . J3, 4

E (hartree)

C(4s4p)/[2s2p],H(3s3p)/[1s1p]

− 78.53139(6)

VMC aDZc VMChyb c MRSD-CId FCIb Expt.

O(4s3p1d),H(2s1p) O(4s3p1d),H(2s1p)

− 76.40588(5) − 76.40741(2) − 76.3963 − 76.431(3) − 76.438(3)e

VMC aDZ

(3s2p) (3s2p1d) (3s2p) (3s2p1d)

Dipole (D)

Qxx (D · Å)

Qyy (D · Å)

Qzz (D · Å)

C2 H4 VMCa Expt.b

1.6983(3) 1.67

1.6031(3) 1.48

− 3.3014(3) − 3.16

1.8832(2) 1.8894(1) 1.870

2.593(1) 2.5875(1) 2.5556

− 0.160(1) − 0.1466(1)

− 2.433(1) − 2.4409(1)

1.8546(6)f

2.63(2)g

− 0.13(3)g

− 2.50(2)g

− 4.024(1) − 4.169(1) − 4.852(1) − 4.363(1)

2.012(1) 2.085(1) 2.426(1) 2.181(1)

2.012(1) 2.085(1) 2.426(1) 2.181(1)

− 0.3078(1) − 0.3054(1) − 0.3064(1)

− 0.3078(1) − 0.3054(1) − 0.3064(1)

H2 O

VMC aTZ

VMC DZ VMC TZ VMC aTZ CISD-CBSj Expt.k

(2s1p) (2s1p) (2s1p)

VMC TZ VMC aTZ

(2s1p) (2s1p)

Be2 (RBe2 = 4.6108 bohrs)h − 29.31971(2) − 29.32155(2) − 29.32091(2) − 29.32213(2) H2 (RHH = 1.401 bohrs)i − 1.174297(2) − 1.174438(1) − 1.174463(1) − 1.1742

0.6156(1) 0.6108(1) 0.6127(1) 0.6123 0.619(28)

H2 (RHH = 3.000 bohrs) − 1.057268(1) − 1.057309(1)

a

Results from the AGPn* wave function optimized in Ref. 5. From Ref. 76. c From Ref. 58. d From Ref. 77. e Reference energy considered exact, from semi-empirical calculations of Ref. 78. f From Ref. 79. g From Ref. 80. h From Ref. 73. i From Ref. 74. j Extrapolations to the CBS limit from Ref. 81. k From Ref. 75. b

results.67, 69 For the heavier Ne atom, although we observe a larger discrepancy with respect to the experimental results, our VMC energies are again in good agreement with accurate VMC results from Refs. 63, 67, and 69. Also compared to a full CISD calculation70 in the complete basis set limit, our energy differs only by 0.025 hartree, despite the fact that we used a relatively small basis set. From our calculation we observe that the most relevant improvement of the wave function comes from the use of at least two p-orbitals in the Jastrow factor. The inclusion of the d-orbital in the Jastrow part gives a smaller improvement in the total energy. For the calculation of the densities we have used the largest basis sets shown in Table I, confident of the good quality of the adopted wave functions based on the comparison with the reported experimental values and accurate calculations present in the literature. The same procedure has been applied for the H2 and Be2 molecules in order to obtain accurate QMC densities following the scheme reported in Sec. III A. In Table II the values of the total energies and of the electronic properties are reported for each molecule with different basis sets. On the basis of these results we have chosen to apply the ZMP procedure to the Quantum Monte Carlo densities obtained from the largest basis sets reported in Table II for each molecule (Sec. III A).

For the ethylene molecule we have used the densities obtained from the biggest optimized wave function used in Ref. 5, while for the water molecule we have used two basis sets, both proposed in Ref. 58 and described in Sec. III A, in order to test the effect of the quality of the wave function on the calculated density. The electronic properties of the water monomer obtained with these two wave functions are reported in Table II are very similar and in particular the dipole moment, presents a discrepancy of 0.03[D] from the experimental values. Further considerations are reported in Ref. 58. If we compare the quadrupole moments, we can see that the basis set build with hybrid orbitals shows a better compatibility with the experimental values. This is due to the bigger basis set used to build this hybrid AGP wave function. B. Kohn-Sham orbitals from radial atomic densities

The first application of our implementation has been the evaluation of KS quantities (orbital energies, exchangecorrelation energies, and potentials) of the three closed shell atoms (He, Be, Ne) from QMC electronic densities, calculated on a radial logarithmic grid defined in Sec. III B. Radial logarithmic grids have a very fine spacing near the nuclear cusp, and the best approach to obtain smooth


054102-8


J. Chem. Phys. 140, 054102 (2014)

TABLE III. Exchange-Correlation and Kohn-Sham orbital energies from radial atomic densities. All energies reported in hartree. ρ(r)

Exc

Refs. CIc CISD ACS

− 1.0685a − 1.068 − 1.066 − 1.067

Refs. CIc CISD ACS

− 2.7614a − 2.776 − 2.764 − 2.766

1s

2s

2p

Helium atom − 0.9037b − 0.9039 − 0.9039 − 0.9052 Beryllium atom − 4.2142 − 4.2091 − 4.1546

− 0.3424b − 0.3384 − 0.3413 − 0.3453

Neon atom Refs. CIc CISD BDd ACS

− 12.498a − 12.49 − 12.40 − 12.436

− 30.82 − 30.824 − 30.881 − 30.708

− 1.654 − 1.650 − 1.6277 − 1.650

− 0.7923b − 0.797 − 0.794 − 0.7683 − 0.792

a

From Ref. 71: exchange-correlation energies from accurate ab initio calculations and HF theory. b Spectroscopic results for Ref. 65. c From Refs. 27 and 72; KS quantities extracted from CI densities. d From Ref. 33: Brueckner doubles calculated with a reduced cc-pVQZ basis set with the f and g orbitals removed. = 900.

QMC densities in this region, is to apply the ACS estimator described in Sec. II C. This estimator has the great advantage of yielding accurate descriptions of the charge density near the nuclear cusp, provided that accurate wave functions are employed. The good quality of the used wave functions can be seen comparing the contact densities, i.e., the values of the electronic densities at the nuclear cusps point, obtained through the ACS density estimator, with the ones calculated at CISD level. Thanks to the presence of the onebody Jastrow factor in the JAGP the real behaviour of the wave function near the nuclear cusp is ensured and we find contact densities in very good agreement with those obtained through the zero-variance zero-bias Diffusion Monte Carlo (DMC) estimator:59 for example, for the Ne atom the contact density was estimated to be 620.550(69) au compatible with our ACS value of 620.44(18) a.u., while the CISD value is of 614.3 a.u. The KS eigenvalues and the exchange correlation energies for the three atoms, calculated from the ACS target densities, are reported in Table III and compared with those obtained from our CISD densities and with other results present in the literature. Experimental ionization energies and exchange correlation energies calculated through high level quantum chemistry methods are also reported as a reference. The ACS densities of the three atoms, to which we have applied the ZMP procedure, were obtained through a VMC sampling on the JAGP wave functions with the largest basis sets (see Table I). The CISD and QMC total energies where used to estimate the exchange-correlation energy through Eq. (21). In order to extract orbitals and potentials we have solved Eq. (19) for different values of and the exchangecorrelation potential (Eq. (20)) was then extrapolated using a polynomial fit. We have also verified that the orbital eigenval-

ues and exchange-correlation energies obtained form the extrapolated vxc (r) potential provided the same values of those extrapolated at different ’s, proof of the robustness of our results. The eigenvalues reported in Table III obtained from the ACS densities are in good agreement with those obtained form our CISD calculations and from those reported in the pioneering works of Zhao, Morrison, and Parr.27, 72 In particular for the last occupied KS eigenvalues, obtained using the ACS densities, a good agreement with the experimental ionization energy for the He, Be, and Ne atoms is observed, with a discrepancy of no more than 0.003 hartree. This result is promising, and shows that even at the variational level, the QMC densities calculated from optimized wave functions have the same quality as those obtained by CISD calculations. The ACS and CISD densities for the three atoms are shown in the left panel of Fig. 1. In the inset the densities are plotted in logarithmic scale to highlight the error bars and the small fluctuations of the ACS densities. A delicate aspect of the ZMP method, that we accurately investigated, is the finiteness of the basis set used in Eq. (19) to build the KS density (Eq. (14)) that targets the accurate ab initio one. The KS density tends to the target density in the limit → ∞, and to achieve this limit the basis set used in the solution of Eq. (19) must be sufficiently large as to correctly reproduce the target ρ 0 (r).33 When dealing with CISD target densities we have verified that the basis set must be at least the same used to build the ab initio ρ 0 (r). For the Ne atom at the CISD level, using the UGBS1V basis set in the ZMP procedure, the maximum difference ρ(r) − ρ 0 (r) in the limit of → ∞ is equal to 0.04 a.u. at the cusp where the CISD density is 614.3 a.u. i.e., an error of the 0.006%. For the ACS density of Ne, the extrapolated difference ρ(r) − ρ 0 (r) grows up to 0.26 a.u. at the cusp, which is of 620.44(18) a.u. with an error of about the 0.04%: this behaviour is expected because the nuclear cusp described through the one body Jastrow factor is difficult to reproduce with a basis set built of primitive Gaussians. Similar energy differences are found for both the He and Be atoms: for the former the CISD and ACS differences are less then 0.01% while for the latter they are, respectively, of 0.02% and 0.05%. Through this approach we have obtained the vxc (r) potential curves, shown in the right column of Fig. 1: in the inset the vxc (r) multiplied by r is also shown to highlight the correct asymptotic behaviour (−1/r) of the exchange correlation potential. The small discrepancies that appear in the potentials obtained from the ACS densities, with respect to those from CISD ones, are due to basis set incompleteness. As grows, the ZMP procedure reduces the density difference with respect the target density in the cusp region, spoiling the middle range regions in which the density is already of the order of 10−5 a.u. This behaviour was not registered when using large basis sets for the CISD densities.The quality of our results can be also seen by comparing the exchange correlation energies (Eq. (21)) obtained through the ZMP procedure, and the values estimated from direct ab initio calculations from Ref. 71. The exchange-correlation energies used as reference values in Table III and reported in Ref. 71 are obtained from accurate Hartree-Fock and CISD total energy calculations used to identify the Correlation and


054102-9


J. Chem. Phys. 140, 054102 (2014)

FIG. 1. Densities (left panels) and exchange-correlation potentials (right panels) for He, Be, and Ne atoms. The densities are plotted multiplied by the square of the radius r2 ρ(r), while in the insets ρ(r) is plotted in logarithmic scale. The errors of the ACS densities have been reported only in the insets. In the insets of the right panels the rvxc (r) is shown in order to highlight the correct long-range asymptotic behaviour of the potentials.

Exchange components which where afterwards summed together. The results obtained for the ACS densities seem to be in better agreement with this reference, at least for the He and Be atoms if compared with the values obtained from CISD densities, with an error for the latter of about 5 mhartree. We

must point out that these values are only indicative, as a matter of fact through Eq. (21) we can see that they directly depend on the total energy value corresponding to that particular density which may be quite distant from convergence, especially for CISD (Table I).


054102-10


J. Chem. Phys. 140, 054102 (2014)

TABLE IV. Exchange-Correlation and Kohn-Sham orbital energies from cubic atomic densities. All energies reported in hartree. The atomic densities used for these calculations are the same used for those on radial logarithmic grids (Table III), allocated on uniform 3d cubes. Exc

ρ(r)

1s

2s

CISD ACS

− 1.063 − 1.067

Helium atom − 0.9044 − 0.9049

CISD ACS

− 2.761 − 2.764

Beryllium atom − 4.204 − 4.167

− 0.3411 − 0.3502

CISD ACS

− 12.392 − 12.449

Neon atom − 30.798 − 30.687

− 1.649 − 1.652

2p

− 0.793 − 0.793

C. Kohn-Sham orbitals and potentials from densities on uniform cubic grids

To investigate the possibility of applying the ZMP procedure to molecular densities calculated through QMC methods we have repeated the calculations of the three atomic systems using densities distributed on uniform 3d cubic grids. The exchange correlation energies and the KS eigenvalues obtained for the ACS and CISD densities on these cubic grids are reported in Table IV. If we compare these results with those presented in Table III and obtained from linear radial densities, we can see that the maximum deviation is about 5 mhartree on both eigenvalues and exchange correlation energies, validating this approach and opening the possibility of the application to molecular systems.

D. Kohn-Sham orbitals and potentials for molecular systems

The successful application of the ZMP procedure on the 3d atomic densities obtained by QMC, gives us the possibility to further investigate its applicability on densities of small molecules. For the first time we use the ZMP method to extract KS quantities from many-body QMC molecular densities of particularly interesting chemical systems. First, we have studied the triatomic H2 O water molecule, the C2 H4 molecule, and the Be2 dimer, a prototype of weakly bonded molecules. The obtained results have been compared with experimental ionization energies and with other ZMP calculations on densities obtained by the Brueckner Doubles (BD) method33 in the case of water. Next we have considered the H2 molecule at both the equilibrium distance of R0 = 1.401 bohrs and in the dissociation range (R = 3.000 bohrs).83, 84 The eigenvalues obtained for the C2 H4 , H2 O, and Be2 molecules are reported in Table V, compared, when possible, with the experimental ionization potentials (IP), CISD and DFT calculations, the latter with the local PW92 potential, and with other ZMP calculations present in literature.33 First of all let us consider the results for the H2 O molecule, for which two different basis sets have been used. The traditional VMC estimator for the density has been used only for the smaller AGP basis set so that we refer to it as VMC. We can see that the eigenvalues obtained through the traditional VMC estimator are more compatible with the ZMP values obtained from BD densities33 (indicated as ZMPBD ), with respect to our CISD calculations. For the latter the discrepancy between the highest eigenvalue and the ionization energy is of 1.3 eV while for VMC this difference reduces to 0.7 eV

TABLE V. Kohn-Sham orbital energies from molecular densities on cubic grids. All energies are reported in eV. 1

2

3

4

ZMPCISD ZMPV MC ZMPV MC+ACS PW92b

− 275.497 − 272.245 − 276.626 − 265.996

− 275.497 − 272.245 − 276.626 − 265.996

C2 H4 (IP = −10.68a ) − 22.951 − 18.317 − 23.466 − 18.863 − 23.187 − 18.648 − 18.659 − 14.133

ZMPBD c ZMPCISD ZMPV MC ZMPV MC+ACS hyb ZMPV MC+ACS PW92b

− 518.104 − 517.207 − 513.805 − 522.567 − 522.720 − 506.273

− 30.580 − 32.137 − 31.320 − 31.238 − 30.579 − 25.183

H2 O (IP = −12.62a ) − 18.161 − 14.427 − 19.728 − 16.000 − 19.211 − 15.429 − 19.021 − 15.238 − 18.350 − 14.562 − 13.246 − 9.374

ZMPCISD ZMPV MC ZMPV MC+ACS PW92b

− 112.836 − 113.553 − 114.477 − 104.699

− 112.836 − 113.553 − 114.477 − 104.699

Be2 (IP = −7.418d ) − 10.664 − 7.992 − 10.572 − 7.898 − 10.563 − 7.891 − 7.239 − 4.537

5

6

7

8

− 15.519 − 15.983 − 15.878 − 11.510

− 14.213 − 14.638 − 14.550 − 10.195

− 12.588 − 13.053 − 12.950 − 8.454

− 10.639 − 10.982 − 10.981 − 6.867

− 12.365 − 13.905 − 13.333 − 13.089 − 12.499 − 7.366

a

From Ref. 76. This work: DFT calculations with aug-cc-pV5Z basis sets. Extracted from Brueckner doubles relaxed density with = 900 from Ref. 33. d From Ref. 82. b c


054102-11


even though the tail of the charge distribution presents error fluctuations. These discrepancies in the higher eigenvalues obtained from the CISD and BD densities can be ascribed to the fact that the aug-cc-pVTZ basis set is not enough to obtain a fully converged density in the tails. Larger discrepancies with ZMPBD are found for the VMC results, in the eigenvalues associated with the more internal core orbitals, which seem on the contrary well described by the CISD densities. This is due to the bad behaviour of the traditional Monte Carlo estimator for the electronic density in proximity of the nucleus where the density varies rapidly, as discussed in Sec. II C. To cure this behaviour we have built up a new density treating the region near the cusps through the ACS estimator, so to correctly reproduce this small region with a relatively small computational cost. In the case of water, the ACS estimator has been applied in a region of 0.2 bohrs around the nuclei. Repeating the ZMP calculations on these new densities we have obtained the results labelled with ZMPV MC+ACS in Table V for the smaller basis set. We can see in this case a big change in the innermost eigenvalue. The core correction of the VMC density also corrects, in a slighter way, the estimation of the valence eigenvalues: the discrepancy between the experimental ionization potential and the ZMP estimations, for example, reduces to 0.5 eV. Moreover, the biggest hybrid basis set was used to obtain a mixed VMC+ACShyb density and the ZMP eigenvalues, labelled hyb with ZMPV MC+ACS , are reported in Table V. We can see that this more accurate wave function leads to better ionization potential with respect to those extracted from CISD and BD densities, with a discrepancy from the experimental value of less then 0.1 eV. Second, we have applied the same procedure to the C2 H4 and Be2 molecules, calculating the KS quantities on both the pure VMC densities obtained with the old estimator, and the mixed VMC+ACS densities. In line with CISD calculations, for the Be2 molecule we found a difference of 0.5 eV between the highest ZMP occupied eigenvalue and the first ionization potential which is the only experimental reference found for this molecule, reported in Table V. For the C2 H4 molecule the ZMPV MC+ACS values of the ionization potential differs of about 0.3 eV from the experimental value, while the ZMPCISD result seems to be nearly exact. The inaccurate description of the density near the nuclear cusp for the VMC traditional estimator, leads to non-physical vxc (r) potentials which reflect in the wrong estimation of the core orbitals. In Fig. 2 the vxc (r) potentials extracted from VMC and VMC+ACS densities along the principal molecular axes are shown, and we can see that the non-physical behaviour due to the underestimation of the density in the case of the VMC calculations is corrected when the ACS density is used near the cusp. In Fig. 3, from left to right, we report the contour plots of the converged vxc (r) potentials extracted from VMC+ACS, CISD densities and the LDA (PW92) potentials of the C2 H4 (upper panels), the H2 O molecule (middle panels), and of the Be2 molecule (lower panels) along the principal molecular plane. Comparing these results we observe similar features between the potentials coming from the CISD and VMC+ACS densities, that are also reflected in similar Kohn-Sham

J. Chem. Phys. 140, 054102 (2014)

FIG. 2. The linear profiles of the vxc (x, y, z) potentials are reported along the principal axes of symmetry for the C2 H4 molecule (upper panel), for the H2 O molecule (centre panel), and for the Be2 molecule along the molecular bond (lower panel). In the insets the region around one of the C atom, for the C2 H4 molecule, the O atom, for the H2 O molecule, and around one of the Be atoms for the Be2 dimer are enlarged. For the C2 H4 and the Be2 molecules we report the vxc profiles extracted from the VMC and VMC+ACS densities. For the H2 O molecule we report the vxc profiles extracted from the VMC and VMC+ACS densities obtained from the first basis set, respectively, ZMPV MC and ZMPV MC+ACS , and that obtained from the VMC+ACS denhyb sity built from the hybrid wave function, ZMPV MC+ACS . Also, for all the molecules, the vxc extracted from CISD densities, and the local PW92 functional are reported for comparison.


054102-12


J. Chem. Phys. 140, 054102 (2014)

FIG. 3. vxc (r) potentials along the molecular planes, for the C2 H4 molecule (top panels), H2 O (middle panel), and for the Be2 molecule (bottom panels). For the two molecules the vxc (r) are extracted with the ZMP procedure from VMC+ACS (left panels) and CISD (centre panels) densities, on the right panel the PW92 local functional is reported. The associated eigenvalues are reported in Table V.


054102-13


FIG. 4. Molecular orbitals obtained with the ZMP procedure applied on the VMC+ACS densities of the C2 H4 (left columns), the H2 O (centre column), and the Be2 (right column) molecules.

J. Chem. Phys. 140, 054102 (2014)

eigenvalues obtained with the two potentials reported in Table V. Substantial differences are found, as expected, with respect to the LDA potential which presents a less structured profile around the nuclei, beside the wrong long range tail. This is also shown in Fig. 2, where the same potentials are drawn along the principal molecular axes, and where the regions with small differences between the potentials from CISD and VMC+ACS densities are highlighted around the nuclei. In Fig. 4 the molecular orbitals, corresponding to the eigenvalues reported in Table V, obtained solving Eq. (19) using the VMC+ACS densities, are shown for the C2 H4 , H2 O, and the Be2 molecules. From the bottom to the top the orbitals are displayed from the most internal to the most external, showing the correct order and symmetry. The last molecular system we have investigated is the H2 molecule at both equilibrium distance and with a stretched bond length of 3 bohrs inspired by the work of Buijse et al.83 which studied the effect of electron correlation in the bond formation of this molecule. The stretched bond length case is

FIG. 5. Density (top panel), effective potential vhxc (r) (middle panel), and exchange-correlation potential (bottom panel) along the bond axis for the H2 molecule at equilibrium distance (left) and at 3 bohrs bonding distance (right). Results from VMC+ACS, CISD, and B3LYP densities are reported. LDA (PW92) potentials and densities are also plotted. Eigenvalues and exchange correlation energies are reported in Table VI.


054102-14


J. Chem. Phys. 140, 054102 (2014)

TABLE VI. Kohn-Sham orbital energy from H2 molecular densities. All energies are reported in eV. H2 (RHH = 1.401 bohrs)

H2 (RHH = 3.0 bohrs)

(IP = −15.42a ) − 16.41 − 16.54 − 16.43 − 16.43 − 11.83 − 10.23

− 13.03 − 13.25 − 13.06 − 13.06 − 8.98 − 7.85

CIb ZMPCISD ZMPV MC ZMPV MC+ACS B3LYPc PW92c a

From Ref. 76. CI calculations from Ref. 83 with a contracted 5s2p1d basis set. c This work: DFT calculations with aug-cc-pV5Z basis sets. b

of particular interest as in this regime the Hartree-Fock approximation becomes inadequate and the correlation part of the exchange-correlation potential starts to play an important role. In Fig. 5 we plot the densities, the vxc (r) and also the total effective potential the vhxc (r) = vh (r) + vxc (r) in order to examine the barrier formed in the KS potential at the centre of the molecular bond.32, 84 The accurate CISD and VMC+ACS derived KS potentials are then compared with the LDA potential (PW92) and with a potential extracted, through the ZMP procedure, from a DFT density obtained with the B3LYP hybrid functional, which contains a fraction of Hartree-Fock exchange. In the top panel of Fig. 5 we can see that at equilibrium distance we found similar profiles for the B3LYP, CISD, and VMC+ACS densities, while the LDA one presents a depletion of the charge in the middle of the bond. The situation is different for the stretched molecule where the B3LYP density noticeably differs from the CI and VMC+ACS density and resembles the LDA one, presenting a greater delocalization in the bond. When looking at the potentials we observe that ZMP ones, extracted from VMC+ACS and accurate CISD densities, have a remarkable correspondence, and the obtained eigenvalues, reported in Table VI, are in good agreement with those calculated by Buijse et al.83 through CISD total energy differences. The ZMP potential based on the B3LYP target density, even if it displays evident discrepancies with respect the potential obtained from VMC+ACS and CISD densities, it exhibits a comparable barrier height at equilibrium distance, while in the stretched case, where the correlation is more important, the agreement is lost, the discrepancy becomes larger, and it presents a more similar shape to the LDA beside the different tail behaviour. V. CONCLUSIONS

In conclusion, we have performed accurate Variational Monte Carlo calculations of atomic densities (He, Be, Ne), and molecular systems: the ethylene, the triatomic water molecule, the Be2 dimer, prototype of weakly bonded molecules, and H2 at bonding and stretched distances as an example of correlated system. We have shown that the total energy calculations are in good agreement with reference results and with those obtained through other high level quantum chemistry calculations. The employment of a compact

JAGP wave function combined with the optimization of all variational parameters including exponents and coefficients of the Gaussian primitives has guaranteed a fast convergence with the basis sets of both the energies and the calculated densities. Using such accurate densities we have shown the possibility to extract Kohn-Sham orbitals and exchange-correlation potentials using the ZMP technique, which is non-trivial for QMC densities because of the presence of the stochastic error. Our results demonstrate the validity of this procedure and, thanks to the fair computational scaling of QMC calculations with the number of electrons, can be extended to larger molecules, where traditional post-Hartree-Fock methods with large basis sets might be computationally too demanding. The possibility to have accurate potentials from many-body densities for extended molecular systems opens different perspectives in computational chemistry and multi-scale modelling. First, it can contribute to the testing and development of new exchange and correlation potentials, in addition, the definition of Kohn-Sham quantities corresponding to QMC densities is the first step for the development of embedding multi-scale methods based on QMC and DFT.

ACKNOWLEDGMENTS

The authors thank Sandro Sorella, Marcella Iannuzzi, Paola Gori-Giorgi, Emanuele Coccia, Andrea Zen, and Daniele Bovi for stimulating discussions. The authors acknowledge funding provided by the European Research Council Project No. 240624 within the VII Framework Program of the European Union. Computational resources were supplied by CINECA, PRACE infrastructure, and the Caliban-HPC centre at the University of L’Aquila. 1 M.

Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001); J. Kolorenˇc, and L. Mitas, Rep. Prog. Phys. 74, 026502 (2011). 2 M. Caffarel, M. Rérat, and C. Pouchan, Phys. Rev. A 47, 3704 (1993). 3 E. Coccia, O. Chernomor, M. Barborini, S. Sorella, and L. Guidoni, J. Chem. Theory Comput. 8, 1952 (2012). 4 O. Valsson and C. Filippi, J. Chem. Theory Comput. 6, 1275 (2010); O. Valsson, C. Angeli, and C. Filippi, Phys. Chem. Chem. Phys. 14, 11015 (2012); C. Filippi, F. Buda, L. Guidoni, and A. Sinicropi, J. Chem. Theory Comput. 8, 112 (2012). 5 M. Barborini, S. Sorella, and L. Guidoni, J. Chem. Theory Comput. 8, 1260 (2012). 6 R. Assaraf and M. Caffarel, J. Chem. Phys. 119, 10536 (2003). 7 S. Chiesa, D. M. Ceperley, and S. Zhang, Phys. Rev. Lett. 94, 036404 (2005). 8 M. Marchi, S. Azadi, C. Casula, and S. Sorella, J. Chem. Phys. 131, 154116 (2009). 9 M. Casula, M. Marchi, S. Azadi, and S. Sorella, Chem. Phys. Lett. 477, 255 (2009). 10 A. Zen, D. Zhelyazov, and L. Guidoni, J. Chem. Theory Comput. 8, 4204 (2012). 11 C. J. Umrigar, Int. J. Quantum Chem. 36, 217 (1989). 12 C. Filippi and C. J. Umrigar, Phys. Rev. B 61, R16291 (2000). 13 L. K. Wagner and J. C. Grossman, Phys. Rev. Lett. 104, 210201 (2010). 14 S. Sorella and S. Capriotti, J. Chem. Phys. 133, 234111 (2010). 15 E. Coccia and L. Guidoni, J. Comput. Chem. 33, 2332 (2012). 16 E. Coccia, D. Varsano, and L. Guidoni, J. Chem. Theory Comput. 9, 8 (2013). 17 M. Barborini and L. Guidoni, J. Chem. Phys. 137, 224309 (2012). 18 S. Saccani, C. Filippi, and S. Moroni, J. Chem. Phys. 138, 084109 (2013). 19 G. Mazzola, A. Zen, and S. Sorella, J. Chem. Phys. 137, 134112 (2012).


054102-15 20 P.


Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964); W. Kohn and L. Sham, ibid. 140, A1133 (1965). 21 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 22 C. Attaccalite, S. Moroni, P. Gori-Giorgi, and G. Bachelet, Phys. Rev. Lett. 88, 256601 (2002); D. Varsano, S. Moroni, and G. Senatore, Europhys. Lett. 53, 348 (2001). 23 M. Levy, Proc. Natl. Acad. Sci. U.S.A. 76, 6062 (1979). 24 C. O. Almbladh and A. C. Pedroza, Phys. Rev. A 29, 2322 (1984). 25 A. Görling, Phys. Rev. A 46, 3753 (1992). 26 Y. Wang and R. G. Parr, Phys. Rev. A 47, R1591 (1993). 27 Q. Zhao, R. C. Morrison, and R. G. Parr, Phys. Rev. A 50, 2138 (1994). 28 R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994). 29 Q. Wu and W. Yang, J. Chem. Phys. 118, 2498 (2003). 30 K. Peirs, D. Van Neck, and M. Waroquier, Phys. Rev. A 67, 012505 (2003). 31 I. G. Ryabinkin and V. N. Staroverov, J. Chem. Phys. 137, 164113 (2012). 32 O. V. Gritsenko, R. van Leeuwen, and E. J. Baerends, Phys. Rev. A 52, 1870 (1995); J. Chem. Phys. 104, 8535 (1996); P. R. T. Schipper, O. V. Gritsenko, and E. J. Baerends, Theor. Chem. Acc. 98, 16 (1997); 99, 329 (1998); J. Chem. Phys. 111, 4056 (1999). 33 D. J. Tozer, V. E. Ingamells, and N. C. Handy, J. Chem. Phys. 105, 9200 (1996); V. E. Ingamells and N. C. Handy, Chem. Phys. Lett. 248, 373 (1996). 34 M. E. Mura, P. J. Knowles, and C. A. Reynolds, J. Chem. Phys. 106, 9659 (1997). 35 M. J. Allen and D. J. Tozer, J. Chem. Phys. 117, 11113 (2002). 36 C. R. Jacob, J. Chem. Phys. 135, 244102 (2011). 37 K. Boguslawski, C. R. Jacob, and M. Reiher, J. Chem. Phys. 138, 044111 (2013). 38 C. Filippi, X. Gonze, and C. J. Umrigar, “Generalized-gradient approximation to density-functional theory: comparison with exact results,” in Recent Developments and Applications of Density Functional Theory, edited by J. M. Seminario (Elsevier, Amsterdam, 1996). 39 M. Casula and S. Sorella, J. Chem. Phys. 119, 6500 (2003); M. Casula, C. Attaccalite, and S. Sorella, ibid. 121, 7110 (2004). 40 R. Assaraf, M. Caffarel, and A. Scemama, Phys. Rev. E 75, 035701 (2007). 41 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103 (2005). 42 O. Roncero, M. P. de Lara-Castells, P. Villarreal, F. Flores, J. Ortega, M. Paniagua, and A. Aguado, J. Chem. Phys. 129, 184104 (2008); O. Roncero, A. Zanchet, P. Villarreal, and A. Aguado, ibid. 131, 234110 (2009); J. D. Goodpaster, N. Ananth, F. R. Manby, and T. F. Miller III, ibid. 133, 084103 (2010); S. Fux, C. R. Jacob, J. Neugebauer, L. Visscher, and M. Reiher, ibid. 132, 164101 (2010). 43 N. Govind, Y. A. Wang, and E. A. Carter, J. Chem. Phys. 110, 7677 (1999); P. Huang and E. A. Carter, ibid. 125, 084102 (2006); C. Huang, M. Pavone, and E. A. Carter, ibid. 134, 154110 (2011); A. Severo Pereira Gomes and C. R. Jacob, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 108, 222 (2012). 44 S. Sorella, M. Casula, and D. Rocca, J. Chem. Phys. 127, 014105 (2007). 45 C. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and H. Rhenning, Phys. Rev. Lett. 98, 110201 (2007). 46 A. G. Anderson and W. A. Goddard, J. Chem. Phys. 132, 164110 (2010). 47 B. Braïda, J. Toulouse, M. Caffarel, and C. J. Umrigar, J. Chem. Phys. 134, 084108 (2011). 48 F. Fracchia, C. Filippi, and C. Amovilli, J. Chem. Theory Comput. 8, 1943– 1951 (2012). 49 L. Pauling, The Nature of the Chemical Bond, 3rd ed. (Cornell University Press, Itaca, New York, 1960), pp. 230–240.

J. Chem. Phys. 140, 054102 (2014) 50 J.

A. Pople, Proc. R. Soc. London, Ser. A 202, 323 (1950); A. C. Hurley, J. E. Lennard-Jones, and J. A. Pople, ibid. 220, 446 (1953). 51 M. Marchi, S. Azadi, and S. Sorella, Phys. Rev. Lett. 107, 086807 (2011). 52 F. Sterpone, L. Spanu, L. Ferraro, S. Sorella, and L. Guidoni, J. Chem. Theory Comput. 4, 1428 (2008). 53 Q. Zhao and R. G. Parr, J. Chem. Phys. 98, 543 (1993). 54 M. Levy, and J. P. Perdew, in Density Functional Methods in Physics, NATO ASI Series B123, edited by R. M. Dreizler and J. da Providencia (Plenum, New York, 1985), pp. 11–30. 55 E. Fermi and E. Amaldi, Mem. R. Ital. Acc. 6, 119 (1934). 56 R. G. Parr and S. K. Ghosh, Phys. Rev. A 51, 3564 (1995). 57 S. Sorella, “Turborvb Quantum Monte Carlo Package” (accessed date December 12, 2012). 58 A. Zen, Y. Luo, S. Sorella, and L. Guidoni, J. Chem. Theory Comput. 9, 4332 (2013). 59 M. C. Per, I. K. Snook, and S. P. Russo, J. Chem. Phys. 135, 134112 (2011). 60 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian 09, Revision A.1, Gaussian, Inc., Wallingford, CT, 2009. 61 See www.cp2k.org for “CP2K Open Source - Molecular Dynamics” (accessed date January 8, 2014). 62 G. Lippert, J. Hutter, and M. Parrinello, Mol. Phys. 92, 477 (1997). 63 A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 71, 066704 (2005). 64 E. R. Davidson, S. A. Hagstrom, S. J. Chakravorty, V. M. Umar, and C. F. Fischer, Phys. Rev. A 44, 7071 (1991). 65 S. J. Chakravorty, S. R. Gwaltney, E. R. Davidson, F. A. Parpia, and C. F. Fischer, Phys. Rev. A 47, 3649 (1993). 66 C. E. Moore, Natl. Stand. Ref. Data Ser. (U.S. Natl. Bur. Stand.) 1, 34 (1970). 67 C.-J. Huang, C. J. Umrigar, and M. P. Nightingale, J. Chem. Phys. 107, 3007 (1997). 68 P. Maldonado, A. Sarsa, E. Buendía, and F. J. Gálvez, J. Chem. Phys. 133, 064102 (2010). 69 J. Toulouse and C. J. Umrigar, J. Chem. Phys. 128, 174101 (2008). 70 D. E. Woon and J. Thom H. Dunning, J. Chem. Phys. 103, 4572 (1995). 71 C. Lee and R. G. Parr, Phys. Rev. A 42, 193 (1990). 72 R. C. Morrison and Q. Zhao, Phys. Rev. A 51, 1980 (1995). 73 J. M. Martin, Chem. Phys. Lett. 303, 399 (1999). 74 K. Huber and G. Herzberg, Molecular Spectra and Molecular Structure. IV. Constants of Diatomic Molecules (Van Nostrand Reinhold Co., 1979). 75 A. Buckingham and J. Cordle, Mol. Phys. 28, 1037 (1974). 76 NIST, National Institute of Standards and Technology, Computational Chemistry Comparison and Benchmark Database. 77 D. Feller, C. M. Boyle, and E. R. Davidson, J. Chem. Phys. 86, 3424 (1987). 78 A. Lüchow, J. B. Anderson, and D. Feller, J. Chem. Phys. 106, 7706 (1997). 79 S. A. Clough, Y. Beers, G. P. Klein, and L. S. Rothman, J. Chem. Phys. 59, 2254 (1973). 80 S. G. Kukolich, J. Chem. Phys. 50, 3751 (1969). 81 D. B. Lawson and J. F. Harrison, J. Phys. Chem. A 101, 4781 (1997). 82 J. M. Merritt, A. L. Kaledin, V. E. Bondybey, and M. C. Heaven, Phys. Chem. Chem. Phys. 10, 4006 (2008). 83 M. A. Buijse, E. J. Baerends, and J. G. Snijders, Phys. Rev. A 40, 4190 (1989). 84 N. Helbig, I. V. Tokatly, and A. Rubio, J. Chem. Phys. 131, 224105 (2009).


Quantum Gibbs ensemble Monte Carlo.

Quantum speedup of Monte Carlo methods.

Hyperon puzzle: hints from quantum Monte Carlo calculations.

A pure-sampling quantum Monte Carlo algorithm.

The liquid structure of tetrachloroethene: molecular dynamics simulations and reverse Monte Carlo modeling with interatomic potentials.

MM) Monte Carlo Techniques for Modeling Condensed-Phase Reactions.

Quantum Monte Carlo Study of the Protonated Water Dimer.

Quantum Monte Carlo study of a resonant Bose-Fermi mixture.

Valence-bond quantum Monte Carlo algorithms defined on trees.

Ab initio molecular dynamics simulation of liquid water by quantum Monte Carlo.

Unbiased reduced density matrices and electronic properties from full configuration interaction quantum Monte Carlo.

Structural stability and defect energetics of ZnO from diffusion quantum Monte Carlo.

Monte Carlo modeling of the dual-mode regime in quantum-well and quantum-dot semiconductor lasers.

Linear-scaling quantum calculations using non-orthogonal localized molecular orbitals.

Ab initio molecular dynamics with noisy forces: validating the quantum Monte Carlo approach with benchmark calculations of molecular vibrational properties.

Excited states from quantum Monte Carlo in the basis of Slater determinants.

Efficiency in nonequilibrium molecular dynamics Monte Carlo simulations.

Wormhole Hamiltonian Monte Carlo.

Molecular Mechanics.

Learning reduced kinetic Monte Carlo models of complex chemistry from molecular dynamics.

π-Conjugation in trans-1,3-butadiene: static and dynamical electronic correlations described through quantum Monte Carlo.

Communication: fixed-node errors in quantum Monte Carlo: interplay of electron density and node nonlinearities.

Concurrent Monte Carlo transport and fluence optimization with fluence adjusting scalable transport Monte Carlo.

Optimized molecular reconstruction procedure combining hybrid reverse Monte Carlo and molecular dynamics.