Box Size Expansion and Gaps between Periodic Images in NPT Simulation

User discussions
1

GROMACS version: 2021.4
GROMACS modification: Yes/No
Hi everyone,

I’m running a simulation of a surface with a water box on top that contains about 160 sodium ions and 160 chloride ions. I’m using an INDUS simulation to remove the water above the surface by applying a bias potential. However, when I apply the bias, the box appears to expand in the x and y directions, and I’m noticing gaps between periodic images in the trajectory. I’m also using the NPT ensemble for this simulation. My guess is that when I remove the water molecules above the surface, the remaining water accumulates in the smaller region of the system (above the created cavity), which increases the local pressure. This pressure increase causes a significant expansion in the x and y dimensions, leading to gaps between periodic images.My guess is that when I remove the water molecules above the surface, the remaining water accumulates in the smaller region of the system (above the created cavity), which increases the local pressure. This pressure increase causes a significant expansion in the x and y dimensions, leading to gaps between periodic images. Could someone please help me understand how to avoid these gaps? Additionally, I should mention that I simulated my system using AMBER and then used acpype to generate a topology file from the final frame of my AMBER simulation for INDUS simulation in Gromacs. Below, I’ve attached my .mdp file and the output I received:

title = OPLS Lysozyme NPT equilibration
; Run parameters
define = -DPOSRES ; Prevent S and C6 from moving
refcoord_scaling = No ; fix restraint
integrator = md ; leap-frog integrator
nsteps = 4500000 ; 2 * 4500000 = 9000000 ps (9 ns)
dt = 0.002 ; 2 fs
; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save compressed coordinates every 1.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 0 ; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off

2

I just realized that using constant surface tension in my simulation (with these commands:
pcoupltype = surface-tension
ref_p = 5 1.0 )
prevents the surfaces from separating. However, the x and y dimensions of my surface (surface area) still increase significantly during the simulation. I previously used a surface tension of 10 dyne/cm in my AMBER simulations, but I’m not sure how to convert it to bar·nm. Any thoughts?

3

I think 1 \text{dyne/cm} = 10 \text{bar}\times\text{nm}, so 10 \text{dyne/cm} = 100 \text{bar}\times\text{nm}.

More in depth, first change to SI knowing that a dyne is 10^{-5} N:

\begin{equation} \begin{split} \frac{\text{dyne}}{\text{cm}} & = \frac{10^{-5}\times\text{N}}{10^{-2}\times\text{m}} \\ & = 10^{-3}\times\text{m}\times\frac{\text{N}}{\text{m²}} \\ & = 10^{-3}\times\text{m}\times\text{Pa} \\ \end{split} \end{equation}

And then you know that a bar is 10^{5} Pa so:

\begin{equation} \begin{split} 10^{-3}\times\text{m}\times\text{Pa} & = 10^{-3}\times10^{-9}\text{m}\times10^{-5}\text{bar} \\ & = 10^{-3}\times10^{9}\text{nm}\times10^{-5}\text{bar} \\ & = 10\times\text{nm}\times\text{bar} \\ \end{split} \end{equation}
1 Like