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

Overview

This post documents the setup of a molecular dynamics simulation containing two Pd₃ palladium nanoclusters solvated in OPC water, where the clusters are placed 45 Å apart along the 3D box diagonal. It is a direct continuation of the single Pd₃ OPC water post and focuses on the differences introduced by having multiple solute molecules.

The key new elements are:

• Diagonal placement geometry — computing exact cluster coordinates for a given separation distance along the space diagonal
• Multi-structure pack.inp — two fixed cluster blocks with two independent outside sphere water exclusions
• Single moleculetype, count = 2 — GROMACS handles both clusters from one topology definition

Everything else — the OPC 4-site correction, mixing rule fix, MW insertion script, MDP files, and SLURM script — is identical to the single-cluster system and reused without modification.

System: 2 Pd₃ clusters + 7,175 OPC water molecules, 60 × 60 × 60 Å cubic box, 298.15 K, 1 bar — total 28,706 atoms.

---

1. What Changes vs Single Pd₃

Property	1× Pd₃	2× Pd₃
Cluster position	(30, 30, 30) Å	(17.01, 17.01, 17.01) and (42.99, 42.99, 42.99) Å
Inter-cluster distance	N/A	45 Å (3D diagonal)
Water count	7,205	7,175 (two exclusion spheres)
Total atoms	28,823	28,706
pack.inp structures	2 (1 cluster + water)	3 (cluster 1 + cluster 2 + water)
[ molecules ]	pd3 1	pd3 2
Moleculetype definitions	1	still 1 (count handles the rest)
OPC block	unchanged	identical
MDP files	unchanged	identical
SLURM script	unchanged	identical

---

2. Diagonal Placement Geometry

The constraint

Two clusters must be separated by exactly 45 Å, symmetric about the box centre (30, 30, 30 Å), inside a 60 × 60 × 60 Å box with 1.5 Å PBC buffer and 5.5 Å exclusion spheres around each cluster.

Derivation

For a 3D diagonal displacement the separation vector is (Δx, Δy, Δz) with Δx = Δy = Δz. The magnitude condition gives:

\[\sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2} = 45\,\text{Å} \implies \Delta x = \frac{45}{\sqrt{3}} = 25.9808\,\text{Å}\]

The half-offset from box centre:

\[\delta = \frac{25.9808}{2} = 12.9904\,\text{Å}\]

Cluster positions:

\[\text{Cluster 1} = (30 - 12.9904,\; 30 - 12.9904,\; 30 - 12.9904) = (17.0096,\; 17.0096,\; 17.0096)\,\text{Å}\] \[\text{Cluster 2} = (30 + 12.9904,\; 30 + 12.9904,\; 30 + 12.9904) = (42.9904,\; 42.9904,\; 42.9904)\,\text{Å}\]

Feasibility checks

Check	Value	Requirement	Status
Inter-cluster distance	45.0000 Å	= 45 Å	✓
Face clearance (cluster 1)	17.01 − 5.5 = 11.51 Å	> 1.5 Å PBC buffer	✓
Face clearance (cluster 2)	60 − 42.99 − 5.5 = 11.51 Å	> 1.5 Å PBC buffer	✓
Water gap between clusters	45 − 2 × 5.5 = 34 Å	> 0 (no overlap)	✓
Cluster–cluster image distance	60 − 45 = 15 Å	> 2 × 5.5 = 11 Å	✓

The last check is the periodic image condition: the nearest image of cluster 2 seen by cluster 1 (through the box face) is 60 − 45 = 15 Å away, which is greater than two exclusion radii. The clusters do not interact with each other’s periodic images.

---

3. Packing the Box

Step 1 — Generate xyz files

# fftool needs the total count of each species for density reporting
# 2 pd3 molecules, 7175 OPC waters
fftool 2 pd3.zmat 7175 opc.zmat --box 60.0

# fftool overwrites pack.inp — restore the custom one
cp pack.inp pack.inp.fftool
cp pack.inp.mine pack.inp

Step 2 — Custom pack.inp

The critical difference from the single-cluster system is the three structure blocks and two outside sphere exclusions for water:

tolerance 2.0
filetype xyz
output simbox.xyz

# Cluster 1: lower diagonal corner
structure pd3_pack.xyz
  number 1
  center
  fixed 17.0096 17.0096 17.0096 0.0 0.0 0.0
end structure

# Cluster 2: upper diagonal corner
structure pd3_pack.xyz
  number 1
  center
  fixed 42.9904 42.9904 42.9904 0.0 0.0 0.0
end structure

# OPC water: fills box interior, excluded from both cluster spheres
structure opc_pack.xyz
  number 7175
  inside box 1.5000 1.5000 1.5000 58.5000 58.5000 58.5000
  outside sphere 17.0096 17.0096 17.0096 5.5
  outside sphere 42.9904 42.9904 42.9904 5.5
end structure

Both cluster blocks reference the same pd3_pack.xyz file — Packmol reads the geometry from the file and applies the fixed transform independently for each block.

Water count rationale

The single-cluster system used 7,205 waters. Each 5.5 Å exclusion sphere displaces approximately:

\[N_{\text{displaced}} = \rho_{\text{water}} \times \frac{4}{3}\pi r^3 = 0.0334\,\text{mol/Å}^3 \times 696.9\,\text{Å}^3 \approx 23\,\text{waters}\]

Two spheres displace ~46 waters. Using 7,175 (= 7,205 − 30, with a small buffer) gives approximately the correct bulk density outside the exclusion zones. The NPT barostat adjusts the box to the exact equilibrium density regardless.

Step 3 — Run Packmol

packmol < pack.inp
# On success: writes simbox.xyz and prints "Success!"

# Verify: 2×3 Pd + 7175×3 OPC sites = 21531 atoms
head -2 simbox.xyz

Step 4 — Generate GROMACS files

fftool 2 pd3.zmat 7175 opc.zmat --box 60.0 --gmx
# Writes: config.pdb  run.mdp  field.top (raw — replace with corrected version)

---

4. Topology (field.top)

The topology is structurally identical to the single-Pd₃ system. The only line that changes is in [ molecules ]:

; Single cluster     →     Two clusters
pd3       1                pd3       2
OPC    7205                OPC    7175

GROMACS uses one pd3 [ moleculetype ] definition and stamps it out twice, applying the per-copy atom offset internally. No second copy of the moleculetype block is needed or valid.

The full corrected topology (corrections identical to the single-cluster post — comb-rule 2, 4-site OPC, MW virtual site):

; field.top  --  2× Pd3 clusters in OPC water
; Total atoms: 2×3 + 7175×4 = 28706

[ defaults ]
; nbfunc   comb-rule   gen-pairs   fudgeLJ   fudgeQQ
  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
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
pd3       3

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

[ constraints ]
;  ai   aj   func    b0(nm)
    2    1     1     0.26760
    3    1     1     0.26760
    2    3     1     0.26760

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

[ 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 ]
2x Pd3 clusters in OPC water

[ molecules ]
pd3       2
OPC    7175

---

5. MW Insertion

The add_mw.py script is reused without modification — it keys on residue name OPC and ignores the solute:

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

# Verify atom order in first water block
grep -m4 ' OPC ' config_mw.pdb
# Must be: O, H1, H2, MW

non-water atoms = 6 confirms both clusters (2 × 3 Pd) are present.

---

6. Index File

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

[ pd3 ] will contain 6 atoms — both clusters combined under one residue name. This is correct for the thermostat: both clusters are identical Pd and are coupled to the same temperature group. If you need to analyse them separately (e.g. track the inter-cluster distance over time), create named subgroups at the make_ndx prompt:

r 1
name 4 pd3_1
r 2
name 5 pd3_2
q

---

7. MDP Files

The three MDP files (01_em.mdp, 02_nvt.mdp, 03_npt.mdp) are identical to the single-cluster system. The tc-grps line references pd3 — which now matches a 6-atom group (both clusters) rather than a 3-atom group. This is physically correct: both clusters thermostat together.

tc-grps  = pd3  OPC
tau-t    = 2.0  2.0
ref-t    = 298.15  298.15

---

8. SLURM Submission

The submission script is identical to the single-cluster system. Update only WORKDIR and --job-name:

#!/bin/bash
#SBATCH --partition=medium
#SBATCH --job-name=2pd3_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

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

WORKDIR="/scratch/20ch91r15/GROMACS/two_pd3_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 has no thread-MPI layer. -ntmpi and -ntomp are invalid. All parallelism is controlled by mpirun -np 40.

---

9. grompp Validation

Always run grompp for EM before submitting:

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

Expected output lines to verify:

Generated 28 of the 28 non-bonded parameter combinations
Excluding 3 bonded neighbours molecule type 'pd3'
Excluding 2 bonded neighbours molecule type 'OPC'
Cleaning up constraints and constant bonded interactions with virtual sites

Then validate the NVT coupling groups are found:

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"

Expected:

Number of degrees of freedom in T-Coupling group pd3 is 12.00
Number of degrees of freedom in T-Coupling group OPC is 85788.00

pd3 DOF = 12 because 2 clusters × 3 atoms × (3 translational − 3 constraint DOF per rigid body) are replaced by the rigid body DOF. If this line instead shows T-Coupling group rest, the index group name does not match tc-grps in the mdp — rename the group in make_ndx or update the mdp to match.

---

10. Post-Run Sanity Checks

Update three variables in check_equilibration.sh:

EXPECTED_TOTAL=28706    # was 28823 for single Pd3
EXPECTED_PD=6           # was 3 for single Pd3  (2 clusters × 3 atoms)
REF_PD_BOND=0.26760     # unchanged

The Pd–Pd distance check (gmx distance -select 'atomnr 1 2') still measures the Pd1–Pd2 bond within the first cluster (atoms 1 and 2 in the coordinate file). This confirms cluster 1 is rigid. To additionally verify cluster 2, add:

gmx_CPUIMPI distance \
    -f npt_center.xtc -s npt.tpr -n index.ndx \
    -select 'atomnr 4 5' \
    -oall check_pd_distance_cluster2.xvg -xvg none

Atoms 4 and 5 are Pd1 and Pd2 of the second cluster (offset by 3 from the first). The distance should be an identical flat line at 0.2676 nm.

Inter-cluster distance tracking

To verify the two clusters maintain their separation and do not drift together (or apart under PBC), measure the COM-to-COM distance over the NPT trajectory:

# Create a group containing both cluster COMs
echo -e "2\n2\n" | gmx_CPUIMPI distance \
    -f npt_center.xtc \
    -s npt.tpr \
    -n index.ndx \
    -select 'com of group "pd3_1" plus com of group "pd3_2"' \
    -oall inter_cluster_distance.xvg -xvg none

This requires the pd3_1 and pd3_2 subgroups created in the make_ndx step. The starting distance is 45 Å; NPT box fluctuations will cause small variations (±1–2 Å) but no systematic drift in a well-equilibrated system.

---

11. Reusability Pattern for N Clusters

To generalise to N clusters of the same type in OPC water:

pack.inp: add one structure block per cluster with a unique fixed position. Add one outside sphere line per cluster to the water block. The pd3_pack.xyz file is reused for every block.

field.top: change only [ molecules ] count. One moleculetype definition serves all copies.

check_equilibration.sh: update EXPECTED_PD = N × 3 and EXPECTED_TOTAL = N×3 + n_water×4.

MDP files: tc-grps = pd3 OPC remains correct — the pd3 group grows with N but the coupling logic is unchanged.

---

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/