Hyper-parallel tempering Monte Carlo: Application to the Lennard-Jones
fluid and the restricted primitive model
Qiliang Yan and Juan J. de Pablo
Department of Chemical Engineering, University of Wisconsin-Madison, Madison, Wisconsin 53706
共Received 14 July 1999; accepted 7 September 1999兲
A new generalized hyper-parallel tempering Monte Carlo simulation method is presented. The
method is particularly useful for simulation of many-molecule complex systems, where rough
energy landscapes and inherently long characteristic relaxation times can pose formidable obstacles
to effective sampling of relevant regions of configuration space. In this paper, we demonstrate the
effectiveness of the new method by implementing it in a grand canonical ensemble for the
Lennard-Jones fluid and the restricted primitive model. Coexistence curves and critical behavior
have been explored by the new method. Our numerical results indicate that the new algorithm can
be orders of magnitude more efficient than previously available techniques. © 1999 American
Institute of Physics. 关S0021-9606共99兲50945-8兴
I. INTRODUCTION
Molecular simulations of complex systems are difficult,
particularly at low temperatures, because configurations can
get easily trapped in local energy minima, thereby precluding
sampling of other, relevant regions of phase space. The in-
herently long characteristic relaxation times that often pre-
vail in complex fluids further aggravate matters. Over the
last few years several, powerful methods have been proposed
to overcome the first difficulty; some examples are provided
by multicanonical sampling,
1,2
1/k sampling,
3
simple
tempering,
4,5
expanded ensembles,
6
J-walking,
7,8
and parallel
tempering.
9–11
While these methods are relatively effective
at overcoming high-energy barriers, they do little to ‘‘accel-
erate’’ the slow relaxation of complex, many-molecule sys-
tems at low temperatures.
It has been increasingly recognized that open ensembles
provide an effective means for overcoming slow-relaxation
problems; molecules can get in and out of a system, thereby
circumventing diffusional bottlenecks. In this paper we draw
elements from simple and parallel tempering, J-walking, ex-
panded ensembles, and histogram reweighting to propose a
new, powerful Monte Carlo method for simulation of many-
molecule systems. In doing so we combine the proven ben-
efits of tempering or extended ensemble techniques with
those of open-ensemble simulations. In this new method,
configurations can hop simultaneously along several inten-
sive variables. Furthermore, our formulation is well-suited
for further combinations of expanded ensemble 共or ‘‘simple
tempering’’兲 with parallel tempering. To distinguish our
more general formulation from conventional 共one-variable兲
parallel tempering, throughout this manuscript we refer to
the new method as ‘‘hyper-parallel tempering Monte Carlo’’
共HPTMC兲. As shown in this manuscript, hyper-parallel tem-
pering is significantly more efficient than existing methods
for simulation of phase transitions. The new algorithm is
relatively simple, and has the added benefit of being easily
incorporated into molecular-dynamics or Monte Carlo simu-
lation programs with minor modifications to existing codes.
In order to demonstrate the usefulness of HPTMC, in
this paper, we have chosen to work with two systems: The
Lennard-Jones fluid and the restricted primitive model of
electrolyte solutions. On the one hand, the Lennard-Jones
fluid has generally been adopted as a benchmark system on
which to test new algorithms. Extensive studies of its phase
behavior facilitate significantly comparisons to literature data
and to existing algorithms. On the other hand, the restricted
primitive model has attracted considerable interest over this
past decade; its liquid–liquid phase transition at low tem-
peratures has posed a challenge to theoreticians, and litera-
ture data have often been called into question. We use
HPTMC to study that phase transition.
The paper is organized as follows. We begin by describ-
ing in detail how hyper-parallel tempering works. We then
define the models and give a brief description of method-
ological technicalities used in our work, such as histogram
reweighting and finite-size scaling techniques. We then
present the results for the Lennard-Jones fluid and the re-
stricted primitive model, and we compare these to available
data. We conclude the paper by discussing the advantages
and disadvantages of the new method, and by outlining some
of its possible future applications.
II. HYPER-PARALLEL TEMPERING MONTE CARLO
For the sake of generality, we consider an arbitrary en-
semble whose partition function is given by
Z⫽
兺
x
⍀
共
x
兲
w
共
x
兲
兿
j
exp
共
f
j
q
j
共
x
兲兲
, 共1兲
where x denotes the state of the system, ⍀(x) is the density
of states, w(x) is an arbitrary weighting function for state x,
f
j
’s are generalized forces or potentials, and the q
j
’s are the
corresponding conjugate generalized coordinates of the sys-
tem. A grand canonical ensemble can thus be recovered by
setting
w
共
x
兲
⫽ exp
共
⫺

U
共
x
兲兲
, f⫽

, q
共
x
兲
⫽ N
共
x
兲
, 共2兲
JOURNAL OF CHEMICAL PHYSICS VOLUME 111, NUMBER 21 1 DECEMBER 1999
95090021-9606/99/111(21)/9509/8/$15.00 © 1999 American Institute of Physics
Downloaded 08 Mar 2007 to 128.104.198.190. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
where

⫽ 1/k
B
T, T is temperature, U(x) is the potential en-
ergy corresponding to configuration x,
is the specified
chemical potential, and N(x) is the number of particles in
configuration x. Alternatively, an isobaric-isothermal NPT
ensemble can be recovered by setting
w
共
x
兲
⫽ exp
共
⫺

U
共
x
兲兲
, f⫽⫺

P, q
共
x
兲
⫽ V
共
x
兲
, 共3兲
where V(x) is the volume that corresponds to configuration
x.
We now consider a composite ensemble consisting of M,
noninteracting replicas of the above mentioned ensemble,
each at a different set of generalized forces and weighting
function. The complete state of the composite ensemble is
specified through x⫽ (x
1
, x
2
, ...,x
M
)
T
, where x
i
denotes
the state of the ith replica. We define the partition function of
the composite ensemble according to
Z
c
⫽
兿
i⫽ 1
M
Z
i
. 共4兲
The un-normalized probability density of the composite state
x is given by
p
共
x
兲
⫽
兿
i⫽ 1
M
w
i
共
x
i
兲
兿
j
exp
共
f
j,i
q
j
共
x
i
兲兲
. 共5兲
To sample configurations from the composite ensemble,
a Markov chain is constructed in such a way as to asymp-
totically generate configurations according to the limiting
distribution function appearing in Eq. 共5兲. Two types of trial
moves are used to realize that Markov chain:
共1兲 Standard Monte Carlo trial moves are used to locally
update each of the replicas of the system. Since replicas
do not interact with each other, standard Metropolis
acceptance–rejection criteria 共for the underlying en-
semble兲 are employed within each replica.
共2兲 Configuration swaps are proposed between pairs of rep-
licas i and i⫹ 1, so that
x
i
new
⫽ x
i⫹ 1
old
,
共6兲
x
i⫹ 1
new
⫽ x
i
old
.
To enforce a detailed-balance condition, the pair of replicas
to be swapped is selected at random, and the trial swap is
accepted with probability:
p
acc
共
x
i
↔x
i⫹ 1
兲
⫽ min
冋
1,
w
i
共
x
i⫹ 1
兲
w
i⫹ 1
共
x
i
兲
w
i
共
x
i
兲
w
i⫹ 1
共
x
i⫹ 1
兲
⫻
兿
j
exp
共
⫺ ⌬ f
j
⌬q
j
兲
册
, 共7兲
where ⌬ f
j
⫽ f
j,i⫹ 1
⫺ f
j,i
is the difference in generalized force
f
j
between the two replicas, and ⌬q
j
⫽ q
j
(x
i⫹ 1
)⫺ q
j
(x
i
)is
the difference of the corresponding conjugate generalized co-
ordinates.
So far tempering has been discussed for the case in
which several fields are sampled at the same time. A further
generalization is to now couple the technique to that of ex-
panded ensembles 共or ‘‘simple tempering’’兲, so one can
simulate long polymer chains at high densities more effi-
ciently. In this latter generalization, the whole simulation
system consists of several expanded grand canonical en-
sembles. Each replica includes a tagged chain, whose length
fluctuates during the simulation. When a swap trial move is
attempted, these tagged chains are also switched. The corre-
sponding acceptance criterion can be obtained by following a
development analogous to the one presented above. More
details of the generalization and its application are given
elsewhere.
27
III. MODELS AND SIMULATION METHODS
A. Models
In this work, we apply the hyper-parallel tempering
method to two different molecular models. The first is the
Lennard-Jones fluid. The potential energy between a pair of
particle i and j, separated by a distance r
ij
, is defined by
U
ij
⫽ 4
⑀
冋
冉
r
ij
冊
12
⫺
冉
r
ij
冊
6
册
, 共8兲
where
⑀
is the depth of the attractive well, and
is a param-
eter that controls the size of the particles. In this work, the
Lennard-Jones potential-energy function is truncated at a
cutoff distance r
c
. To be consistent with previous studies,
we choose r
c
⫽ 2.5
, and the potential is left unshifted.
The restricted primitive model 共RPM兲 of an ionic solu-
tion consists of 2N hard spheres of diameter
, half of them
carrying a negative charge and the other half carrying a posi-
tive charge. A solvent is not modeled explicitly; instead, it is
simply considered as a dielectric continuum with dielectric
constant D. The charged spheres interact via a Coulomb po-
tential,
U
ij
⫽
再
⫹ ⬁, r
ij
⭓
e
2
4
DD
0
z
i
z
j
r
ij
, r
ij
⬍
, 共9兲
where z
i
e and z
j
e are the charges carried by ions i and j,
respectively, e is the charge of the electron e⫽ 1.602
⫻ 10
⫺ 19
C, and D
0
is the dielectric permeability of vacuum,
D
0
⫽ 8.85⫻ 10
⫺ 12
C
2
N
⫺ 1
m
⫺ 2
. Note that the restricted
primitive model requires that
兩
z
i
兩
⫽
兩
z
j
兩
.
Throughout this work results are reported in reduced
units. For the Lennard-Jones fluid, the reduced temperature
is defined as T
*
⫽ k
B
T/
⑀
, where k
B
is Boltzmann’s constant.
The reduced density is defined as
*
⫽ N
3
/V, with N being
the number of particles and V the volume of the simulation
box. For the restricted primitive model, the reduced tempera-
ture is defined as T
*
⫽ 4
DD
0
k
B
T/e
2
, and the reduced
density is given by
*
⫽ 2N
3
/V.
B. Simulation details
In this work, the hyper-parallel tempering method is
implemented in the grand canonical ensemble. By substitut-
ing Eq. 共2兲 into Eq. 共7兲, we arrive at the following accep-
tance criteria for swapping two replicas:
p
acc
共
x
i
↔x
i⫹ 1
兲
⫽ min
关
1,exp
共
⌬

⌬U⫺ ⌬
共

兲
⌬N
兲
兴
共10兲
9510 J. Chem. Phys., Vol. 111, No. 21, 1 December 1999 Q. Yan and J. J. de Pablo
Downloaded 08 Mar 2007 to 128.104.198.190. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
where ⌬

⫽

i⫹ 1
⫺

i
, ⌬U⫽ U(x
i⫹ 1
)⫺ U(x
i
), ⌬(

)
⫽

i⫹ 1
i⫹ 1
⫺

i
i
, and ⌬N⫽ N(x
i⫹ 1
)⫺ N(x
i
).
Two kinds of trial moves are employed, namely, trial
displacements of particles and trial creations or destructions
of particles. The relative frequency of these two kinds of trial
moves is set to be unity. For the restricted primitive model,
the conventional Ewald sum method is used to calculate
long-range contributions to the energy arising from the slow-
decaying Coulomb potential. Conducting boundary condi-
tions are employed in our calculations; It has been pointed
out in the literature that such a boundary condition is essen-
tial for simulations of ionic systems.
20,21
Particle creations
and destructions are conducted in pairs to preserve electro-
neutrality. A distance-bias Monte Carlo scheme is also em-
ployed to further facilitate particle transfer moves.
23
During the simulation, the number of particles N and the
total potential-energy U of each replica is recorded. The joint
distributions p(N,U) are accumulated in the form of a his-
togram. Note that for the restricted primitive model, N is the
number of particle pairs 共the chemical potential is defined as
the total chemical potential per ion pair兲.
C. Histogram reweighting technique
The histogram reweighting technique
12–14
is designed to
extract a maximal amount of information from the results of
a set of molecular simulations. It has been widely used to
determine the near-critical and sub-critical coexistence prop-
erties of fluids. For completeness, in this work we provide a
brief description of its implementation; for additional details,
readers are referred to the original publications. In a grand
canonical simulation, the relevant thermodynamic variables
are the temperature T and the chemical potential
. A simple
grand canonical simulation can be conducted at T⫽ T
0
and
⫽
0
, and a two-dimensional histogram of energy and
number of molecules can be constructed. The entries to such
a histogram H(N,U) represent the number of times that the
system is observed with N particles and potential-energy U.
If the total number of realizations of the system is denoted by
K, the probability P(N,U,T
0
,
0
) that the system has N par-
ticles and energy U at T
0
and
0
is given by
P
共
N,U,T
0
,
0
兲
⫽ H
共
N,U
兲
/K . 共11兲
The probability distribution for a grand canonical ensemble
is given by
P
共
N,U,T,
兲
⫽
⍀
共
N,V,U
兲
exp
共
⫺

U⫹ N

兲
⌶
共
,V,T
兲
, 共12兲
where ⍀(N,V,U) is the microcanonical partition function
共density of states at N and U), and ⌶ is the grand partition
function, given by
⌶
共
,V,T
兲
⫽
兺
N
兺
U
⍀
共
N,V,U
兲
exp
共
⫺

U⫹ N

兲
.
共13兲
Combination of Eqs. 共11兲 and 共12兲 provides a Monte
Carlo estimate of ⍀(N,V,U), given by
⍀
共
N,V,U
兲
⫽ wH
共
N,U
兲
exp
共

0
U⫺ N

0
0
兲
, 共14兲
where w⫽ ⌶(
,V,T)/K is a proportionality constant. Be-
cause ⍀(N,V,U) is independent of T and
, the probability
that the system has N particles and energy U at a different
state point T and
can be estimated according to
P
共
N,U,T,
兲
⫽
H
共
N,U
兲
exp
关
⫺
共

⫺

0
兲
U⫹ N
共

⫺

0
0
兲
兴
兺
N
兺
U
H
共
N,U
兲
exp
关
⫺
共

⫺

0
兲
U⫹ N
共

⫺

0
0
兲
兴
.
共15兲
The average value of any function of N and U can therefore
be estimated from
具
A
典
⫽
兺
N
兺
U
A
共
N,U
兲
P
共
N,U,T,
兲
. 共16兲
In practice, the extrapolation scheme provided by Eq. 15 is
only useful if the distance between (T
0
,
p
) and (T,
) is not
too large. The range of energy and number of particles cov-
ered by a single histogram is limited, and extrapolations be-
yond that range are unfounded.
In order to expand the applicability of the histogram re-
weighting technique, several simulations can be carried out
at different state points. Such simulations have traditionally
been conducted independently of each other, i.e., in series. In
this work, such simulations are conducted in parallel accord-
ing to our tempering scheme; each replica corresponds to a
different state point, and the resulting histograms are com-
bined to generate thermodynamic predictions over a wide
range of conditions.
The underlying bridge to combine multiple histograms is
the fact that the microcanonical partition function
⍀(N,V,U) is independent of temperature and chemical po-
tential. Estimates of ⍀(N,V,U) by Eq. 共14兲 from different
simulations should be consistent with each other 共within the
statistical uncertainty of the calculations兲. If histograms from
different runs overlap sufficiently, it is possible to find a set
of proportionality constants 兵w其 so that this requirement is
fulfilled. An optimal estimate of ⍀(N,V,U) can then be de-
termined as a weighted average of ⍀’s extracted from dif-
ferent runs. Ferrenberg and Swendsen proposed an efficient
method to determine the sought-after proportionality con-
stants and the weighting factors: for R simulation runs, the
optimal 共unnormalized兲 probability distribution is given by
P
共
N,U;
,

兲
⫽
兺
n⫽ 1
R
H
n
共
N,U
兲
exp
共
⫺

U⫹ N

兲
兺
m⫽ 1
R
K
m
exp
共
⫺

m
U⫹ N

m
m
⫺ C
m
兲
,
共17兲
where K
m
is the total number of realizations in run m. The
constants C
m
must be determined self-consistently by the
iterative relationship
exp
共
C
m
兲
⫽
兺
N
兺
U
P
共
N,U;
m
,

m
兲
. 共18兲
In this work, we use a quasi-Newton method to solve the
nonlinear equations 关Eqs. 共17兲 and 共18兲兴.
At sub-critical temperatures, the grand canonical density
distribution is characterized by a double peaked structure,
provided the chemical potential is close to the coexistence
9511J. Chem. Phys., Vol. 111, No. 21, 1 December 1999 Hyper-parallel tempering Monte Carlo
Downloaded 08 Mar 2007 to 128.104.198.190. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
value. If the so-called ‘‘two-state’’ approximation is
adopted,
19
the criteria for phase equilibrium are equivalent to
requiring that the areas under the two peaks be equal. The
determination of the precise location of a coexistence point
can, therefore, be achieved by simply tuning the chemical
potential at any given temperature, until the areas under the
two peaks become the same; the coexisting densities of the
two coexisting phases then correspond to the mean densities
under these two peaks.
D. Mixed-field finite-size scaling
The periodic boundary conditions generally employed in
molecular simulations suppress long-range fluctuations; near
a critical point, thermodynamic properties calculated by
simulation are, therefore, different from those corresponding
to the thermodynamic limit. However, it is possible to esti-
mate with high accuracy the coordinates of the critical point
by resorting to finite-size scaling techniques.
15,16
For the systems studied here, the relevant scaling fields
consist of combinations of the temperature and chemical po-
tential and are given by
⫽

c
⫺

⫹ s
共
⫺
c
兲
, 共19兲
h⫽
⫺
c
⫹ r
共

c
⫺

兲
,
where
is the thermal scaling field, h is the ordering scaling
field, and subscript c serves to denote a quantity evaluated at
the critical point. Parameters s and r, which control the de-
gree or extent of field mixing, are system specific.
Conjugate to the two scaling fields are two scaling op-
erators, namely, the ordering operator M and an energy-like
operator E. They consist of linear combinations of the par-
ticle density
and the energy density u⫽ U/V
M⫽
1
1⫺ sr
共
⫺ su
兲
, 共20兲
E⫽
1
1⫺ sr
共
u⫺ r
兲
.
For the simple case of models with Ising symmetry 共where
s⫽ r⫽ 0), M is simply the magnetization and E is just the
energy density.
In the critical region, the probability distributions of M
and E exhibit a scaling behavior of the form
P
L
共
M
兲
⫽ AL

/
P
M
*
共
AL

/
␦
M
兲
,
共21兲
P
L
共
E
兲
⫽ BL
(1⫺
␣
)/
P
E
*
共
BL
(1⫺
␣
)/
␦
E
兲
where
␣
,

, and
are universal critical point exponents. the
fixed point limiting operator distributions P
M
*
and P
E
*
are
also universal for all systems of a given universality class.
The quantities
␦
M and
␦
E measure deviations from critical-
ity;
␦
M⫽ M⫺ M
c
,
␦
E⫽ E⫺ E
c
; A and B are system-
specific constants.
For the Ising universality class, the universal distribution
functions mentioned above can be estimated from high-
resolution Monte Carlo simulations.
17
To estimate the critical
point of an Ising-class fluid, the temperature, chemical po-
tential and the mixing parameters are tuned so that the re-
sulting distributions P
L
(M) and P
L
(E) collapse onto those
of the Ising model. Values corresponding to the thermody-
namic limit can be estimated by extrapolating finite-size val-
ues according to the scaling behavior
T
c
共
L
兲
⫺ T
c
共
⬁
兲
⬀L
⫺ (
⫹ 1)/
, 共22兲
c
共
L
兲
⫺
c
共
⬁
兲
⬀L
⫺ (1⫺
␣
)/
, 共23兲
where
␣
,
, and
are universal exponents. For the 3D
共three-dimensional兲 Ising universality class,
␣
⬇0.11,
⬇0.54,

⬇0.312, and
⬇0.629.
IV. RESULTS AND DISCUSSION
A. Results for Lennard-Jones fluids
In our calculations for Lennard-Jones fluids, we used 18
replicas, each having a box size L⫽ 7
. The temperature and
chemical potential of each replica are reported in Table I.
These values were selected to guarantee frequent swaps be-
tween replicas. Configuration swaps were attempted every
10 Monte Carlo steps; a total of 10
6
Monte Carlo steps were
performed to generate the energy and density histograms re-
quired to construct joint (N,U) distributions for each replica.
During a simulation, we keep track of ‘‘physical’’ repli-
cas as well as ‘‘logical’’ replicas. A physical replica is an
actual collection of atoms that we follow throughout the
course of the simulation. A logical replica is whatever con-
figuration happens to visit a specific box at some specified
conditions of temperature and chemical potential 共e.g., box i,
at
i
and T
i
). Figure 1 shows the evolution of a logical
replica at T
*
⫽ 0.73 and

⫽⫺5.30, as a function of Monte
Carlo steps; the ordinate axis indicates which physical rep-
lica happens to be visiting the logical replica at T
*
⫽ 0.73 at
any given time during the simulation. This serves to illustrate
how configurations are swapped during the course of the
simulation. As can be inferred from Fig. 1, a physical replica
visits each logical replica relatively frequently and uni-
formly. After a successful configuration swap, two logical
replicas adopt a completely new configuration. Further, con-
figurations corresponding to logical replicas at low tempera-
tures and high densities are ‘‘passed’’ over to logical replicas
at higher temperatures and lower densities; such configura-
tions can then relax more expeditiously and then be passed
back to lower temperature replicas. This process greatly ac-
TABLE I. Temperatures and chemical potentials used for the simulation of
the truncated Lennard-Jones fluid, where T
*
⫽ k
B
T/, with k
B
being Boltz-
mann’s constant.
Replica T
*

Replica T
*

1 1.20 ⫺2.74 10 0.82 ⫺4.52
2 1.17 ⫺2.83 11 0.79 ⫺4.76
3 1.11 ⫺3.02 12 0.76 ⫺5.02
4 1.07 ⫺3.18 13 0.73 ⫺5.30
5 1.03 ⫺3.35 14 0.70 ⫺5.63
6 0.99 ⫺3.51 15 0.67 ⫺5.99
7 0.94 ⫺3.75 16 0.64 ⫺6.38
8 0.90 ⫺3.98 17 0.62 ⫺6.66
9 0.86 ⫺4.23 18 0.60 ⫺6.97
9512 J. Chem. Phys., Vol. 111, No. 21, 1 December 1999 Q. Yan and J. J. de Pablo
Downloaded 08 Mar 2007 to 128.104.198.190. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
celerates the relaxation of the global system and facilitates
sampling of phase space under adverse conditions of tem-
perature and chemical potential.
Figure 2 shows the probability density of the marginal
distribution for the number of particles in the logical replica
at T
*
⫽ 0.73 and

⫽⫺5.30. The chemical potential of that
replica corresponds to a thermodynamic state slightly off the
saturated liquid line. We would like to emphasize the occur-
rence of two distinct peaks in that density distribution. Given
the fact that the temperature is well below the critical tem-
perature, the observed ‘‘tunneling’’ behavior between a con-
densed phase and its vapor indicates that the large free-
energy barrier associated with the vapor–liquid phase
transition can be overcome by the hyper-parallel tempering
method. Such a crossing of high-energy barriers greatly fa-
cilitates the simulation of phase transitions. In contrast to a
multicanonical sampling method, this tunneling is achieved
by simple configuration swaps, rather than by artificial, trial-
and-error flattening of free energy barriers; intermediate
samples in unstable regimes are not necessary in this algo-
rithm.
To demonstrate how accurately the new method can
sample the relative weights of the two peaks, we also show
in the figure the ‘‘real’’ density probability distribution cor-
responding the same thermodynamic conditions; the latter
was calculated from histogram reweighting of data obtained
from all 18 boxes. For clarity, that distribution is shifted by
0.05.
Figure 3 shows the phase diagram calculated from his-
tograms corresponding to all 18 logical replicas. Also shown
are literature data for the same fluid.
18
As expected from a
correct algorithm, the agreement between the two sets of data
is satisfactory. The slight discrepancies at high temperatures
are due to different definitions of the equilibrium saturated
density. In this work we regard the mean density correspond-
ing to a peak of the distribution as the equilibrium value;
Wilding defines it as the peak value of the distribution. As
can be seen in Fig. 3, the proposed method is able to gener-
ate phase equilibrium data at temperatures and densities in
the near vicinity of the triple point of the truncated Lennard-
Jones fluid 共e.g., T
*
⫽ 0.60 and
*
⫽ 0.86). A simple Gibbs
ensemble method or conventional grand canonical simula-
tions would be difficult and unreasonably demanding under
such conditions.
FIG. 4. Autocorrelation function for the average potential energy per atom
for a Lennard-Jones fluid, corresponding to the lowest temperature replica
(T
*
⫽ 0.60,

⫽⫺6.97). The solid line corresponds to results of hyper-
parallel tempering simulations. The dashed line shows results from a con-
ventional grand canonical Monte Carlo simulation.
FIG. 1. Replica number as a function of Monte Carlo steps, for T
*
⫽ 0.73
and

⫽⫺5.30.
FIG. 2. Probability density of the marginal distribution of number of par-
ticles for T
*
⫽ 0.73,

⫽⫺5.30. The solid line represents the original his-
togram for the system at these conditions. The dashed line depicts the com-
bined histogram that results from a full histogram reweighting analysis of all
results at all conditions studied here 共i.e., it is the ‘‘true’’ histogram兲.
FIG. 3. Phase diagram 共vapor–liquid equilibria兲 for a truncated Lennard-
Jones fluid. The squares correspond the results of this work, and the tri-
angles show results reported by Wilding 共Ref. 18兲. Statistical errors are
smaller than the symbol size. The solid line is an Ising form fit to the
simulation data.
9513J. Chem. Phys., Vol. 111, No. 21, 1 December 1999 Hyper-parallel tempering Monte Carlo
Downloaded 08 Mar 2007 to 128.104.198.190. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
