The path integral molecular dynamics (PIMD) method provides a convenient way to compute the quantum mechanical structural and thermodynamic properties of condensed phase systems at the expense of introducing an additional set of high frequency normal modes on top of the physical vibrations of the system. Efficiently sampling such a wide range of frequencies provides a considerable thermostatting challenge. Here we introduce a simple stochastic path integral Langevin equation (PILE) thermostat which exploits an analytic knowledge of the free path integral normal mode frequencies. We also apply a recently developed colored noise thermostat based on a generalized Langevin equation (GLE), which automatically achieves a similar, frequency-optimized sampling. The sampling efficiencies of these thermostats are compared with that of the more conventional Nosé–Hoover chain (NHC) thermostat for a number of physically relevant properties of the liquid water and hydrogen-in-palladium systems. In nearly every case, the new PILE thermostat is found to perform just as well as the NHC thermostat while allowing for a computationally more efficient implementation. The GLE thermostat also proves to be very robust delivering a near-optimum sampling efficiency in all of the cases considered. We suspect that these simple stochastic thermostats will therefore find useful application in many future PIMD simulations.

1.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
2.
R. P.
Feynman
,
Statistical Mechanics
(
Addison-Wesley
,
New York
,
1972
).
3.
D.
Chandler
and
P. G.
Wolynes
,
J. Chem. Phys.
74
,
4078
(
1981
).
4.
J. A.
Barker
,
J. Chem. Phys.
70
,
2914
(
1979
).
5.
M.
Parrinello
and
A.
Rahman
, in
Monte Carlo Methods in Quantum Problems
, edited by
M. H.
Kalos
(
Reidel
,
Dordrecht
,
1984
), p.
105
.
6.
M.
Parrinello
and
A.
Rahman
,
J. Chem. Phys.
80
,
860
(
1984
).
7.
R. W.
Hall
and
B. J.
Berne
,
J. Chem. Phys.
81
,
3641
(
1984
).
8.
H. C.
Andersen
,
J. Chem. Phys.
72
,
2384
(
1980
).
9.
B. J.
Berne
and
D.
Thirulamai
,
Annu. Rev. Phys. Chem.
37
,
401
(
1986
).
10.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
11.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
12.
G. J.
Martyna
,
M. L.
Klein
, and
M. E.
Tuckerman
,
J. Chem. Phys.
97
,
2635
(
1992
).
13.
M. E.
Tuckerman
,
B. J.
Berne
,
G. J.
Martyna
, and
M. L.
Klein
,
J. Chem. Phys.
99
,
2796
(
1993
).
14.
G.
Bussi
and
M.
Parrinello
,
Phys. Rev. E
75
,
056707
(
2007
).
15.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
,
J. Chem. Phys.
126
,
014101
(
2007
).
16.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
,
Phys. Rev. Lett.
102
,
020601
(
2009
).
17.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
,
J. Chem. Theory Comput.
6
,
1170
(
2010
).
18.
H.
Trotter
,
Proc. Am. Math. Soc.
10
,
545
(
1959
).
19.
L. S.
Schulman
,
Techniques and Applications of Path Integration
(
Dover
,
New York
,
2005
).
20.
M. F.
Herman
,
E. J.
Bruskin
, and
B. J.
Berne
,
J. Chem. Phys.
76
,
5150
(
1982
).
21.
D. M.
Ceperley
,
Rev. Mod. Phys.
67
,
279
(
1995
).
22.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
23.
B. J.
Braams
and
D. E.
Manolopoulos
,
J. Chem. Phys.
125
,
124105
(
2006
).
24.
J. M.
Sanz-Serna
and
M. P.
Calvo
,
Numerical Hamiltonian Problems
(
Chapman and Hall
,
London
,
1994
).
25.
B.
Leimkuhler
and
S.
Reich
,
Simulating Hamiltonian Dynamics
(
Cambridge University Press
,
Cambridge
,
2004
).
27.
The number of Fourier transforms per time step can be reduced from four to three per degree of freedom in the absence of a thermostat by keeping the momenta in the normal mode representation and transforming the forces into this representation for the updates in Eqs. (21)–(25); this also eliminates the need for the Fourier transforms in the PILE thermostat in Eqs. (27) and (29).
28.
A.
Ricci
and
G.
Ciccotti
,
Mol. Phys.
101
,
1927
(
2003
).
29.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
Oxford
,
2001
).
30.
C. W.
Gardiner
,
Handbook of Stochastic Methods
, 3rd ed. (
Springer
,
Berlin
,
2003
).
31.
M.
Ceriotti
,
G.
Bussi
, and
M.
Parrinello
,
Phys. Rev. Lett.
103
,
030603
(
2009
).
32.
G. J.
Martyna
,
M. E.
Tuckerman
,
D. J.
Tobias
, and
M. L.
Klein
,
Mol. Phys.
87
,
1117
(
1996
).
33.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
107
,
9514
(
1997
).
34.
G.
Bussi
and
M.
Parrinello
,
Comput. Phys. Commun.
179
,
26
(
2008
).
35.
R. A.
Kuharski
and
P. J.
Rossky
,
J. Chem. Phys.
82
,
5164
(
1985
).
36.
A.
Wallqvist
and
B. J.
Berne
,
Chem. Phys. Lett.
117
,
214
(
1985
).
37.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
131
,
024501
(
2009
).
38.
T. E.
Markland
and
D. E.
Manolopoulos
,
J. Chem. Phys.
129
,
024105
(
2008
).
39.
T. E.
Markland
and
D. E.
Manolopoulos
,
Chem. Phys. Lett.
464
,
256
(
2008
).
40.
S. L.
Adler
,
Phys. Rev. D
23
,
2901
(
1981
).
41.
X. W.
Zhou
,
J. A.
Zimmerman
,
B. M.
Wong
, and
J. J.
Hoyt
,
J. Mater. Res.
23
,
704
(
2008
).
42.
A suitable choice of τ0 can be determined automatically for the GLE thermostat in advance of the calculation, as follows. The maximum internal frequency of the free ring polymer is at ωmax=2/(βn). Since the GLE we have used has been fit to a frequency range of 4 orders of magnitude (see Fig. 1), this implies that we should take ω0 in Eq. (41) to be ω0=102ωmax=0.02/(βn), giving τ0=1/(2ω0)=25βn. For liquid water at 298 K with n=32 ring polymer beads, this gives τ0=20fs, and for a hydrogen atom in palladium at 350 K with ten beads it gives τ0=55fs; inspection of Figs. 2–5 shows that these values of τ0 do indeed give nearly optimal sampling for the GLE in every case.
43.
In our simulations of liquid water the use of the PILE thermostat was found to introduce an overhead of 5% over an NVE simulation. The overhead for the GLE thermostat (which requires more random numbers to be generated and matrix-vector multiplications to be performed) was about 30%. Numerical integration of the NHC equations requires the computation of several exponentials for each degree of freedom and chain variable. This effort is then further magnified by the necessity of a multiple time step integration to avoid any drift in the conserved quantity, which we found to result in a 250% overall overhead in our simulations. We should however stress that such a large thermostatting overhead will only be obtained in cases where the physical forces acting on each bead of the ring polymer can be evaluated extremely cheaply, as they can for an empirical water model using ring polymer contraction techniques. In ab initio PIMD simulations the NHC overhead will typically be a tiny fraction of the overall cost of the calculation.
You do not currently have access to this content.