Longer NPT equilibration

User discussions
1

GROMACS version: 2021.4
GROMACS modification: Yes/No
Here post your question : I prepared an equilibrated box of cyclohexane. But the desired density was achieved after an NPT equilibration of 20 ns. Is this okay or advisable for so long NPT equilibration? Or am I wrong at some initial step?

npt.mdp is as follows-
title = charmm36 CYHE system NPT equilibration

; Run parameters

integrator = md ; leap-frog integrator

nsteps = 10000000 ; 2 * 10000000 = 20000 ps

dt = 0.002 ; 2 fs

; Output control

nstenergy = 500 ; save energies every 1.0 ps

nstlog = 500 ; update log file every 1.0 ps

nstxout-compressed = 500 ; save coordinates every 1.0 ps

; Bond parameters

continuation = yes ; continuing from NVT

constraint_algorithm = lincs ; holonomic constraints

constraints = h-bonds ; bonds to H are constrained

lincs_iter = 1 ; accuracy of LINCS

lincs_order = 4 ; also related to accuracy

; Neighbor searching and vdW

cutoff-scheme = Verlet

ns_type = grid ; search neighboring grid cells

nstlist = 20 ; largely irrelevant with Verlet

rlist = 1.2

vdwtype = cutoff

vdw-modifier = force-switch

rvdw-switch = 1.0

rvdw = 1.2 ; short-range van der Waals cutoff (in nm)

; Electrostatics

coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics

rcoulomb = 1.2

pme_order = 4 ; cubic interpolation

fourierspacing = 0.16 ; grid spacing for FFT

; Temperature coupling

tcoupl = V-rescale ; modified Berendsen thermostat

tc-grps = system ; two coupling groups - more accurate

tau_t = 1.0 ; time constant, in ps

ref_t = 300 ; reference temperature, one for each group, in K

; Pressure coupling

pcoupl = C-rescale ; pressure coupling is on for NPT

pcoupltype = isotropic ; uniform scaling of box vectors

tau_p = 5.0 ; time constant, in ps

ref_p = 1.0 ; reference pressure, in bar

compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1

refcoord_scaling = com

; Periodic boundary conditions

pbc = xyz ; 3-D PBC

; Dispersion correction is not used for proteins with the C36 additive FF

DispCorr = no

; Velocity generation

gen_vel = no ; velocity generation off after NVT

2

Dear @SP06 ,

What is the problem? You think the equilibration is too long? Why?

I would say that, in principle, an equilibration cannot be too long, but only too short. Clearly, if the equilibration is too long this means that you are ‘wasting’ computational time, but at least you are equilibrated, so it’s better to have a somewhat ‘too’ long equilibration than a ‘too’ short one.

1 Like
3

While preparing an equilibrated box of cyclohexane, after NVT equilibration I got the required temperature. The first NPT equilibration was for 6 ns, for which I got the density plot as follows-

IMG_20240603_082915
IMG_20240603_0829151920×2550 695 KB

As the average density was not close to the density of cyclohexane, so I increased the NPT equilibration, so finally at 20 ns, I got the desired density. But the average Pressure went down to -9 bar. SO I wanted to know, for NPT equilibration, should I always go on increasing the time till density of cyclohexane is achieved?

4

It seems strange that the density plot would flatten like that unless you are in equilibrium. Or are you using restraints? Restraints may be required for early equilibration stages (often just a single NVT equilibration), but after that they may do more harm than good.

5

Sir
I wanted to know, for a protein-water, protein-cyclohexane or just any system involving a protein, how much approx. temperature deviation from reference temperature (after NVT) and pressure deviation from reference pressure (after NPT) is permissible ? Is there any such permissible deviation range?

6

It is very difficult to give a permissible deviation range. It will vary from system to system, but when the system is equilibrated you want the temperature and pressure to fluctuate (how large the fluctuations are will depend on your system and thermostat/barostat settings) with an average close (how close may depend on your system and thermostat/barostat settings) to the specified temperature and pressure.

When deciding if your system is equilibrated, I would also recommend checking that the potential energy has reached a plateau.

7

Ok. I understand. I wanted to know even if I achieved the desired density but Pressure is far away from reference Pressure of 1 bar(after NPT equilibration), should I continue extending the NPT equilibration till reference pressure reached? For my system, after NPT I got -2 bar pressure while my reference is 1 bar. I already equilibrated it for 40 ns.

8

How have you measured the pressure (when getting -2 bar)? Is that what you observe when it has reached a plateau? If you use, e.g., gmx energy with the -b option to start from the beginning of the plateau, what is the uncertainty of the average?

By the way, above you said that the average density after 6 ns was not close to the desired value. The value I found was 0.77 g/cm3. The plot you showed is very close to that from 3 ns and onwards. Make sure that you start measuring from your stable state when you calculate an average during equilibration.

Some systems may need to be equilibrated for more than 100 ns, but it’s best to first make sure that you really need to extend it that much. I don’t think a box of cyclohexane should take very long to equilibrate.

9

Yes I used gmx energy to calculate average pressure and density. After 6 ns, the average density by gmx energy (gmx energy -f npt.edr -o density.xvg) showed 716 kg/m3, but average pressure (gmx energy -f npt.edr -o pressure.xvg) was -26 bar. So I extended it to 20 ns, after which average density and pressure (by gmx energy) shows as 756 kg/m3 and -2 bar. So, I wanted to know if I should extend the equilibration ? Or should I leave equilibration at that if in the plot it reaches a plateau (as in the attached plot), without caring about the average density and pressure calculation, which was calculated on the total npt.edr and not only for the plateau region?

With that, can you please explain me this part of pressure from the lysozyme md tutorial-
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 7.5 ± 160.5 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 large root-mean-square fluctuation (160.5 bar), so statistically speaking, one cannot distinguish a difference between the obtained average (7.5 ± 160.5 bar) and the target/reference value (1 bar).

According to this is it okay to have larger fluctuations in pressure upto 7-8 bar, even when reference pressure is set to 1 bar?

10

As I said before, do not calculate the average values over the whole equilibration simulation. The values before reaching equilibrium will affect the average value. What you are interested in is the average value after reaching equilibrium, which seems to be after 3 ns according to the plot above. Try gmx energy -b 3000 and see what values you get.

11

Should I check for pressure too from 3 ns , where density reaches a plateau?

12

@SP06 How did you parameterise the cyclohexane? and what is the compressibility used?

13

Yes, you should make sure that it reaches a plateau and measure its value from that point and onwards. Otherwise it can be heavily influenced by the starting conditions.

14

As pressure is fluctuating constantly, so plateau is difficult to observe, but from density plot where density reaches a plateau (3 ns) from there I used gmx energy with -b 3000 for pressure calculation, I got average pressure to be -0.47 bar. Is it acceptable w.r.t to reference pressure of 1, or should I equilibrate more?

15

What is the error bar? Pressure fluctuates like crazy during a simulation, on the order of sqrt(N).

16

Average pressure = -0.472165
Err.Est. = 0.51
RMSD = 106.715
Tot-Drift = 1.85977
(in bar)

17

-0.47 +/- 0.51 might be a bit too far from the target. The drift is quite large, so the pressure has probably not reached it’s equilibrium when you start measuring. So, that is from 3-6 ns? Or is that 3-40 ns? If it is 3-6 ns, I’d extend it. If it’s 3-40 ns, I’d measure in the range 20-40 ns instead.

18

I got this pressure in 3-20 ns. Okay I will check for 20-40 ns range to see how the average pressure is.

With this, can you please explain me the following from lysozyme md tutorial -

" 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 7.5 ± 160.5 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 large root-mean-square fluctuation (160.5 bar), so statistically speaking, one cannot distinguish a difference between the obtained average (7.5 ± 160.5 bar) and the target/reference value (1 bar)."

Generally what should be the err.est., RMSD and tot-drift during average pressure and density calculation. Is it okay for the RMSD for avg. density (~600) to be so high than in avg.pressure?

19

I’m sorry, I thought you said above that you had extended it to 40 ns already. That’s why I suggested to analyse 20-40 ns.

If you don’t have data after 20 ns, you can check 10-20 ns or 15-20 ns instead. Just select a period that is stable and long enough for you to be sure that the simulation is in equilbrium. If there is no such period you might need to extend it - to 30 or 40 ns.

You would not expect a large drift, of e.g. temperature, pressure, density and potential energy, when the simulation is in equilibrium. There are no hard rules, but if the average pressure is not according to the target (within 2 err. est.) and the drift is indicating that it is approaching the target, then you should probably extend (or start measuring the observable at a later time point). Pressure RMSD will usually be quite high, and it depends on the size of your system - large systems usually have smaller pressure fluctuations. However, that does not mean that it is better to run large systems.

20

Now I have extended it to 50 ns, from 20-50 ns, the average pressure shows to be 1.9 bar, but as I change the range from 30-50 ns, 40-50 ns, the pressure changes to -3 and -4 bar respectively. Actually I had the idea that a box of cyclohexane (5 * 5 * 5), won’t take so long for equilibrating (as in biphasic tutorial it took 10 ns for equilibration). But seeing the average pressure, it seems, I need to extend the npt equilibration again. As I am new to MD simulations, I am worried if this is normal to have as long as more than 50 ns to equilibrate a box of cyclohexane.