MBWR
MBWR
MBWR
Corresponding author.
E-mail address: may@fbm.fh-darmstadt.de (H.-O. May).
0378-3812/$ see front matter 2003 Elsevier B.V. All rights reserved.
doi:10.1016/S0378-3812(03)00283-8
2 P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19
feature in the interaction of these associating systems is the occurrence of strong, directional bonds, in case
of water hydrogen bonds. These bonds are responsible for highly orientational-dependent interactions.
Very surprisingly, the prospect of a second CP and related density anomaly can also be found for
orientational-independent interactions of simple liquids, when their potentials have a region of negative
curvature in their repulsive core. These so-called core-softened (CS) potentials are used in many ap-
plications of condensed matter physics. CS models are relevant to a number of liquid metals [6] as well
as for colloids [7].
First studies on CS potentials were published by Hemmer and Stell more than 30 years ago [8]. A very
simple CSpotential, the Gaussian core model (GCM), was introduced more than 25 years ago by Stillinger
[9]. It could be shown by many researchers that water-like anomalies exist in this model liquid [10,11]. In
recent years, the GCMbecame very popular in the eld of soft condensed matter physics [12]. Aso-called
CSramp potential was introduced by Jagla [13] and analysed in three dimensions (3D) by a Monte Carlo
(MC) simulation. A stable second CP could be found as well as a density anomaly. Head-Gordon and
Stillinger [14] showed that the inversion of an experimentally found oxygenoxygen radial distribution
function of water leads again to a CS potential. It has not been tested until now if this experimental
CS potential leads to anomalies and even a second CP. Instead of this, Stanley and co-workers [1518]
devoted signicant effort to a similar CS potential recently, which we will call (following Wilding and
Magee [19]) CS shoulder (CSS) model. Stanley and co-workers argued that such a potential is able
to mimic the effect of bonding and that it could be therefore seen as a zero-order approximation for
water. In the case of 1D and for the discrete form of the CSS potential, the phase behaviour could be
calculated analytically [16]. Molecular dynamics (MD) simulations in 2Dwere used for the simulation of
the discrete and smooth form of the potential [18]. These studies reveal a phase diagram similar to that of
water. The model promotes a local correlation between low-density and low-energy states succeeding in
reproducing anomalies for 1Dand 2Dwhich were attributed to the existence of a metastable liquidliquid
CP [15].
These results were seen in perspective in a recent paper of Wilding and Magee [19]. They argue that the
thermodynamic anomalies in the 2D CSS model are connected to pseudo-critical uctuations associated
with the 2D freezing transition. Moreover, they say that this interpretation is consistent with the recently
discussed results of a discrete so-called generic CS model [20]. Stanley and co-workers show [21] that
two rst-order uiduid phase transitions exist for this model, but that a density anomaly is absent. This
is related to the presence of only one stable crystal structure.
All these ndings suggest immediately the great inuence of the potentials shape and parameters
on the phase behaviour and the appearance of liquid-state anomalies. Moreover, the anomalous phase
behaviour of a CS systemis inuenced by the dimension of the uid. This can be seen fromsimple general
thermodynamic arguments given from Debenedetti et al. [22]. They show that the relation (/T)
<
2Dk
B
(where k
B
is Boltzmann constant and D the dimensionality of the system) must be fullled for
the virial if the thermal expansion coefcient is negative which is possible for a CS system.
A detailed study of the gasliquid phase behaviour for the 3D case of the CSS potential is still lacking
and seems to be very interesting. The aim of this paper is the lling of this gap. Using intensive MC
simulations for a N, V, T-ensemble and a Gibbs-ensemble (GE) simulation [23] for the determination of
the gasliquid coexistence line (CL) and of the gasliquid CP, we develop an equation of state (EOS)
by means of a modied Benedict Webb Rubin (MBWR) equation with 33 adjustable parameters. The
MBWR-equation has been used successfully by Nicolas et al. [24] and by Johnson et al. [25] for tting
the thermodynamic data of a LennardJones uid. The most data Johnson used belong to the supercritical
P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19 3
region. It is an interesting question how far this equation giving a complete thermodynamic description
could be applied to a more sophisticated system like a CSS model, especially for the metastable region.
In the following, we describe a MC simulation for a N, V, T-ensemble and a GE, and the method
for the evaluation of a MBWR-EOS for the CSS model. The formulas of the MBWR-equation and the
parameters resulting from the regression will be given. We close with a discussion of our results.
2. Monte Carlo simulation
Our CSS potential is obtained by adding a Gaussian well to the LennardJones potential:
v(r) = 4
_
_
r
_
12
r
_
6
_
exp
_
_
w
_
r
r
o
__
2
_
(1)
where and are used as units of length and energy. This leads to reduced units v
= v/, r
= r/, p
3
p/,
=
3
, T
= k
B
T/ for potential energy, distance between two particles, pressure, density and
temperature (k
B
is Boltzmann constant). The free parameters are taken from [18] with = 1.7, w = 5.0,
and r
o
= 1.5. Notice that our denition for the reduced density is different from [18]. These settings
are usually applied for a LennardJones potential and a direct comparison with this potential is therefore
possible.
To develop an EOS for the CSS liquid in 3D we made standard N, V, T MC simulations. In order to
examine our simulation code we performed MC calculations in two dimensions and compared our results
with 2D MD simulations taken from [18]. Our results showed exactly the same anomalous behaviour as
described in [18].
For nding the CP in 3D we started our 3D simulation with the calculation of the gasliquid CL by
means of a GE simulation. We performed this simulation on two different systems with N = 512 and
1024 particles. The resulting gasliquid CL is shown in the T
c
= 3.8 so that the truncation error for the calculation
of pressure and energy could assumed to be small. The rst temperature was T
= 0.025 until
= 0.2 until T
projection of the
4 P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19
Fig. 1. Pressuredensity isotherms, gasliquid coexistence line and spinodal line from MC simulations as well as the corre-
sponding MBWR-EOS t. (a) Several isotherms (top to bottom) for T
projection of the gasliquid coexistence line obtained from GE simulations for N = 512 (squares) and N = 1024 (circles)
particle simulation run. CP is the gasliquid critical point. The thick solid line is a t of the simulation data by means of a law of
rectilinear diameters where the critical exponent was 0.32. The thin solid line represents a constant Z = 1 line, where Z denotes
the compressibility factor. (b) Enlarged view around CP, symbols as in (a).
P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19 5
Fig. 2. Congurational energydensity isotherms for T
= 1.6, . . . , 3.2. Symbols are as in Fig. 1. The thin solid lines represent
the MBWR-EOS t of the reduced energy excluding the unstable region. The different energy isotherms are shifted by 10(T
c
)
each in order to avoid overlaps.
phase diagram in an enlarged view of Fig. 1b. The reduced congurational energy u
obtained in our MC
simulation is given in Fig. 2. Each curve has been shifted by a factor of 10(T
c
) so that they are clearly
distinguishable.
3. Equation of state
The MBWR-equation used here is the same as described by Nicolas et al. [24] and Johnson et al. [25]
with 32 linear parameters x
i
and one non-linear parameter :
p
, T
) =
+
2
(x
1
T
+x
2
T
1/2
+x
3
+x
4
T
1
+x
5
T
2
)
+
3
(x
6
T
+x
7
+x
8
T
1
+x
9
T
2
) +
4
(x
10
T
+x
11
+ x
12
T
1
)
+
5
(x
13
) +
6
(x
14
T
1
+x
15
T
2
) +
7
(x
16
T
1
) +
8
(x
17
T
1
+x
18
T
2
)
+
9
(x
19
T
2
) +exp(
2
) [
3
(x
20
T
2
+x
21
T
3
) +
5
(x
22
T
2
+x
23
T
4
)
+
7
(x
24
T
2
+x
25
T
3
) +
9
(x
26
T
2
+x
27
T
4
) +
11
(x
28
T
2
+x
29
T
3
)
+
13
(x
30
T
2
+x
31
T
3
+x
32
T
4
)] (2)
The reduced congurational energy u
, T
) =
_
0
d
2
_
p
_
p
_
(3)
6 P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19
Table 1
Parameters for the MBWR-EOS regressed from the simulation data
x
1
= 1.25160592890 x
2
= 18.3765471239 x
3
= 47.2084083590
x
4
= 27.2586402440 x
5
= 26.6048675417 x
6
= 8.30436023454
x
7
= 123.260018043 x
8
= 797.059526607 x
9
= 4.66312138808E6
x
10
= 16.2253941449 x
11
= 352.620691527 x
12
= 2262.02875340
x
13
= 122.732866058 x
14
= 6133.71601777 x
15
= 357977.540376
x
16
= 9961.12284242 x
17
= 5212.04976785 x
18
= 2.56816784576E6
x
19
= 1.55152549900E6 x
20
= 4.66453560123E6 x
21
= 735.689215836
x
22
= 1.39396961868E7 x
23
= 2640.58426419 x
24
= 1.95597531234E7
x
25
= 3994.79185494 x
26
= 1.59651718894E7 x
27
= 29731.5751161
x
28
= 8.68274937307E6 x
29
= 90698.7713399 x
30
= 3.24193893121E6
x
31
= 121064.734813 x
32
= 45568.1751952
and the reduced second virial coefcient is:
B
2
=
3
2
(x
1
+x
2
T
1/2
+x
3
T
1
+x
4
T
2
+x
5
T
3
) (4)
Moreover, the conditions for the critical point are:
_
p
_
T
c
,T
c
= 0 and
_
2
p
2
_
T
c
,T
c
= 0 (5)
Our procedure of calculating the unknown parameters was very similar to Johnson et al. [25]. First, we
have generated exact virial coefcients at 440 temperatures with Mayers f
12
function and used these
data to determine the rst ve parameters x
1
x
5
of Eq. (4) by regression. Then, we used the MC data
to regress the remaining parameters. We estimated the slope of the spinodal (SP) line and we omitted
the values belonging to the unstable region. In order to avoid unphysical regions we only used the
simulation data in the temperature range from T
and u
, T
) and in
Fig. 2 for the reduced energy u
, T
projection, we
used Maxwells construction.
P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19 7
As shown in Figs. 1 and 2, the agreement between simulated results and MBWR-EOS for the pressure
and the reduced energy is reasonably good within the range of temperature and density we have studied.
The slope of pressure in Fig. 1 indicates clearly that the CSS model promotes no density anomaly in 3D
because the pressure isotherms do not cross each other, which is an important difference to the 1Dand 2D
case [18]. The compressibility factor Z = p
_
T
(6)
and the thermal expansion coefcient
p
= K
T
_
p
(7)
Fig. 3. Isothermal compressibility K
T
(a) and thermal expansion coefcient
p
(b) along various isobars for p
= 0.6, . . . , 4.0.
K
T
and
p
are calculated using the MBWR-EOS t.
8 P. Mausbach, H.-O. May / Fluid Phase Equilibria 214 (2003) 19
can be calculated directly from the MBWR-EOS. In water, the isothermal compressibility is expected to
increase upon isobaric cooling while the thermal expansion coefcient becomes negative. We show K
T
and
p
for the CSS model along various isobars in Fig. 3, where the pressure is varied fromp
= 0.6 to
4. K
T
always increases upon increasing T