Molecular Dynamics - Temperature

Contents

Introduction

Temperature is a state variable that specifies the thermodynamic state of the system and is also an important concept in dynamics simulations. This macroscopic quantity is related to the microscopic description of molecular simulations through the kinetic energy, which is calculated from the atomic velocities. The temperature and the distribution of atomic velocities in a system are related through the Maxwell-Boltzmann equation:

Eq. 5-3:
This well known formula expresses the probability f(v) that a molecule of mass m has a velocity of v when it is at temperature T. Figure 5-6 plots this distribution at various temperatures. The x, y, z component of the velocities, on the other hand, has a Gaussian distribution:

Eq. 5-4:
In the Discover program, the initial velocities are generated from the Gaussian distribution of vx, vy, and vz. The Gaussian distribution is generated from the random number generator and a random number seed which can be defined by the user or generated automatically.

How Temperature Is Calculated

Temperature is a thermodynamic quantity, which is only meaningful at equilibrium. It can be related to the average kinetic energy of the system through the equipartition principle. This principle states that every degree of freedom (either in momenta or in coordinates), which appears as a squared term in the Hamiltonian, has an average energy of kT/2 associated with it. This is the case for momenta pi which appear as pi2/2m in the Hamiltonian. Hence we have:

Eq. 5-5:
The left side of Eq. 5-5 is also called the average kinetic energy of the system, Nf is the number of degrees of freedom, and T is the thermodynamic temperature. In an unrestricted system with N atoms, Nf is 3N because each atom has three velocity components, i.e., vx, vy, and vz.

Because of this, it is convenient to define an instantaneous kinetic temperature function whose average is the thermodynamic temperature T:

Eq. 5-6:
The average of the instantaneous temperature is the thermodynamic temperature.

Temperature is calculated from the total kinetic energy and the total number of degrees of freedom. Therefore, for a nonperiodic system:

Eq. 5-7:
Six degrees of freedom are subtracted because both the translation and rotation of the center of mass are ignored.

And for a periodic system:

Eq. 5-8:
Only the three degrees of freedom corresponding to translational motion can be ignored, since rotation of a central cell imposes a torque on its neighboring cells.

How Temperature Is Controlled

Although the initial velocities are generated so as to maintain a Maxwell-Boltzmann distribution at the desired temperature, the distribution does not remain constant as the simulation continues. This is especially true when the system does not start at an equilibrium configuration. This often occurs, since the system is minimized to eliminate the hot spots. During dynamics, kinetic energy is changed to potential energy as the minimized structure changes to the equilibrium structure, and the temperature also changes.

To maintain the correct temperature, the computed velocities have to be adjusted appropriately. Besides getting the temperature to the right target, the temperature-control mechanism also must produce the correct statistical ensembles. This means that the probability of occurrence of a certain configuration obeys the laws of statistical mechanics. For constant-temperature, constant-volume dynamics to generate the canonical ensemble, this means that P(E), the probability that a configuration with energy E will occur, is proportional to exp(-E/kBT), also called the Boltzmann factor.

Several methods for temperature control are used in the Discover program. In version 2.9.5, these are:

In version 95.0, the temperature-control methods from which you can choose (during either stage of the simulation) are:

Direct Velocity Scaling

Direct velocity scaling is a drastic way to change the velocities of the atoms so that the target temperature can be exactly matched whenever the system temperature is higher or lower than the target by some user-defined amount. The velocities of all atoms are scaled uniformly as follows:

Eq. 5-9:
This adds (or subtracts) energy from the system efficiently, but it is important to recognize that the fundamental limitation to achieving equilibrium is how rapidly energy can be transferred to, from, and among the various internal degrees of freedom of the molecule. The speed of this process depends on the potential energy expression, the parameters, and the nature of the coupling between the vibrational, rotational, and translational modes. It also depends directly on the size of the system, larger systems taking longer to equilibrate.

Berendsen's Method of Temperature-Bath Coupling

After equilibration, a more gentle exchange of thermal energy between the system and a heat bath. can be introduced through the Berendsen et al. (1984) method, in which each velocity is multiplied by a factor given by:

Eq. 5-10:
where t is the timestep size, is a characteristic relaxation time, T0 is the target temperature, and T the instantaneous temperature. To a good approximation, this treatment gives a constant-temperature ensemble that can be controlled, both by adjusting the target temperature T0 and by changing the relaxation time, which is equivalent to in units of picoseconds.

Nosé Dynamics

Nosé dynamics (Nosé 1984a, 1984b, 1991) is a relatively new method for performing constant-temperature dynamics, which produces true canonical ensembles both in coordinate space and in momentum space. The actual formalism implemented is based on a simplified reformulation by Hoover (1985). Thus, the method is also called the Nosé-Hoover thermostat (or Nosé-Hoover dynamics).

The main idea behind Nosé-Hoover dynamics is that an additional (fictitious) degree of freedom is added to the real physical system. This fictitious mass is given a mass Q. The equations of motion for the extended (i.e., real plus fictitious) system are solved. If the potential chosen for that degree of freedom is correct, the constant-energy dynamics (or the microcanonical dynamics, NVE) of the extended system produces the canonical ensemble (NVT) of the real physical system.

The Hamiltonian H* of the extended system is:

Eq. 5-11:
Equations of motion for the real-atom coordinates q and moments p, as well as for the fictitious coordinates S and momentum (where is the interaction potential) are:

Eq. 5-12:
Eq. 5-13:
Eq. 5-14:

where Q = the user-defined q_ratio x a constant x g x T;
g = number of degrees of freedom;
T = temperature.

The choice of the fictitious mass Q of that additional degree of freedom is arbitrary but is critical to the success of a run. If Q is too small, the frequency of the harmonic motion of the extended degree of frequency is too high. This forces a smaller timestep to be used in integration. However, if Q is too large, the thermalization process is not efficient--as Q approaches infinity, there is no energy exchange between the heat bath and the real system. (The Nosé method also requires an accurate integrator, for energy to be conserved. Thus, the default integrator for Nosé dynamics is ABM4 with 0.5 fs as the timestep.)

Q should be different for different systems--Nosé (1991) suggests that Q should be proportional to gkBT, where g is the number of degrees of freedom in the system, kB is the Boltzmann constant, and T is the temperature. To determine the proportionality constant, studies with a box of liquid argon containing 343 atoms at 87 K were done. Setting Q to 2.5 10-5 kcal mol-1 fs-2 for this system yielded good results. This proportionality constant, together with the gkBT term, was then used to generate Q for a box of water and an amorphous cell of polypropylene, which also yielded satisfactory results.

Tests on polypropylene indicate that Nosé dynamics needs a timestep of 0.5 fs with the velocity Verlet method in order to approach within 3% of the target temperature of 298 K. In contrast, 1.0 fs is sufficient for the direct velocity scaling method. For the Nosé method, the smaller the timestep, the closer it approaches the target temperature. Apparently, the need for a smaller timestep to achieve accuracy is unrelated to Q, as long as Q is in the appropriate range (because the error in reaching the target temperature remains the same when Q is increased or decreased by a factor of 2).

Andersen Method

There are two versions of the Andersen method of temperature control. One version is to randomize the velocities of all atoms at a predefined frequency. The other version is to pick one atom at each time step and change it according to the Boltzmann distribution. The Discover 95.0 program implements the first version. The predefined frequency is proportional to N2/3, where N is the number of atoms in the system. Although this frequency is calculated by the program, you can change it.

Main access page Theory/Methodology access.

Dynamics access Dynamics - Statistical Ensembles Dynamics - Pressure & Stress

Copyright Biosym/MSI