Barostats in Molecular Dynamics

Molecular Dynamics (MD) is a computational technique that allows us to virtually simulate the motion of a molecule. This technique involves a series of approximations that might result in non-realistic simulations.

To assess the quality of our results we may decide to compare them to a real experiment. If the properties of our simulation get close to experimental data we can be satisfied.

Otherwise, we probably made some errors during the procedure, or we simply were not able to correctly reproduce the properties of the real molecule through this approximation. This might also be due to the low quality of the available starting models of our system.

In the laboratory, it is common practice to perform chemical reactions at constant pressure (P). Reproducing this type of data required computational chemists to develop mathematical tools to simulate a system evolving in time at constant pressure.

We can think of these algorithms as implementing a barostat in our experiment, allowing the instantaneous pressure to fluctuate while the average over a certain period of time remains fixed.

The most widely sampled ensemble in Molecular Dynamics is the NPT ensemble, where the number of atoms (N), the pressure (P), and the temperature (T) is conserved. For this reason, a barostat algorithm is mostly implemented together with a thermostat so that both P and T are conserved along the course of our simulation.

Berendsen barostat

Much of the background, concerning barostat algorithms is analogous to the one needed for thermostat algorithms, which were treated in detail in another article.

In that section, we saw that one possible approach to keep the temperature of the system constant is to scale the velocities of atoms. If they move faster we have higher temperatures and vice versa.

The same idea can be applied to implement a barostat in a given MD simulation. However, to adjust the P we need to scale the volume (V) of the simulation box instead of the velocities. Decreasing the size of the box leads to atoms being more compressed, and therefore the overall pressure of the system increases. We have the opposite situation when we increase the size, and atoms can easily expand in a larger volume.

From a practical point of view, this is achieved by scaling the coordinates of each atom depending on our necessities.

What we saw until now is the main idea behind the Berendsen barostat. As already mentioned, this algorithm is conceptually similar to the analogous thermostat. We conceive our system as weakly coupled to a pressure bath scaling the volume by a certain scaling factor ($\lambda$). At each time step we have that:

$$ V_i^{new} = V_i^{old} \cdot \lambda $$

As mentioned above, in practice we are scaling the positions of atoms instead of the volume. Scaling the volume by $\lambda$ corresponds to scaling the atomic coordinates by $\lambda^{1/3}$. So we write:

$$ r_i^{new} = r_i^{old} \lambda^{1/3} $$

The change of $P$ in time is proportional to the difference in pressure between the bath ($P_{bath}$), and the system ($P(t)$).

$$ \frac{\partial P(t)}{\partial t} = \frac{1}{\tau}(P_{bath}- P(t)) $$

where $\tau$ represents the coupling constant between the two subsystems, and it needs to be increased if we desire a weaker coupling. The explicit expression for the scaling factor in the Berendsen barostat is:

$$ \lambda = \left( 1- \frac{k\delta t}{\tau}(P(t)- P_{bath})\right) $$

where $k$ is a constant specifying the isothermal compressibility of the substance and is larger for easily compressible ones. As we can see, the larger the pressure difference, the larger will be the resulting scaling factor. This means that in cases of high differences the volume is scaled to a larger extent than for smaller differences.

Andersen barostat

Contrary to what we saw for the Berendsen barostat, the Andersen one does not behave as the analogous thermostat.

The Andersen barostat is an extended system method, just like the Nose-Hoover thermostat. Also in this case the system is coupled to a pressure bath, and the coupling is achieved by including an additional degree of freedom in the equations of motion.

This mathematical artifact mimics the action of a piston scaling the volume of the system. The piston has mass $Q$, and has its own kinetic and potential energy. Low mass $Q$ results in rapid box oscillations, and vice versa for large $Q$ values.

The coordinates of the system are thus scaled as follows:

$$ r_i^{new}=r_i^{old} \cdot V^{1/3} $$

Parrinello-Rahman barostat

The idea behind the Andersen barostat was further expanded by Parrinello and Rahman, who developed the analogous barostat.

The main difference with the previous algorithm is that in the present case we are allowed to change the shape of the box, as well as the size. We call this an anisotropic deformation of the box.