I’m trying to simulate a confined system with two walls and have a question about whether NPT equilibration prior to NVT production is necessary.
Here’s my system setup: Two CaCO3 slabs are positioned at the bottom and top (z=0, z=L), with their orientations flipped; i.e., the top slab is generated by flipping the bottom slab. Between the two slabs, there are water molecules, ions, and polymers.
Initially, I attempted both NVT and NPT equilibrations. However, the NPT equilibration failed because the box’s shrinking speed was faster than the top CaCO3 slab could adjust. I have noticed that many slab simulations skip NPT equilibration. Is it truly necessary in this case?
No I didn’t restrain the slabs. I’m using the C-rescale with semiisotropic pressure coupling.
This is mdp setting.
integrator = md ; “steep” for energy minimization and “md” for a leap frog algorithm for integrating Newton’s eq of motion
dt = 0.001 ; 1 femtosecond
nsteps = 5000000 ; 5 nanoseconds
comm-mode = Linear ; “Linear” (remove center for mass translation) “Angular” ( >> >> and rotation around the center of mass) “None” (no restric)
nstcomm = 100 ; [steps] frequency for center of mass motion removal
comm-grps = System ; groups for which the center of mass motion is removed
;;; Output control ;;;
nstxout = 0 ; [steps] number of steps that elapse between writing coordinates to the output trajectory file (trr)
nstvout = 0 ; [steps] number of steps that elapse between writing velocities to the output trajectory file (trr)
nstfout = 0 ; [steps] number of steps that elapse between writing forces to the output trajectory file (trr)
nstcalcenergy = 100 ; [steps] number of steps that elapse between calculating the energies, 0 is never
nstenergy = 1000 ; [steps] number of steps that elapse between writing energies to energy file
nstlog = 5000 ; [steps] number of steps that elapse between writing energies to the log file
nstxout-compressed = 5000 ; [steps] number of steps that elapse between writing to .xtc file
;;; Neighbor searching ;;;
cutoff-scheme = Verlet ; algorithm to generate neigbor list
nstlist = 10 ; [steps] frequency to update the neighbor list
pbc = xy
rlist = 1.2000 ; [nm] cut-off distance for the short-range neighbor list
;;; Electrostatics ;;;
coulombtype = PME ; algorithm to generate electrostatics
rcoulomb = 1.2000 ; [nm] The distance for the Coulomb cut-off
;;; Van der Waals ;;;
vdw-type = cut-off ; algorithm to generate Van der Waals
rvdw = 1.2000 ; [nm] distance for the LJ or Buckingham cut-off
DispCorr = AllEnerPres ; corrections to apply for long-ranged energy and/or pressure
;;; Ewald ;;;
fourierspacing = 0.12 ; [nm] for ordinary Ewald, the ratio of the box dimensions and the spacing determines a lower bound for the number of wave vectors to use in each (signed) direction
pme-order = 4 ; interpolation order for PME. 4 equals cubic interpolation.
ewald-rtol = 1.0e-5 ; the relative strength of the Ewald-shifted direct potential at rcoulomb
ewald-geometry = 3dc
;;; Temperature coupling ;;;
Tcoupl = nose-hoover ; Temp coupling using velocity rescaling. Temperature coupling using a Nose-Hoover extended ensemble. The reference temperature and coupling groups are selected as above, but in this case tau-t controls the period of the temperature fluctuations at equilibrium, which is slightly different from a relaxation time. For NVT simulations the conserved energy quantity is written to the energy and log files.
tc-grps = System ; groups to couple separately to Temp bath
tau-t = 0.50 ; [ps] time constant for coupling (one for each group in tc-grps)
ref-t = 300
;;; Pressure Coupling ;;;
Pcoupl = C-rescale ; No pressure coupling, This means fixed box size. Extended-ensemble pressure coupling where the box vectors are subject to an equation of motion. The equation of motion for the atoms is coupled to this. No instantaneous scaling takes place.
Pcoupltype = semiisotropic ; specifies the kind of isotropy of the pressure coupling used.
tau-p = 5.0 ; [ps] The time constant for pressure coupling (one value for all directions).
ref-p = 1 1
compressibility = 4.5e-15 4.5e-5 ; [bar^{-1}] for water at 1 atm and 300K the compressibility is 4.5e-5 bar-1
refcoord-scaling = all ; The reference coordinates are scaled with the scaling matrix of the pressure coupling.
;;; Velocity generation ;;;
gen-vel = no ; “no” do not generate velocities “yes” generate velocities in grompp at temp gen-temp = 300
gen-temp = 300
;;; Bonds ;;;
constraints = h-bonds ; controls which bonds in the topology will be converted to rigid holonomic constraints.
constraint-algorithm = lincs ; chooses which solver satisfies any non-SETTLE holonomic constraints. “lincs” is faster but does not work with angle constraints
lincs-order = 4 ; accuracy of lincs: the number of matrices in the expansion for the matrix inversion
lincs-iter = 1 ; number of iterative corrections to matrix inversion to compensate for lengthening due to rotation
lincs-warnangle = 30 ; print warning to log file and stderr if bond rotations be more than this angle
;;; Vacuum parameters ;;;
nwall = 2
wall-type = 10-4 ; type of wall potential
wall-atomtype = WR WR ; the atom type name in the force field for each wall.
wall-density = 50 50 ; [nm-3] / [nm-2] the number density of the atoms for each wall
wall atom WR is a fictitious particle whose mass, charge, LJ parameters are 0.
I used 3dc Ewald-sum geometry and xy pbc due to the net dipole of the slabs
The system moves downward even though I removed the center of mass motion and the two slabs merged at the end of NPT equlibration.
Have you tried to equilibrate the system at constant pressure with pbc=xyz and comm-mode = None first? These are the only two options that I believe can cause this issue. You can then change them for NVT and production runs.
Also, you can set the compressibility to 0 in the x/y direction (not that 4.5e-15 should do anything…).