Molecular Dynamics - Pressure and Stress

Introduction
Units and Sign Conventions
How Pressure is Calculated
How Pressure is Controlled

Introduction

Pressure is another basic thermodynamic variable that defines the state of the system. It is a familiar concept, defined as the force per unit area. Standard atmospheric pressure is 1.013 bar, where 1 bar = 10⁵ N m^-2. A single number for the pressure implies that pressure is a scalar quantity, but in fact, pressure is a tensor of the more general form (McQuarrie 1976):

Eq. 5-15:

Each element of the tensor is the force that acts on the surface of an infinitesimal cubic volume that has edges parallel to the x, y, and z axes. The first subscript refers to the direction of the normal to the plane on which the force acts, and the second subscript refers to the direction of the force.

In an isotropic situation, where forces are the same in all directions and there is no viscous force, the pressure tensor is diagonal. With the same diagonal elements, the tensor can be written as:

Eq. 5-16:

where the scalar quantity p is the pressure (or hydrostatic pressure), as usually referred to.

Sometimes, (especially in materials science studies), the pressure tensor is also referred to as the stress tensor or simply the stress. The diagonal elements are known as the tensile stress, and the non-diagonal elements are the shear stress.

Units and Sign Conventions

Pressure and stress can be expressed in many different units. The most common ones are bar and GPa. In SI units, 1 bar = 10⁵ N m ^-2 and 1 GPa = 10⁹ N m ^-2. Hence, 1 GPa = 10⁴ bar. Pressure is usually expressed in bars, but in materials science, stress is often expressed in terms of GPa.

Since the Discover 2.9.5 program is geared towards controlling pressure, the unit used is bar. In the Discover 95.0 program, constant-stress dynamics has GPa as the default unit.

Also note that pressure and stress, although defined in the same physical quantities, have different sign conventions. Positive pressure implies a compressive force pushing the system inward, but positive stress means a force acting outward to pull the system apart.

How Pressure is Calculated

In the Discover program, pressure is calculated through the use of the virial theorem (Allen and Tildesley 1987). Like temperature, pressure is a thermodynamic quantity and is, strictly speaking, meaningful only at equilibrium.

Thermodynamic pressure, thermodynamic temperature, volume, and internal virial can be related in the following way:

Eq. 5-17:

And W is defined as:

Eq. 5-18:

Note that pressure is defined only when the system is placed in a container having a definite volume. In a computer simulation, the unit cell under periodic boundary conditions is viewed as the container. Volume, pressure, and density are calculated only when the system is created with periodic systems.

Analogous to the temperature, an instantaneous pressure function P can be defined so that thermodynamic pressure is the average of the instantaneous values:

Eq. 5-19:

where P is the instantaneous pressure and T is the instantaneous kinetic temperature, which is related to the instantaneous kinetic energy K of the system as:

Eq. 5-20:

The instantaneous pressure function can be written as:

Eq. 5-21:

As mentioned above, pressure is a tensor and its components can also be expressed in tensorial form. Eq. 5-17 can be recast in the form of:

Eq. 5-22:

In detail, the two terms on the right-hand side of the equation are:

Eq. 5-23:
Eq. 5-24:

where r_ix, v_ix, and f_ix indicate the x components of the position, velocity, and force vectors of the i^th atom, respectively. From the definition of the instantaneous pressure tensor, the instantaneous hydrostatic pressure is calculated as 1/3 the trace of the pressure tensor (P_xx + P_yy + P_zz).

When periodic boundary conditions are used, atoms in the unit cell interact not only with the other atoms in the unit cell but also with their translated images. Forces on the images in the virial W must be included correctly, especially for implicit images. If the interaction is pairwise, using Newton's third law, W can be written as:

Eq. 5-25:

instead of:

Eq. 5-26:

Berendsen et al. (1984) uses the r_ij · f_ij formalism by evaluating the virial and the kinetic energy tensor based on the centers of mass, which is valid because the internal contribution to the virial is canceled (on the average) by the internal kinetic energy. Because of the way forces are evaluated, rescaling of coordinates is also done on the basis of the centers of mass of the molecules.

The Discover program, however, calculates pressure on the atomic basis and performs atomic scaling. The major advantage of using atomic scaling of the coordinates is that atom overlapping can be avoided. Such overlap can occur if centers of mass are moved instead of individual atoms. In addition, for large molecules having internal flexibility, atomic scaling has been found to yield a smoother response to pressure changes.

Because the Discover program uses explicit images (also called ghost atoms) in periodic boundary conditions, the r_ij · f_ij formalism is used with ghost atoms (explicit images) included. When the option of using the minimum image (implicit image) is used under periodic boundary conditions, the formalism of r_ij · f_ij is used so that forces between real and image (ghost) atoms can be properly accounted for.

How Pressure is Controlled

Versions 2.9.5 and 95.0 of the Discover program differ in the methods used to control the pressure: version 2.9.5 uses Berendsen's method of pressure-bath coupling, and version 95.0 uses the Parrinello-Rahman method. With the Berendsen method there is no change in the shape of a system, and only the pressure is controlled. With the Parrinello-Rahman method, the system's shape can change, and therefore both pressure and stress can be controlled.

Berendsen Method of Pressure Control

Pressure changes can be accomplished by changing the coordinates of the particles and the size of the unit cell in periodic boundary conditions.

During a molecular dynamics simulation, the Discover program uses Berendsen's method (Berendsen et al. 1984), which couples the system to a pressure ``bath'' to maintain the pressure at a certain target. The strength of coupling is determined both by the compressibility of the system (using a user-defined variable ) and by a relaxation time constant (a user-defined variable ). At each step, the x, y, and z coordinates of each atom are scaled by the factor:

Eq. 5-27:

where

t is the time step, P is the instantaneous pressure, and P₀ is the target pressure. The Cartesian components of the unit cell vectors are scaled by the same factor

Note that this method changes the cell uniformly, so that the size of the cell is changed, but not its shape. In cases such as phase transitions of crystals, where both the cell size and shape are expected to change, this method is not appropriate.

The pressure fluctuation has been observed to be large during test runs with the Discover program and in other studies in the literature (Brown and Clark 1984). Negative pressures sometimes occur because the virial can be negative, even though this defies the usual sense that pressure is a positive number. The calculated pressure depends on the cutoff distances used in the simulation. To compensate for the missing long-range part of the potential contributions to the energy, pressure at r > r_c can be estimated by assuming the radial distribution g(r) ~ 1 (uniform distribution) in that region. With this assumption and the known form of nonbond potentials, the correction can be estimated analytically (Allen and Tildesley 1987). However, this correction is not built into the calculation. If the cutoff distance is too short, the calculated pressure may be wrong. Therefore in practice, you should test the effect of cutoff on pressure by gradually increasing the cutoff and choosing the appropriate cutoff accordingly.

Andersen Method of Pressure Control

In the Discover 9.50 program, the Andersen method (1980) is used to control the pressure instead of the Berendsen method (above). The volume of the cell can change, but its shape is preserved by allowing the cell to change isotropically. The Anderson method is useful for liquid simulations since the box could become quite elongated in the absence of restoring forces if the shape of the cell were allowed to change. A constant shape also makes the dynamics analysis easier. However, this method is not very useful for studying materials under nonisotropic stress or phase transitions, which involve changes in both cell lengths and cell angles (for these conditions, the Parrinello-Rahman method, below, should be used).

The basic idea of the method is to treat the volume V of the cell as a dynamic variable in the system. The Lagrangian of the system is modified so that it contains a term in the kinetic energy with a user-defined mass M and a potential term which is the pV potential derived from an external pressure P_ext acting on volume V of the system.

Parrinello-Rahman Method of Pressure and Stress Control

The Parrinello-Rahman method of pressure and stress control can allow simulation of a molecular system under externally applied stress. This is useful for studying the stress-strain relationship of materials. Both the shape and the volume of the cell can change, so that the internal stress of the system can match the externally applied stress.

The method was presented in detail by Parrinello and Rahman (1981) and is only summarized here. The Lagrangians of the system are modified such that a term representing the kinetic energy of the cell depends on a user-defined mass W. An elastic energy term p is related to the pressure and volume or the stress and strain of the system. The equations of motion for the atoms and cell vectors can be derived from this Lagrangian. The motion of the cell vectors, which determines the cell shape and size, is driven by the difference between the target and internal stress.

With hydrostatic pressure only:

Eq. 5-28:

where h = {a · b · c} is the cell vector matrix, G = h'h, r_i = hs_i,

is the interaction potential, and

is the volume of the cell. The dots above some symbols indicate time derivatives and the primes indicate matrix transposition. Tr is the trace of a matrix.

With stress, the elastic term p is replaced by:

Eq. 5-29:

where S is the applied stress,

₀ is the initial volume, and

is the strain;

is also a tensor, defined as (h₀'^-1 Gh₀^-1 - 1)/2.

Note the user-defined variable W. A large W means a heavy, slow cell. In the limiting case, infinite W reverts to constant-volume dynamics. A small W means fast motion of the cell vectors. Although that may mean that the target stress can be reached faster, there may not be enough time for equilibration. In tests, too small a W also induced artificial periodic motions of the cell. A value of 20 seems to give satisfactory results for a cell with 76mers of polyethylene, and hence is used as the default value for W.

Main access page

Theory/Methodology access.

Dynamics access Dynamics - Temperature Dynamics - Constraints