GROMACS Equilibration and Production Protocol for Solvated Nanoclusters

Overview

This post documents the generic equilibration and production protocol used for all solvated nanocluster simulations in this series. The three-phase equilibration pipeline (EM → NVT → NPT) and the subsequent production run are system-independent — the same MDP files, SLURM scripts, and sanity-check workflow apply to any system once the topology and coordinate files are correctly prepared.

This post assumes the following inputs already exist:

File	Source
`field.top`	Corrected topology (see system-specific posts)
`config_mw.pdb`	Packed + MW-inserted coordinate file
`index.ndx`	Generated by `gmx make_ndx`

1. The Equilibration Pipeline

Three sequential phases bring the system from a Packmol-packed initial geometry to a thermally and mechanically equilibrated NPT state.

Packmol box  →  EM  →  NVT (1–10 ns)  →  NPT (1–10 ns)  →  Production
(steric overlaps)  (geometry)  (temperature)  (density/pressure)

The order is mandatory. Skipping directly to NPT from a packed box causes pressure spikes that destabilise the Parrinello-Rahman barostat. Each phase has a specific physical purpose:

Phase	Box vectors	Temperature	Purpose
EM	Fixed	None	Remove steric clashes
NVT	Fixed	Coupled	Equilibrate temperature
NPT	Free	Coupled	Equilibrate density and pressure
Production	Free	Coupled	Sample the equilibrium ensemble

2. Consistent Design Decisions Across All Phases

These choices apply to every phase and every system in this series:

constraints = all-bonds everywhere — constrains all bonds including metal–metal distances (via LINCS) and water O–H (via SETTLE). This allows a 2 fs timestep safely. Using h-bonds instead would leave metal–metal bonds as flexible harmonic springs requiring a sub-femtosecond timestep.

coulombtype = PME everywhere including EM — the OPC water model places a point charge (MW = −1.3582 e) on a virtual site off the oxygen. A plain cutoff severely mis-treats this discrete off-atom charge. PME handles it correctly in all phases.

DispCorr = No in EM only — the long-range dispersion tail correction assumes a converged homogeneous liquid density. A freshly Packmol-packed box has voids. Applying the correction during EM gives a nonsense pressure value. It is switched on from NVT onward.

integrator = md for dynamics — leap-frog integration with explicit thermostat and barostat. Not sd (stochastic dynamics), which is its own thermostat and makes tcoupl irrelevant. Using md keeps the thermostat and barostat clearly defined and separable in the methods section.

tcoupl = V-rescale with separate groups — the Bussi et al. 2007 velocity-rescaling thermostat correctly samples the canonical ensemble. Separate coupling groups for solute and solvent prevent the “hot solute / cold solvent” artefact where a small metal cluster drains or injects energy into the large water bath.

pcoupl = Parrinello-Rahman for NPT — produces a correct isobaric-isothermal ensemble. The GROMACS manual explicitly states Berendsen barostat does not give a correct NPT ensemble because it suppresses volume fluctuations. Berendsen may be used for a short pre-equilibration if the starting pressure is very far from 1 bar, but all production-quality runs use Parrinello-Rahman.

lincs-order = 6, lincs-iter = 2 — the default LINCS order of 4 is insufficient for closed constraint networks (triangles, tetrahedra). Order 6 with 2 iterative corrections handles these topologies correctly.

3. MDP Files

3.1 Energy Minimisation (`01_em.mdp`)

  
; =============================================================
;  01_em.mdp  --  Steepest-descent energy minimisation
;  Purpose: remove steric clashes from the Packmol-packed box
; =============================================================

integrator            = steep
emtol                 = 1000.0        ; stop when Fmax < 1000 kJ/mol/nm
emstep                = 0.01          ; initial step size [nm]
nsteps                = 50000         ; max steps (early exit expected)

nstlog                = 500
nstenergy             = 500

cutoff-scheme         = Verlet
pbc                   = xyz

coulombtype           = PME
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = No            ; off: packed box density not converged

constraints           = all-bonds
constraint-algorithm  = LINCS
lincs-order           = 6
lincs-iter            = 2

Key parameters explained:

emtol = 1000 — the minimiser stops when the largest force on any atom falls below this threshold. 1000 kJ/mol/nm is deliberately loose: the goal is to remove catastrophic overlaps, not to find the true energy minimum. Tighter tolerances waste time on a configuration that will be immediately heated in NVT.

nsteps = 50000 — a ceiling, not a target. For a ~29k-atom water box convergence typically occurs in 500–5000 steps. If the full 50000 steps run without converging, the starting configuration has severe overlaps — inspect Packmol’s output and the exclusion sphere radii.

3.2 NVT Equilibration (`02_nvt.mdp`)

  
; =============================================================
;  02_nvt.mdp  --  NVT equilibration
;  Purpose: heat the system to target T at fixed volume
;  Duration: 1–10 ns depending on system
; =============================================================

integrator            = md
dt                    = 0.002         ; 2 fs
nsteps                = 5000000       ; 10 ns  (adjust as needed)

nstlog                = 5000          ; every 10 ps
nstenergy             = 5000
nstxout-compressed    = 5000
nstxout               = 0
nstvout               = 0

cutoff-scheme         = Verlet
verlet-buffer-tolerance = 0.005
pbc                   = xyz

coulombtype           = PME
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = EnerPres      ; tail correction on from NVT onward

; --- Temperature coupling ---
; V-rescale: correctly samples canonical ensemble (Bussi et al. 2007)
; Separate groups prevent hot-solute/cold-solvent artefact
; GROUP NAMES must match groups in index.ndx
tcoupl                = V-rescale
tc-grps               = <solute>  <solvent>   ; e.g. pd3  OPC
tau-t                 = 2.0       2.0
ref-t                 = 298.15    298.15

; --- No pressure coupling in NVT ---
; Box vectors fixed: equilibrate T before releasing box
pcoupl                = no

; --- Fresh start after EM ---
gen-vel               = yes
gen-temp              = 298.15
gen-seed              = -1

constraints           = all-bonds
constraint-algorithm  = LINCS
lincs-order           = 6
lincs-iter            = 2
continuation          = no

comm-mode             = Linear
nstcomm               = 100

Key parameters explained:

tc-grps — replace <solute> and <solvent> with the actual group names from your index.ndx (e.g. pd3 OPC, pd4 OPC). These must match exactly — a mismatch causes grompp to fall back to a single rest group, thermostating everything together.

gen-vel = yes — generate velocities from a Maxwell-Boltzmann distribution at gen-temp. This is a fresh start after the static minimisation — the system has no velocities yet.

pcoupl = no — box vectors are completely fixed during NVT. Releasing the box before temperature is equilibrated causes pressure spikes that can crash the Parrinello-Rahman barostat in the subsequent NPT phase.

3.3 NPT Equilibration (`03_npt.mdp`)

  
; =============================================================
;  03_npt.mdp  --  NPT equilibration
;  Purpose: equilibrate density and pressure at target T, P
;  Duration: 1–10 ns depending on system
; =============================================================

integrator            = md
dt                    = 0.002
nsteps                = 5000000       ; 10 ns

nstlog                = 5000
nstenergy             = 5000
nstxout-compressed    = 5000
nstxout               = 0
nstvout               = 0

cutoff-scheme         = Verlet
verlet-buffer-tolerance = 0.005
pbc                   = xyz

coulombtype           = PME
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = EnerPres

tcoupl                = V-rescale
tc-grps               = <solute>  <solvent>
tau-t                 = 2.0       2.0
ref-t                 = 298.15    298.15

; --- Parrinello-Rahman barostat ---
; Produces correct NPT ensemble; Berendsen does not.
; tau-p must be > tau-t to avoid coupling resonance.
pcoupl                = Parrinello-Rahman
pcoupltype            = isotropic     ; cubic box stays cubic
tau-p                 = 2.0           ; [ps]
ref-p                 = 1.0           ; [bar]
compressibility       = 4.5e-5        ; [bar-1] liquid water at 298 K

; --- Continue from NVT checkpoint ---
; gen-vel = no  : velocities come from nvt.cpt via -t flag in grompp
; continuation  : tells grompp not to check for missing velocities
gen-vel               = no
continuation          = yes

constraints           = all-bonds
constraint-algorithm  = LINCS
lincs-order           = 6
lincs-iter            = 2

comm-mode             = Linear
nstcomm               = 100

Key parameters explained:

tau-p = 2.0 ps — the coupling time constant for pressure. Must be significantly larger than tau-t (also 2.0 ps here — at the minimum safe ratio). For systems where pressure oscillations appear in the NPT trajectory, increase to 5.0 ps.

compressibility = 4.5e-5 bar-1 — the isothermal compressibility of liquid water at 298 K. This is the standard value for any aqueous simulation and should not be changed unless simulating at very different conditions (high pressure, supercritical water, etc.).

continuation = yes with gen-vel = no — together these tell grompp that velocities come from an external checkpoint (-t nvt.cpt in the grompp command). Without -t nvt.cpt, this combination would give zero velocities and a dead simulation.

4. SLURM Submission — Equilibration

All three phases run inside a single SLURM job to avoid queue re-entry delays and ensure each phase inherits the exact checkpoint from the previous one. The set -e directive halts execution immediately if any phase fails — NVT never starts if EM crashes.

  
#!/bin/bash
#SBATCH --partition=medium
#SBATCH --job-name=sys_equil
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=40
#SBATCH --time=03-00:00:00
#SBATCH --output=equil_%j.out
#SBATCH --error=equil_%j.err

# Source GMXRC BEFORE set -euo pipefail
# GMXRC references variables that may be unbound on a clean login shell;
# the -u flag would abort immediately if sourced after set -euo pipefail.
source /home/apps/gromacs/gromacs-2022.2/installCPUIMPI/bin/GMXRC.bash

set -e

WORKDIR="/scratch/${USER}/GROMACS/<system>/equilibration"
cd "${WORKDIR}" || { echo "ERROR: cannot cd to ${WORKDIR}"; exit 1; }

MPIRUN="mpirun -np 40"
GMX="gmx_CPUIMPI"

# Pre-flight: verify all inputs exist before consuming queue allocation
for f in field.top config_mw.pdb index.ndx 01_em.mdp 02_nvt.mdp 03_npt.mdp; do
    [[ -f "$f" ]] || { echo "ERROR: missing $f"; exit 1; }
done

# ── Phase 1: Energy Minimisation ──────────────────────────────────────────
# grompp: pre-process inputs into a portable run input (.tpr)
# -maxwarn 1: allow the one expected CONECT warning from the PDB format
${GMX} grompp -f 01_em.mdp -c config_mw.pdb -p field.top \
              -n index.ndx -o em.tpr -maxwarn 1
${MPIRUN} ${GMX} mdrun -v -deffnm em

# ── Phase 2: NVT Equilibration ─────────────────────────────────────────────
${GMX} grompp -f 02_nvt.mdp -c em.gro -p field.top \
              -n index.ndx -o nvt.tpr -maxwarn 1
${MPIRUN} ${GMX} mdrun -v -deffnm nvt

# ── Phase 3: NPT Equilibration ─────────────────────────────────────────────
# -t nvt.cpt: CRITICAL — passes NVT velocities into NPT.
# Without this flag the NPT run starts from zero velocities,
# discarding the entire NVT equilibration.
${GMX} grompp -f 03_npt.mdp -c nvt.gro -t nvt.cpt \
              -p field.top -n index.ndx -o npt.tpr -maxwarn 1
${MPIRUN} ${GMX} mdrun -v -deffnm npt

gmx_CPUIMPI parallelism: this binary is compiled with real Intel MPI only — no thread-MPI layer. The flags -ntmpi and -ntomp are invalid and crash all 40 MPI ranks with the error “Setting the number of thread-MPI ranks is only supported with thread-MPI”. All parallelism is controlled exclusively by mpirun -np 40.

PARAM Shakti filesystem policy: all job I/O must be on /scratch/$USER. The /home path is not on the fast parallel filesystem. Never run mdrun from /home.

5. grompp Validation Before Every Submission

Always run grompp interactively before submitting any job. It costs zero queue time and catches errors that would otherwise waste hours of allocation:

  
# Validate EM
gmx_CPUIMPI grompp -f 01_em.mdp -c config_mw.pdb \
    -p field.top -n index.ndx -o em.tpr -maxwarn 1

# Validate NVT — check T-coupling groups specifically
gmx_CPUIMPI grompp -f 02_nvt.mdp -c config_mw.pdb \
    -p field.top -n index.ndx -o nvt_test.tpr -maxwarn 1 \
    2>&1 | grep "T-Coupling\|ERROR"

rm -f nvt_test.tpr

What to look for in grompp output:

Excluding 3 bonded neighbours molecule type '<solute>'   ← topology read
Excluding 2 bonded neighbours molecule type 'OPC'        ← water read
Number of degrees of freedom in T-Coupling group <solute> is X.00
Number of degrees of freedom in T-Coupling group OPC is Y.00

If instead you see T-Coupling group rest — the index group names do not match tc-grps in the mdp. Either rename the groups in make_ndx or update the mdp tc-grps line to match the actual group names.

6. Post-Equilibration Sanity Checks

After the equilibration job finishes, extract and verify all key observables. First find the correct energy term indices for your system (they vary between EM, NVT, and NPT edr files):

  
# List all available terms in each edr file
echo "0" | gmx_CPUIMPI energy -f em.edr  2>&1 | grep -A40 "Select the terms"
echo "0" | gmx_CPUIMPI energy -f nvt.edr 2>&1 | grep -A40 "Select the terms"
echo "0" | gmx_CPUIMPI energy -f npt.edr 2>&1 | grep -A40 "Select the terms"

For a typical solvated nanocluster system the indices are:

Term	EM edr	NVT edr	NPT edr
Potential	4	5	5
Temperature	—	9	9
Pressure	5	11	11
Volume	—	—	16
Density	—	—	17

Always verify indices from your own edr — they shift when energy groups change. Never use hardcoded indices from another system.

Extract and check each observable:

  
# EM: potential energy descent
echo "4 0" | gmx_CPUIMPI energy -f em.edr -o em_potential.xvg -xvg none

# NVT: temperature
echo "9 0" | gmx_CPUIMPI energy -f nvt.edr -o nvt_temp.xvg -xvg none

# NPT: temperature, density, pressure, volume
echo "9 0"  | gmx_CPUIMPI energy -f npt.edr -o npt_temp.xvg     -xvg none
echo "17 0" | gmx_CPUIMPI energy -f npt.edr -o npt_density.xvg  -xvg none
echo "11 0" | gmx_CPUIMPI energy -f npt.edr -o npt_pressure.xvg -xvg none
echo "16 0" | gmx_CPUIMPI energy -f npt.edr -o npt_volume.xvg   -xvg none

Check averages over the last 50% of each trajectory (the steady-state half):

  
avg_last_half() {
    local file=$1
    awk 'NF==2 && !/^[@#]/{lines[NR]=$2; total++}
         END{start=int(total/2)+1; sum=0; n=0;
             for(i=start;i<=total;i++){sum+=lines[i];n++}
             printf "%.3f\n", sum/n}' "$file"
}

echo "NVT temperature  : $(avg_last_half nvt_temp.xvg) K     (target 298.15)"
echo "NPT temperature  : $(avg_last_half npt_temp.xvg) K     (target 298.15)"
echo "NPT density      : $(avg_last_half npt_density.xvg) kg/m3  (target ~997)"
echo "NPT pressure     : $(avg_last_half npt_pressure.xvg) bar   (target 1)"

Expected convergence criteria:

Observable	Target	Acceptable range
Temperature	298.15 K	± 5 K
Density	~997 kg/m³	± 10 kg/m³
Pressure (average)	1 bar	± 300 bar
Pd–Pd distance	FF value nm	± 0.001 nm

If density is still drifting at the end of the NPT run, extend by another 1–10 ns before starting production.

7. PBC Correction

Raw GROMACS trajectories use the minimum-image convention — molecules that straddle a box face appear split across opposite walls. Two trjconv calls fix this in the correct order.

Order is critical: making molecules whole first, then centring. If centring is applied before making molecules whole, it can re-break molecules that were just fixed.

  
# Step A: make all molecules whole (un-wrap across box faces)
# Select group 0 (System) for output
echo "0" | gmx_CPUIMPI trjconv \
    -f npt.xtc -s npt.tpr -n index.ndx \
    -o npt_pbc.xtc -pbc mol -ur compact

# Step B: centre the solute cluster in the box
# Select solute group for centering (e.g. 2 = pd3), then 0 (System) for output
echo -e "2\n0" | gmx_CPUIMPI trjconv \
    -f npt_pbc.xtc -s npt.tpr -n index.ndx \
    -o npt_center.xtc -pbc mol -center

# Apply the same correction to the final frame
echo -e "2\n0" | gmx_CPUIMPI trjconv \
    -f npt.gro -s npt.tpr -n index.ndx \
    -o npt_center.gro -pbc mol -center

npt_center.xtc and npt_center.gro are the files to use for all structural analysis and as input to the production run.

8. Production MDP (`production.mdp`)

  
; =============================================================
;  production.mdp  --  NPT production run
;  Duration: 100 ns  (adjust nsteps as required)
; =============================================================

integrator            = md
dt                    = 0.002
nsteps                = 50000000      ; 100 ns

; Write every 10 ps → 10,000 frames over 100 ns
nstlog                = 5000
nstenergy             = 5000
nstxout-compressed    = 5000          ; .xtc trajectory
nstxout               = 0
nstvout               = 0

cutoff-scheme         = Verlet
verlet-buffer-tolerance = 0.005
pbc                   = xyz

coulombtype           = PME
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = EnerPres

tcoupl                = V-rescale
tc-grps               = <solute>  <solvent>
tau-t                 = 2.0       2.0
ref-t                 = 298.15    298.15

pcoupl                = Parrinello-Rahman
pcoupltype            = isotropic
tau-p                 = 2.0
ref-p                 = 1.0
compressibility       = 4.5e-5

; Continue from NPT checkpoint — do NOT set gen-vel = yes
gen-vel               = no
continuation          = yes

constraints           = all-bonds
constraint-algorithm  = LINCS
lincs-order           = 6
lincs-iter            = 2

comm-mode             = Linear
nstcomm               = 100

9. Production grompp

Always validate grompp interactively before submitting:

  
gmx_CPUIMPI grompp \
    -f production.mdp \
    -c ../equilibration/npt_center.gro \
    -t ../equilibration/npt.cpt \
    -p field.top \
    -n index.ndx \
    -o production.tpr \
    -maxwarn 1

Three lines to verify in the output:

Number of degrees of freedom in T-Coupling group <solute> is X.00
Number of degrees of freedom in T-Coupling group OPC is Y.00
There were 0 warnings

Why npt_center.gro and npt.cpt together:

-c npt_center.gro — PBC-corrected coordinates with the cluster centred in the box. Cluster starts at maximum distance from all faces.
-t npt.cpt — carries the full phase-space state (positions + velocities + thermostat/barostat state) from the last NPT frame. grompp takes coordinates from -c and velocities from -t, discarding the checkpoint’s positions. You get clean geometry with physically correct thermal velocities.

Starting time for run is 0 ps — grompp resets the production clock to zero regardless of the NPT end time. This is correct and expected.

10. SLURM Submission — Production

  
#!/bin/bash
#SBATCH --partition=large
#SBATCH --job-name=sys_prod
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=40
#SBATCH --time=03-00:00:00
#SBATCH --output=prod_%j.out
#SBATCH --error=prod_%j.err

source /home/apps/gromacs/gromacs-2022.2/installCPUIMPI/bin/GMXRC.bash

set -e

PROD_DIR="/scratch/${USER}/GROMACS/<system>/production"
EQUIL_DIR="/scratch/${USER}/GROMACS/<system>/equilibration"

cd "${PROD_DIR}"

MPIRUN="mpirun -np 40"
GMX="gmx_CPUIMPI"
DEFFNM="production"

# Pre-flight checks
for f in production.mdp field.top index.ndx; do
    [[ -f "$f" ]] || { echo "ERROR: missing $f"; exit 1; }
done
for f in "${EQUIL_DIR}/npt_center.gro" "${EQUIL_DIR}/npt.cpt"; do
    [[ -f "$f" ]] || { echo "ERROR: missing $f"; exit 1; }
done

# grompp
${GMX} grompp \
    -f  production.mdp \
    -c  "${EQUIL_DIR}/npt_center.gro" \
    -t  "${EQUIL_DIR}/npt.cpt" \
    -p  field.top \
    -n  index.ndx \
    -o  "${DEFFNM}.tpr" \
    -maxwarn 1

# mdrun
# -cpi  : checkpoint input — enables clean restart if job hits walltime
# -append : append to existing output files on restart
${MPIRUN} ${GMX} mdrun \
    -v \
    -deffnm  "${DEFFNM}" \
    -cpi     "${DEFFNM}.cpt" \
    -append

Checkpoint-based restart

If the job hits the walltime limit before completing 100 ns, resubmit the identical sbatch production.sh command without any changes. mdrun finds production.cpt, appends to the existing output files, and continues from the last saved frame. No manual intervention needed.

  
# Check how far the run got
tail -5 production.log   # shows last simulation time
grep "Performance" production.log | tail -1   # shows ns/day

# Resubmit if incomplete
sbatch production.sh

Partition choice

Use large (7-day max) rather than medium (3-day max) for 100 ns production runs. At ~100 ns/day on 40 cores the run needs ~24 h, but leaving headroom for queue wait time and slower-than-expected throughput is good practice. If the large queue is congested, switch to medium — 100 ns at 100 ns/day fits within 3 days.

11. Common Errors and Fixes

`Invalid directive found in batch script: ;`

SBATCH directives do not allow inline comments with ;. Every #SBATCH line must end at the value — no trailing comments:

  
# Wrong
#SBATCH --time=03-00:00:00    ; 3 days

# Correct
#SBATCH --time=03-00:00:00

Fix all occurrences:

  
sed -i 's/#SBATCH \(.*\);.*/#SBATCH \1/' production.sh

`Setting the number of thread-MPI ranks is only supported with thread-MPI`

The -ntmpi flag was passed to gmx_CPUIMPI, which is a real MPI build with no thread-MPI layer. Remove -ntmpi and -ntomp from every mdrun call. All parallelism is controlled by mpirun -np 40.

`T-Coupling group rest` instead of named groups

The group names in tc-grps in the mdp do not match the group names in index.ndx. Check actual group names:

grep '^\[' index.ndx

Then update the mdp to match:

  
sed -i 's/tc-grps\s*=.*/tc-grps               = <actual_solute_name>  OPC/' \
    02_nvt.mdp 03_npt.mdp production.mdp

`File 'field.top' does not exist`

Files were not copied to the production directory. Always copy the topology and index file before running grompp from a new directory:

cp ../equilibration/field.top ../equilibration/index.ndx .

`Unknown left-hand 'nstxout-compressed-precision'`

This option was introduced in GROMACS 2023+. The PARAM Shakti installation (2022.2) does not recognise it. Remove from the mdp:

  
sed -i '/nstxout-compressed-precision/d' production.mdp

Energy term index gives wrong values

Energy term indices are not fixed — they differ between EM, NVT, and NPT edr files and between systems. Always list available terms first:

  
echo "0" | gmx_CPUIMPI energy -f npt.edr 2>&1 | grep -A50 "Select the terms"

Then use the correct index for your specific edr file.

12. Files Required for Each Submission

Equilibration

equilibration/
├── config_mw.pdb      ← packed + MW-inserted coordinates
├── field.top          ← corrected topology
├── index.ndx          ← from gmx make_ndx
├── 01_em.mdp
├── 02_nvt.mdp
├── 03_npt.mdp
└── equilibration.sh

Production

production/
├── production.mdp
├── field.top          ← copied from equilibration/
├── index.ndx          ← copied from equilibration/
└── production.sh

Production inputs from equilibration (referenced by path, not copied):

equilibration/npt_center.gro   ← PBC-corrected final frame
equilibration/npt.cpt          ← equilibrated velocities

References

Bussi, G., Donadio, D. & Parrinello, M. (2007). Canonical sampling through velocity rescaling. J. Chem. Phys. 126, 014101.
Parrinello, M. & Rahman, A. (1981). Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 52, 7182.
Hess, B. (2008). P-LINCS: A parallel linear constraint solver for molecular simulation. J. Chem. Theory Comput. 4, 116–122.
GROMACS Reference Manual 2022.2. https://manual.gromacs.org/2022.2/reference-manual/
GROMACS mdp options reference. https://manual.gromacs.org/current/user-guide/mdp-options.html

GROMACS Equilibration and Production Protocol for Solvated Nanoclusters

Overview

1. The Equilibration Pipeline

2. Consistent Design Decisions Across All Phases

3. MDP Files

3.1 Energy Minimisation (01_em.mdp)

3.2 NVT Equilibration (02_nvt.mdp)

3.3 NPT Equilibration (03_npt.mdp)

4. SLURM Submission — Equilibration

5. grompp Validation Before Every Submission

6. Post-Equilibration Sanity Checks

7. PBC Correction

8. Production MDP (production.mdp)

9. Production grompp

10. SLURM Submission — Production

Checkpoint-based restart

Partition choice

11. Common Errors and Fixes

Invalid directive found in batch script: ;

Setting the number of thread-MPI ranks is only supported with thread-MPI

T-Coupling group rest instead of named groups

File 'field.top' does not exist

Unknown left-hand 'nstxout-compressed-precision'

Energy term index gives wrong values

12. Files Required for Each Submission

Equilibration

Production

References

Further Reading

Setting Up a Pd₃ Nanocluster in OPC Water for GROMACS MD Simulation

Setting Up a Pd₄ Nanocluster in OPC Water for GROMACS MD Simulation

Two Pd₃ Nanoclusters in OPC Water — Diagonal Placement and GROMACS Setup

3.1 Energy Minimisation (`01_em.mdp`)

3.2 NVT Equilibration (`02_nvt.mdp`)

3.3 NPT Equilibration (`03_npt.mdp`)

8. Production MDP (`production.mdp`)

`Invalid directive found in batch script: ;`

`Setting the number of thread-MPI ranks is only supported with thread-MPI`

`T-Coupling group rest` instead of named groups

`File 'field.top' does not exist`

`Unknown left-hand 'nstxout-compressed-precision'`