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

Overview

This post documents the setup of a molecular dynamics simulation of a Pd₄ palladium nanocluster solvated in OPC water using GROMACS on the PARAM Shakti HPC cluster at IIT Kharagpur. It is a direct continuation of the Pd₃ OPC water post and focuses on the differences introduced by going from a 3-atom triangular cluster to a 4-atom tetrahedral one.

The key new element is the rigid tetrahedron constraint network: whereas Pd₃ needed 3 constraints to rigidify a triangle, Pd₄ needs all 6 edges of the tetrahedron. Everything else — the OPC 4-site correction, the mixing rule fix, the MW insertion script, and the equilibration protocol — is identical to the Pd₃ system and reused without modification.

System: 1 Pd₄ cluster + 7,205 OPC water molecules, 60 × 60 × 60 Å cubic box, 298.15 K, 1 bar — total 28,824 atoms.

---

1. What Changes vs Pd₃

Property	Pd₃	Pd₄
Cluster geometry	Equilateral triangle (D₃h)	Regular tetrahedron (Tₐ)
Atoms	3	4
Constraints needed	3 (all edges of triangle)	6 (all edges of tetrahedron)
Pd–Pd–Pd angle	60°	109.471° = arccos(−⅓)
Total system atoms	28,823	28,824
[ molecules ] count	pd3 1	pd4 1
OPC block	unchanged	identical
Mixing rule fix	comb-rule 2	comb-rule 2
MW virtual site	identical	identical

Everything not in this table is a direct copy from the Pd₃ workflow.

---

2. Force Field and Z-matrix

pd4.ff — key parameters

The LJ parameters are identical to Pd₃ (same element, same INTERFACE FF reference):

Parameter	Value
σ	2.5114 Å
ε	25.7316 kJ/mol
Charge	0.0 e
Pd–Pd bond (DFT equilibrium)	2.6754 Å
Cluster geometry	Regular tetrahedron (Tₐ)

The pd4.ff BONDS section lists all 6 edges explicitly with cons (rigid constraint, zero force constant):

BONDS
Pd1  Pd2  cons  2.6754  0.0000
Pd1  Pd3  cons  2.6754  0.0000
Pd1  Pd4  cons  2.6754  0.0000
Pd2  Pd3  cons  2.6754  0.0000
Pd2  Pd4  cons  2.6754  0.0000
Pd3  Pd4  cons  2.6754  0.0000

A note on the fftool warnings

Running fftool 1 pd4.zmat 7205 opc.zmat --box 60.0 produces many warnings. All are benign:

Angle warnings (pd4 angle ... removed) — fftool found Pd–Pd–Pd angle triplets in the z-matrix and looked for harmonic angle parameters. Finding none, it removed them. This is correct: a rigid tetrahedron held by 6 distance constraints has no internal degrees of freedom — harmonic angles are physically meaningless and are correctly absent.

Dihedral and improper warnings — same logic. fftool searched for torsion parameters, found none, and skipped them. With 6 fixed-length edges the dihedral angles are geometrically determined and cannot change.

The one number that confirms the z-matrix is correct:

pd4.zmat   pd4   1   pd4.ff   4   6   file   +0.0000

natom 4 and nbond 6 = all edges of C(4,2) = 6 pairs. ✓

Z-matrix geometry note

The z-matrix bond lengths (2.5622–2.5628 Å) differ from the DFT equilibrium target (2.6754 Å). This is normal: z-matrix internal coordinates are a parameterisation tool, not the simulation geometry. LINCS snaps all six Pd–Pd distances to 0.26754 nm at the first constrained step.

---

3. Packing the Box

# Generate pd4_pack.xyz and opc_pack.xyz
cp pack.inp pack.inp.mine
fftool 1 pd4.zmat 7205 opc.zmat --box 60.0
cp pack.inp.mine pack.inp

# Custom pack.inp: centre Pd4 at (30,30,30), 5.5 Å exclusion sphere
packmol < pack.inp

# Verify: 4 Pd + 7205×3 = 21619 atoms
head -2 simbox.xyz

# Generate GROMACS files
fftool 1 pd4.zmat 7205 opc.zmat --box 60.0 --gmx

The pack.inp is identical in structure to the Pd₃ version — only the structure file name changes from pd3_pack.xyz to pd4_pack.xyz.

---

4. Correcting the Topology (field.top)

fftool produces the same two errors as for Pd₃, plus one additional issue specific to Pd₄: the angle warnings caused fftool to emit a non-empty [ angles ] section with zero force constants, which must be removed.

The full list of corrections:

1. comb-rule 3 → 2 (Lorentz-Berthelot for Pd–Ow)
2. ! created by fftool → ; (LAMMPS comment char, invalid in GROMACS)
3. OPC: charge moved from Ow to MW virtual site
4. OPC nrexcl 3 → 2
5. MW atomtype added (zero LJ, zero mass)
6. Stale Ow/Hw atomtype charges zeroed
7. OPC [ angles ] harmonic entry removed, replaced with [ settles ]
8. [ virtual_sites3 ] added for OPC MW site
9. OPC [ exclusions ] added (mandatory for SETTLE)
10. Pd₄ [ exclusions ] added for all 6 pairs (closed polyhedron)

The corrected topology in full:

; field.top  --  Pd4 cluster solvated in OPC water  (GROMACS topology)
; Total atoms: 4 + 7205×4 = 28824
; Ref: GROMACS Reference Manual 2022.2

[ defaults ]
; nbfunc   comb-rule   gen-pairs   fudgeLJ   fudgeQQ
; comb-rule 2 = Lorentz-Berthelot (arithmetic σ, geometric ε)
; fftool emits comb-rule 3 (geometric σ) → wrong for INTERFACE FF + OPC
  1        2           yes         0.5       0.5

[ atomtypes ]
; name  at.num     mass    charge  ptype   sigma(nm)     epsilon(kJ/mol)
Pd1        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Pd2        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Pd3        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Pd4        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Ow          8   15.9990   0.0000     A     3.16655e-01   8.90360e-01
Hw          1    1.0080   0.0000     A     1.00000e-01   0.00000e+00
MW          0    0.0000   0.0000     A     0.00000e+00   0.00000e+00

[ moleculetype ]
; name    nrexcl
pd4       3

[ atoms ]
;  nr   type  resnr  residue  atom  cgnr     charge
    1   Pd1      1   pd4      Pd1     1     0.000000
    2   Pd2      1   pd4      Pd2     2     0.000000
    3   Pd3      1   pd4      Pd3     3     0.000000
    4   Pd4      1   pd4      Pd4     4     0.000000

[ constraints ]
;  ai   aj   func    b0(nm)
; All 6 edges of the tetrahedron: C(4,2) = 6 pairs
; b0 = 0.26754 nm from pd4.ff DFT geometry
; 6 constraints on 4 atoms = fully rigid tetrahedron
    1    2     1     0.26754
    1    3     1     0.26754
    1    4     1     0.26754
    2    3     1     0.26754
    2    4     1     0.26754
    3    4     1     0.26754

[ exclusions ]
; All 6 Pd-Pd pairs mutually excluded.
; nrexcl=3 + func-1 constraints auto-generates these,
; but listed explicitly for the closed polyhedral network.
    1    2    3    4
    2    1    3    4
    3    1    2    4
    4    1    2    3

[ moleculetype ]
; name    nrexcl
OPC       2

[ atoms ]
;  nr   type  resnr  residue  atom  cgnr     charge       mass
    1   Ow       1   OPC      O       1     0.000000   15.99900
    2   Hw       1   OPC      H1      1     0.679100    1.00800
    3   Hw       1   OPC      H2      1     0.679100    1.00800
    4   MW       1   OPC      MW      1    -1.358200    0.00000

[ settles ]
;  OW   func    doh(nm)    dhh(nm)
    1    1      0.08724    0.13713

[ exclusions ]
    1    2    3    4
    2    1    3    4
    3    1    2    4
    4    1    2    3

[ virtual_sites3 ]
; site  ai   aj   ak  func      a           b
    4    1    2    3    1      0.147722    0.147722

[ system ]
Pd4 cluster in OPC water

[ molecules ]
pd4       1
OPC       7205

---

5. The Rigid Tetrahedron — Why 6 Constraints

The number of constraints required to fully rigidify N points in 3D is the number of edges in the complete graph K_N, which equals C(N, 2) = N(N−1)/2.

Cluster	Atoms	Edges (constraints)	Geometry
Pd₂	2	1	dimer
Pd₃	3	3	triangle
Pd₄	4	6	tetrahedron
Pd₅	5	10	trigonal bipyramid

For Pd₄ the six pairs are: Pd1–Pd2, Pd1–Pd3, Pd1–Pd4, Pd2–Pd3, Pd2–Pd4, Pd3–Pd4. Once all six distances are fixed at 0.26754 nm, the 3D shape is uniquely determined. No angle terms, no dihedral terms, no improper terms are needed or valid.

Degrees of freedom

DOF	Status
All 6 Pd–Pd bond lengths	Frozen at 0.26754 nm (LINCS)
All Pd–Pd–Pd angles	Frozen at 109.471° (implicit from 6 fixed edges)
Internal vibration	Zero (no internal DOF remain)
Cluster translation	Free — COM moves under net solvent force
Cluster rotation	Free — torque from solvent rotates the body

The true tetrahedral Pd–Pd–Pd angle is:

\[\theta = \arccos\!\left(-\frac{1}{3}\right) = 109.4712°\]

This is the value listed in the ANGLES section of pd4.ff. It differs from the ~60° and ~69° values in the z-matrix, which are internal coordinate construction variables — not the physical bond angle in 3D Cartesian space.

LINCS order recommendation

A tetrahedral 4-atom constraint network is more complex than the Pd₃ triangle. The closed polyhedral topology requires a higher LINCS expansion order for accurate constraint satisfaction:

lincs-order           = 6    ; raise to 8 if LINCS warnings appear
lincs-iter            = 2

---

6. MW Insertion

The add_mw.py script from the Pd₃ post is reused without any modification. It keys on residue name OPC and ignores the solute entirely:

python3 add_mw.py config.pdb config_mw.pdb
# Expected:
# water residues     : 7205   (one MW added to each)
# non-water atoms    : 4
# total atoms written: 28824  (expected 28824)

grep -m4 ' OPC ' config_mw.pdb
# Must be: O, H1, H2, MW in that order

The only difference from Pd₃: non-water atoms = 4 (not 3) and total atoms written = 28824 (not 28823). The script itself is unchanged.

---

7. Index File

gmx_CPUIMPI make_ndx -f config_mw.pdb -o index.ndx
# At the prompt: q
grep '^\[' index.ndx
# Expected: [ System ]  [ Other ]  [ pd4 ]  [ OPC ]

---

8. MDP Files

The three MDP files are identical to the Pd₃ system. The only relevant note is on LINCS, which is already set correctly in those files:

lincs-order           = 6
lincs-iter            = 2

The tetrahedral constraint network is handled without any other changes to the equilibration protocol.

---

9. SLURM Submission

The submission script is identical to the Pd₃ system. Copy it into the Pd₄ working directory on PARAM Shakti and update WORKDIR:

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

set -e
set -x

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

WORKDIR="/scratch/20ch91r15/GROMACS/single_pd4_water/equilibration"
cd "${WORKDIR}"

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

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

${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

${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

${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

Reminder: gmx_CPUIMPI is a real Intel MPI build with no thread-MPI layer. -ntmpi and -ntomp are invalid flags that crash all 40 MPI ranks. Parallelism is controlled exclusively by mpirun -np 40.

---

10. Post-Run Sanity Checks

The check_equilibration.sh script from the Pd₃ post runs without modification. Update two values at the top of the script:

EXPECTED_TOTAL=28824    # was 28823 for Pd3
EXPECTED_PD=4           # was 3 for Pd3
REF_PD_BOND=0.26754     # was 0.26760 for Pd3

All seven checks, the PBC two-step correction, and the Pd–Pd rigidity distance check apply identically. The expected outputs remain:

Observable	Expected	Normal fluctuation
Temperature	298.15 K	± 2–3 K
Pressure	1 bar	± 100–300 bar
Density	~997 kg/m³	± 2–3 kg/m³
Pd–Pd distance	0.26754 nm	< 0.001 nm

---

11. Reusability Pattern

This workflow establishes a general pattern for any rigid Pd_N cluster solvated in OPC water:

What changes with the cluster:

• [ moleculetype ] name and atom count
• [ constraints ]: C(N,2) entries, one per unique pair, with the DFT equilibrium distance from the .ff file
• [ exclusions ]: same C(N,2) pairs
• [ atomtypes ]: add new type names if used (Pd5, Pd6, etc.) with the same INTERFACE FF σ and ε
• [ molecules ]: update name and count
• add_mw.py output atom count: N_Pd + 7205 × 4

What never changes:

• [ defaults ]: always comb-rule 2
• The entire OPC [ moleculetype ] block (SETTLE + virtual_sites3 + exclusions)
• The three MDP files (EM, NVT, NPT)
• The SLURM script
• The add_mw.py script
• The check_equilibration.sh script (update three variables only)

---

References

1. Izadi, S. & Onufriev, A. V. (2016). Accuracy limit of rigid 3-point water models. J. Chem. Phys. 145, 074501.
2. Heinz, H. et al. (2008). Accurate Simulation of Surfaces and Interfaces of Face-Centered Cubic Metals Using 12-6 and 9-6 Lennard-Jones Potentials. J. Phys. Chem. C 112, 17281.
3. Bussi, G., Donadio, D. & Parrinello, M. (2007). Canonical sampling through velocity rescaling. J. Chem. Phys. 126, 014101.
4. Páduai, A. A. H. fftool. https://github.com/paduagroup/fftool
5. GROMACS Reference Manual 2022.2. https://manual.gromacs.org/2022.2/reference-manual/