Molecular dynamics 243
be modified to include a thermostat on each particle or on each degree of freedom
(“massive” thermostatting), as discussed in Section 4.10.
To illustrate the use of eqns. (5.9.5), consider the simple example of aparticle of
massmmoving in a one-dimensional box with lengthLsubject to periodic poten-
tial. Letpandqbe the momentum and coordinate of the particle, respectively. The
potential is given by
U(q,L) =
mω^2 L^2
4 π^2
[
1 −cos
(
2 πq
L
)]
, (5.9.8)
whereωis a parameter having units of inverse time. Such a potential could beused,
for example, as a simple model for the motion of particles through a nanowire. We will
use eqns. (5.9.5) to determine the position and box-length distributions for a given
pressurePand temperatureT. These distributions are given by
P(q)∝
∫∞
q
dLexp [−βPL] exp
{
−β
mω^2 L^2
4 π^2
[
1 −cos
(
2 πq
L
)]}
P(L)∝exp [−βPL]
∫L
0
dqexp
{
−β
mω^2 L^2
4 π^2
[
1 −cos
(
2 πq
L
)]}
∝Lexp [−βPL]
∫ 1
0
dsexp
{
−β
mω^2 L^2
4 π^2
[1−cos (2πs)]
}
, (5.9.9)
where the last line is obtained by introducing the scaled coordinates=q/L. The one-
dimensional integrals can be performed using a standard numericalquadrature scheme,
yielding “analytical” distributions that can be compared to the simulated ones. The
simulations are carried out using a numerical integration scheme that we will present
in Section 5.12. Fig. 5.3 shows the comparison for the specific case thatω= 1,m= 1,
kT= 1, andP= 1. The parameters of the simulation are:W= 18,M= 4,Qj= 1,
Q′j= 1, and ∆t= 0.05. It can be seen that the simulated and analytical distributions
match extremely well, indicating that eqns. (5.9.5) generate the correct phase space
distribution.
5.10 Molecular dynamics in the isothermal-isobaric ensemble II:
Anisotropic cell fluctuations
Suppose we wish to map out the space of stable crystal structures for a given substance.
We can only do this within a molecular dynamics framework if we can sample different
cell shapes. For this reason, the development of molecular dynamics approaches with
a fully flexible cell or box is an extremely important problem. We have already laid
the groundwork in the isotropic scheme developed above and in our derivation of the
pressure tensor estimator in eqn. (5.7.1). The key modification we need here is that
the nine components of the box matrixhmust be treated as dynamical variables with
nine corresponding momenta. Moreover, we must devise a set of equations of motion
whose compressibility leads to the metric factor
√
g= [det(h)]^1 −dexp
[
dNη 1 +ηc+d^2 ξ 1 +ξc