Problem while applying pressure coupling on a system containing frozen slab

GROMACS version: VERSION 2024.2-Homebrew
GROMACS modification: Yes/No
I have a simulation box of dimensions 70x70x600 Ang^3, containing, a solid slab (70x70x100) in middle of box (whose atoms are frozen) and there is water (liquid) on both of the outer surfaces. I am performing NPAzT simulations. While applying npt, I want to compress the box from both sides in z direction but the box size is changing in only in one direction (extreme right box boundary is changing), however I did’nt see any changes in the left box boundary. I understand it could be because of shifting of box origin to 0 during pressure scaling however the coordinates of the surface is not changing (because of freeze command that I have applied). After running simulations, the right side water molecules which were initially adsorb onto the slab are getting desorbed and moving away from the surface while the water remains adsorbed on the slab on the left side. I want to keep the surface frozen and want to applying pressure scaling in z direction. What should I do? I want to ensure that the pressure should apply in both direction.
Here are the mdp parameters which I have given:-
;NPT input file
;------------------------------------
integrator = md ;
dt = 0.002 ; in ps
nsteps =5000 ;

;------Names of freeze group
freezegrps = CCT ; freezing Calcite
freezedim = Y Y Y ; freezing in all directions

;------- Continuation
;continuation = no ; by default
tinit = 0 ; starting t value in ps,
;------- Output printing
nstcalcenergy = 100 ; calculating energy freq
nstenergy = 500 ; energy output frequency in edr
nstlog = 500 ; energy output frequency in log
nstxout-compressed = 500 ; trajectory output frequency (output is .xtc file)

nstcomm = 100 ; nstcomm = nstcalcenergy (gmx note)

;------- Generate velocities is off
gen_vel = yes
gen_temp = 300

nstlist = 100 ;default

;---------bond constraints
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all h containing bonds (even heavy atom-H bonds)
lincs_iter = 4 ; accuracy of LINCS (default)
lincs_order = 8 ; also related to accuracy (default)
lincs_warnangle = 90 ;

;------Wall Parameters
;nwall = 2
;wall-atomtype = mmsl_900 mmsl_900
;wall-type = 12-6
;ewald-geometry = 3dc

;------LJ and electrostatic interaction
cutoff-scheme = Verlet ; integrator
verlet-buffer-tolerance = -1

vdwtype = Cut-off ; cut-off by default
rvdw = 1.2 ; vdw cut-off
vdw-modifier = none

coulombtype = PME ; Ewald type
pme_order = 4 ; cubic interpolation (bydefault)
rcoulomb = 1.2 ; cutoff
coulomb-modifier = none

fourierspacing = 0.12 ; grid spacing for FFT (bydefault)
ewald_rtol = 1e-6 ; (bydefault)
rlist = 1.4
;----------dispcorrection
DispCorr = Ener

;----------PBC
pbc = xyz ; pbc only in two direction

;----------------Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = CCT HEP TOL ; multiple coupling groups - more accurate
tau_t = 0.1 0.1 0.1 ; time constant, in ps
ref_t = 300 300 300 ; reference temperature, one for each group, in K

;----------------Presure Coupling
pcoupl = Parrinello-Rahman
pcoupltype = anisotropic

;xx, yy, zz, xy/yx, xz/zx and yz/zy
compressibility = 0.0 0.0 4.5e-11 0.0 0.0 0.0
ref-p = 0.0 0.0 1.0 0.0 0.0 0.0
tau-p = 2.0

With freeze groups such issues can not be avoided completely, but you can translate your system such that the slab is just above z=0. That will minimize the artifacts.

Another solution could be to use position restraints with refcoord-scaling=com.

1 Like