Liquid structure 207Note,{s(∆t/2),s ̇(∆t/2)}are evaluated by substituting the present momenta and
forces into eqns. (4.12.20) witht= ∆t/2. The symbol, “←−,” indicates that on the
computer, the values on the left-hand side are overwritten in memory by the values
on the right-hand side.
The isokinetic ensemble method has recently been shown to be a useful method
for generating a canonical coordinate distribution. First, it is a remarkably stable
method, allowing very long time steps to be used, particularly when combined with
the RESPA scheme. Unfortunately, the isokinetic approach suffers from some of the
pathologies of the Nos ́e–Hoover approach so some care is neededwhen applying it.
Minaryet al.(2004b) showed that such problems can be circumvented by combining the
isokinetic and Nos ́e–Hoover chain approaches (Minaryet al., 2004b). Such approaches
also eliminate resonance problems in the RESPA algorithm (see Problem3.6).
4.13 Applying canonical molecular dynamics: Liquid structure
Figures 4.2 and 4.3 showed radial distributions functions for liquid argon and water,
respectively. The importance of the radial distribution function in understanding the
structure of liquids and approximating their thermodynamic properties was discussed
in Section 4.6. In this section, we will describe how these plots can be extracted from a
molecular dynamics trajectory. Since the radial distribution function is an equilibrium
property, it is appropriate to employ a canonical sampling method such as Nos ́e–Hoover
chains or the isokinetic ensemble for this purpose.
The argon system represented in Fig. 4.2 was simulated using the “massive” Nos ́e–
Hoover chain approach on the argon system described in Section 3.14.2. The ther-
mostats maintained the system at a temperature of 300 K by controlling the kinetic
energy fluctuations. The system was integrated for a total of 10^5 steps using a time
step of 10.0 fs. Each Cartesian degree of freedom of each particlewas coupled to its own
Nose ́e–Hoover chain thermostat withM= 4, usingnsy= 7 andn= 4 in the Suzuki–
Yoshida integration scheme of eqn. (4.11.15). The parameterτused to determine the
value ofQ 1 ,...,Q 4 was taken as 200.0 fs.
The water system represented in Fig. 4.3 was simulated using, once again, the
“massive” Nos ́e–Hoover chain approach on a system of 64 water molecules in a cubic
box of length 12.4164 ̊A subject to periodic boundary conditions. The forces were
obtained directly from density functional theory electronic structure calculations per-
formed at each molecular dynamics step via the Car–Parrinello approach (Car and
Parrinello, 1985; Tuckerman, 2002; Marx and Hutter, 2009). Thesystem was main-
tained at a temperature of 300 K using a time step of 0.1 fs. For the thermostats, the
following parameters were used:nsy= 7,n= 4, andτ= 20 fs. The system was run
for a total of 60 ps.
After the molecular dynamics calculation has been performed, the trajectory is
subsequently used to compute the radial distribution function using the following
algorithm:
- Divide the radial interval betweenr= 0 andr=rmax, wherermaxis some radial
value beyond which no significant structure exists, intoNrintervals of length ∆r.
It is important to note that the largest valuermaxcan have is half the length of