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.



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.


The main advantage of the Berendsen barostat is that we can adjust the strength of the coupling between the system and the pressure bath. However, it also shows major disadvantages, as it does not correctly sample from the NPT ensemble.

Therefore, it is generally used in the equilibration phase while it is not good for the production run. If you are not comfortable with these terms check the article on how to set up a Molecular Dynamics simulation.



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} $$



Contrarily to what we saw for the previous barostat, the Andersen barostat does sample from the NPT ensemble. So it can be coupled to a thermostat to perform an NPT production run.

However, this scheme only allows the isotropic deformation of the unit cell. In other words, the box is uniformly modified in all orientations.



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.


The Parrinello-Rahman barostat is generally linked to the Nose-Hoover thermostat to correctly reproduce the correct statistical ensemble. Using these two methods we can carry a production run within the NPT ensemble.

Overall, the Parrinello-Rahman is currently the most employed barostat algorithm in MD. You can set it in a GROMACS simulation by using the Parrinello-Rahman keyword in the mdp file as showed here