Constraints: The ROLL algorithm 2550 10 20 30 40 50 60 70
t (ps)36
38
40
42
44
46
48
L (Å)
0 2 4 6 8 10
r (Å)0
0.5
1
1.5
2
g(r)P = 0.5 kbar
P = 1.0 kbar
P = 1.5 kbar(a) (b)
Fig. 5.4(a) Box-length fluctuations at pressures ofP= 0.5 kbar (top curve),P = 1. 0
kbar (middle curve), andP= 1.5 kbar (bottom curve), respectively. (b) Radial distribution
functions at each of the three pressures.
5.12.1 Example: Liquid argon at constant pressure
As an illustrative example of molecular dynamics in the isothermal-isobaric ensemble,
we consider first the argon system of Section 3.14.2. Three simulations at applied
external pressures of 0.5 kbar, 1.0 kbar, and 1.5 kbar and a temperature of 300 K
are carried out, and the radial distribution functions computed ateach pressure. The
parameters of the Lennard-Jones potential are described in Section 3.14.2, together
with the integration time step used in eqn. (5.12.4). Temperature control is achieved
using the “massive” Nos ́e-Hoover chain scheme of Section 4.10. The values ofτfor the
particle and barostat Nos ́e-Hoover chains are 100.0 fs and 1000.0fs, respectively, while
τb= 500.0 fs. Nos ́e-Hoover chains of lengthM= 4 are employed usingnsy= 7 and
n= 4 in eqn. (4.11.16). Each simulation is 75 ps in length and carried out in acubic
box with periodic boundary conditions subject only to isotropic volume fluctuations.
In Fig. 5.4(a), we show the fluctuations in the box length at each pressure, while in
Fig. 5.4(b), we show the radial distribution functions obtained at each pressure. Both
panels exemplify the expected behavior of the system. As the pressure increases, the
box length decreases. Similarly, as the pressure increases, the liquid becomes more
structured, and the first and second peaks in the radial distribution function become
sharper. Fig. 5.5 shows the density distribution (in reduced units) obtained form the
simulation atP = 0.5 kbar (P∗=Pσ^3 /ǫ= 1.279). The solid and dashed curves
correspond toτbvalues of 500.0 fs and 5000.0 fs, respectively. It can be seen that the
distribution is fairly sharply peaked in both cases around a density valueρ∗≈ 0 .704,
and that the distribution is only sensitive to the value ofτbnear the peak. Interestingly,
the distribution can be fit very accurately to a Gaussian form,
PG(ρ∗) =1
√
2 πσ^2e−(ρ∗−ρ 0 ) (^2) / 2 σ 2
(5.12.29)
with a widthσ= 0.01596 and averageρ 0 = 0.7038. Such a fit is shown in circles on
the solid curve in Fig. 5.5.