A Formulation of the Ring Polymer Molecular Dynamics

0.24 MB·

by  Atsushi Horikoshi

A Formulation of the Ring Polymer Molecular Dynamics
Atsushi Horikoshi∗
Department of Natural Sciences, Faculty of Knowledge Engineering,
Tokyo City University, Tamazutsumi, Setagaya-ku, Tokyo 158-8557, Japan
The exact formulation of the path integral centroid dynamics is extended to include composites
of the position and momentum operators. We present the generalized centroid dynamics (GCD),
which provides a basis to calculate Kubo-transformed correlation functions by means of classical
averages. We define various types of approximate GCD, one of which is equivalent to the ring
polymer molecular dynamics (RPMD). The RPMD and another approximate GCD are tested in
4
one-dimensional harmonic system, and it is shown that the RPMD works better in the short time
1
region.
0
2

n
a
J
I. INTRODUCTION
4

] There has been great interest in revealing quantum dynamical aspects of condensed-phase molecular systems
h
p suchastheliquidhydrogenandliquidhelium. Thepathintegralformulationofquantummechanics[1]ismostsuited
-
t for numerical analysis of complex molecular systems to provide the basis of the path integral Monte Carlo (PIMC)
n
a and path integral molecular dynamics (PIMD) [2–4]. Most of the static equilibrium properties of finite temperature
u
quantum systems can be computed by means of the PIMC/PIMD techniques. However, it is difficult to apply the
q
[ PIMC/PIMD directly to compute dynamical properties such as real time quantum correlation functions. This is

1 becausetheimaginarytime pathintegralformalismisusedinthe PIMC/PIMD,inwhichweexploittheisomorphism
v betweenthe imaginarytime pathintegralrepresentationofthe quantumpartitionfunctionandthe classicalpartition
0
9 functionofafictitiousringpolymer[2–4]. InthePIMDsimulations,physicalquantitiesareevaluatedbytimeaverages
7 of the dynamical variables as in ordinary classical molecular dynamics simulations. However, the time evolution of
0
. each bead qj in the ring polymer is just fictitious, and, therefore, we cannot directly evaluate real time-dependent
1
0 properties by means of the PIMD method.
4 A number of quantum dynamics methods to calculate real time quantum correlation functions have been pro-
1
: posedso far [5–7]. In this article, we focus onquantum dynamics methods to calculate Kubo-transformedcorrelation
v
functions [8] Bˆ(0)Aˆ(t) K by means of the PIMD technique. Kubo-transformed correlation functions are important
Xi h iβ
quantities that characterize dynamical effects in quantum mechanical systems, and play a central role in the linear
r
a response theory [8]. Full quantum correlation functions Bˆ(0)Aˆ(t) can be reproduced from the Kubo-transformed
β
h i
correlation functions using the relation in the frequency space: BˆAˆ (ω)=E(ω) BˆAˆ K(ω), where E(ω) is a known
h iβ h iβ
factor [8].
Two methods have been proposed so far. One is the centroid molecular dynamics (CMD) [9] and the other is the
ring polymer molecular dynamics (RPMD) [10]. The CMD is a classical dynamics of the centroid variables (q ,p )
c c
on the effective classical potential surface [11, 12]. The name “CMD” comes from the fact that the variable q cor-
c
responds to the “centroid” of the ring polymer beads: q = (1/P) P q . The effective classical potential is a
0 j=1 j
quantum mechanically correctedpotential, and can be evaluated by PIMC, PIMD, or other methods. In the most of
P
CMD simulations, however, the “on the fly” integration scheme [13] is employed, and, therefore, the CMD is usually
∗Electronicaddress: [email protected]
2
implemented by means of the PIMD technique. The CMD can be derived from a more fundamental dynamics, the
centroiddynamics(CD),whichisdefinedbythe quasi-densityoperator(QDO)formalismforoperators(qˆ,pˆ)[14,15].
It has been shown that a correlation function given by the CD is equivalent to the corresponding Kubo-transformed
correlationfunction if the operatorBˆ is linear in qˆand pˆ[14]. Therefore,the CD and CMD are not applicable to the
casethat Bˆ is nonlinear inthem. This is the nonlinearoperatorprobleminthe CD. This difficulty canbe avoidedby
considering an improved CD, correlation functions given by which correspond to the higher order Kubo-transformed
correlation functions [16]. However, it is in general difficult to convert those correlation functions to Bˆ(0)Aˆ(t) or
β
h i
Bˆ(0)Aˆ(t) K.
h iβ
On the other hand, the RPMD is a quite simple method. We just identify the fictitious time evolution of the ring
polymer beads in the PIMD as the real time evolution, and then calculate correlation functions B (0)A (t) RPMD,
h 0 0 iβ
where A andB are the centroidsof position-dependentoperatorsA(qˆ) and B(qˆ), respectively [10]. The RPMD has
0 0
someadvantages: wecantreatoperatorsnonlinearinqˆ,anditworkswellintheshorttimeregiontobeexactatt=0
[17]. These good properties are ensured by the fact that an identity between a static Kubo-transformed correlation
function and the corresponding quantity defined by the path integral, Bˆ(0)Aˆ(0) K= B (0)A (0) PI, holds for any
h iβ h 0 0 iβ
position-dependent operators. The RPMD looks promising, however, this is just a model and has not been derived
from more fundamental theories so far.
Recently, Hone et al. have pointed out that for the case of Aˆ = Bˆ = qˆ, the RPMD is equivalent to the CMD
approximated by the instantaneous force approximation [18]. Their observation is quite suggestive. This correspon-
dence implies that the RPMD can be formulated as an approximate dynamics of a more fundamental dynamics, just
like the CMD is formulated as an approximate CD. However, there is a problem here: The correspondence found by
them is limited for the case of Aˆ = Bˆ = qˆ. This is because the CD and CMD have the nonlinear operator problem.
Therefore, the QDO formalism should be extended to include operators nonlinear in qˆ.
In this article, we develop a new QDO formalism which includes composite operators of qˆ and pˆ. A new CD,
the generalized centroid dynamics (GCD), is defined by means of this formalism. We then show that one of the
approximate GCD is equivalent to the RPMD.
This article is organized as follows. In Sec. 2, we introduce the effective classical potential and density for com-
posite operators. In Sec. 3, the QDO formalism is extended to include composite operators, and the GCD is then
defined. In Sec. 4, severalapproximate dynamics are presented and tested in a simple system. Conclusions are given
in Sec. 5.
II. EFFECTIVE CLASSICAL POTENTIALS AND EFFECTIVE CLASSICAL DENSITIES FOR
COMPOSITE OPERATORS
A. Effective classical potentials for composite operators
Consider a quantum system, the Hamiltonian of which is given by
1
Hˆ = pˆ2+V(qˆ). (1)
2m
For simplicity, we use one-dimensionalnotation in this article, but the multidimensional generalizationis straightfor-
ward. Thedensityoperatorofthissystemisρˆ =e−βHˆ andthequantummechanicalpartitionfunctionisgivenbythe
β
trace of it, Z =Tr ρˆ ; here, β =1/(k T) is the inverse temperature. The phase space path integral representation
β β B
3
of Z is written as [19]
β
∞ q(β~)=q
Z = dq q p e−S[q,p]/~ (2)
β
D D
Z−∞ ZZq(0)=q
P
dq dp
= Pl→im∞ 2jπ~j e−βPHPps, (3)
j=1ZZ
Y
where the imaginarytime interval[0,β~]is discretizedinto P slices, q andp arethe positionandmomentum atthe
j j
j-th time slice, respectively, and β =β/P. The action S[q,p] is defined by
P
β~
S[q,p]= dτ(H[q,p] ipq˙), (4)
−
Z0
and the Hamiltonian of 2P variables q ,p is defined by
j j
{ }
P p2 p (q q )
Hps( q ,p )= j +V(q ) i j j − j−1 . (5)
P { j j} "2m j − βP~ #(cid:12)
Xj=1 (cid:12)q0=qP
(cid:12)
(cid:12)
(cid:12)
Consider position-dependent Hermitian operators Oˆ(i) = O(i)(qˆ) (i = 1 N) which are given by operator
{ } ∼
products of qˆ or sums of them. We refer to these operators as composite operators. For each operator Oˆ(i), there
exists a corresponding “centroid” O(i) in the path integral representation,
0
1 β~ 1 P
O(i) = dτ O(i)(q(τ))= lim O(i)(q ). (6)
0 β~Z0 P→∞P j=1 j
X
Inserting N identities 1= ∞ dO(i)δ(O(i) O(i)) into the partition function Z (Eq. (2)), we rewrite the partition
−∞ c 0 − c β
function as
R
N ∞
Zβ =C dOc(i) e−βUβc({Oc(i)}), (7)
i=1Z−∞
Y
whereO(i) isthestaticcentroidvariablethatcorrespondstothecentroidO(i),andUc( O(i) )istheeffectiveclassical
c 0 β { c }
potential for composite operators
1 ∞ q(β~)=q N
Uc( O(i) )= log dq q p δ(O(i) O(i))e−S[q,p]/~ +C. (8)
β { c } −β "Z−∞ ZZq(0)=q D D i=1 0 − c #
Y
Here and hereafter, we represent irrelevant constant factors by a symbol C. This is the constraint effective potential
with many constraints on composite operators, which has been introduced by Fukuda and Kyriakopoulos in the
context of quantum field theories [20].
B. Effective classical densities for composite operators
Next, consider Hermitian composite operators Oˆ(i) = O(i)(qˆ,pˆ) (i = 1 N) which are given by operator
{ } ∼
products of qˆ, operator products of pˆ, or sums of them. It should be noted that operators which are products of qˆ
and pˆ, such as (qˆpˆ+pˆqˆ)/2, are excluded here. The corresponding centroids are given as
1 β~ 1 P
O(i) = dτ O(i)(q(τ),p(τ)) = lim O(i)(q ,p ). (9)
0 β~Z0 P→∞P j=1 j j
X
4
InsertingN identities1= ∞ dO(i)δ(O(i) O(i))intoEq. (2),weobtainanotherexpressionofthequantumpartition
−∞ c 0 − c
function
R
N ∞
Z =C dO(i) ρc( O(i) ), (10)
β c β { c }
i=1Z−∞
Y
where
1 ∞ q(β~)=q N
ρc( O(i) ) = dq q p δ(O(i) O(i))e−S[q,p]/~ (11)
β { c } C D D 0 − c
Z−∞ ZZq(0)=q i=1
Y
= e−βHβc({Oc(i)}) (12)
is the effective classicaldensity for composite operators,and Hc( O(i) ) is the effective classicalHamiltonian. These
β { c }
effective classical quantities enable a classical description of static quantum properties. If we choose N = 2 and
(Oˆ(1),Oˆ(2))=(qˆ,pˆ), the effective classical Hamiltonian becomes the ordinary one [14]
p2
Hc(q ,p )= c +Vc(q ), (13)
β c c 2m β c
where Vc(q ) is the ordinary effective classical potential [11, 12].
β c
III. EXACT CENTROID DYNAMICS FOR COMPOSITE OPERATORS
Inthis section, we extend the QDO formalism[14]to include Hermitiancomposite operators Oˆ(i) =O(i)(qˆ,pˆ)
{ }
introduced in the preceding section. Then, we define an exact dynamics of centroid variables and show the exact
correspondence between centroid correlation functions and Kubo-transformed correlationfunctions.
A. QDO for composite operators
The canonical density operator for the Hamiltonian (Eq. (1)) can be decomposed as
N ∞
ρˆ =C dO(i) ϕˆc( O(i) ), (14)
β c β { c }
i=1Z−∞
Y
where we introduced an operator
1 N ∞ dη
ϕˆcβ({Oc(i)})= C 2πi e−βHˆ+iPNi=1ηi(Oˆ(i)−Oc(i)). (15)
i=1Z−∞
Y
One can see the validity of this decomposition as follows. A position space matrix element of this operator has its
phase space path integral representation,
hqb|ϕˆcβ({Oc(i)})|qai= C1 q(β~)=qbDqDp N ∞ d2ηπi eiηi(O0(i)−Oc(i))e−S[q,p]/~. (16)
ZZq(0)=qa i=1Z−∞
Y
Taking the trace of this matrix element and using the integral expression of the δ function, (2π)δ(O O ) =
0 c
−
∞ dη eiη(O0−Oc), we reproduce the effective classical density for composite operators (Eq. (11))
−∞
R ρc( O(i) )=Tr ϕˆc( O(i) ) . (17)
β { c } β { c }
h i
The quantum partition function Z is then reproduced using the expression (Eq. (10)). We define the QDO for
β
composite operators by normalizing the operator ϕˆc,
β
ϕˆc( O(i) )
δˆc( O(i) )= β { c } . (18)
β { c } ρc( O(i) )
β { c }
5
B. Generalized centroid dynamics
We define an exact time evolution of the QDO as
δˆc(t; O(i) )=e−iHˆt/~ δˆc( O(i) ) eiHˆt/~. (19)
β { c } β { c }
A dynamical centroid variable O(i)(t) is then defined by
c
O(i)(t)=Tr δˆc(t; O(i) )Oˆ(i) =Tr δˆc( O(i) )Oˆ(i)(t) . (20)
c β { c } β { c }
h i h i
This is the GCD, an exact CD for composite operators. In the case of N = 2 and (Oˆ(1),Oˆ(2)) = (qˆ,pˆ), the GCD is
reduced to the original CD proposed by Jang and Voth [14].
Next, consider a Hermitian operatorAˆ=A(qˆ,pˆ) which is givenby operatorproducts of qˆ, operator products of pˆ,
or sums of them. In terms of the QDO, a classical counterpart to the operator Aˆ is defined as
Ac( O(i) )=Tr δˆc( O(i) )Aˆ . (21)
β { c } β { c }
h i
ThisistheeffectiveclassicaloperatorforAˆ,whichisafunctionofstaticcentroidvariables O(i) . Thetimedependent
c
{ }
effective classical operator is given by using time dependent QDO (Eq. (19)),
Ac(t)=Ac(t; O(i) )=Tr δˆc(t; O(i) )Aˆ =Tr δˆc( O(i) )Aˆ(t) . (22)
β β { c } β { c } β { c }
h i h i
IftheoperatorAˆbelongstotheoperatorset Oˆ(i) whichischosentodefinetheQDO,theeffectiveclassicaloperator
{ }
is equal to the static centroid variable, Ac( O(i) ) = A , and the time-dependent effective classical operator equals
β { c } c
the dynamical centroid variable Ac(t, O(i) )=A (t).
β { c } c
C. Correlation functions
ConsideracoupleofHermitiancompositeoperators,Aˆ=A(qˆ,pˆ)andBˆ =B(qˆ,pˆ). Intermsofthecorresponding
effective classical operators,we define the centroid correlation function by means of a “classical” ensemble average,
1 N ∞
Bc(0)Ac(t) CD = C dO(i) ρc( O(i) ) BcAc(t). (23)
h β β iβ Z c β { c } β β
β i=1Z−∞
Y
If the operator Bˆ belongs to the operator set Oˆ(i) , the centroid correlation function (Eq. (23)) is identical to the
{ }
Kubo-transformed correlationfunction [8, 14]:
Bc(0)Ac(t) CD = B (0)Ac(t) CD
h β β iβ h c β iβ
= 1 βdλ Tr e−(β−λ)HˆBˆ e−λHˆAˆ(t)
Z β
β Z0 h i
= Bˆ(0)Aˆ(t) K, (24)
h iβ
This identity gives us an effective classical way to calculate Kubo-transformed correlation functions. The GCD
procedure to calculate Bˆ(0)Aˆ(t) K is summarized as follows: First, we choose an operator set Oˆ(i) which includes
h iβ { }
the operator Bˆ, and calculate the effective classical density ρc( O(i) ). Second, we compute the GCD trajectories
β { c }
Ac(t) by evolving Ac( O(i) ), where the initial values of the centroid variables O(i) are given by the distribution
β β { c } { c }
ρc( O(i) ). Averaging over the GCD trajectories, we obtain the centroid correlation function (Eq. (23)), which is
β { c }
identical to the Kubo-transformed correlationfunction.
6
IV. APPROXIMATE DYNAMICS: CLASSICAL CD
Although the GCD identity (Eq. (24)) is exact for the case Bˆ Oˆi , the exact computation of the GCD
∈ { }
trajectoriesAc(t)isingeneraldifficulteveninsimpleone-dimensionalsystems. Therefore,approximationstotheGCD
β
are necessary, and we can define various types of approximations. Among them, two approximations, the decoupled
centroid approximation and the classical approximation, are especially important for practical applications. The
former was introduced by Jang and Voth [15], and severaleffective classicalmolecular dynamics of centroidvariables
canbederivedviathisapproximation. Adetaileddescriptionofsuchapproximatedynamicswillbeprovidedelsewhere
[21]. On the other hand, the latter, which is presented in this section, reduces the exact CD to various classical
dynamics of the path integral beads.
A. Classical approximation
The classical approximationconsists of two steps. We first impose an assumption
Ac(t; O(i) ) Ac( O(i)(t) ), (25)
β { c } ≃ β { c }
which means that the time-dependent effective classical operator is assumed to be a function of dynamical centroid
variables O(i)(t) . Then, we approximate dynamical centroid variables by the corresponding centroids defined by
c
{ }
the dynamical path integral beads,
P
1
O(i)(t) O(i)(t)= lim O(i)(q (t),p (t)), (26)
c ≃ 0 P→∞P j j
j=1
X
where the variables q (t),p (t) are assumed to evolve in a classical fashion. This is the classical centroid dynamics
j j
{ }
(CCD),anapproximatedynamicsoftheGCD.Varioustypesofclassicaltimeevolutionsofthepathintegralbeadscan
beintroducedtodefinevarioustypesofCCD.IntheCCD,thecentroidcorrelationfunctions(Eq. (23))areevaluated
by means of the molecular dynamics techniques; if the dynamics is ergodic, ensemble averages can be obtained by
takingtimeaveragesoftrajectoriesgivenbysolvingasetofclassicalequationsofmotion. Inthefollowingsubsections,
we present several variations of the CCD and resulting correlation functions.
Here,wegivesomecommentsonsymmetriesofKubo-transformedcorrelationfunctions. Inthecasethat(Aˆ,Bˆ)=
(A(qˆ),B(qˆ)), the Kubo-transformed correlation function has the same symmetries as the corresponding correlation
functions defined in classical statistical mechanics [10]. Therefore, in this case, the centroid correlationfunction (Eq.
(23)) might be calculated by means of the classical molecular dynamics techniques. However, in case of (Aˆ,Bˆ) =
(A(pˆ),B(pˆ)), the time-reversal symmetry of the Kubo-transformed correlation function holds only if Aˆ = Bˆ or the
potential V(qˆ) is parity symmetric. Furthermore, in the case that (Aˆ,Bˆ) = (A(qˆ),B(pˆ)) or (A(pˆ),B(qˆ)), the time-
reversal symmetry is in general broken. These facts suggest that in those cases, molecular dynamical calculations of
centroid correlationfunctions do not work very much.
B. Phase space CCD
We first consider a classical dynamics where dynamical variables q (t),p (t) evolve according to Hamilton’s
j j
{ }
equations of motion with the Hamiltonian Hps (Eq. (5)),
P
d ∂Hps
q (t) = P , (27)
j
dt ∂p
j
d ∂Hps
p (t) = P . (28)
j
dt − ∂q
j
7
If the system is ergodic, the correlation function can be evaluated as
P
hBβc(0)Acβ(t)iCβD ≃Pl→im∞Z1β j=1ZZ dq2jπd~pj e−βPHPpsBβc({O0(i)(0)})Acβ({O0(i)(t)}). (29)
Y
Werefertothisdynamicsasthephasespaceclassicalcentroiddynamics(PS-CCD)becausetheclassicaltimeevolution
is governed by the Hamiltonian Hps in the phase space path integral representation.
P
It shouldbe noted here that the PS-CCDis notsuited for practical calculationsbecause the imaginaryterm in the
Hamiltonian, ip (q q )/(β ~), makes the dynamics ill-defined.
j j j−1 P
− −
C. Fictitious momenta and masses: summary of the PIMD
The PS-CCD has the problem originatingin the fact that the Hamiltonian Hps is complex. On the other hand,
P
in ordinaryPIMD simulations,we never encounter this kind of difficulty because we use the configurationspace path
integral representation of the quantum partition function,
P
m
Zβ = Pl→im∞j=1r2πβP~2 Z dqj e−βPVPcs. (30)
Y
This expression can be derived from the phase space path integral expression (Eq. (3)) by integrating out the
momentum variables p . The potential Vcs is given by
{ j} P
P
1
Vcs( q )= k (q q )2+V(q ) , (31)
P { j} 2 P j − j−1 j
Xj=1 (cid:20) (cid:21)(cid:12)(cid:12)q0=qP
(cid:12)
where kP =m/(βP2~2). Inserting P identities 1= βP/(2πm˜j) −∞∞dp˜j e−β(cid:12)Pp˜2j/(2m˜j) into Zβ (Eq. (30)), we obtain
Pp R
m dq dp˜
Zβ =Pl→im∞j=1rm˜j ZZ 2jπ~j e−βPHPcs, (32)
Y
where the Hamiltonian is given by
P p˜2
Hcs = j +Vcs. (33)
P 2m˜ P
j
j=1
X
New momenta p˜ and masses m˜ are introduced as fictitious momenta and masses, respectively. This is the
j j
{ } { }
partition function used in the configuration space PIMD simulations [3, 4].
In most of PIMD simulations, we often use a transformationfrom the bead variables q to normal modes x ,
j n
{ } { }
P
x = U q , (34)
n nj j
j=1
X
where U is a unitary matrix. If we transform the variables q in Eq. (30) using an orthogonal matrix
j
{ }
1 2πnj 2πnj
U = cos sin , (35)
nj
√P P − P
(cid:18) (cid:19)
and insert P identities 1 = βP/(2πm˜n) −∞∞dp˜n e−βPp˜2n/(2m˜n) into Eq. (30), we obtain another expression of the
quantum partition function,
p R
P
m dx dp˜
Zβ =Pl→im∞n=1rm˜n ZZ 2nπ~n e−βPHPNM, (36)
Y
8
where the Hamiltonian is given by
P p˜2 πn
HNM( x ,p˜ )= n + 2k sin2 x2 +V˜( x ) . (37)
P { n n} 2m˜ P P n { n}
nX=1(cid:20) n (cid:16) (cid:17) (cid:21)
This kind of expression of the partition function is used in the normal mode PIMD simulations [3, 4].
D. Configuration space CCD
Here we introduce another type of CCD, the configuration space classical centroid dynamics (CS-CCD). There
are two expressions of the CS-CCD. One is a classical dynamics governed by Hamilton’s equations of motion for the
bead variables q (t),p˜ (t) ,
j j
{ }
d ∂Hcs
q (t) = P , (38)
j
dt ∂p˜
j
d ∂Hcs
p˜ (t) = P , (39)
j
dt − ∂q
j
where Hcs is the Hamiltonian (Eq. (33)) used in the configuration space PIMD. If the dynamics is ergodic, the
P
correlation function is given by
P
1 m dq dp˜
Bc(0)Ac(t) CD lim j j e−βPHPcsBc( O(i)(0) )Ac( O(i)(t) ). (40)
h β β iβ ≃P→∞Zβ j=1rm˜j ZZ 2π~ β { 0 } β { 0 }
Y
We refer to this dynamics as the ring polymer CS-CCD because the variables q (t) form a ring polymer where the
j
{ }
nearest neighbors are connected with the spring constant k .
P
The other expression of the CS-CCD is a classical dynamics of the normal mode variables x (t),p˜ (t) . This is
n n
{ }
the normal mode CS-CCD, which is governed by Hamilton’s equations of motion with the Hamiltonian HNM (Eq.
P
(37)),
d ∂HNM
x (t) = P , (41)
n
dt ∂p˜
n
d ∂HNM
p˜ (t) = P . (42)
n
dt − ∂x
n
If the ergodicity of the dynamics is satisfied, the correlationfunction is given by
P
1 m dx dp˜
Bc(0)Ac(t) CD lim n n e−βPHPNMBc( O(i)(0) )Ac( O(i)(t) ). (43)
h β β iβ ≃P→∞Zβ n=1rm˜n ZZ 2π~ β { 0 } β { 0 }
Y
Here,we givetwocomments onthe CS-CCD:(1)Ingeneralsystems,the CS-CCDis exactonlyfor the calculation
of the static correlation functions Bc(0)Ac(0) for position-dependent operators (A(qˆ),B(qˆ)). In the case of t = 0
h β β i 6
or momentum-dependent operators, the fictitious momenta p˜ (or p˜ ) enlarge the gap between the CS-CCD
j n
{ } { }
correlationfunctionandtheKubo-transformedcorrelationfunction. Therefore,thevalidityoftheCS-CCDexpressions
(Eqs. (40)or(43))shouldbe checkeddepending onthe situation. (2)The choiceofthe fictitiousmassesiscrucial. In
the normal mode CS-CCD, if the fictitious mass of P-th normal mode is chosen as m˜ =m, then the CS-CCD gives
P
the exact correlation function Bc(0)Ac(0) for Aˆ = Bˆ = pˆ Oˆ(i) . This is because the choice m˜ = m connects
h β β i ∈ { } P
the fictitious momentum p˜ withthe physicalmomentumcentroidp via the centroidexpressionof the density (Eqs.
P 0
(12) and (13)). In the ring polymer CS-CCD, this condition can be satisfied by setting m˜ = m , i.e. setting all
j
{ }
fictitious masses equal to the physical mass.
9
E. Ring polymer molecular dynamics
Finally,wepresentaspecialcaseoftheringpolymerCS-CCD.Iftheoperators(Aˆ,Bˆ)belongtotheoperatorset
Oˆ(i) ,andallfictitious masses arechosento be the physicalmass, m˜ =m ,the ringpolymer CS-CCDcorrelation
j
{ } { }
function (Eq. (40)) becomes
P
1 dq dp˜
hBβc(0)Acβ(t)iCβD ≃Pl→im∞Zβ j=1ZZ 2jπ~j e−βPHPcsB0(0)A0(t). (44)
Y
This dynamics is equivalent to the RPMD proposed by Craig and Manolopoulos [10]. As we mentioned in the
preceding subsection, the RPMD reproduces the exact correlation functions at t = 0 for (Aˆ,Bˆ) = (A(qˆ),B(qˆ)) or
(pˆ,pˆ). Recently, it has been shown that for the calculations of qˆ(0)qˆ(t) K and pˆ(0)pˆ(t) K, the RPMD is correct up
h iβ h iβ
to O(t6) and O(t4), respectively [17].
F. Simple example: a harmonic system
Here, we consider the case Aˆ = Bˆ = qˆ2 Oˆ(i) , and show the results of two types of approximate CD: the
∈ { }
RPMD and the normal mode CS-CCD. We consider a simple system with a harmonic potential V(qˆ) = (mω2/2)qˆ2.
In this system, the exact Kubo-transformed “position squared” autocorrelationfunction is given by
~2 2 β~ω β~ω
qˆ2(0)qˆ2(t) K = coth cos2ωt+2coth2 1 . (45)
h iβ 4m2ω2 β~ω 2 2 −
(cid:20) (cid:21)
The corresponding RPMD correlation function is calculated as
P P P
1 1 1
qˆ2(0)qˆ2(t) RPMD = lim (cos2ω t+1)+ , (46)
h iβ P→∞β2m2 "n=1ωn4 n n=1l=1 ωn2ωl2#
X XX
where ω = ω2+(4k /m)sin2(πn/P). The corresponding normal mode CS-CCD correlation function is also
n P
obtained as q
P P P
1 1 1
qˆ2(0)qˆ2(t) NM = lim (cos2Ω t+1)+ , (47)
h iβ P→∞β2m2 "n=1ωn4 n n=1l=1 ωn2ωl2#
X XX
where Ω =ω m/m˜ .
n n n
Figure 1 shows the plot of the three correlation functions, Eqs. (45)–(47), at two different temperatures with the
p
parameters ~ = k = m = ω = 1. We set the number of beads P as 1000, which is so large as to make the results
B
convergedsufficiently. Thefictitiousmassesofthenormalmodesarechosenas m˜ =m+(4m/β2~2ω2)sin2(πn/P) .
{ n P }
This choice makes each frequency 2Ω equal to 2ω. At the higher temperature β = 1 (Fig. 1 (a)), the RPMD
n
correlation function and the normal mode CS-CCD correlation function are almost identical, and they coincide with
theexactKubo-transformedcorrelationfunctionatt=0. However,asthetimetincreases,theyslightlydeviatefrom
the exact correlation function. These deviations become more significant at the lower temperature β = 10 (Fig. 1
(b)). As is observedclearlyinFig. 1 (b), the RPMDcorrelationfunctiondampswith time. This is becausethe mode
summation P causesadephasingeffect. Thelowerthetemperaturefalls,themoremodeswithdifferentfrequencies
n=1
arerelevantinthe summation P . Therefore,the dephasingeffectbecomesstrongeratlowertemperature. Onthe
P n=1
other hand, the normalmode CS-CCD correlationfunction is free from such a dephasing effect, thanks to the special
P
choiceoffictitiousmasses m˜ . However,asonecanalsoseeinFig. 1(b), theshorttime behaviorof qˆ2(0)qˆ2(t) NM
{ n} h iβ
is worse than qˆ2(0)qˆ2(t) RPMD. This is because the mass choice adopted in the RPMD, m˜ =m , is the optimized
h iβ { j }
10
one to reproduce the short time behaviors of exact Kubo-transformed correlation functions [17]. These observations
in the simple harmonic system show the importance of the mass choice. If you focus on the short time behaviors of
correlationfunctions,youshouldchoose m˜ =m . Onthe otherhand, ifyourespectthe dynamicalquantumeffects
j
{ }
such as coherent oscillations of correlationfunctions, another choice might be better.
V. CONCLUDING REMARKS
Inthis work,wehavedevelopedanexactCDforcompositeoperatorsofqˆandpˆ, the GCD,whichgivesanexact
identity between a centroid correlation function and a Kubo-transformed correlation function (Eq. (24)). We have
then proposed the classical approximationand the corresponding approximate GCD, the CCD (Eqs. (25) and (26)).
Introducing severaltypes of classical time evolutions, we have defined severaltypes of CCD, the PS-CCD (Eqs. (27)
and(28)),the ringpolymerCS-CCD(Eqs. (38)and(39)),andthe normalmodeCS-CCD(Eqs. (41)and(42)). Ifwe
consideroperatorsAˆandBˆ whichbelongtotheoperatorset Oˆ(i) ,andsetthefictitiousmassesequaltothephysical
{ }
mass, m˜ =m , the ring polymer CS-CCD becomes equivalent to the RPMD proposed by Craig and Manolopoulos
j
{ }
(Eq. (44)) [10]. A schematic diagram of various approximate dynamics is given in Fig. 2. The results of simple
calculations in a harmonic system have shown that the choice of fictitious masses is crucial in CS-CCD calculations.
The RPMD might be the best CS-CCD method to calculate correlation functions in the short time region.
We haveshownthat the RPMD canbe formulatedas anapproximateGCD. However,the classicalapproximation
employed there is rather crude. The physical meaning of the approximation should be clarified, and systematic
schemes to improve the approximation should be developed. Recently, a relationship between the RPMD and the
semiclassicalinstantontheoryhasbeendiscussedinthedeeptunnelingregime,andithasbeenshownthattheRPMD
can be systematically improved in that regime [22].
Finally, we briefly mentionanother approximationscheme. Applying the decoupledcentroidapproximationto the
GCD, we obtain the generalized centroid molecular dynamics (GCMD) [21]. In the case of N =2 and (Oˆ(1),Oˆ(2))=
(qˆ,pˆ), the GCMD is reduced to the CMD proposed by Cao and Voth [9]. The relations between these methods are
summarized in Fig. 2. The GCMD is a dynamics of centroid variables O(i)(t) rather than centroids O(i)(t)
{ c } { 0 }
defined by the summation of the path integral beads. Therefore, the GCMD correlation functions are expected to
be free from the summation-induced dephasing effect seen in the RPMD correlation functions (Fig. 1). However, in
general, it is not so easy to formulate the GCMD and to implement it for practical calculations. This is because the
GCMD isa non-Hamiltoniandynamicswhichshouldbe formulatedinanextendedphasespacespannedbycanonical
multiplets [21]. A properformulationof the GCMD might be givenby a generalizedHamiltoniandynamics proposed
by Nambu [23].
[1] R.P. Feynman and A.R. Hibbs,Quantum Mechanics and Path integrals, McGraw-Hill, New York,1965.
[2] D.M. Ceperley, Path integrals in the theory of condensed helium,Rev.Mod. Phys.67 (1995), pp.279–355.
[3] B.J.BerneandD.Thirumalai, On the Simulation of Quantum Systems: Path Integral Methods, Annu.Rev.Phys.Chem.
37 (1986), pp. 401–424.
[4] B.J.Berne,G.Ciccotti,andD.F.Coker(ed.),Classical andQuantum DynamicsinCondensed Phase Simulations,World
Scientific,Singapore, 1998.
[5] W.H.Miller,ChemicalTheory andComputationSpecial Feature: Quantum dynamicsofcomplexmolecularsystems,Proc.
Natl. Acad.Sci. U.S.A.102 (2005), pp.6660–6664.
[6] E. Rabani, D. R. Reichman, G. Krilov, and B. J. Berne, The calculation of transport properties in quantum liquids using
the maximum entropy numerical analytic continuation method: Application to liquidpara-hydrogen, Proc.Natl.Acad.Sci.