Lecture Notes

Consider the integral of interest

                                  (1)

which can be written as

                         (2)

?(x) = an arbitrary probability function; for example, a uniform distribution in the range of (x₁, x₂) may be used:

                                 (3)

Imagine we perform a total of T trials, each consisting of choosing a random number R_i from the distribution ?(x) in the range (x₁, x₂). Then:

     (4)

The integral of equation (1) can thus be estimated by a numerical sampling procedure (eq 4).

Example: The shaded area of a circle with a unit radius, which is p, is determined by the integral

Following equation (4)

     (5)

Exercise: Write a simple Fortran program to estimate the value of p by carrying out 10⁷ uniform sampling.

For classical systems, the mechanical (i.e., microscopic) state of the system is characterized by points in phase space

v(r^N, p^N) = (r₁, r₂, ..., r_N, p₁, p₂, ..., p_N)

where r_i and p_i are the coordinates and conjugate momenta, respectively, for particle I. The flow in this space is determined by the time integration of Newton’s equation of motion. A basic concept in statistical mechanics is that if we wait long enough, the system will eventually travel trough all the microscopic states with constraints we imposed to control the system. The assembly of all possible microstates - all states consistent with these constraints (characterizing the system macroscopically) is an ensemble. For example, the microcanonical ensemble is the collection of all states with fixed total E, fixed volume V, and number of molecules N. Ensemble averages are determined by

                  (6)

It would be very difficult, and impractical, to estimate the numerator and denominator separately by using the uniform sampling method, i.e., performing T independent measurements and hoping to observe all states and G. However, non-uniform sampling methods can be designed and used to carry out the simulation.

Configurations (microscopic states) are generated (sampled) in such a way that the regions of configurational space that make the largest contributions to the integral <G>_NVT are also regions that are sampled most frequently. IF ?(v) is the probability of choosing a configuration v(r^N, p^N)

                                         (7)

A natural choice of ?(v) is to sample the system on the Boltzmann distribution itself:

                                                                 (8)

Then, the canonical ensemble average is therefore obtained as an unweighted average over configurations sampled:

                                                                 (9)

A trajectory is a chronological sequence of configurations for a system. The sequence may or may not be time dependent. A Monte Carlo trajectory is generated by performing a random walk through configuration space.

Metropolis Monte Carlo Algorithm: (N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, J. Chem. Phys., 21, 1087 (1953)).

In the case of N particles (for example, N Ar atoms) in a box, and E_v = E(r₁, r₂, ..., r_N),

                            I = N*R(0,1) + 1

                            v^new ® v^old: E_old(r₁, r₂, ..., r_I, ..., r_N) ® E_new(r₁, r₂, ..., r_I^new, ..., r_N)

                            If ?E > 0 and exp(-ß ?E) ³ 0, accept the move and new configuration.

otherwise, reject the move.

Let

w_v®v’ = probability per unit time for a transition of the system from state v to state v’.

P_v = the probability that the trajectory resides at a given time in state v.

If we follow a first-order kinetics, the rate of transition is expressed by the master equation:

                         (11)

At equilibrium in the canonical ensemble, _v = 0, and

                                 (12)

Consequently,

                                         (13)

Metropolis Monte Carlo generates a trajectory that samples configurations in accord with the canonical Boltzmann distribution of configurations.

Exercise. Write a simple Monte Carlo (Metropolis) program to compute the magnetization <M> of a square lattice of 20x20 = 400 spins, where M = Sµs_i with µ the magnetic moment, and s_i = ±1, at temperatures J/k_B, 2J/k_B, 4J/k_B, and 8J/k_B. Start with an initial configuration of phase separation (200 spins up and 200 spins down) along with periodic boundary conditions.

The Metropolis Monte Carlo method is designed to sample low energy regions of the system. Thus, it is impractical to compute statistical properties, such as the free energy A and entropy S, which requires the knowledge of the partition function Z itself. On the other hand, the computation of a free energy difference is somewhat less demanding. Consider two fluids characterized by potentials E₁(r) and E₀(r), the free energy difference between the two fluids can be determined by (R. W. Zwanzig, J. Chem. Phys., 22, 1420 (1954)):

                     (14)

where ?E = E₁(r) - E₀(r), and the ensemble < ... >₀ is taken in the reference system E₀(r). Use of equation (14) in computer simulations to compute free energy differences is known as free energy perturbation (FEP) method. An informative reference for estimating the difference in free energies of solvation can be found in J. Chem. Phys.,83, 3050 (1985), by W. L. Jorgensen and C. Ravimohan; a recent review is given by P. A. Kollman, Chem. Rev., 93, 2395 (1993).

A general solution to sample regions of large negative ?E, which can not be efficiently sampled using equation (14), is on the basis of a general non-Boltzmann distribution (Torrie and Valleau):

                         (15)

where W(r) is a positive-valued weighting function. The average of any property in the original ensemble, <Q>₀, can be related to averages taken over Monte Carlo trials in the weighted ensemble:

                                                         (16)

Thus, both averages <Q/W>_W and <1/W>_W are needed to obtain the final result. This method is also termed umbrella sampling, and it has been used in numerous simulation studies. Notably is the computation of conformation distribution of liquid butane and dilute aqueous solution of butane.

The first applications of molecular dynamics for molecular simulations were made to simple fluids of hard spheres by Alder and Wainwright in 1957. Molecular fluid simulations were carried out by Rahman and Stillinger on liquid water, and proteins by McCammon, Glen, and Karplus. In contrast to Monte Carlo procedures, molecular dynamics is a deterministic procedure. Atoms are moved according to the laws of Newtonian mechanics by integrating equations of motion.

                                                     (17)

where r_i(t) is the position of atom I at time t. The forces on atom I is determined by

                                                     (18)

The most widely used numerical integration scheme was that first used by Verlet, which is derived by truncating the Taylor expansion of r_i(t+?t) at r_i(t):

             (19)

The velocities were computer one time-step behind positions:

                                     (20)

An equivalent expression, though perhaps numerically more advantageous, is the velocity Verlet algorithm:

                             (21)

The temperature of the system is then computed from the relation

                                             (22)

It is typically desirable to carry out constant temperature and pressure molecular dynamics simulations. However, these are not as conveniently controlled and precisely defined as in Monte Carlo method. Nose in 1984 described an implementation by treating the dynamics of a system in contact with a thermal reservoir and including an additional degree of freedom to represent the reservoir. Energy is allowed to flow dynamically from and into the reservoir. The technique is like controlling the volume by using a piston. Hoover extended Nose’s method in 1985, and derived a set of equations:

                                                     (23)

                                                             (24)

where Q (unit: energy.time²) is the thermal inertial parameter and f is the number of degrees of freedom. The Nose-Hoover method provides a much gentler control of the temperature than the Berendsen method.

For constant pressure molecular dynamics, Andersen proposed a method in 1980, which involves coupling the system to an external variable, the volume of the simulation box. The coupling mimics the action of a piston on a real system, with a "mass" Q (unit: mass.(length)^-4) associated with a kinetic energy, K = QV’²/2. The potential energy of V is E(V) = PV, where P is the specified pressure.

An excellent textbook on computer simulations is Allen & Tildesley, "Computer Simulation of Liquids", Oxford University Press, New York, 1987.