We show that the two phase points considered in the recent simulations of liquid para hydrogen by Hone and Voth lie in the liquid-vapor coexistence region of a purely classical molecular dynamics simulation. By contrast, their phase point for ortho deuterium was in the one-phase liquid region for both classical and quantum simulations. These observations are used to account for their report that quantum mechanical effects enhance the diffusion in liquid para hydrogen and decrease it in ortho deuterium

Konieczny, Martin

Madden, Paul A.

Manolopolous, David E.

Miller, Thomas F., III

Oberhofer, Harald

English

Caltech Authors - Main

Comment on ‘‘A centroid molecular dynamics study of liquid para
hydrogen and ortho deuterium’’ J. Chem. Phys. 121, 6412 2004
Thomas F. Miller III, David E. Manolopoulos and Paul A. Madden
Physical and Theoretical Chemistry Laboratory, Oxford University, South Parks Road, Oxford OX1 3QZ,
United Kingdom
Martin Konieczny
Institut fu¨r Theoretische Physik II, Heinrich-Heine-Universita¨t Du¨sseldorf, Universita¨tsstraße 1,
D-40225 Du¨sseldorf, Germany
Harald Oberhofer
Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna, Austria
~Received 11 October 2004; accepted 5 November 2004; published online 18 January 2005!
@DOI: 10.1063/1.1839867#
Hone and Voth1 have recently reported both classical and
centroid molecular dynamics ~CMD! simulations of liquid
para hydrogen at two different state points (T525 K, r
50.0190 Å23 and T514 K, r50.0235 Å23), and of ortho
deuterium at one (T520.7 K, r50.0254 Å23). A striking
result of their study was that quantum mechanical effects
were found to increase the diffusion in liquid para hydrogen,
but to decrease it in ortho deuterium.1 Here we provide the
explanation for this result.
In search of this explanation, we have investigated the
system-size dependence of the classically computed diffu-
sion constants. Our simulations were performed at the same
three thermodynamic state points, and used the same isotro-
pic ~Silvera-Goldman2! pair potential as Hone and Voth for
the interaction between J50 hydrogen ~deuterium! mol-
ecules. In each (NVE) simulation, we equilibrated the sys-
tem for 20 ps, and then calculated the velocity autocorrela-
tion function for 2 ps by averaging over 100 consecutive 4 ps
trajectories with a time step of 0.5 fs. The temperature was
controlled by resampling momenta from the Maxwell distri-
bution between each trajectory. This procedure was repeated
five times to obtain an average value for ~and a standard
error in! the diffusion constant of each simulation.
The resulting diffusion constants are compared with
those of Hone and Voth1 in Table I. It can be seen from this
table that our classical results agree well with theirs for the
one classical system size ~343 molecules! that they consid-
ered. However, it is also clear that, while it suffices for ortho
deuterium, this system is too small to give a converged result
for para hydrogen. The reason for this becomes apparent
when one considers Fig. 1, which shows a typical configu-
ration from our simulation at the 25 K state point of a system
with 500 molecules (53 face-centered cubic unit cells!. A
large bubble has formed in the middle of the simulation cell,
and we have found that a similar bubble also forms at the 14
K state point. The spontaneous formation of these bubbles is
a clear indication that both state points lie in the liquid-vapor
coexistence region of the classical phase diagram.
The densities that Hone and Voth chose for their two
para hydrogen state points were obtained from earlier
isothermal-isobaric path integral Monte Carlo calculations at
zero pressure.3 These quantum mechanical densities are
about 2/3 of the zero-pressure liquid densities obtained in
purely classical Monte Carlo simulations.3 It is not therefore
surprising that these state points lie in the classical liquid-
vapor coexistence region and lead to phase separation.
The effect of the phase separation on the calculated dif-
fusion constant can be understood from the following argu-
ment. The diffusion of the molecules at the surface of the
bubble is faster than that of the molecules in the bulk liquid,
because the surface of the bubble has a lower density. Since
the fraction of molecules at the surface is proportional to
N21/3, where N is the number of molecules in the simulation
cell, one would expect the calculated diffusion constant to
satisfy
D~N !5D‘1bN21/3, ~1!
where D‘ is the diffusion constant in the high-density liquid
phase and b is a positive constant.
This expectation is confirmed in Fig. 2, which plots the
classical data for para hydrogen in Table I as a function of
N21/3. The resulting values of D‘ are 0.56(2) Å2/ps at the
25 K state point and 0.02(1) Å2/ps at 14 K, where the num-
bers in parentheses are the standard errors in the last digit as
obtained from linear regression. To the extent that one ac-
cepts the extrapolations to N21/350 in Fig. 2, these values
are the correct classical molecular dynamics results for the
two para hydrogen state points. @For reference, earlier clas-
sical molecular dynamics simulations at densities corre-
sponding to zero external pressure have been found to give
D50.5 and 0 Å2/ps ~solid! at 25 and 14 K, respectively.4
Strictly speaking, the pressure at which D‘ is obtained from
the extrapolations in Fig. 2 is not zero, but rather the classi-
cal liquid vapor pressure at each temperature.#
The implications of all this for quantum diffusion in liq-
uid para hydrogen and ortho deuterium are as follows. The
classical and quantum CMD results for para hydrogen in the
paper by Hone and Voth1 pertain to different situations. The
quantum fluid is in the liquid phase at both para hydrogen
state points, whereas the classical fluid is in the liquid-vapor
THE JOURNAL OF CHEMICAL PHYSICS 122, 057101 ~2005!
122, 057101-10021-9606/2005/122(5)/057101/2/$22.50 © 2005 American Institute of Physics
Downloaded 24 Jun 2008 to 131.215.225.137. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp
coexistence region. This does have an indirect effect on the
diffusion coefficient, but it is a rather artificial one. The clas-
sical fluid phase separates into a low-density vapor region
and a high-density liquid region, the latter of which deter-
mines the diffusion constant as N→‘ . Since this liquid re-
gion has a significantly higher density than the correspond-
ing quantum liquid, the classical diffusion constant is found
to be lower, but then this involves a comparison between two
very dissimilar systems.
On the other hand, it is clear from the system-size inde-
pendence of our classical results for ortho deuterium in Table
I that the state point considered by Hone and Voth for this
molecule is in the stable branch of the classical liquid. In this
case, the comparison between CMD and classical molecular
dynamics does provide a valid way to assess the importance
of quantum mechanical effects in liquid diffusion, and one
sees from Table I that these effects lead to a reduction in the
diffusion constant. This is consistent with the fact that quan-
tum fluctuations cause molecules to ‘‘swell,’’ as was illus-
trated for a hard sphere fluid by Chandler and Wolynes.5
When comparing quantum and classical systems with the
same ~uniform! number density, one would clearly expect the
swollen molecules of the quantum system to diffuse more
slowly.
Finally, we note that diffusion constants at the present
three state points calculated using the ring polymer molecu-
lar dynamics method6 corroborate the CMD results of Hone
and Voth and actually increase slightly rather than decrease
when the system size is increased beyond 180 molecules.7
This indicates that all three state points do indeed lie in the
quantum liquid region, as we have assumed throughout our
discussion.
1 T. D. Hone and G. A. Voth, J. Chem. Phys. 121, 6412 ~2004!.
2 I. F. Silvera and V. V. Goldman, J. Chem. Phys. 69, 4209 ~1978!.
3 D. Scharf, G. J. Martyna, and M. L. Klein, Low Temp. Phys. 19, 364
~1993!.
4 J. Cao and G. J. Martyna, J. Chem. Phys. 104, 2028 ~1996!.
5 D. Chandler and P. G. Wolynes, J. Chem. Phys. 74, 4078 ~1981!.
6 I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 121, 3368 ~2004!.
7 T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. ~to be published!.
TABLE I. Calculated diffusion constants for para hydrogen and ortho deu-
terium.
Calculationa System size
Diffusion constant (Å2/ps)b
H2 ~25 K! H2 ~14 K! D2 ~20.7 K!
Hone and Voth
CMD 180 1.52~8! 0.35~5! 0.40~6!
Classical 343 1.20~6! 0.18~2! 0.52~4!
Present work
hcp 4@5332# 180 1.40~1! 0.230~2! 0.500~2!
fcc 4@43# 256 1.31~1! 0.206~2! 0.500~2!
sc 1@73# 343 1.23~1! 0.183~2! fl
fcc 4@53# 500 1.16~1! 0.168~2! fl
fcc 4@63# 864 1.06~1! 0.145~2! fl
aThe CMD and classical results of Hone and Voth are taken from Ref. 1.
bThe number in parentheses after each entry is the standard error in the last
digit. Note that the present ortho deuterium results are already well con-
verged with a system size of 180 molecules.
FIG. 1. A typical configuration from our 500 molecule classical molecular
dynamics simulation of para hydrogen at the T525 K and r
50.0190 Å23 state point. Note the bubble.
FIG. 2. System-size scaling of our computed classical diffusion constants
for para hydrogen ~filled circles! at the T525 K, r50.0190 Å23 and T
514 K, r50.0235 Å23 state points. The empty diamonds are the results
reported by Hone and Voth ~Ref. 1!.
057101-2 Miller et al. J. Chem. Phys. 122, 057101 (2005)
Downloaded 24 Jun 2008 to 131.215.225.137. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp


Comment on "A centroid molecular dynamics study of liquid para hydrogen and ortho deuterium" [J. Chem. Phys. 121, 6412 (2004)]

https://authors.library.caltech.edu/11039/1/MILjcp05a.pdf

Comment on "A centroid molecular dynamics study of liquid para hydrogen and ortho deuterium" [J. Chem. Phys. 121, 6412 (2004)]

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main