For its first 30 years, molecular dynamics was limited to the microcanonical ensemble. We have already seen how the development of thermostats has enabled simulations to be performed in the isochoric, canonical and isokinetic ensembles. We will now describe molecular dynamics algorithms for performing simulations at constant pressure or constant enthalpy. The technique used to make the pressure rather than the volume, the independent state defining variable, uses essentially the same ideas as those employed in §6.6 to design Norton ensemble algorithms. The methods we describe now are of use for both equilibrium and nonequilibrium simulations.
It is often advantageous, particularly in studies of phase transitions, to work within the isobaric ensemble. It is possible to stabilise the pressure in a number of ways: we will describe the Gaussian method (Evans and Morriss, 1983b) since it was both the first deterministic isobaric technique to be developed and it is conceptually simpler than the corresponding Nosé-Hoover (Hoover, 1985) and Rahman-Parrinello (1980a,b, 1981) schemes. Although it may be slightly more difficult to write the computer programs, once written they are certainly easier to use. The Gaussian method has the distinct advantage that the pressure is a rigorous constant of the motion whereas the Nosé based schemes (Nosé, 1984) and those of Parrinello and Rahman allow fluctuations in both the pressure and the volume.
If one makes a poor initial guess for the density, Nosé-Hoover isobaric algorithms induce sharp density changes in an attempt to correct the density, to that appropriate for the specified mean pressure. Because bulk oscillations damp quite slowly, Nosé-Hoover methods can easily result in the system exploding - a situation that cannot be reversed due to the finite range of the interaction potentials. Gaussian isobaric algorithms are free of these instabilities.
Consider the SLLOD equations of motion where the strain rate tensor ∇u is isotropic . The equations of motion become
(6.83)
(6.84)
Now if the system was cold (pi=0 for all i), and non-interacting (Φij=0, these equations would reduce to
(6.85)
Since this equation is true for all particles i, it describes a uniform dilation or contraction of the system. This dilation or contraction is the same in each coordinate direction, so if the system initially occupied a cube of volume V, then the volume would satisfy the following equation of motion.
(6.86)
For warm, interacting systems, the equation of motion for qi shows that the canonical momentum pi is in fact peculiar with respect to the streaming velocity 
. The dissipation for the system (6.83 & 6.84) is
(6.87)
Since H0 is the internal energy of the system we can combine (6.87) with the equation of motion for the volume to obtain the first law of thermodynamics for adiabatic compression,
(6.88)
It is worth noting that these equations are true instantaneously. One does not need to employ any ensemble averaging
to obtain equation (6.88). By choosing the dilation rate 
to be a sinusoidal function of time, these equations of motion
can be used to calculate the bulk viscosity. Our purposes are however to use the
dilation rate as a multiplier to maintain the system at a constant hydrostatic
pressure. Before we do this however, we will introduce a Gaussian thermostat into
the equations of motion;
(6.89)
(6.90)
The form for the thermostat multiplier is determined by the fact that the momenta in (6.89 & 6.90) are peculiar with respect to the dilating coordinate frame. By taking the moment of (6.90) with respect to pi, and setting the time derivative of the peculiar kinetic energy to zero we observe that,
(6.91)
Differentiating the product pV, (6.87) with respect to time gives,
(6.92)
The first term on the LHS is zero because the pressure is constant, and the first
term on the RHS is zero because the peculiar kinetic energy is constant.
Substituting the equations of motion for 
and 
, and we can solve for the dilation rate.
(6.93)
Combining this equation with (6.91) gives a closed expression for the thermostat multiplier α.
In summary our isothermal/isobaric molecular dynamics algorithm involves solving 6N+1 first order equations of motion (equations (6.86, 6.89 & 6.90)). There are two subtleties to be aware of before implementing this method. Firstly the pressure is sensitive to the long range tail of the interaction potential. In order to obtain good pressure stability the long range truncation of the potential needs to be handled carefully. Secondly, if a Gear predictor corrector scheme is used to integrate the equations of motion then some care must be taken in handling the higher order derivatives of the coordinates and momenta under periodic boundary conditions. More details are given in Evans and Morriss (1983b) and (1984a).
For the adiabatic constant pressure equations of motion we have already shown that the first law of thermodynamics for compression is satisfied
(6.94)
It is now easy to construct equations of motion for which the enthalpy I=H0+pV, is a constant of the motion. The constraint we wish to impose is that
(6.95)
Combining these two equations we see that for our adiabatic constant pressure equations of motion the rate of change of enthalpy is simply
(6.96)
This equation says that if our adiabatic equations preserve the pressure then the enthalpy is automatically constant. The isobaric-isoenthalpic equations of motion are simply obtained from the isothermal-isobaric equations by dropping the constant temperature constraint. The isoenthalpic dilation rate can be shown to be (Evans and Morriss, 1984a),
(6.97)