Surface tension calculation

GROMACS version: 2024.3
GROMACS modification: No

Hi everyone,

I am trying to compute the surface tension of an aqueous solution. I built the simulation box in CHARMM-GUI and got the .gro file - a relatively small system, just to get started. I then modified the .gro file dimensions, doubling the size of the z-axis (from 3.0 to 6.0). I don’t know if this is the correct approach. EM and NVT equilibration run smoothly; size of the box remains the same. However, during NPT equilibration, sizes become: 1.27914 1.27914 13.66222. So for the production run I get the (understandable) error:

ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist.

Am I doing it wrong? Should I increase the size of the box?

My NPT equilibration .mdp file is:

; Run control
integrator = md
tinit = 0
dt = 0.001
nsteps = 20000000 ; 20 ns
nstcomm = 100

; Output control
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 20000
nstenergy = 20000
nstxout-compressed = 20000

; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.2

; Electrostatics
coulombtype = PME
rcoulomb = 1.2

; van der Waals
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2

; Apply long range dispersion corrections for Energy and Pressure
DispCorr = no

; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12

; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
epsilon_surface = 0

; Temperature coupling
tcoupl = v-rescale
tc_grps = system
tau_t = 1.0
ref_t = 298

; Pressure coupling is on for NPT
pcoupl = C-rescale
pcoupltype = surface-tension
tau_p = 1.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
refcoord_scaling = com

; Generate velocities
gen_vel = yes
gen_temp = 298
gen_seed = -1

; options for bonds
constraints = h-bonds

; Type of constraint algorithm
constraint-algorithm = lincs

; Constrain the starting configuration
; since we are continuing from NVT
continuation = no

; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4

Thank you in advance!