Fast NPT equilibration

User discussions
1

GROMACS version: 2021.4
GROMACS modification: Yes/No
Here post your question: I’ve created a simulation box including 1128 EC solvents, 37 magnesiums, 74 TFSI anions and 16 12-crown-4 in a 5*5*5 nm³ simulation box applying OPLSAA FF. However, after 10 ns of NPT equilibration, the GRO/coordinate file reveals that my simulation box expanded from 5 nm to 5.3 nm in three directions. After generating the density graph, the initial density is >1600 kg/m³, the average density is 1400 kg/m³, and the density of my EC solvent is 1320 kg/m³. However, I started my pressure at 1 atm, but after equilibration, the pressure ended at 4 bar. I think my equilibration is too fast, and I really don’t know how I can control the pressure of my system. I also want to know if the expansion of my box is fine. Should I proceed with annealing based on the size of this box, or do I need to make it 5 nm in all three directions, as I initially specified for the box size? Can anyone please help me know where I’m going wrong? I’m new to this. npt.mdp is as follows:

; Run parameters

integrator = md
nsteps = 10000000 ;
dt = 0.001 ; 1 fs

; Output control
nstxout-compressed = 10000 ; Save .xtc every 10 ps
nstxout = 1000 ; Save coordinates every 1 ps
nstvout = 1000 ; Save velocities every 1 ps
nstenergy = 1000 ; Save energies every 1 ps
nstlog = 1000 ; Update log every 1 ps

; Bond parameters
constraint_algorithm = lincs
constraints = none
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

; Neighbor searching and vdW
cutoff-scheme = Verlet
nstlist = 10 ; largely irrelevant with Verlet
ns_type = grid
rlist = 1.2
rvdw = 1 ; short-range van der Waals cutoff (in nm)
vdw-type = Cut-off

; Electrostatics
coulombtype = pme
rcoulomb = 1.2

; Ewald
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
ewald_rtol = 1e-05
epsilon_surface = 0
pme-order = 4

; Temperature coupling
tcoupl = v-rescale
tc-grps = System
tau_t = 1
ref_t = 303K

; Pressure coupling
pcoupl = Berendsen
pcoupltype = isotropic
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 ; 1 atm
compressibility = 4.5e-5

; Periodic boundary conditions
pbc = xyz

; Dispersion correction
DispCorr = EnerPres

; Velocity generation
gen_vel = no

2

Hi @penelope,

Pressure coupling in molecular simulations works by changing the size of the simulation box, so the change you are seeing is entirely expected. Also, the instantaneous pressure is expected to fluctuate by quite a bit, so it’s important to check the average values rather than the ones the simulation ended on.

Looking at your mdp settings, you’re using the Berendsen barostat, which is outdated and should not be used. I would try c-rescale instead. Also, 4.5e-5 bar^-1 is the compressibility of water, so if you’re using a different solvent you should set the appropriate value for the compressibility.

1 Like
3

Hello Imullender,

Thank you very much for the reply.