scieee Science in your language
[en] (orig)
Isabel Nitzke, Jadran Vrabec
Numerical Discrimination of Thermodynamic Monte
Carlo Simulations in All Eight Statistical Ensembles
Open Access via institutional repository of Technische Universität Berlin
Document type
Journal article | Accepted version
(i. e. final author-created version that incorporates referee comments and is the version accepted for
publication; also known as: Author’s Accepted Manuscript (AAM), Final Draft, Postprint)
This version is available at
https://doi.org/10.14279/depositonce-18524
Citation details
Nitzke, I., & Vrabec, J. (2023). Numerical Discrimination of Thermodynamic Monte Carlo Simulations in All
Eight Statistical Ensembles. In Journal of Chemical Theory and Computation (Vol. 19, Issue 12, pp.
3460–3468). American Chemical Society (ACS). https://doi.org/10.1021/acs.jctc.3c00252.
This document is the Accepted Manuscript version of a Published Work that appeared in final form in Journal
of Chemical Theory and Computation, copyright © 2023 The Authors. after peer review and technical editing
by the publisher. To access the final edited and published work see https://doi.org/10.1021/acs.jctc.3c00252.
Terms of use
This work is protected by copyright and/or related rights. You are free to use this work in any way permitted by
the copyright and related rights legislation that applies to your usage. For other uses, you must obtain
permission from the rights-holder(s).
Numerical discrimination of thermodynamic
Monte Carlo simulations in all eight statistical
ensembles
Isabel Nitzke and Jadran Vrabec∗
Thermodynamics, Technische Universtität Berlin, 10587 Berlin, Germany
E-mail: vrabec@tu-berlin.de
Abstract
Generalized expressions for thermodynamic properties in terms of ensemble averages
are discussed for adiabatic and isothermal ensembles. They are implemented in the sim-
ulation code ms2 and are validated by Monte Carlo simulations for the Lennard-Jones
fluid. A comparison of the eight statistical ensembles regarding size scaling behavior,
convergence and stability is provided for state points throughout the homogeneous fluid
region. The resulting data are in good agreement, but differ in their statistical distri-
butions. In closed systems, the statistical quality of the data is better than in open
systems. Overall, the microcanonical ensemble performs best.
Introduction
Various properties of matter are accessible on the basis of interaction potentials via the
statistical mechanics route in eight natural ensembles. An ensemble comprises a large col-
lection of microstates of a mechanical system subject to a specific macrostate. For a simple
1
one-component system, it is specified by three independent variables which determine a ther-
modynamic potential. Gibbs, who coined the term ensemble in 1902,1was concerned with
systems of fixed volume. He defined the microcanonical ensemble (NV E), the canonical
ensemble (NV T)and the grand-canonical ensemble (µV T), where Ndenotes the number of
particles, Vthe volume, Ethe energy, Tthe temperature and µthe chemical potential.
In molecular simulation, closed systems with fixed Nare much easier to realize than open
systems, where insertion and deletion of particles are required. Accordingly, the NV T en-
semble was used for the first Monte Carlo (MC) simulations by Metropolis et al. in 1953.2Six
years later Alder and Wainwright3laid the foundation for molecular dynamics simulations in
the NV E ensemble. The other closed ensembles, which are widely in use for simulations up
to the present day, are the isobaric-isothermal ensemble (NpT)4and the isobaric-isoenthalpic
ensemble (NpH)5introduced by Guggenheim and Byers Brown,6respectively. Here, pde-
notes the pressure and H:= E+pV the enthalpy.
In the open systems, instead of the particle number, the chemical potential is kept con-
stant. Guggenheim also introduced the generalized ensemble (µpT)that he called completely
open.4Technically, this ensemble violates the requirement for independence of its three in-
tensive variables. However, for consistent (µ, p, T)triplets it can be treated as the other
ensembles,7although the system size remains undetermined. Therefore, it was hardly ever
used in molecular simulations.
The six physical quantities number of particles, temperature, energy, pressure, volume
and chemical potential may be placed on the vertices of an octahedron, defining the fixed
macroscopic constraints of eight ensembles as its faces, cf. Fig. 1. The variable Estands for
either of the thermodynamic combinations of energy E,H,L:= E−µN or R:= E+pV −µN,
the latter being specified by Hill8and Ray9and associated to the open systems with constant
energy and volume (µV L)or pressure (µpR), respectively.
Graben and Ray10 provided a joined perception of these eight ensembles by showing how
they mutually arise from one another by successive Legendre-Laplace transformations. For
2
Figure 1: Thermodynamic octahedron showing triplets of macroscopic constraints of eight
ensembles.10 North of the NpµV -plane the four ensembles are adiabatic, south they are
isothermal. Similar for TpEV-plane. To each face is associated a thermodynamic potential
and a statistical mechanical partition function.
each ensemble, derivatives of its thermodynamic potential with respect to the independent
variables yield all time-independent thermodynamic properties.
A rigorous methodology for the derivatives of the thermodynamic potential up to arbi-
trary order has been developed in the seminal work of Lustig11–14 for the NV E ensemble.
As a result, properties like isochoric and isobaric heat capacities, speed of sound, isothermal
compressibility, thermal expansion coefficient or Joule-Thomson coefficient can be expressed
in terms of kinetic and potential energy and volume derivatives of the potential energy.
Lustig employed the same methodology to the NV T ensemble,15,16 yielding expressions for
the Helmholtz energy up to arbitrary order and the thermodynamic properties.
A decade passed until Ströker et al. revisited this formalism and applied it to the NpT
ensemble17 such that derivatives of the Gibbs energy up to arbitrary order became avail-
able. Ströker and Meier continued this work and then provided analogous expressions first
for the µV T ensemble18 and shortly thereafter for the NpH ensemble.19 Recently, the un-
3
derlying procedure has been transferred to the remaining ensembles, i.e. µV L,µpR and
µpT.20 Consequently, the aforementioned properties are now available for all eight triplets
of macroscopic constraints.
In the present study, these eight ensembles are thoroughly compared using MC simula-
tions with rigorous acceptance criteria.14–20
Expressions for thermodynamic properties
Expressions for thermodynamic properties in terms of derivatives of the partition function in
all ensembles have been derived and discussed in the meticulous works of Lustig and Ströker
et al.11–20 Therefore, it is not intended to repeat the method in all details. Instead, an
attempt is made to provide a brief summary that offers generalized expressions applicable
to all ensembles.
In equilibrium, a homogeneous system of pure matter can be characterized by the ther-
modynamic potential S(E, V, N)as
dS=1
TdE+p
TdV−µ
TdN . (1a)
Manipulation involving series of one to three simultaneous Legendre transformations
with respect to the conjugate pairs of variables (E, 1/T),(V, p/T)and (N, −µ/T)yields a
thermodynamic potential as
dS=1
TdH−V
Tdp−µ
TdN , (1b)
dS=1
TdL+p
TdV+N
Tdµ , (1c)
dS=1
TdR−V
Tdp+N
Tdµ , (1d)
for the adiabatic ensembles in Fig. 1. For isothermal ensembles, the thermodynamic poten-
4
tials Helmholtz energy A:= E−TS, Gibbs energy G:= E+pV −TS, grand potential
J:= E−µN −TS and Null potential Z:= E+pV −µN −TS can be expressed in a reduced
form as21
−dA
T=−Ed1
T+p
TdV−µ
TdN , (1e)
−dG
T=−Ed1
T−Vdp
T−µ
TdN , (1f)
−dJ
T=−Ed1
T+p
TdV+Ndµ
T,(1g)
−dZ
T=−Ed1
T−Vdp
T+Ndµ
T.(1h)
These reduced thermodynamic potentials equate to the natural logarithm of a partition
function generated by corresponding Laplace transformations.
For a detailed discussion of the relations between the ensembles in terms of statistical
mechanics and the according Legendre-Laplace mapping scheme, the reader is referred to
Graben and Ray.10
To develop a generalized formulation of thermodynamic properties, an ensemble is de-
noted by (x, y, z)with y∈{p, V },z∈{µ, N}and x∈{E, T}, wherein E∈{E, H, L, R}. If
necessary, functional arguments identify the affiliation with the respective ensemble.
Adiabatic ensembles
First, the adiabatic cases are considered, i.e. x=E. In these four ensembles, each thermo-
dynamic potential corresponds to a different expression of the entropy S(E, y, z)according
to Eqs. (1a) to (1d). Entropy can be defined by either the phase space volume Ωor by the
5
phase space density ω
S(E, y, z)
kB
:= ln Ω(E, y, z),(2)
S(E, y, z)
kB
:= ln ω(E, y, z),(3)
where kBis the Boltzmann constant. Note that ωis given by the partial derivative of the
phase space volume with respect to energy
ω(E, y, z) = ∂Ω(E, y, z)
∂E.(4)
It has been shown that both entropy definitions (2) and (3) are equivalent in the ther-
modynamic limit.22,23 Since they lead to different expressions for thermodynamic properties,
the derivations in this work are restricted to the entropy definition (2). The second definition
is addressed in Refs.11,14,19,20
So-called phase space functions that are applicable to all four adiabatic ensembles arise
from partial derivatives of the phase space volume
Ωmno := 1
ω
∂m+n+oΩ
∂Em∂yn∂zo, m, n, o = 0,1,2. . . , (5)
and allow for the construction of the entire set of thermodynamic properties of a system, as
was first shown by Lustig for the NV E ensemble.11 Functional arguments are omitted for
notational convenience.
The total differentials of Sin Eqs. (1a) to (1d) yield first derivative thermodynamic
properties in terms of phase space functions. The inverse temperature is given by
1
kBT=!∂S/kB
∂E"y,z
=!∂ln Ω
∂E"y,z
=ω
Ω,(6)
6
such that
kBT=Ω000 .(7)
For por Vderivatives, complementary properties to yare defined as
yc:= #
$
$
%
$
$
&
−Vif y=p ,
pif y=V .
(8)
Then, with Eqs. (1a) to (1d), ycis obtained by
yc
kBT=!∂S/kB
∂y"E,z
=!∂ln Ω
∂y"E,z
=Ω010
Ω000
.(9)
Similarly, Narises in open ensembles, where z=µ, as
N
kBT=!∂S/kB
∂µ"E,y
=!∂ln Ω
∂µ"E,y
=Ω001
Ω000
.(10)
For the construction of other thermodynamic properties, derivatives of T,ycand Nwith
respect to E,yand µare needed. Conveniently, higher order derivatives of the phase space
functions Ωmno are recursive
∂Ωmno
∂E=Ωm+1,no −ΩmnoΩ200 ,(11)
∂Ωmno
∂y=Ωm,n+1,o −ΩmnoΩ110 ,(12)
∂Ωmno
∂µ=Ωmn,o+1 −ΩmnoΩ101 .(13)
In this work, only derivatives up to second order are considered. Explicit results are
summarized in Table S1 of the Supporting Information.
7
Isothermal ensembles
In the isothermal ensembles, where x=T, the thermodynamic potentials obtained by Legen-
dre transformation correspond to the Helmholtz energy A(T, V, N), Gibbs energy G(T, p, N),
grand potential J(T, V, µ)or Null potential Z(T, p, µ)according to Eqs. (1e) to (1h). The
latter is termed Null potential because both Z= dZ= 0 according to Euler’s theorem in
thermodynamics.21 However, this peculiarity is irrelevant for further derivations.
To allow for generalization, Φ∈{−A/(kBT),−G/(kBT),−J/(kBT),−Z/(kBT)}denotes
any thermodynamic potential reduced as in Eqs. (1e) to (1h). It is related to the according
partition function Qby
Φ(β, y, z) = ln Q(β, y, z),(14)
with inverse temperature β= (kBT)−1.
In analogy to their adiabatic counterparts Ωmno, the phase space functions have the
structure
Qmno := 1
Q
∂m+n+oQ
∂βm∂yn∂zo, m, n, o = 0,1,2. . . . (15)
While entropy cannot be related to Qmno, internal energy can according to Eqs. (1e) to
(1h)
E=−!∂Φ
∂β "y,z
=−!∂ln Q
∂β "y,z
=−Q100 .(16)
With the same notation as introduced above, pressure or volume arise by
βyc=!∂Φ
∂y"β,z
=!∂ln Q
∂y"β,z
=Q010 ,(17)
and in open systems, the number of particles is given by
βN=!∂Φ
∂µ"β,y
=!∂ln Q
∂µ"β,y
=Q001 .(18)
8
Figure 2: From any ensemble to any thermodynamic property.
Similar to the adiabatic case, higher order derivatives can be obtained recursively
∂Qmno
∂β =Qm+1,no −QmnoQ100 ,(19)
∂Qmno
∂y=Qm,n+1,o −QmnoQ010 ,(20)
∂Qmno
∂µ=Qmn,o+1 −QmnoQ001 .(21)
Generalized expressions for derivatives up to second order in terms of averages of any isother-
mal ensemble are listed in Table S2 of the Supporting Information.
As summarized in Fig. 2, the phase space functions Ωmno and Qmno are the building
9
blocks for the construction of any thermodynamic equilibrium property, like isochoric and
isobaric heat capacity cvand cp, thermal expansion coefficient αp, isothermal compressibility
βTand thermal pressure coefficient γvdefined as
cv=!∂E
∂T"V,N
, cp=!∂H
∂T"p,N
,αp=1
V!∂V
∂T"p,N
,
βT=−1
V!∂V
∂p"T,N
,γv=!∂p
∂T"V,N
.
(22)
These properties are related by
cv=cp−V TβTγ2
v,αp=βTγv,(23)
such that is is sufficient to compute only three of them if they are suitably chosen. Moreover,
the Joule-Thomson coefficient µJT and the speed of sound ware defined by
µJT =V(Tαp−1)
cp
, w2=V cp
NMβTcv
.(24)
The expressions in Eq. (22) can always be rewritten in terms of derivatives of the ther-
modynamic potential up to second order introduced above with the help of Jacobian deter-
minants21 to match any desired ensemble. To illustrate the procedure, expressions for cvare
given in Table S3 of the Supporting Information. Further specific results can be found in
Refs.15–20
Simulations
The statistical mechanical expressions for thermodynamic properties outlined in the previous
section were implemented into the simulation tool ms224 and were evaluated by means of
10
MC simulations for the classical Lennard-Jones fluid25,26
uij(rij) = 4ε'!σ
rij "12
−!σ
rij "6(.(25)
All reported data are reduced by the energy parameter εand the size parameter σsuch
that T∗=TkB/ε,ρ∗=ρσ3,p∗=pσ3/ε,E∗=E/ε,H∗=H/ε,c∗
v=cv/kB,c∗
p=cp/kB,
α∗
p=αpkB/ε,β∗
T=βTσ3/ε,µ∗
JT =µJT kB/σ3and w∗=w(MN/ε)1/2with molar mass M.
In the following, asterisks are omitted for brevity.
Throughout the homogeneous fluid region with 0.7< T < 6and ρ<1.25, a total of 167
state points were sampled with N= 1372 and N= 4000 particles, cf. Fig. 3. From specified
temperature-density pairs, the input parameters for the other ensembles were determined
from NV T simulation data. The cut-offradius was set to rc= 5σ, beyond which analytical
Figure 3: Investigated state points. For each state, MC simulations were carried out for two
system sizes (N= 1372 and N= 4000) in each of the eight ensembles. Crosses represent
states in Table 1 which were analyzed in more detail. A subdivision of the homogeneous
fluid region, as well as the vapor-liquid equilibrium curves and the freezing line are depicted
for orientation.
cut-offcorrections27 were applied. The chemical potential in the NV T ensemble was sampled
with Widom’s test particle method.28 For each ensemble, the simulations started with an
equilibration of 105cycles in the NV T ensemble and 3·105cycles in the respective ensemble
under study, followed by a production run of 106cycles. In one cycle, the ratio between
attempted translations and possible other attempted moves, like volume changes and particle
11
insertions or deletions, was on average Nto 1. To shorten the time-to-solution, one overall
simulation was divided into four independent copies of the system such that 2.5·105cycles
were performed in each.
Eight state points were selected for closer examination. Thus, their input parameters
were obtained from the equation of state (EOS) of Thol et al.,29 as listed in Table 1.
Table 1: Data set used as an input for selected state points obtained from the EOS of Thol
et al.29
#Tρp E/N H/N µ
1 Cold saturated liquid 0.7 0.844 0.022 -5.055 -5.026 -4.357
2 Warm saturated vapor 0.95 0.02 0.017 1.230 2.072 -3.937
3 Warm saturated liquid 1 0.75 0.400 -3.721 -3.187 -3.321
4 Dense hot gas 1.2 0.05 0.049 1.380 2.359 -4.044
5 Near-critical state 1.33 0.32 0.136 -0.330 0.110 -3.562
6 Hot supercritical liquid-like 2 0.75 4.008 -1.542 3.800 2.584
7 Very hot supercritical 4 0.50 3.522 3.396 10.438 1.792
8 Extremely hot and dense 6 1.25 84.901 11.607 79.514 86.394
Results and discussion
Detailed analysis of selected state points
For the eight selected state points listed in Table 1, besides average values at the end of a
simulation, coarse grained averages over 103cycles were taken, resulting in 103samples to
be histogrammed for distributions.
Since the eight ensembles rest on different sets of specified properties, their comparison
would be pointless. One property allowing for a mutual assessment of all ensembles is the
residual internal energy ur=)i<j uij, that is also generally sampled with great statistical
accuracy. The distributions for two state points 1 and 4 listed in Table 1 are shown in
Fig. 4. Other state points are depicted in section 2 of the Supporting Information. Each
plot provides a representation of all ensembles for two system sizes.
All simulations are in very good agreement with the EOS of Thol et al.29 Fig. 4 shows
12
Figure 4: Residual internal energy ur. Light color corresponds to N= 1372 and dark color
to N= 4000. The dashed lines inside the distribution indicate the median of the entire
dataset, as well as to the upper and lower quartile. The solid line represents the EOS of
Thol et al.29
13
that some internal energy distributions in the open ensembles µV T and µV L are bimodal.
Apparently, the systems attempt to occupy at least two distinct regions in phase space. Since
the state point 1 is close to the triple point, it can not be excluded that any additional phase
attempts to develop. The effect is more pronounced for the larger system size.
In Fig. 4, the distribution of µpT data at the dense hot gas state 4 is very broad for
the smaller system (N= 1372). As shown in Fig. 5, fluctuations of Nat this state point
are significantly larger than in the other ensembles. This might be expected, as its entirely
intensive ensemble constraints do not specify the system size. For dZ= 0, Eq. (1h) is
the Gibbs-Duhem equation, relating points on a thermodynamic equilibrium surface f(Z=
0, µ, p, T) = 0.21 If the triplet (µ, p, T )offered to a molecular simulation is consistent, Z= 0,
dZ= 0 and f= 0 for a stable output 〈N〉. If (µ, p, T)is (slightly) inconsistent, the
underlying Euler equation of thermodynamics is necessarily compromised to Z∕= 0.21 In
this situation, no stable 〈N〉can be maintained and a drift must be expected. However,
slight inconsistencies in (µ, p, T)have little effect on other thermodynamic properties.
Figure 5: Evolution of the number of particles Nin the open ensembles at the hot gas state
4. The total number of MC cycles was distributed between four independent copies of the
system that are represented by individual lines. Bold lines correspond to large systems,
initiated with N= 4000, and faint lines to small systems, initiated with N= 1372.
14
Parallel sampling of four independent systems highlights the element of randomness in
MC as they traverse very different paths in phase space. In some µpT simulation runs,
successive deletions of particles may have eventually culminated in the complete eradication
of the system. This was avoided by the introduction of an artificial lower bound for that
ensemble beyond which no further deletions were allowed. In the large systems that were
initiated with N= 4000, the number of particles fluctuated approximately by the same
amount as in those initiated with N= 1372. However, as collated in section 2 of the
Supporting Information, the fluctuations are similar for all open ensembles at the other
state points.
Statistical uncertainties were estimated with a block averaging method by Flyvbjerg
and Peterson.30 To quantify system size effects in all ensembles at the state points listed
in Table 1, ratios of statistical uncertainties for N= 4000 and 1372 were evaluated and
visualized for urin Fig. 6. Those ratios are mostly between 0.5 and 0.75, regardless of
ensemble and state point which corresponds to a size scaling that is inversely proportional
to √N.31 Results for all state points, density and pressure show a consistent behavior, cf.
Fig. S1 of the Supporting Information.
Figure 6: Ratio of statistical uncertainties of urfor N= 4000 and 1372. In case of the open
systems, Nstands for the initial particle number.
15
The extreme differences observed for gas states (state points 2 and 4) in the µpT ensemble
are outliers for possible reasons outlined above. The only other case where a similar result
was obtained, is the hot supercritical liquid-like state (state point 6) in the µpR ensemble.
The evolution of Nfor this particular state point reveals that for N= 1372, by chance two
of the four independent copies of the system evolved in opposite directions, cf. section 2.6
of the Supporting Information. Apparently, this behavior is attributable to details of one
particular numerical experiment.
The inverted ratio for the µV L ensemble at the near-critical state (state point 5) can be
similarly explained, cf. section 2.5 of the Supporting Information. While the smaller system
showed a tendency to increase Nby several hundred particles, the larger system exhibited the
opposite effect. Since the volume was constant, this led to a significant underestimation of
urby 47%in the first case and an overestimation by 20%, although with larger uncertainties.
Properties from first order derivatives
The overall performance of the eight ensembles is rather similar, such that the deviations
between the EOS of Thol et al.29 and simulation data are within statistical uncertainties, cf.
Fig. 7. For better perception, the homogeneous fluid region was divided into five sections,
cf. Fig. 3: the subcritical gas (I) and liquid (II) regions and the supercritical regions with
low (III), medium (IV) and high density (V).
Figure 7: Relative deviation of the residual internal energy urin % between simulation data
(N= 4000) and the EOS of Thol et al.29 .
16
Further insights are gained by the analysis of the relative statistical uncertainties in %,
i.e. δξ/ξfor any property ξ. These data are shown for the residual internal energy in Fig. 8.
For individual state points, they are equivalent to the width of distributions in Fig. 4.
Figure 8: Statistical uncertainties of the residual internal energy urin %(N= 4000).
Supercritical (SC) region: LD: ρ/ρc<0.6; MD: 0.6≤ρ/ρc≤1.5; HD: 1.5<ρ/ρc.
Clearly, the NV E ensemble is superior throughout the entire fluid region because its
unique feature that all ensemble constraints are extensive apparently promotes stability.
Opposed to that is the µV T ensemble, where the uncertainties are 10 to 40 times larger.
As expected, simulations in closed systems are more stable than in open systems, where the
variation of the particle number imposes major interferences.
Mostly, the adiabatic ensembles show lower uncertainties than their isothermal counter-
parts. A noteworthy exception is the µpT ensemble, which not only clearly outperforms all
other open ensembles, but even compares to the closed NpT ensemble.
For isochoric ensembles, an evaluation of the pressure is desirable. As can be seen in
Fig. S2 of the Supporting Information, the differences between open and closed systems are
insignificant. Closed systems mainly show deviations from the EOS of less than 3% in the
subcritical and less than 0.4% in the supercritical regions, whereas for the open systems,
they are below 4% and 0.5%, respectively. The average statistical uncertainties determined
for the residual pressure pr=p−ρkBTwere obtained by the error propagation law, cf.
Fig. S3. In the NV E and NV T ensembles, they are very similar and range from 0.02%
in the supercritical fluid with high density to 0.32%in the liquid. Density fluctuations in
17
the open systems lead to considerably larger uncertainties of up to 5%in the liquid. In
the supercritical fluid with low density, the uncertainty of 29%in the µV T ensemble is
inconclusive due to residual pressure values close to zero.
Conjugate to the pressure, the volume, or density as its inverse, can be assessed in the
isobaric systems, cf. Figs. S2 and S3 of the Supporting Information. The ensembles show
very similar results with up to 0.3%and 0.1%in the subcritical and supercritical regions,
respectively. The less stable µpR simulations deviate a little more in the supercritical region.
Statistical uncertainties are smallest in the NpH and slightly larger in the NpT and µpT
ensembles which again yield very similar results. The µpR ensemble performs poorest with
uncertainties up to 10 times larger for low densities than in the other ensembles.
Properties from second order derivatives
Properties related to second order derivatives are discussed using the isochoric heat capacity
cvas an example. Results for the other properties can be found in section 3 of the Supporting
Information. Fig. 9 shows the deviation of simulation results from the EOS of Thol et al.29
for the state points listed in Table 1. The near-critical state was excluded because the results
were inconclusive.
Overall, the MC results are rather similar. They deviate by less than 2%from the EOS
for w,4% for cv,5% for αpand βT, and up to 10%for cpand µJT . For properties discussed
above, the EOS is reliable enough to trace deviations of such magnitude to flaws in the
course of the simulations, but this is not the case here. Except for occasional outliers, the
results show little dependence on ensembles or state points. This is even more evident for
the statistical uncertainties in relative terms. The averages for the formerly introduced parts
of the fluid region are given in Fig. 10. With rising density, uncertainties increase from 0.2%
to 0.8%. The pattern is broken by the µpR ensemble which slightly underperforms for cv.
The speed of sound has the smallest uncertainties, mostly below 0.3%. In the gas region,
isochoric ensembles yield better results for the other properties, as opposed to liquid states.
18
Figure 9: Relative deviation of the isochoric heat capacity cvbetween simulation data and
EOS of Thol et al.29 (baseline) for N= 1372 (light symbols) and 4000 (dark symbols). In
case of the open systems, Nstands for the initial particle number.
Figure 10: Statistical uncertainties for the isochoric heat capacity cvin %(N= 4000).
Supercritical (SC) region: LD: ρ/ρc<0.6; MD: 0.6≤ρ/ρc≤1.5; HD: 1.5<ρ/ρc.
19
Obtaining more precise results by exploiting size effects proves futile, cf. Fig. 11. As
Fig. 9 suggests, improvement is only achieved in individual cases and size effects are usually
negligible. As shown in the Supporting Information, similar statements hold for all properties
of second order.
Figure 11: Ratio of statistical uncertainties of cvfor N= 1372 and 4000. In case of the open
systems, Nstands for the initial particle number.
Conclusion
Thermodynamic equilibrium properties in eight statistical ensembles were obtained from
derivatives of the respective thermodynamic potential. The construction of fluctuation prop-
erties, like isochoric and isobaric heat capacity, thermal expansion coefficient, isothermal
compressibility, thermal pressure coefficient, speed of sound and Joule-Thomson coefficient
involves derivatives up to the second order. A generalized methodology based on the formal-
ism developed by Lustig was outlined to express these derivatives in terms of phase space
functions in adiabatic and isothermal ensembles.
Rigorous expression as combinations of ensemble averages were implemented into the
simulation code ms2 and validated for a large range of states in the homogenous region of
the Lennard-Jones fluid with Monte Carlo simulation for two system sizes. For eight selected
20
state points, the distribution of the obtained data throughout a simulation run was assessed.
The resulting data were analyzed regarding size scaling behavior, convergence and stability.
In general, a good agreement with the EOS of Thol et al. was observed for all ensembles
for systems with N= 4000. The size scaling was as expected, except for the open isobaric
ensembles, which experienced convergence issues in different fluid regions for N= 1372.
Overall, the NV E ensemble with entirely extensive constraints yields data with the best
statistical quality.
In the µpT ensemble, inconsistency of the input triplet may cause a random walk to
become unstable. Such inconsistencies may be accidentally compensated for by varying
initial system size. For second order derivatives, possible (µ, p, T)-inconsistencies could not
be detected. For consistent input triplets, the µpT ensemble clearly outperforms the other
open ensembles and is comparative to the NpT ensemble regarding the resulting data.
Acknowledgement
Molecular simulations were performed either on the Atos BullSequana XH2000 system Noc-
tua 2 at the Paderborn Center for Parallel Computing (PC2) or on the HPE Apollo system
Hawk at the High Performance Computing Centre Stuttgart (HLRS) contributing to the
project MMHBF2.
Supporting Information Available
The following files are available free of charge.
• SI_phase_space_funs.pdf: expressions for phase space functions in terms of ensemble
averages and expressions for cvin terms of phase space functions
• SI_additional_plots.pdf: deviation from EOS, average statistical uncertainties and
size effects for properties from first and second order derivatives, distribution of ur,T,
21
pand ρ, evolution of Nfor selected states, deviation from EOS
References
(1) Gibbs, J. Elementary Principles in Statistical Mechanics: Developed with Especial Ref-
erence to the Rational Foundations of Thermodynamics; Elementary Principles in Sta-
tistical Mechanics: Developed with Especial Reference to the Rational Foundation of
Thermodynamics; C. Scribner’s sons, 1902.
(2) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N. Equation of State Calculations
by Fast Computing Machines. 1953,21, 1087.
(3) Alder, B. J.; Wainwright, T. E. Studies in Molecular Dynamics. I. General Method.
The Journal of Chemical Physics 1959,31, 459–466.
(4) Guggenheim, E. A. Grand Partition Functions and So-Called ”Thermodynamic Prob-
ability”. The Journal of Chemical Physics 1939,7, 103–107.
(5) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temper-
ature. The Journal of Chemical Physics 1980,72, 2384–2393.
(6) Brown, W. B. Constant pressure ensembles in statistical mechanics. Molecular Physics
1958,1, 68–82.
(7) Sack, R. Pressure-dependent partition functions. Molecular Physics 1959,2, 8–22.
(8) Hill, T. An Introduction to Statistical Thermodynamics; Addison-Wesley series in chem-
istry; Dover Publications, 1986.
(9) Ray, J. R.; Graben, H. W. Fourth adiabatic ensemble. The Journal of Chemical Physics
1990,93, 4296–4298.
22
(10) Graben, H.; Ray, J. R. Eight physical systems of thermodynamics, statistical mechanics,
and computer simulations. Molecular Physics 1993,80, 1183–1193.
(11) Lustig, R. Statistical thermodynamics in the classical molecular dynamics ensemble. I.
Fundamentals. The Journal of Chemical Physics 1994,100, 3048–3059.
(12) Lustig, R. Statistical thermodynamics in the classical molecular dynamics ensemble.
II. Application to computer simulation. The Journal of Chemical Physics 1994,100,
3060–3067.
(13) Lustig, R. Statistical thermodynamics in the classical molecular dynamics ensemble.
III. Numerical results. The Journal of Chemical Physics 1994,100, 3068–3078.
(14) Lustig, R. Microcanonical Monte Carlo simulation of thermodynamic properties. The
Journal of Chemical Physics 1998,109, 8816–8828.
(15) Lustig, R. Direct molecular NVT simulation of the isobaric heat capacity, speed of
sound and Joule-Thomson coefficient. Molecular Simulation 2011,37, 457–465.
(16) Statistical analogues for fundamental equation of state derivatives. 110 .
(17) Ströker, P.; Hellmann, R.; Meier, K. Systematic formulation of thermodynamic prop-
erties in the NpT ensemble. Physical Review E 2021,103, 023305.
(18) Ströker, P.; Meier, K. Classical statistical mechanics in the grand canonical ensemble.
Physical Review E 2021,104, 014117.
(19) Ströker, P.; Meier, K. Rigorous expressions for thermodynamic properties in the NpH
ensemble. Physical Review E 2022,105, 035301.
(20) Ströker, P. Bestimmung thermodynamischer Eigenschaften von Fluiden mit einer
weiterentwickelten molekularen Simulationsmethodik und hochgenauen ab initio-
Potentialen. Ph.D. thesis, Helmut Schmidt University, 2023.
23
(21) Münster, A. Classical thermodynamics; Wiley-Interscience: Chichester, New York,
1970.
(22) Becker, R. Theory of Heat; Springer Berlin Heidelberg: Berlin, Heidelberg, 1967.
(23) Pearson, E. M.; Halicioglu, T.; Tiller, W. A. Laplace-transform technique for deriving
thermodynamic equations from the classical microcanonical ensemble. Physical Review
A1985,32, 3030–3039.
(24) Fingerhut, R.; Guevara-Carrion, G.; Nitzke, I.; Saric, D.; Marx, J.; Langenbach, K.;
Prokopev, S.; Celný, D.; Bernreuther, M.; Stephan, S.; Kohns, M.; Hasse, H.; Vrabec, J.
ms2: A molecular simulation tool for thermodynamic properties, release 4.0. Computer
Physics Communications 2021,262, 107860.
(25) Jones, J. E. On the determination of molecular fields.-I. From the variation of the
viscosity of a gas with temperature. 1924,106, 441–462.
(26) Jones, J. E. On the determination of molecular fields.-II. From the equation of state of
a gas. 1924,106, 463–477.
(27) Allen, P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press, 1987.
(28) Widom, B. Some Topics in the Theory of Fluids. The Journal of Chemical Physics
1963,39, 2808–2812.
(29) Thol, M.; Rutkai, G.; Köster, A.; Lustig, R.; Span, R.; Vrabec, J. Equation of State for
the Lennard-Jones Fluid. Journal of Physical and Chemical Reference Data 2016,45,
023101.
(30) Flyvbjerg, H.; Petersen, H. G. Error estimates on averages of correlated data. The
Journal of Chemical Physics 1989,91, 461–466.
24
(31) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Ap-
plications, 2nd ed.; Computational Science Series; Academic Press: San Diego, 2002;
Vol. 1.
25
Graphical TOC Entry
26
SUPPORTING INFORMATION
Numerical discrimination of thermodynamic Monte Carlo
simulations in all eight statistical ensembles
Isabel Nitzkeaand Jadran Vrabeca
aThermodynamics, Technische Universtit¨at Berlin, 10587 Berlin, Germany
Phase space functions in terms of ensemble averages
To lighten notation, K=3N
2−1(E − ˆ
E)−1is introduced. Therein, Eis the fixed
energy of the respective ensemble and ˆ
Eis its configurational counterpart, wherein E
is replaced by the potential energy Uobtained from simulation.
Table S1. Phase space functions in terms of ensemble averages for
adiabatic systems.
(E, y, z) (E, p, z) (E, V, z)
Ω000 D2
3N(E − ˆ
E)E
Ω100 1
Ω200 ⟨K⟩
Ω010 −⟨V⟩D2
3V(E − ˆ
E)E−D∂U
∂V E
Ω020 KV 2D2(N−1)
3V2(E − ˆ
E)E−2DN
V
∂U
∂V E
−D∂2U
∂V 2E+K∂U
∂V 2
Ω001 ⟨N⟩
Ω002 KN2
Ω110 − ⟨KV ⟩ − DK∂U
∂V E+DN
VE
Ω101 ⟨KN⟩
Ω011 − ⟨KNV ⟩ − DK∂U
∂V NE+DN2
VE
Table S2. Phase space functions in terms of ensemble averages for isothermal systems.
(T, y, z) (T, p, z) (T, V, z)
Q000 1
Q100 −β−1D3N
2E− ⟨ ˆ
E⟩
Q200
β−2D3N
23N
2−1E
+3β−1⟨ˆ
EN⟩+⟨ˆ
E2⟩
Q010 β⟨V⟩DN
VE−βD∂U
∂V E
Q020 β2⟨V2⟩
⟨N2⟩−⟨N⟩
V2−2β
VDN∂U
∂V E
−βD∂2U
∂V 2E+β2∂U
∂V 2
Q001 β⟨N⟩
Q002 β2⟨N⟩2
Q110 D3N
2−1VE+β⟨ˆ
EV⟩−2
3β−1DN2
VE−1
V⟨ˆ
EN⟩
+βDˆ
E∂U
∂V E+D3N
2−1∂U
∂V E
Q101 D1−3N
2NE−β⟨ˆ
EN⟩
Q011 −β2⟨V N⟩βDN2
VE−β2DN∂U
∂V E
Table S3. Isochoric heat capacity cvin terms of phase space functions ins all ensembles. In closed ensembles, index
o= 0 is omitted for brevity.
NV E cv=kB(1 −Ω00Ω20)−1
NV T cv=kBβ2(Q20 −Q2
01)
NpH cv=kB(1 −Ω00Ω20)−1"1 + Ω00(Ω11−Ω20Ω01)2
(1−Ω00Ω20)(Ω11Ω01−Ω02)−(Ω01−Ω00Ω11)(Ω20Ω01−Ω11)#
NpT cv=kB"β2Q20 −Q2
10−Q01−β(Q11−Q10Q01)2
Q02−Q2
01 #
µV L cv=kBΩ002−2Ω001Ω101+Ω200Ω2
001
(1−Ω000Ω200)(Ω002−Ω001Ω101)−(Ω101−Ω001Ω200)(Ω001−Ω000Ω101)
µV T cv=kB"β2Q200 −Q2
100−Q001−β(Q101−Q100Q001)2
Q002−Q2
001 #
µpR cv=kBhΩ002−2Ω001Ω101+Ω200Ω2
001
(1−Ω000Ω200)(Ω002−Ω001Ω101)−(Ω101−Ω001Ω200)(Ω001−Ω000Ω101)
+"(1−Ω000Ω200)(Ω002−Ω001Ω101)−(Ω101−Ω001Ω200)(Ω001−Ω000Ω101)
Ω000(Ω110−Ω010Ω200)(Ω101Ω001−Ω002)−(Ω101−Ω001Ω200)(Ω010Ω101−Ω011)2
×(Ω010Ω110 −Ω020)(1 −Ω000Ω200)(Ω002 −Ω001Ω101)−(Ω101 −Ω001Ω200)(Ω001 −Ω000Ω101)
−(Ω010Ω200 −Ω110)(Ω010 −Ω000Ω110)(Ω002 −Ω001Ω101)−(Ω011 −Ω001Ω110)(Ω001 −Ω000Ω101)
+(Ω010Ω101 −Ω011)(Ω010 −Ω000Ω110)(Ω101 −Ω001Ω200)−(Ω011 −Ω001Ω110)(1 −Ω000Ω200)i#−1#
µpT cv=kB"β2Q200 −Q2
100−Q001−β(Q101−Q100Q001)2
Q002−Q2
001
−hQ010−β(Q110−Q100Q010)−(Q011−Q010Q001)Q001−β(Q101−Q100Q001)
Q002−Q2
001 i2
Q020−Q2
010−(Q011−Q010Q001)2
Q002−Q2
001 #
2
SUPPORTING INFORMATION
Numerical discrimination of thermodynamic Monte Carlo
simulations in all eight statistical ensembles
Isabel Nitzkeaand Jadran Vrabeca
aThermodynamics, Technische Universtit¨at Berlin, 10587 Berlin, Germany
1 Size scaling, deviation from EOS and statistical uncertainties
Figure S1. Ratio of statistical uncertainties for reduced internal energy ur(top), pressure p(left) and density
⇢(right) in for N= 4000 and 1372. In case of the open systems, Nindicates the initial particle number.
Figure S2. Relative deviation from EOS for pressure p(left) and density ⇢(right) in %. Supercritical region:
LD: ⇢/⇢c<0.6; MD: 0.6⇢/⇢c1.5; HD: 1.5<⇢/⇢
c.
Figure S3. Statistical uncertainties for residual pressure pr(left) and density ⇢(right) in %. Supercritical
region: LD: ⇢/⇢c<0.6; MD: 0.6⇢/⇢c1.5; HD: 1.5<⇢/⇢
c
2