Accurate, efficient and simple forces with Quantum Monte Carlo methods

0.17 MB·English

by  Simone Chiesa

Accurate, efficient and simple forces with Quantum Monte Carlo methods
Simone Chiesa∗ and D. M. Ceperley†
Dept. of Physics, University of Illinois Urbana-Champaign, Urbana, IL 61801
Shiwei Zhang‡
Dept. of Physics, College of William and Mary, Williamsburg, VA 23187
5
0 Computation ofionic forces usingquantumMonteCarlo methodshaslong beena challenge. We
0 introduce a simple procedure, based on known properties of physical electronic densities, to make
2 thevarianceoftheHellmann-Feynmanestimatorfinite. Weobtainveryaccurategeometriesforthe
n molecules H2, LiH, CH4, NH3, H2O and HF, with a Slater-Jastrow trial wave function. Harmonic
frequencies for diatomics are also in good agreement with experiment. An antithetical sampling
a
J method is also discussed for additional reduction of thevariance.

3
1
The optimization of molecular geometries and crys-

tal structures and ab initio molecular dynamics simula-
]
h tionsareamongthemostsignificantachievementsofsin-
p gle particle theories. These accomplishments were both 0.06
- possible thanks to the possibility of readily computing
p
m fOoprcpeesnhoenimthere aiopnpsrowxiitmhaintiotnh.e fTrahmeeawpoprrkoxoifmtahtee Btroerant-- ensity0.04
o ment of electron interactions typical of these approaches ce d
c or
can,however,leadtoquantitatively,andsometimesqual- f
.
s itatively, wrongresults. This fact, together witha favor- 0.02
c
i able scaling of the computational cost with respect to
s
the number of particles, has spurred the development of
y
0
h stochastictechniques,i.e. quantumMonteCarlo(QMC) 0 0.5 1 1.5
distance from H (a.u.)
p methods. Despite the higher accuracy achievable for
[ many physicalproperties, the lack of an efficient estima-
FIG.1: Forcedensityalongthez-direction fortheHatomin
torforforceshasprevented,untilrecently[1,2,3],theuse
2 LiH.Thebondisalongthez-axis,withalengthof3.316Bohr.
of QMC methods to predict even the simplest molecular
v The continuous black curve is calculated from the Hartree-
7 geometry. The chief problem is to have a Monte Carlo Fock orbitals. The dashed line is the estimate of fz using
8 (MC) estimator for the force with sufficiently smallvari- thebareestimator. CirclesareobtainedinanidenticalQMC
0 ance. Forexample,inall-electroncalculations,astraight- simulation using the antithetic sampling technique outlined
9 forward application of MC sampling of the Hellmann- in the text.
0
Feynman estimator has infinite variance. This can be
4
easily seen from the definition of the force. For a nu-
0
/ cleus of charge Z at the origin, the force can be written,
s together with its variance, as a function of the charge the origin. In a QMC calculation based in configura-
c
i density ρ(r) as tion space, the charge density is a sum of delta func-
s tions: ρ(r) ∝ δ(r −r′), where the sum is over all
y r′
h F=Z drρ(r) r ; σ2 =Z2 drρ(r) 1 −F2. (1) sNeepaerlaectetrlyonthpeoesPlieticotnrosnasnwditahllinMaCdsisatmanpcleesR. Wofetchoenastidoemr
p r3 r4
: Z Z and those outside. The contribution to the force from
v charges outside, FO, can be calculated directly with the
i Since the electronic density is finite at the origin, the z
X Hellmann-Feynman estimator in Eq. (1). The contribu-
variance integral diverges.
tion from inside the sphere is responsible for the large
r In this paper, we propose a modified form for the
a variances in the direct estimator. It is convenient to in-
force estimator which has finite variance. This estima-
troducea“forcedensity”definedastheforcearisingfrom
tor is then used to calculate forces and predict equilib-
electron charges at a distance r from the origin:
rium geometry and vibrational frequencies for a set of
small molecules. Without loss of generality we will con-
sider only the z-component of the force on an atom at f (r)=Z dΩρ(r,θ,φ)cosθ (2)
z
Z
Then the force is given as:
∗Electronicaddress: [email protected]
R
†Electronicaddress: [email protected] F =FO+ f (r)dr. (3)
‡Electronicaddress: [email protected] z z z
Z0
2
Theforcedensityisasmoothfunctionofrthattendsto0
linearlyas r approachesthe origin. The force density for 0.02
HinaLiHmoleculecomputedwithHartree-Fockandtwo
0.01
differentQMCestimatorsisshowninFig.1. Asexpected
the bare force estimator fluctuates wildly at small r. ) 0.00
u.
Becausethe forcedensityis asmoothfunction, wecan a.-0.01
(
represent it in the interval (0,R) with a polynomial e
c
or-0.02 1.10 R
M F eq
1.05 R
f (r)= a rk (4) -0.03 eq
z k R
eq
kX=1 -0.04 0.95 R
e eq
and determine the coefficients, a , by minimizing 0.90 R
k -0.05 eq
R 2 -1 0 1 2 3 4 5 6 7 8 9 10
χ2 = drrm fz(r)−fz(r) (5) projection time (Hartree-1)
Z0 h i
whererm isaweightfactorusedtobaelancecontributions
FIG.3: Projection oftheforceinLiHusingforward walking.
from different values of r. Thepoints at negativeimaginary timegive theVMC values.
Sincethe relationbetweentheforceandtheforceden- Valuesat 0 are the mixed estimates of the DMC simulation.
sity is linear, and the relation between the fitting coeffi-
cientsandtheelectronicdensityislinear,wecandirectly
writetheforceasaveragesovermomentsoftheforceden-
Note that for the bare estimator g(r) = θ(R−r). Be-
sity. After some manipulations we arrive at:
cause of the restriction on the basis, the variance of the
Ne z newestimatorisfiniteaslongasm>−1/2. Wehavenu-
F =FO+Z g(r ) i , (6)
z z i r3 merically found that the weighting factor m = 2, where
*Xi=1 i +MC eachvolumeelementisweightedequally,givesthelowest
where the new estimator function is: variance estimate of the force.
To derive the estimator we have used the fact that
M
f (r) goes linearly at small r. [17] This is the crucial
g(r)=θ(R−r) c rk+m. (7) z
k property that allows to filter out the s-wave component
k=1
X of the density responsible for the variance divergence.
The coefficients c ’s are determined by c= S−1h where Theoriginalestimatoriscorrectforanyarbitrarycharge
k
the Hilbert matrix S and the residual vector h are density while the new filtered one uses physical proper-
ties of the charge density to reduce the variance. The
Rm+k+j+1 Rj+1
S = , h = . (8) variance depends on the fitting radius R and on the ba-
kj j
m+k+j+1 j+1
sis set size M. As R increases,the size of the basis must
increase, which increases the variance. Charge densities
corresponding to low energy states must be smooth and
0.053 we typically find that only 2 or 3 basis functions are
needed. The size of the basis can be reduced by using
0.02
0.052
more appropriate basis sets. For example, in all calcula-
0.051 tions reported below we used the expansion
0.00
a.u.) 0.050 M
orce (-0.02 Lithium 0.049 Hydrogen fz(r)=fzSD(r) akrk, (9)
F k=0
X
0.048 improved e
-0.04 0.047 polynomial where fzSD is the force density of a single determinant
wave function, which can be readily computed from the
-0.06 0.046 orbitals. The improved basis allows a smaller polyno-
1 2 3 4 5 1 2 3 4 5
mial set and a reduction of the variance. In Fig. 2 the
Number of basis functions
dependence of the bias on the basis set type and size is
shown for the case of a variational Monte Carlo (VMC)
FIG.2: DependenceoftheVMCforceontheexpansionbasis,
for LiH with a bond length of 3.316 Bohr. The fitting radius simulation on LiH at a bond length of 3.316 Bohr.
R =0.6 Bohr. The definitions of the basis functions are in The trial wave functions ΨT used in all cases were of
Eq.’s(4)and(9). TheforcesonHandLiaredifferentbecause theSlater-Jastrowform. Theorbitalswereobtainedfrom
ofthelackoffulloptimizationoftheVMCwavefunction(see a Hartree-Fock calculation using CRYSTAL98 [4]. The
text). electron-electronandelectron-protonJastrowfactorshad
the form of exp(ar/(1 + br)), with a and b optimized
3
TABLE I: Equilibrium distances in ˚A. Experimental, TABLE II: Harmonic frequencies in cm−1. Experimental,
CCSD(T) and B3LYP values were taken from Ref[6]. The CCSD(T) and B3LYP values were taken from Ref[6]. The
CCSD(T)andtheB3LYPresultswereobtainedusingthecc- CCSD(T)andtheB3LYPresultswereobtainedusingthecc-
pVTZbasisset withtheexceptionofLiHwherethe6-311G* pVTZbasissetwith theexceptionofLiHwherethe6-311G*
set was used. PBE results [7] were all obtained using the set was used. PBE results [14] were obtained using ad hoc
aug-cc-pVTZ basis set. gaussian basis sets.
QMC Exp. CCSD(T) B3LYP PBE QMC Exp. CCSD(T) B3LYP PBE
H2 0.7419(4) 0.741 0.743 0.743 0.751 H2 4464(18) 4410 4420 4401 4323
LiH 1.592(4) 1.596 1.618 1.595 1.606 LiH 1445(20) 1369 1414 1405 1380
CH4 1.091(1) 1.094 1.089 1.088 1.096 HF 4032(266) 4181 4085 4138 4001
NH3 (N-H) 1.009(2) 1.012 1.014 1.014 1.023
NH3 (H-H) 1.624(2) 1.624 1.616 1.624 1.634
H2O (O-H) 0.959(2) 0.956 0.959 0.961 0.971 that cannot be simply addressed is the fixed-node er-
H2O (H-H) 1.519(3) 1.517 1.508 1.520 1.531 ror. Infixed-nodeDMC,therandomwalkisforbiddento
HF 0.919(1) 0.918 0.917 0.923 0.932 crossthe nodes of the trial wavefunctionin order to pre-
ventthe loss of efficiency due to the fermion antisymme-
try. Ifthe nodesareaccurate,sois the QMC energyand
electronic density; hence the force. For incorrect nodes,
by minimizing |E −hEi| [5] over points sampled from
loc
|Ψ |2. ThetimestepinthediffusionMonteCarlo(DMC) theenergyisanupperboundtothetrueenergy,butsuch
T
cannotbesaidfortheforce. Itisalsonotnecessarilythe
simulations was chosen to give an acceptance ratio of
case that the forces obtained from Eq. (1) are equal to
98%, a value for which the time step bias on forces was
the gradient of the fixed-node energy[11, 12, 13]: this is
within the statistical error bars.
onlyguaranteedinthelimitofexactnodalsurfaces. The
Since the exact density is needed for the Hellmann-
highqualityofthegeometriesandvibrationalfrequencies
Feynmantheorem,forwardwalking[8]oroneofthevaria-
suggests that these errors, at least for the cases treated
tionalpathintegralalgorithms[9,10]isneededinorderto
inthis paper,are negligible. This is perhapsnot surpris-
evaluate the force estimator. An example of the conver-
ing, since the electronic density is a 1-electron property,
genceofforwardwalkingisshowninFig.3. Theforceas
while the nodal error is a many-body effect.
afunctionoftheforward-walkingprojectiontimequickly
reaches a plateau corresponding to the exact value. In We have also tested another method to further reduce
the variance of the Hellmann-Feynman estimator. The
this example, the variational forces are far from correct.
Thisdiscrepancyresultsfromthelackoffulloptimization filtered estimator performs well on the hydrogen atom
but for heavier nuclei the error bar grows and seems to
ofthetrialwavefunctionmadeoflocalizedbasisorbitals
and atom centered Jastrow factors, and can be reduced scale as Z3. In those cases the new method can poten-
tially be very useful, with error bars scaling between Z
somewhatbyincludingPulayterms[2]. InDMC,forward
walking eliminates the need for the Pulay corrections. and Z2. The method is based on the observation that,
while electronsinthe corecauselargefluctuations inthe
The equilibrium geometries were computed by fitting
force density, they contribute very little to it. A stan-
the QMC forces in the proximity of the equilibrium ge-
dard approach to reduce the variance of a Monte Carlo
ometry to a polynomial with the appropriate symmetry.
estimate is the use of antithetic variates[15]: a positive
Fig. 4 shows the force in hydrogen fluoride in a 2% in-
fluctuation is paired with a negative fluctuation. Sup-
terval around the equilibrium geometry. The equilib-
posetherandomwalkarrivesatamultidimensionalelec-
rium geometries are reported in Table I together with
tronic configuration R, with p (≥ 1) electrons inside a
those given by CCSD(T), DFT using the B3LYP and
radiusR ≤Rofanatomlocatedattheorigin. Weob-
the PBE functional, and experiments. The differences av
tainanantitheticconfigurationR′ byreflectingallpcore
between QMC and experimental values are in all cases
electrons about the origin. We then estimate the force
less than 0.4% and closer to the experiment than the
contributiondue to thep electrons using bothR andR′,
other techniques. For diatomics it is easy to provide an
assigning a weight factor of w(R′) = |ψ(R′)/ψ(R)|2 to
estimateoftheharmonicvibrationalfrequenciesstarting
R′. Their joint contribution to the estimator in Eq. (6)
from the derivative of the force curve at equilibrium ge-
ometry. The QMC frequencies, reported in Table II, are isZ1−w2(R′) ig(ri)zir−3 wherethesumrunsoverthep
ingoodagreementwiththeexperiment,witherrorscom- core electrons. Since w →1 as Rav →0, fluctuations in
P
parabletothatfromCCSD(T)andDFTPBEorB3LYP. the core are much reduced.
This suggests that forces computed within our approach WithinVMCthisschemecanbeimplementedexactly,
areaccuratealsoawayfromtheequilibriumandcouldbe leadingtoadramaticreductionofthevarianceascanbe
used in molecular dynamics calculations or to optimize noticed from Fig. 1. However this estimator is non-local
molecular geometries. and, in DMC, suffers from the same problems as non-
The onlysourceofsystematicerrorinourcalculations localpseudopotentials,making anunbiasedimplementa-
4
0.03 theHellmann-Feynmanestimatorandcanbeunderstood
in the framework of this paper: one can prove that it
0.02 correspondsto filtering out thes-wavecomponentofthe
density leaving the force density unchanged. The semi-
localcharacterofthe“zero-variance”estimatormakesits
a.u.) 0.01 DMCimplementationtrickier. Toovercomethisproblem
e ( there have been attempts[2, 16] to use correction terms
c
or 0.00 similarinnaturetothePulaytermsinsingle-particleap-
F
proaches. In practice, this scheme requires extensive op-
timizationand,althoughpromising,itisunclearifitwill
-0.01
be viable for more complicated cases. In addition, the
value of the force is very sensitive to small errors[16] in
-0.02 thechargedensityandtheoptimizationwithinastochas-
0 1 2 3 4 5 6 7 0.99 1 1.01
-1 H-F Distance (R ) tic technique is probably not sufficiently stable to elimi-
projection time (Hartree ) eq
nate these errors.
FIG. 4: DMC force in hydrogen fluoride. Left panel: evo-
In conclusion, we have developed a simple method for
lution of the force over forward-walking time. Right panel:
fullyprojectedforcesasafunctionofnucleardistance. Req is computing forceswithinquantumMonte Carloandused
theexperimental equilibrium distance. ittofindtheequilibriumgeometriesforsmallpolyatomic
molecules. This has been the first time that a QMC
technique is used to predict geometries of molecules be-
yond diatomics. The only overhead in the calculation is
tionnotstraightforward. Wepostponefurtherdiscussion the necessity of determining unbiased estimators, which
of the antithetic method to a future article. requires the use of either forward-walking or reptation
Two other approaches have been introduced recently MC techniques. The new method leads to very accu-
for the computation of forces in QMC. Filippi and Um- rateforcesdespiteerrorsfromthefixed-nodeapproxima-
rigar have computed forces for diatomics by correlating tion and from its contribution to the energy derivatives.
random walks for interatomic separations a and a′. In Extension of the method, including the antithetic esti-
DMC the difficulty associated with the nodal error and mator technique, to heavier atoms and to atoms with
the branching factor was overcome by neglecting some pseudopotentials[18]is under investigation.
types of correlation. The main drawback of a finite dif-
ferencemethodisthedifficultyofcalculatingallthecom- Thismaterialisbaseduponworksupportedinpartby
ponents of the force simultaneously; for a system of N theU.SArmyResearchOfficeunderDAAD19-02-1-0176.
atoms this method would require 3N separate force cal- Computational support was provided by the Materials
culations. Computational Center and the National Center for Su-
The otherapproach,introducedinRef. [3], is closerto percomputing Applications at the University of Illinois.
our method. It is based on a “zero-variance” version of S.Z. acknowledges support from NSF.
[1] C. Filippi and C. J. Umrigar, Phys. Rev. B 61, R16291 [9] D. Ceperley, Rev.Mod. Phys. 67, 279 (1995).
(2000). [10] S. Baroni and S. Moroni, Phys. Rev. Lett. 82, 4745
[2] M. Casalegno, M. Mella, and A. M. Rappe, J. Chem. (1999).
Phys.118, 7193 (2003). [11] F. Schautz and H. J. Flad, J. Chem. Phys. 110, 11700
[3] R. Assaraf and M. Caffarel, J. Chem. Phys. 113, 4028 (1999).
(2000). [12] F. Schautz and H. J. Flad, J. Chem. Phys. 112, 4421
[4] V. R. Saunders, R. Dovesi, C. Roetti, M. Causa, N. M. (2000).
Harrison, R.Orlando, and C. M. Zicovich, CRYSTAL98 [13] K. C. Huang, R. J. Needs, and G. Rajagopal, J. Chem.
User’s Manual, University of Torino (1998). Phys. 112, 4419 (2000).
[5] D.Bressanini, G.Morosi, andM. Mella, J. Chem.Phys. [14] D.C.Patton,D.V.Porezag, andM.R.Pederson,Phys.
116, 5345 (2002). Rev. B 55, 7454 (1997).
[6] Computational chemistry comparison and bechmark [15] M. H. Kalos and P. A. Whitlock, Monte Carlo Methods.
database, http://srdata.nist.gov/cccbdb/ NIST Volume I: Basics (John Wiley & Sons, 1986).
Standard Reference Database (2004). [16] R. Assaraf and M. Caffarel, J. Chem. Phys. 119, 10536
[7] X. Xu and W. A. Goddard, J. Chem. Phys. 121, 4068 (2003).
(2004). [17] The force density fz(r)is proportional to thepz compo-
[8] B.L.Hammond,W.A.Lester,andP.J.Reynolds,Monte nentofthedensity.Anon-zerovalueasr→0wouldim-
Carlo Methods in Ab Initio Quantum Chemistry (World plyadiscontinuityofρattheoriginalongthez-direction.
Scientific,1994). [18] Thefilteredestimatorcanbeappliedtoatomswithnon-
5
local pseudopotentials, where the density is replaced by be defined bysumming over partial waves.
a density matrix and a corresponding force density can