The previous step, NVT equilibration, stabilized the temperature of the system. Prior to data collection, we must also stabilize the pressure (and thus also the density) of the system. Equilibration of pressure is conducted under an NPT ensemble, wherein the Number of particles, Pressure, and Temperature are all constant. The ensemble is also called the "isothermal-isobaric" ensemble, and most closely resembles experimental conditions.

The .mdp file used for a 500-ps NPT equilibration can be found here. It is not drastically different from the parameter file used for NVT equilibration, though this phase of equilibration usually takes somewhat longer than NVT equilibration because pressure takes a little longer to relax than temperature. Note the addition of the pressure coupling section, using the C-rescale barostat.

A few other changes:

• continuation = yes: We are continuing the simulation from the NVT equilibration phase
• gen_vel = no: Velocities are read from the trajectory (see below)

We will call grompp and mdrun just as we did for NVT equilibration. Note that we are now including the -t flag to include the checkpoint file from the NVT equilibration; this file contains all the necessary state variables to continue our simulation. To conserve the velocities produced during NVT, we must include this file. The coordinate file (-c) is the final output of the NVT simulation.

gmx grompp -f inputs/npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

gmx mdrun -deffnm npt

Let's analyze the pressure progression, again using energy:

gmx energy -f npt.edr -o pressure.xvg

Type "17 0" at the prompt to select the pressure of the system and exit. The resulting plot should look something like the following:

The pressure value fluctuates widely over the course of the 100-ps equilibration phase, but this behavior is not unexpected. The running average of these data are plotted as the red line in the plot. Over the course of the equilibration, the average value of the pressure is -3 ± 11 bar. Note that the reference pressure was set to 1 bar, so is this outcome acceptable? Pressure is a quantity that fluctuates widely over the course of an MD simulation, as is clear from the error bar (11 bar), so statistically speaking, one cannot distinguish a difference between the obtained average (-3 ± 11 bar) and the target/reference value (1 bar).

Let's take a look at density as well, this time using energy and entering "23 0" at the prompt.

gmx energy -f npt.edr -o density.xvg

As with the pressure, the running average of the density is also plotted in red. The average value over the course of 500 ps is 1025.3 ± 0.5 kg m^-3, higher than the experimental value of 1000 kg m^-3 and the expected density of the TIP3P model of 986 kg m^-3. Is this result reasonable? We have not simulated pure water; instead, there is a protein in our system and some diffuse ions. Therefore, we should not expect an exact match, and during equilibration, we need only to see that the density has stabilized.

Please note: I frequently get questions about why density values obtained do not match my results. Pressure-related terms are slow to converge, and there are elements of randomness in any MD simulation. You will almost certainly not reproduce this outcome exactly.