Table Of ContentFully quantum description of the Zundel ion: combining variational quantum Monte
Carlo with path integral Langevin dynamics
F´elix Mouhat,1 Sandro Sorella,2 Rodolphe Vuilleumier,3 Antonino Marco Saitta,1 and Michele Casula1
1IMPMC, Sorbonne Universit´es, Universit´e Pierre et Marie Curie,
∗
CNRS, IRD, MNHN, 4 place Jussieu, 75252 Paris, France
2International School for Advanced Studies (SISSA), Via Bonomea 26, 34136 Trieste,
Italy and INFM Democritos National Simulation Center, 34151 Trieste, Italy
3PoˆledePhysico-chimieTh´eoriquedel’ENS,LaboratoirePASTEUR-UMRCNRS8640,24rueLhomond,75005Paris,France
7
(Dated: January 13, 2017)
1
0 We introduce a novel approach for a fully quantum description of coupled electron-ion systems
2 from first principles. It combines the variational quantum Monte Carlo solution of the electronic
n partwiththepathintegralformalismforthequantumnucleardynamics. Ontheonehand,thepath
a integral molecular dynamics includes nuclear quantum effects by adding a set of fictitious classical
J particles (beads) aimed at reproducing nuclear quantum fluctuations via a harmonic kinetic term.
1 On the other hand, variational quantum Monte Carlo can provide Born-Oppenheimer potential
1 energy surfaces with a precision comparable to the most advanced post Hartree-Fock approaches,
and with a favorable scaling with the system size. In order to deal with the intrinsic noise due to
] the stochastic nature of quantum Monte Carlo methods, we generalize the path integral molecular
h
dynamics using a Langevin thermostat correlated according to the covariance matrix of quantum
p
MonteCarlonuclearforces. ThevariationalparametersofthequantumMonteCarlowavefunction
-
m areevolvedduringthenucleardynamics,suchthattheBorn-Oppenheimerpotentialenergysurface
is unbiased. Statistical errors on the wave function parameters are reduced by resorting to bead
e grouping average, which we show to be accurate and well controlled. Our general algorithm relies
h
onaTrotterbreakupbetween thedynamicsdrivenbyionicforces and theoneset bytheharmonic
c
interbead couplings. The latter is exactly integrated even in presence of the Langevin thermostat,
.
s thanks to the mapping onto an Ornstein-Uhlenbeck process. This framework turns out to be very
c
efficientalsointhecaseofnoiseless(deterministic)ionicforces. Thenewimplementationisvalidated
i
s on the Zundel ion (H5O+2) by direct comparison with standard path integral Langevin dynamics
y calculations made with a coupled cluster potential energy surface. Nuclear quantum effects are
h
confirmed to be dominant over thermal effects well beyond room temperature giving the excess
p
proton an increased mobility by quantumtunneling.
[
1
v I. INTRODUCTION quired to fully understand the microscopic properties of
2 water. Fromatheoreticalpointofview,liquidwaterhas
9 In spite of about half a century of numerical inves- been studied through Molecular Dynamics (MD) tech-
3 tigations aimed at providing a quantitative and reliable niques such as the empirical force fields methods15–17
3 microscopicdescriptionofliquidwater,thereisstilladif- and ab initio molecular dynamics (AIMD)18–23. Be-
0
. ficultyindescribingaccuratelyprotonhoppinginsidethe cause of its good scalability with the system size which
1 strongandwellstructuredhydrogenbondnetworkgiving is suited to study liquids, Density Functional Theory
0 water its numerous peculiar features. Since water is the (DFT)hasbeenthemostemployedtechniquetoperform
7
solvent of most chemical and biological processes occur- AIMD simulations of bulk water. Nevertheless, DFT-
1
: ring on Earth, one of the long-term goals is to explicitly based calculations performed with Generalized Gradient
v take into account, at a high accuracy level, hydrogen- Approximation (GGA) functionals give rise to an over-
Xi bond driven phenomena in the treatment of species in- structurationofwater: melting point temperature24 and
r teracting with this solvent. The solvation and transport oxygen-oxygen radial distribution function25 are for in-
a properties of hydronium (H O+) and hydroxide (HO−) stance poorly reproduced. Some improvements by in-
3
ions in aqueous environment have enormous impacts on cluding Van der Waals corrections such as the BLYP-
areasrangingfromaqueousacid-basechemistry1–6,enzy- DCACP, BLYP-DG2 or PBE0+TS-vdW(SC) function-
matic protontransfer7,8, as wellas protontransfer inbi- als have been tested and proven to get results closer to
ological channels9, through fuel cell membranes10,11 and experiments26,27 but these approaches suffer from a lack
on ice surfaces facilitating atmospheric reactions12. Re- of generality and transferability. Gillan et al. recently
cent progress in our understanding of proton transport established a numerical scoring system to evaluate the
in water has been reviewed in Refs.13 and 14. qualityofdifferentexchange-correlationfunctionalsfora
Due to a variety of interactions present in water such wide range of water systems28.
as weak and dispersive dipolar interactions that com- As an alternative to DFT-AIMD, very accurate cal-
pete during the formation and the breaking of hydro- culationsusingembedded-fragmentsecond-orderMøller-
gen bonds inside a strongly polar network, an accurate Plessetperturbationtheory(MP2)havebeencarriedout
description of the Potential Energy Surface (PES) is re- and gave satisfactoryresults despite the shortness of the
2
MDtrajectories(6psofequilibrationand5psofproduc- tionsanddistortinghydrogenbondsbecauseofquantum
tion run)29. Perturbative theories such as MP2/MP4 or disorder37–41. In this paper, we propose to extend Zen’s
the most advanced quantum chemistry techniques such pioneeringworkbyincludinganuclearquantumdescrip-
asCoupledCluster(CC)calculationsareindeedsuffering tion within the QMC-driven dynamics. This is achieved
from a poor scalability with the number of particles N within a Path Integral Langevin Dynamics (PILD) ap-
in the system. Thus, the field of interest of these meth- proach which, to the best of our knowledge, has never
ods is typically the study of electronic properties at zero been used in the case of non-deterministic forces.
temperature andstructuralrelaxationsofsmallsystems. On the one hand, we derive an original algorithm to
Quantum Monte Carlo (QMC) techniques constitute integrate the Langevin equations of motion for quantum
a very interesting alternative to study aqueous systems, particles. The designedalgorithmis verygeneralsince it
since the scalability of the method (N3−4) makes the canbeusedtopropagatebothdeterministicandstochas-
simulationsofverylargesystemscomputationallyafford- tic forces. On the other hand, within our framework we
able with a much greater accuracy than DFT in most showthattheinclusionofthequantumnatureofthe nu-
cases. Recently, it has been demonstrated that QMC cleiinthe QMC-drivenMDdoes notslowdownthe sim-
techniques provide results as accurate as the basis set ulation with respect to the classical counterpart. This is
convergedCoupledCluster (CCSD(T)) for small neutral thanks to an effective average of wave function param-
or charged water clusters30,31. Consequently, QMC can eters in the polymer representation of quantum nuclei,
now be used as a benchmark method and certainly ben- which relies upon a bead-grouping technique. There-
efits from its intrisically parallel formulation in modern fore, the cost of the wave function optimization is not
supercomputers. The price to pay here is the systematic increasedby the inclusion of the additional quantumde-
statistical incertitude arising from the stochastic nature grees of freedom.
ofthe QMC approach,whichcanhoweverbe reducedby Totestourmethodologicaldevelopment,weapplyitto
generatinglonger Markovchains as the erroris inversely H O+, namely the Zundel ion42, widely used as bench-
5 2
proportional to the square root of the number of QMC mark system. Indeed, it is the simplest system to ex-
generations. hibit a non trivial proton transfer and its reduced size
Some strategies have thus been developed to incor- makesacomprehensiveandsystematicstudyoftheprob-
porate the intrinsic QMC noise into a classical statisti- lem easier. More importantly, there is a huge amount of
cal mechanics framework of nuclei at finite temperature data43–52 to compare our results with, because this sys-
T. For instance, the Coupled Electron Ion Monte Carlo tem has been massively studied both theoretically and
(CEIMC)method reliesonthe Born-Oppenheimer(BO) experimentally in the last fifty years. We can cite for in-
approximation for treating finite temperature ions cou- stance the extremely accurate Potential Energy Surface
pled with ground state electrons. The Boltzmann distri- (PES)generatedbyBowmanandcoworkersfromalmost
bution function of the ionic degrees of freedom is sam- exactCCSD(T)calculations53. Moreover,becauseofthe
pled at fixed temperature via a Metropolis MC simula- greatimportance of the NQE in this system, it is almost
tionbasedontheelectronicenergiescomputedduringin- always used to test and validate new approaches. For
dependent groundstate QMC calculations32,33. Another example, Ceriotti and coworkers studied this system to
strategy is to resort to a Langevin MD approach to cor- validate their novel multiple time steps integrator to ac-
relate the noise driving the Langevin Dynamics (LD) of celerate MD by including NQE at low cost using a ther-
the system by the QMC forces covariance. In that case, mostat which works with a colored noise54.
the wave function parameters are optimized by Varia-
Thepaperisorganizedasfollows. InSectionII,wein-
tional Monte Carlo (VMC) calculations along the MD
troduce the theoretical framework by first reviewing the
pathsuchthattheelectronicsolutionhasalwaysenergies
algorithms for classical Langevin simulations with QMC
andforcesascloseaspossibletothetrueBOsurface34,35.
forces (Secs. IIA and IIB). In Secs. IIC and IID we
Lastyear,Zenandcoworkersappliedthelattermethod extend the classical description of nuclei by introducing
to perform the very first QMC-based MD simulation of original and efficient algorithms to integrate the equa-
bulk liquid water at ambient conditions using 32 wa- tionsofmotioninthenuclearquantumcasewithorwith-
ter molecules in a cubic cell with periodic boundary outnoisyforces. Astudyofthenumericalstabilityofthe
conditions36. By obtaining a better position of the first new path integral integrators is done in Sec. IIE for de-
peak of the oxygen-oxygen radial distribution function terministic forces, where a direct comparison with exist-
compared with most advanced functionals, they demon- ing PILD algorithms is possible. The inclusion of QMC
strated the ability of QMC to tackle this kind of prob- forces is discussed in Sec. III. The general form of the
lems. However, the authors had to make some approxi- QMC wave function for the electronic part is described
mations: they expanded the wave function over a small in Sec. IIIA, the bead-grouping approximation is intro-
basis set and more importantly a classical description of ducedinSec.IIIB,ashortreviewontheQMCionicforce
the nuclei was adopted. evaluationisprovidedinSec.IIIC,whiletheQMC-noise
Indeed,itiswellknownthatNuclearQuantumEffects correction scheme is detailed in Sec. IIID, which yields
(NQE) play a crucial role in the description of water an unbiased sampling of the canonical quantum parti-
or ice by deeply affecting the radial distribution func- tion function even in the presence of stochastic errors.
3
Sec. IIIE reports on the numerical stability analysis of 3N 3N friction matrix γ and a random force vector
×
the new integrators with noisy QMC forces. In Section η. The latter controls the temperature T of the system
IV, we show the results obtained for the Zundel ion at via the Fluctuation-Dissipation Theorem (FDT)56
low and room temperatures by taking into account, or
α
not, the NQE and demonstrate the excellent agreement γ = , (5)
2k T
between the noisy QMC PES and the very accurate de- B
terministic CCSD(T) one, which validates our method. with k the Boltzmann constant, and
B
Finally, in Section V, we provide some concluding re-
marks and we discuss the possible future applications of αijδ(t t′)= ηi(t)ηj(t′) (6)
− h i
the algorithms developed in this paper.
is the force covariance matrix which simply reduces to a
diagonal form for white noise57,58. The dynamics gen-
erated by the stochastic differential Eqs. (3) and (4) is
II. METHODS
Markovian, as the noise fluctuations are locally corre-
lated in time. In the case the ionic forces are computed
In this Section, we derive the methodology to sam-
byQMCmethods(aswewillseeinSec.IIID),thecovari-
ple the canonical partition function corresponding to
ance matrix has non-trivial off-diagonalelements, and η
Hamiltonians including either classical or quantum nu-
fluctuations are correlatedby the QMC statistics. Thus,
clei. This is achievedusing the generalizedLangevindy-
the Langevin approach can naturally deal with noisy
namics, which has been demonstrated to be compatible
QMC forces, as already shown by Sorella and cowork-
with the path integral approach to include a quantum
ers in the classical MD framework34–36.
description of the nuclei55. Thanks to the stochastic na-
Eq. (3) belongs to a general class of Stochastic Differ-
ture of this category of non-Newtonian dynamics, the
entialEquations(SDEs),whosesolutionisexplicitwhen-
Langevin formulation of the equations of motion can be
ever both the friction matrix and the force are indepen-
extendedto the QMC case,wherethe forcesareintrinsi-
dent of q59
cally noisy34,35.
t′
p(t′)=e−γ(t′−t)p(t)+ dseγ(s−t′)(f +η(s)). (7)
A. Classical Langevin Dynamics: basic algorithm Zt
Starting from this form, Attaccalite et al.34 proposed to
Let us consider a system constituted by N classical use this integration algorithm in a QMC-based MD, be-
particles described by the Hamiltonian causeinshorttime intervalt′ t=δtthe q dependence
− −
of the force and the friction canbe neglectedand can be
3N p2 assumed to remain at their initial values f(t) and γ(t),
H = i +V (q ,...,q ), (1)
1 3N leading to an explicit integration of the LD at discrete
2m
i
Xi=1 equally spaced times tn =nδt:
where pi and qi are the atomic momenta and positions, pn+1 =e−γnδtpn+Γ(fn+η˜) (8)
respectively, and V is the PES in which the system
q =q +p δt, (9)
evolves. Hereafter, we aregoingto use transformedvari- n+1 n n
ables, more convenient to manipulate, achieved by the where the time evolution has been discretized with time
following mass scaling: step δt, and the subscripts refer to the corresponding
time slice, such that p = p(t ), q = q(t ), and f =
q =q0√m n n n n n
i i i f(q(t )). IntheaboveEquations,theΓmatrixisdefined
n
pi =p0i/√mi as
η =η0√m
i i i Γ=γ−1(1 e−γnδt), (10)
f =f0/√m , (2) n −
i i i
and it multiplies the forces f under the approximation
where q0,p0,η0,f0 aretheoriginalcoordinates
i i i i i=1,..,3N thatδtissmallenoughtodisregardthetime dependence
and {q(cid:8)i,pi,ηi,fi}i=1(cid:9),..,3N denote the transformed ones. of f(q) in Eq. (7). In the same interval δt, the time de-
Afterapplyingtheabovetransformation,thecorrespond- pendentnoiseη(t)isintegrated,leadingtothestochastic
ing second order Langevin dynamics reads: variable η˜, whose mean and variance-covarianceare
p˙ = γp+f(q)+η(t) (3) η˜ =0,
− h i
q˙ =p, (4) δt
η˜Tη˜ =k Tγ2coth γ , (11)
h i B n n 2
where f = ∇qV(q) is the applied force deriving from (cid:18) (cid:19)
−
the PES, and the vectors are 3N-dimensional. To en- respectively. These relations have been derived directly
sure that the equilibrium distribution generated by this from the FDT (5), with the hypothesis that the fluctua-
dynamics is canonical, the r.h.s of Eq. (3) includes a tions are local in time (6).
4
The time step δt and the friction matrix γ are algo- Using joint coordinates would be of little use, if we were
rithmic parameters to be defined in order to carry out not be able to evaluate the exponential e−γˆδt in a closed
the classical LD based on this scheme. If the forces analytic form. Indeed, this is possible because the block
are deterministic, γ can take a diagonal form, such that matrix γˆ can be more conveniently rewrittenin terms of
γ = γ δ , with γ a value appropriately chosen to Pauli matrices σ , σ , σ , as follows:
ij BO ij BO x y z
optimize the efficiency of the canonical sampling. Ac-
cording to Eq. (11), this implies a diagonal η˜Tη˜ , i.e. a γ I I γ
non-correlated noise. In the case of noisy fhorces,iη˜ can γˆ = 2 ⊗I− 2 ⊗σx+i2 ⊗σy+ 2 ⊗σz. (16)
be separated into a contribution δf (already included in
Then, the exponentiation can be straightforwardly ob-
the QMC forces f) and η˜ext. The latter contribution
tainedbyusingstandardPaulimatricesalgebra,andthe
has to be explicitly added in Eq. (8), such that its co-
solution can be given in a closed form:
variance added to the external noise variance satisfies:
(η˜ext)Tη˜ext + δfTδf = η˜Tη˜ , where η˜Tη˜ is given p =e−γδtp +Γ(f +η˜) (17)
bhy Eq. (11)iandh δfTδif ish theiQMC forhce coivariance, n+1 n n
which will have inh generial non-zero off-diagonal matrix qn+1 =qn+Γpn+Θ fn+η˜˜ , (18)
elements. Thus,ifthe forcesarenoisyγ isnon-diagonal, (cid:16) (cid:17)
as discussed in Sec. IIID. In that case, one has to rotate with Γ previously defined in Eq. (10), and
the momenta and positions in the frame which diago-
Θ=γ−2( 1+γδt+e−γδt), (19)
nalizes γ, in order to perform numerically the Markov −
iteration in Eqs. (8) and (9). tn+1
η˜=Γ−1 dteγ(t−tn+1)η(t), (20)
Z
tn
B. Classical Langevin Dynamics: Classical
tn+1
Momentum-Position Correlator (CMPC) algorithm
η˜˜=(Θγ)−1 dt(1 eγ(t−tn+1))η(t). (21)
−
Z
The previous basic algorithm has been further devel- tn
oped in Ref. 35 and later in Ref. 36, to improve the in-
Weimmediatelynoticethatthepresentalgorithmdiffers
tegration of the positions q in Eq. (9), by taking into from those developed in literature60–65, since momenta
accountthe evolutionofp(t)during the time intervalδt.
and positions are propagated simultaneously in a single
Here, we provide an alternative derivation, based on the
iteration thanks to the use of momentum-position cor-
joint momentum-position coordinates
relation matrices. In particular, according to Eqs. (17)
and (18), not only the momenta but also the positions
p
X= , (12) are affected by the integrated Langevin noise, η˜ and η˜˜.
q
(cid:18) (cid:19) Their variance-covariancecan be computed by imposing
where X is a 6N-dimensional vector. Analogously, we the FDT (5) to be fulfilled, under the hypothesis that
the α matrix is q-time independent (as we have done in
extend the definition of the ionic and random force vec-
f η Sec. IIA), and by exploiting that [α,γ]=0. One finds:
tors to be F= and E= , respectively. In this
0 0
(cid:18) (cid:19) (cid:18) (cid:19) δt
extended basis, Eqs. (3) and (4) can be rewritten into a α˜ = η˜Tη˜ =k Tγ2coth γ ,
11 B
h i 2
generalized SDE, which reads as (cid:18) (cid:19)
α˜ = η˜˜Tη˜˜ =k T 2Θ Γ2 Θ−2,
X˙ = γˆX+F+E, (13) 22 h i B −
− α˜ = η˜Tη˜˜ =α˜ =(cid:0) η˜˜Tη˜ =(cid:1)k TγΓΘ−1, (22)
12 21 B
where in this notation the matrix γˆ represents a gener- h i h i
alized friction that couples both momenta and positions while the mean of η˜ and η˜˜ is zero. It is interesting to
by: note that the momentum-position correlator formalism
correlatestheintegratednoise,evenwithoutdealingwith
γ 0 non-diagonalαmatrices,aswewilldointheQMCcase.
γˆ = , (14)
I 0 We also remark that in order to derive the integrated
(cid:18)− (cid:19)
equations of motion (17) and (18), we have made the
with the “physical” friction γ being the same 3N 3N same approximation as in Sec.IIA, as far as the ionic
×
matrixintroducedinEq.(3),andIistheidentitymatrix. forcesf areconcerned. Indeed,inEq.(15)wedisregarded
Analogously to Eq. (7), the formal solution of (13) is the q-time dependence of f in the time interval δt. For
provided by deterministic forces, where α is position independent, it
is the only approximation left in the classical Langevin
t′
X(t′)=e−γˆ(t′−t)X(t)+ dseγˆ(s−t′)(F(X(s))+E(s)). dynamics driven by (17) and (18).
As inthe basicalgorithmofSec.IIA, γ =γ δ for
Zt ij BO ij
(15) deterministic f, while the definition of γ is more general
5
for noisy QMC forces (see Sec. III). Once γ and δt have reads as
been set,the numericalevolutionis performedaccording
3N M
1 1 2
to Eqs. (17) and (18) in the frame which diagonalizes γ. H (p,q)= [p(j)]2+ ω˜2 q(j) q(j−1)
M 2 i 2 M i − i
Xi=1Xj=1(cid:18) (cid:16) (cid:17) (cid:19)
M
+ V(q(j),..,q(j)), (24)
1 3N
j=1
C. Path Integral Langevin Dynamics: Path X
Integral Momentum-Position Correlator (PIMPC) written in mass-scaled coordinates (2), with ω˜ =
M
algorithm M/β~, and subjected to the ring boundary condition
(0) (M)
q q . Thus, the PIMD method is a way
i ≡ i
to calculate quantum mechanical properties using the
As we mentioned in the Introduction, it is well known
recipes from classical statistical mechanics with a mod-
that NQE are extremely important, in liquid water and
ified Hamiltonian containing an additional quantum ki-
aqueous species such as the Zundel ion, because of the
netic term. Using the well-known formula for a generic
light mass of the hydrogenatom. NQE are consequently
observable A
essential to ensure a good description of both quantum
and thermal fluctuations of the hydrogen bond network 1
A = Tr e−βHA , (25)
and must be taken into account since they dramatically h i Z
modify the effective energy landscape of the problem.
(cid:2) (cid:3)
it is straightforward to evaluate the potential energy of
SeveralapproachestodealwithNQEhavebeenproposed
so far, such as the Quantum Thermal Bath (QTB)66–68, the system:
but we have chosen here to focus on the Path Integral 1
(PI) framework for the following reasons. First, the for- hVi≃ (2π~)fZ dfpdfqe−τHM(p,q)VM(q), (26)
mulation of the theory is intrinsically parallel and per- Z
fectlycompatiblewiththeQMCmulti-walkerimplemen- M
tation. This constitutes a major advantage to perform where VM(q) = M1 V(q1(j),..,q3(jN)) is the bead-
highly-parallel simulations for High Performance Com- j=1
dependent potential avePraged over all the replicas. The
puting (HPC) applications. Second, the PI relies on ap-
quantum kinetic energy can also be evaluated according
proximations which are systematically improvable and
to
the exactresultcanbe recoveredwithcontrolledextrap-
olations. Third, this approach can be easily combined T 1 dfpdfqe−τHM(p,q)T (q), (27)
withthe LDleadingtothePILDmethod, forwhichvery h i≃ (2π~)fZ M
Z
efficient algorithms already exist69,70.
whereT (q)ismostoftenchosentobethecentroidvirial
M
WefirstrecallthebasicformalismofthePIMDmethod
estimator,
(see Ref. 69 for further details), before introducing our
approach to propagate the corresponding equations of 3N M
N 1
mZo=tionT.r eIn−βPHI t(hweoitrhy,βth=e q1u/aknBtTum) ispaurstuitailolny afupnpcrtoixoin- TM,vir(q)= 2β + 2M Xi=1Xj=1(cid:16)qi(j)−q¯i(cid:17)∂qi(j)V˜, (28)
mated after a standard Trotter factorization71 by
(cid:2) (cid:3) M
with q¯ = 1 q(j), ∂ ∂ and V˜ =
1 i M j=1 i qi(j) ≡ ∂qi(j)
Z ≃ (2π~)f dfpdfqe−τHM(p,q), (23) M V(q1(j),..,q3(jN))P= MVM. Let us notice that it is not
Z j=1
tPhe only way to estimate the quantum kinetic energy.
with f =MN andτ =β/M. We have replacedhere the Indeed, a more noisy and therefore less convenient es-
true quantum particles by fictitious classical ring poly- timator is obtained with the primitive approximation,
mers consisting of M replicas(beads) ofthe systemcon- yielding:
nectedto eachotherbyharmonicspringsevolvinginthe
imaginary time τ. Without loss of generality,we can ex- 3NM 1 M 1 q(j) q(j−1) 2
T (q)= i − i . (29)
tend the definition of classical vectors to the quantum M,pri 2β − M 2 ~τ !
case, by including all M replicas in the string, such that Xj=1
q q(1),...,q(i),...,q(M) T and it becomes 3NM- It is clear from Eq. (29) that the primitive estimator is
≡
dimensional. Bold vectors will automatically include all more difficult to control because of the negative contri-
(cid:0) (cid:1)
replicas in the quantum case. If interpreted classically, butionsthatcanbecomeimportantwhensomepathsare
the partition function Z in Eq. (23) describes a system more spread out than they should. This apparent draw-
at the effective temperature MT. The Hamiltonian H back actually gives precious information about the qual-
M
corresponding to this quantum-to-classical isomorphism ity of the quantum sampling over the paths since one
6
should obtain T = T in fully converged calcula- matrix, which reads as
vir pri
h i h i
tions in the limit δt 0 for any value of M. In practice
thisisalmostneverm→etexactlyanditisparticularlydif- γ K
γˆ = , (30)
ficult to fulfill it for large M. Indeed, two constraints I 0
(cid:18)− (cid:19)
limit the choice of M. On one side, M must be large
enoughto recoverall quantumproperties of the physical where K is 3NM 3NM matrix defined as follows:
×
system; on the other side, M should not be too large,
otherwise the quantum imaginary time τ = β/M would K(j)(k) =ω˜2 δ 2δ(j)(k) δ(j)(k−1) δ(j)(k+1) . (31)
ih M ih − −
become too smallandtoo large fluctuations in the prim-
(cid:16) (cid:17)
itive energy evaluation would appear. Intheabovedefinition,wehaveusedlowerindexestoin-
dicatetheparticlecomponents,theupperones(inparen-
Bymakingthe choiceofParrinelloandRahmantoset thesis)indicatethebeadcomponents,whileδ (δ(i)(j) for
ij
the physical masses for each ring polymer bead, we are the bead indexes) is the usual Kronecker delta. The K
using the Ring Polymer Molecular Dynamics (RPMD) matrixisdiagonalinthe particlesector,asthe harmonic
approximation to real quantum correlation functions72, springsinthefictitious H (24)coupledifferentreplicas
M
although real-time properties are lost when a thermo- onlyforthesameparticlecomponents. Cyclicconditions
stat is applied to the system, as we shall do in the fol- are implicitly applied in Eq. (31) to the matrix bound-
lowing. Various choices are available for the thermo- ariesinthebeadsector(i.e. (0)=(M)),asthepolymers
stat, ranging from the more conventional Nos´e-Hoover are necklaces.
chain73,74 tostochasticthermostattinginbothitssimple As we have done in Eq. (16) for the classical case, we
path integral Langevin equation (PILE) form69, which expand γˆ in Pauli matrices, as follows:
fulfills the FDT, and its Generalized Langevin Equation
(GLE) flavors55,69,75, where the noise is coloredto simu- γ K I K+I γ
γˆ = I+ − σ +i σ + σ . (32)
late quantum effects and accelerate the convergence of 2 ⊗ 2 ⊗ x 2 ⊗ y 2 ⊗ z
the PIMD with the number of beads, at the price of
breaking the validity of the FDT. Here, we would like Moreover, we require that [K,γ] = 0, which is fulfilled
to take the Langevin approach of stochastic thermostat, when the forces f are deterministic. Indeed, γ obeys
fulfilling the FDT, in order to sample the PI canonical the FDT (5) in the extended 3NM dimensional space.
distribution in Eq. (23). This is the so-called “stan- Fordeterministicforces,γ isbeadindependent,whilefor
dard” PILD approach. However, instead of using the noisy forces the matrix α could be correlated according
PILE integration algorithm69, we will extend our clas- to the force covariance δf(j)δf(k) , and γ could acquire
h i h i
sical Momentum-Position Correlator (MPC) algorithm, a bead dependent contribution to fulfill the FDT. For
described in Sec. IIB, to the quantum case. As already the time being, let us assume γ to be bead independent,
mentioned, this framework will allow us to incorporate which implies that K and γ commute. We will see in
noisyQMCforcesintheequationsofmotion,byusingan Sec. IIID how to recover this condition even in the case
appropriately tailored Langevin noise without breaking of QMC forces. Then, Eq. (32) and [K,γ]=0 allow one
the FDT, as it is explained in Sec. III. Moreover, in our toevaluateeγˆδt inaclosedanalyticformforeach(upper
final algorithm, which makes use of a Trotter breakup71 and lower) block component of the SDE formal solution
betweentheharmonicandthephysicalmodes,thequan- in Eq. (15). The related algebra is quite tedious and we
tum harmonic part is integrated exactly together with refer the reader to the Appendix B. The resulting inte-
the Langevin thermostat for the harmonic frequencies. gratedequationsofmotionsleadtothefollowingMarkov
Indeed, the corresponding Langevinequations of motion chain:
arelinear. Itisthenpossibletoexploitthemappingonto
an Ornstein-Uhlenbeck process (OUP), which is exactly pn+1 =Λ1,1pn+Λ1,2qn+Γ(fn+η˜),
integrable, as explained in Sec. IID. q =Λ p +Λ q +Θ f +η˜˜ , (33)
n+1 2,1 n 2,2 n n
We start by applying the same idea as of our classi- (cid:16) (cid:17)
with the integrated 6NM-dimensional noise (E ) and
calalgorithmwithmomentum-positioncorrelatorstothe int
force (F ) which are, respectively,
quantum-to-classical mapping Hamiltonian described in int
Eq.(24). The equations of motion for H canbe rewrit-
M tn+1
ten in the same form as in the Eq. (13), provided X is Γη˜ η
now interpreted as a 6NM-dimensionalvector, the ther- Eint = Θη˜˜ = dteγˆ(t−tn+1) 0 , (34)
mal noise η lives in the 3NM-dimensional space, and (cid:18) (cid:19) tZn (cid:18) (cid:19)
the ionic forces are defined such as f f(j) j=1,...,M, tn+1
where fi(j) = −∂qi(j)V˜. The main differ≡enc{eiwi}thi=r1e,.s..p,3eNct Fint =(cid:18)ΘΓffnn(cid:19)= Z dteγˆ(t−tn+1)(cid:18)f0n(cid:19)
to the classicalformalismis inthe generalizedfrictionγˆ, tn
whichmustberedefinedinordertoincludetheharmonic f
=γˆ−1 I e−γˆδt n . (35)
couplings between the beads. γˆ is now a 6NM 6NM − 0
× (cid:18) (cid:19)
(cid:0) (cid:1)
7
The expressionsfor Λ, Γ and Θ matrices are quite com- γBO (bead independent) and K (particle independent).
plexsotheyaregiveninAppendixB.Usingthisintegra- In this frame γ is diagonal, with eigenvalues:
tion scheme, one is in principle able to work with large
time steps for the dynamics since all the high frequency γ(k) =γBO+γ(k) (39)
i i harm
vibrations are substantially damped in this type of dy-
forthe setofeigenvectorsofbothγBO (lowerindex)and
namics. Similarly to the classical case, we can compute
K (upper index). This uniquely defines the full matrix
the noise correlation matrix by using its definition:
γ, as it can be reconstructed from its eigenvectors and
eigenvalues. For noisy forces the γBO spectrum will be
0 i
ET E = dteγtαˆeγ†t, (36) nontrivial, while for deterministic forces we haveγiBO =
h int inti γBO i. Thus, even in the quantum case which mixes
Z ∀
−δt harmonic with physical modes, the only parameter left
to tune is γ , as the optimal damping for harmonic
BO
where
forces is fully determined by Eq. (38).
To summarize, we have introduced here an algorithm
α 0
αˆ = . (37) which is able to integrate simultaneously both the phys-
0 0
(cid:18) (cid:19) ical molecular and the fictitious ring polymer vibration
modes in a single time step δt. This could be appeal-
Onceagain,wereferthereadertotheAppendixBforthe
ing, but it turns into a disadvantage, as we will see in
analytic expression of the noise correlators in the quan-
the tests reported in Sec. IIE. Indeed, the compromise
tum case.
betweenphysicalandharmonic modes made for the fric-
After setting up the generalformalism, we need to de-
tioninEq.(39)doesnotseemoptimal. Thatchoiceleads
fine the γ matrix in order to fully determine the nu-
to ovedamped dynamics with lower diffusion. Moreover,
merical scheme, as we have done for the classical algo-
we recallthat in Eq. (33) we made the approximationof
rithms (Secs. IIA and IIB). In the quantum case γ is a
constantBOforcesf duringthetimestepδt. Thiscondi-
3NM 3NM matrix. Moreover,wewantthat[K,γ]=0.
tion is hardto meet for fast vibrationalmodes drivenby
×
Let us first introduce γBO as a 3N 3N matrix which
K. The resultingalgorithmturnsouttobe lessstableas
×
is bead independent and defined as in the classical case,
a function of the number of beads M and the time step
namely γBO = γ δ for deterministic forces, with γ
ij BO ij BO δt. We are going to improve it, by resorting to a Trotter
a free parameter, to be optimized. In the case of noisy
breakupin orderto explictly separatethe harmonic part
forces, γBO will have off-diagonal elements, as we have
from the BO one.
anticipatedfortheclassicalalgorithms. Thisisexplained
in Sec. IIID. Therefore, γBO should retainall “optimal”
damping propertiesrelatedto the physics ofeachreplica D. Trotter breakup and Path Integral
(i),setbythe BOforcef(i). We assumethatthe physics Ornstein-Uhlenbeck Dynamics (PIOUD)
within each replica does not change from one bead to
another, thus the optimal γBO can be taken as bead in-
The previous PIMPC algorithm can be improved by
dependent.
employing the Liouvillian approach, that is used to in-
The additional complication of the quantum case is
tegrate the Fokker-Planck differential equation, describ-
represented by the harmonic matrix K, coupling the
ing the evolution of the ensemble probability density
beads to each other, whose eigenvalues are λk = ωk2 for ρ(p,q,t)63. Since the Fokker-Planck is a linear differen-
k =1,...,M, where
tialequation,thisquantityevolvesformallyasρ(t+δt)=
e−Lδtρ(t),with beingtheLiouvilleoperator,definedas
(k 1)π L
ω =2ω˜ sin − for k =1,...,M. (38)
k M 3NM
M
(cid:18) (cid:19) = ∂ +p ∂
L Fi pi i qi
The energy scale set by K can be much larger than the Xi=1
one of the BO potential energy surface. In the case the 3NM
BcaOllyfotrhceesrefleavraentseatuttoocozerrreol,atoinoenctaimnecoτmHp(uutseinagnaHlytais- − Xj=1 γij(cid:18)∂pipj +kBTM∂pi∂pj(cid:19)! (40)
operator in Eq. (48)) of the related Langevin dynamics. in mass-scaled coordinates. is built upon the Hamil-
It has been shown69,76 that the optimal damping which tonian propagation, driven bLy the first two terms, and
minimizesτ isgivenbyγ(k) =2ω ,foreachharmonic the Langevin thermostat, represented by the last two,
H harm k
eigenmode k. deriving from the Fokker-Planck equation. For the sake
We would like γ to be optimal for both BO and har- of readability, i and j run here over all particles and
monic dynamics. A way to define such a matrix is to beads indexed together. In Eq. (40), fBO+fharm
Fi ≡ i i
exploit the commutator [K,γ] = 0, which implies that is the generalized force, comprising the BO and har-
γ and K share a common set of eigenvectors. Thus, monic contributions, i.e. fBO f = ∂ V˜ and
we can work in a rotated frame which diagonalizes both fharm K q , respectiviely.≡ i − qi
i ≡− j ij j
P
8
As we have seen in Sec. IIC, the quantum-to-classical the time dependence of f was the only approximation
n
isomorphismHamiltonianinEq.(24)includesverydiffer- made in the PIMPC case. Here, there is no time de-
entenergy scales. Forinstance, if one workswith a large pendence either in f or in γharm, and thus e−Lharmδt is
n
number of beads M, the ring polymer vibration modes appliedexactly. γharm isdefinedasinSec.IIC,namelyit
can become much stiffer than the molecular vibrations istheoptimaldampingforthedynamicsdrivenbyahar-
arising from the BO potential. For the PIMPC algo- monic H in its normal modes representation. However,
rithm, introduced in Sec. IIC and tested in Sec. IIE, we inthis casethere is nolow-energydamping forsmallfre-
foundthatpropagatingsimultaneouslyallthesemodesis quencies in γharm. Thus, we need to introduce another
notoptimalsincesoftmodeswillbeoverdampedandthe one,γ ,specificto harm,suchthatinthenormalmodes
0
phase space sampling will be less efficient. To overcome representation γharLm is diagonal and reads:
this problem, we would like to split the Liouvillian in
Eq.(40) into just two operators,one containing only the 2ω if 2ω γ
(k) k k 0
γ = ≥ (45)
physical(BO)modes,theotherdependingexclusivelyon harm (γ0 otherwise,
the fictitious (harmonic) modes. To do so, we first sepa-
rate the friction matrix into two contributions, in a way for k = 1,...,M. In other words, the thermostat con-
analogous to Eq. (39): trolledbythefrictionγ0 ismainlyappliedtothecentroid
rototranslationalmodes which wouldnot be thermalized
γ =γBO+γharm. (41) otherwise because their frequency is zero. Contrary to
the other free quantum harmonic modes, the optimal
We can then rewrite the total Liouvillian as the sum of
damping γ is not general and has to be optimized for
0
two terms, = BO+ harm, where
the system under study.
L L L
The equation of motion corresponding to BO is the
L
harm = fharm∂ +p ∂ one in Eq. (8). Indeed, only the momenta are evolved
L i pi i qi in BO. Therefore, no MPC algorithm is needed, and
Xi L
the resulting equations are similar to those of the sim-
γharm ∂ p +k TM∂ ∂ , (42) ple classical Langevin algorithm introduced in Sec. IIA,
− j ij (cid:18) pi j B pi pj(cid:19)! restricted to p, which in the present case is a 3NM-
X dimensional vector. The corresponding γBO is taken
BO = fBO∂ blockdiagonalaccordingtothe beadindex,andpossibly
L i pi beaddependent. Indeed,thereisnoneedofimposingany
i
X commutation relation between γ and K in this case (for
γBO ∂ p +k TM∂ ∂ , (43) the Trotter decomposition in Eq. (44)), at variance with
− j ij (cid:18) pi j B pi pj(cid:19)! the PIMPC algorithmofSec. IIC. Thus, every bead can
X have its own γBO. On the one hand, this is very useful
in such a way that we can break up the evolution op- for QMC noisy forces, as the statistics of the QMC force
erator via Trotter factorization71. In order to preserve covariance matrix is genuinely bead dependent, because
the reversibility of the process and reduce at most the thewavefunctionissampledindependentlyforeachbead
time step error,we approximatethe propagationusing a (see Sec. III). In the QMC case, γBO is fundamental to
symmetric form77,78: correct the intrinsic QMC noise affecting the BO forces.
On the other hand, for deterministic BO forces, the cor-
e−Lδt e−LBOδt/2e−Lharmδte−LBOδt/2. (44) responding γBO can be taken equal to zero, and Eq. (8)
≃
willresult ina simple velocityupdate p =p +δtf .
The equations of motion corresponding to harm are n+1 n 2 n
L In fact, any additional damping in the BO modes will
those in Eq. (13) with the generalized γˆ of Eq. (30)
turn into a slower algorithmic diffusion, as the thermal-
where F = 0. They are linear in both p and q, thus
izationisalreadyguaranteedbytheLangevinthermostat
exactly solvable in an analytic closed form. They de- integratedin harm. Inthislimit, BO willreduceto ,
V
scribe an OUP, and their algebra and properties have L L L
using the same notations as in Ref. 69.
been exploited in the GLE framework, where they are For deteriministic forces, Ceriotti and coworkers69 de-
used to propagate external auxiliary variables together
veloped a different algorithm (dubbed PILE), where the
with the physical momenta in order to generate a col- maindifferencewithrespecttooursisthat harm issplit
orednoisewithacorrespondingeffectivememorykernel. L
into , which propagates the motion of the harmonic
0
In the present algorithm we are going to use the OUP L
oscillatorsin their normalmodes representation,and
to integrate directly the SDE for both p and q, together (the partproportionaltoγharm inEq.(42), being tLhγe
with the Langevin noise, namely without further split- ij L0
remaining part in the first line of the same equation),
ting the Langevin thermostat in harm. The resulting
L whichcorrespondsto the Langevinthermostatactingon
evolution is given by the coupled equations in (33), but
thenormalmodes. ThetotalLiouvillianpropagatorpro-
with the notable difference of f set to zero, as the BO
n posed in Ref. 69 is:
componentofthe totalforce has beenloadedinthe BO
Trotter factors. As we mentioned in Sec. IIC, dropLping e−Lδt e−Lγδt/2e−LVδt/2e−L0δte−LVδt/2e−Lγδt/2. (46)
≃
9
Apart from the additional Trotter breakup, in the above dynamics, which does not constitute a realimprovement
sequence of operators the Langevin thermostat on the since the phase space willbe poorlyvisited. To this pur-
normal modes is the outermost part of the Liouvillian pose, Ceperley introduced the algorithmic diffusion con-
decomposition, for a better control of the target tem- stant D which reads as79
perature. Our algorithm in Eq. (44) performs better
2
than the solution proposed in Eq. (46), since in our case 3N M
(j) (j)
τ q (t+δt) q (t)
Lharm = L0 + Lγ corresponds to an exact propagator " i=1j=1 i − i #
including both the coordinates change and the thermal- D = (cid:16) (cid:17) . (47)
P P
T
izationofthefreequantumringpolymerafterhavingap- * iter +
pliedanormalmodestransformation. As aconsequence,
T isthetotalamountofCPUtimespentforonesingle
one Trotter factorization is saved which in principle en- iter
MDiterationatagiventimestepδt. SinceT isalmost
ablesustoworkwithalargertimestepδtormorequan- iter
equal for the three algorithms under study, the algorith-
tum replicas M. Indeed, the Trotter break-up is obvi-
mic diffusion can be interpreted as the usual diffusion of
ously exact if the two operators involved commute. The
quantum particles in the position space. Furthermore,
splitof harm proposedbyCeriottiisaveryunfortunate
L we will also focus our attention on the potential auto-
choiceasthecommutatorof[ , ]isdominatedbythe
0 γ
term γharmfharm∂ which isLdivLerging as M3 for large correlation time known to be very sensitive to Langevin
ij j pi dampingsincethesoftestvibrationmodesaresometimes
M,asitcanbeeasilyseenbythe exactexpressionofthe
verylongto sample. This useful quantityhas the follow-
harmonic forces and the choice of the optimal friction.
ing expression
Thisargumentsuggeststhatourmethod,involvingcom-
mutators at most diverging as M2, should have a time
Niterδt
step error much better behaved for large M. 1
τ = dt δV(0)δV(t) , (48)
Thesealgorithmicimprovementswillbetestedanddis- V V2 V 2 h i
h i−h i Z
cussed in the next section. 0
with δV = V V the fluctuation of the potential en-
−h i
ergy. τ needs to be minimized to obtain an optimal
V
E. Algorithm stability with deterministic forces production run and this is achieved by working with the
largestpossibletime stepδt. Simultaneously,anoptimal
Wetesttherobustnessofthethreedifferentalgorithms centroid damping γ has to be found to minimize τ , as
0 V
(PIMPC, PILE and PIOUD) detailed in Sec. IIC and we shall detail later.
IID by performing PILD simulations on the Zundel ion
with almost exact deterministic forces which are simply
computedbyfinite differencesofthe CCSD(T)PESpro- 1. Stability with respect to the number of beads M
videdbyHuangandcoworkers53. Sinceouronlyaimhere
istotestthedifferentintegrationschemes,thenumberof We report here the results obtained for the three dif-
MD iterations Niter = 100000 is quite small to visit the ferent integration schemes (PIMPC, PILE and PIOUD)
entirephasespacebutsufficienttoidentifysomepossible at low temperature (T =50 K) in Figure 1 and at room
weakness of the propagators. We recall that PIMPC de- temperature(T =300K)inFigure2asafunctionofthe
notesthe propagatorcorrespondingto the Eq.(33)since number of quantum replicas M. In each case, the four
positions and velocities evolve jointly, whereas the PI- key observables (virial versus primitive kinetic energy,
OUD integrator corresponds to Eq. (44), in which we average temperature, algorithmic diffusion constant and
separatethe physicalfrom the fictitious vibrationmodes potential autocorrelation time) described in the previ-
of the system. ous section are compared. We have decided to perform
Thesealgorithmswillbecomparedtotheveryefficient simulations working with an increasing number of beads
PILEintegratorofEq.(46). AnefficientPILDalgorithm M = 4,8,16,32,64,128and 256, using a reasonable but
must be stable with both largetime stepandlarge num- not-so-small time step δt=0.5 fs.
ber of quantum replicas. Indeed, the collective modes of The quantum kinetic energy estimators are known to
alargeringpolymerbecomeverystiff,somuchmoredif- beusefultoolstodeterminethenumberofbeadsrequired
ficult to control, especially if the integrationtime step is at a given temperature to capture all NQE during the
large. A correct averagetemperature and close values of dynamics. However, other quantities such as the algo-
the virialand the primitive estimators,given in Eq. (28) rithmic diffusion constant and the potential autocorre-
and (29) respectively, point that we are sampling prop- lation time can also be used to diagnose the good con-
erly both the positions and the momenta of all particles. vergence of the quantum kinetic energy. Indeed, these
Beside temperature and kinetic energy, other observ- observableskeepincreasing(quantumkineticenergyand
ables will also be used here to quantify the numerical diffusion) or decreasing (potential autocorrelation time)
efficiency of the tested algorithms. A very stable prop- when the number of beads increases until they reach a
agation of the equations of motion can for instance be plateau value. In the case of the Zundel ion, M = 128
due to an overcautious (too small δt) integration of the quantum replicas seem to be enough to fully recover the
10
0.06
50.4
0.05 T = 50 K T = 50 K
a)
nergy (H 0.04 ature (K) 50.2
etic e 0.03 mper 50
n e
m ki ge t
Quantu 0.02 <<TT<pvTriri>> PP>II MMPIPPLCCE Avera 49.8
0.01 vir
<T<vTir>pr iP>I OPIULDE 49.6 PIMPIPLCE
<Tvir> PIOUD PIOUD
0
0 50 100 150 200 250 0 50 100 150 200 250
Number of beads M Number of beads M
1.8 45
PIMPC
PILE
1.6 u.) 40 PIOUD
a.
onstant 11..24 τme (V 3305
n c n ti
mic diffusio 0 .18 ocorrelatio 2205 T = 50 K
h 0.6 T = 50 K ut 15
Algorit 0.4 ntial a 10
e
PIMPC ot
0.2 PILE P 5
PIOUD
0 0
0 50 100 150 200 250 0 50 100 150 200 250
Number of beads M Number of beads M
FIG. 1. Evolution of the quantum kinetic energy estimators hT i (top left panel), the temperature T (top right panel),
vir/pri
the algorithmic diffusion constant D (bottom left panel) and the potential autocorrelation time τ (bottom right panel) as
V
a function of the number of quantum replicas M evolving at T = 50 K. Solid lines correspond to the virial estimator of the
kineticenergywhereastheprimitiveestimatorcurvesaredashed. Eachalgorithmhasitsowncolorcode: blackforthePIMPC
algorithm, blueforthePILEpropagatorandredforthePIOUDalgorithm. Thetimestepandthefriction arerespectivelyset
to δt=0.5 fs and γ0 =1.4610−3 a.u. (γBO for thePIMPC algorithm).
quantum kinetic energy at low temperature (see Figure ing a single iteration, one has to deal with very differ-
1), whereas only M = 32 beads are required at room ent energy scales, spanning various orders of magnitude.
temperature according to Figure 2. We note that these This leads to an overdamping of the softest intermolec-
values are consistent with those frequently used in the ular modes which are strongly penalized because of the
literature for this system80–82. presence of very high frequency modes. This also ex-
plains why the potential autocorrelation time is slightly
A rapid inspection of Figures 1 and 2 shows that the
larger for this algorithm while the simulation tempera-
PIMPCapproach,inwhichalloperationsarepropagated
ture seems to be rather well-controlled.
into one single block without any mode separation dur-
ing the dynamics, is by far less efficient and stable than The performances of the PIOUD propagator of
the other two. Indeed, the primitive energy estimator Eq.(44)andthePILEalgorithmaremuchclosertoeach
becomes unstable for M = 256 at 50 K and for M 32 other. Looking atthe algorithmic diffusion constantand
≥
at room temperature, showing that this algorithmis not thepotentialautocorrelationtime,thePILEandPIOUD
able to handle the stiffest vibration modes of the ring algorithms have the same computational efficiency. This
polymer. The asymptotic value of the algorithmic dif- is expected because the same normal mode transforma-
fusion constant is also dramatically reduced and even tion is applied in both approaches.
slightly decreases at room temperature for very large
PILE exhibits a built-in stability of the average simu-
number of quantum replicas.
lation temperature with respect to the number of beads.
Indeed,byworkinginamixedspacewereallthephys- On the contrary, our algorithms (both PIMPC and PI-
ical and fictitious vibration modes are propagated dur- OUD)displaya naturalandsmalltime steperroronthe
