Running a Pd3 Cluster in OPC Water: A Complete LAMMPS Simulation Guide

What this post covers: Every error, fix, and design decision encountered while setting up a production MD simulation of a Pd3 palladium nanocluster in OPC water — from force field files to HPC job submission. Written as a living debug log so others can skip the weeks of trial and error.

System Description

A single Pd3 cluster (3 palladium atoms in a rigid triangle) solvated in 7205 OPC water molecules inside a 60 × 60 × 60 Å periodic cubic box.

Property	Value
Total atoms	21,618
Atom types	5 (Pd1, Pd2, Pd3, Ow, Hw)
Box	60 × 60 × 60 Å, periodic
LAMMPS units	`real` (kcal/mol, Å, fs, AMU)
Simulation target	100 ns NPT at 298.15 K, 1 atm

Force Fields

Pd3 — INTERFACE Force Field (Heinz 2008)

The cluster is an asymmetric rigid triangle. Three distinct atom labels (Pd1, Pd2, Pd3) are required by fftool to map the topology correctly — all three carry identical physical parameters.

# pd3.ff
ATOMS
Pd1   106.4200  0.0000  lj  2.48000  3.58509
Pd2   106.4200  0.0000  lj  2.48000  3.58509
Pd3   106.4200  0.0000  lj  2.48000  3.58509

BONDS
Pd1  Pd2  cons  2.9440  0.0
Pd1  Pd3  cons  2.7476  0.0
Pd2  Pd3  cons  2.7476  0.0

ANGLES
Pd1  Pd2  Pd3  harm  57.61  0.0
Pd2  Pd1  Pd3  harm  57.61  0.0
Pd1  Pd3  Pd2  harm  64.79  0.0

Conversion: ε = 15.0 kJ/mol ÷ 4.184 = 3.58509 kcal/mol

Angle sum check: 57.61 + 57.61 + 64.79 = 180.01° ✓ (0.01° is rounding)

OPC Water — Izadi et al. (2014)

OPC is a 4-site rigid TIP4P-family model. In LAMMPS the virtual M-site is not an explicit atom — pair_style lj/cut/tip4p/long computes it on-the-fly at d_OM = 0.1594 Å from Ow.

# opc.ff
ATOMS
Ow   15.9990   0.0000  lj  3.16655  0.21280
Hw    1.0080   0.6791  lj  0.00000  0.00000

BONDS
Ow  Hw  cons  0.8724  0.0

ANGLES
Hw  Ow  Hw  cons  103.60  0.0

Conversion: ε = 0.890359 kJ/mol ÷ 4.184 = 0.21280 kcal/mol

Charge note: fftool writes Ow with q = 0.0. The script must apply set type 4 charge -1.3582 after read_data to restore charge neutrality. q(Ow) = −2 × 0.6791 = −1.3582 e.

Pd–Water Cross Terms (Lorentz–Berthelot)

\[\sigma_{\text{Pd-Ow}} = \frac{2.48000 + 3.16655}{2} = 2.82328\ \text{Å}\] \[\varepsilon_{\text{Pd-Ow}} = \sqrt{3.58509 \times 0.21280} = 0.87352\ \text{kcal/mol}\]

All Pd–Hw cross terms are explicitly zero (Hw has no LJ site). The Ow–Hw cross term must also be set explicitly to zero:

pair_coeff  4 5  0.00000  0.00000   # Ow-Hw: zero, not arithmetic mix

Without this line, pair_modify mix arithmetic auto-generates ε(Ow-Hw) = 0.10640 kcal/mol — physically wrong for TIP4P.

Workflow

fftool → Packmol → fftool -l → data.lmp → Minimisation → NPT Equil → Production

Step 1 — Pack with fftool + Packmol

  
fftool 1 pd3.xyz 7205 water.xyz -b 60
packmol < pack.inp
fftool 1 pd3.xyz 7205 water.xyz -b 60 -l

The Mandatory LAMMPS Parser Order

This is the single most important rule. LAMMPS is a single-pass parser. It reads each line once, in order. If a style is not registered when the parser hits a coefficient section inside the data file, it crashes.

bond_style / angle_style   ← BEFORE read_data  (Bond Coeffs, Angle Coeffs)
pair_style                 ← BEFORE read_data  (Pair Coeffs)
kspace_style               ← BEFORE read_data  (safe to declare early)
read_data                  ← defines box + atom types
pair_coeff                 ← AFTER  read_data  (needs atom types to exist)
set type 4 charge ...      ← AFTER  read_data

Violating any of these produces one of three hard errors:

ERROR: Must define bond_style before Bond Coeffs
ERROR: Must define pair_style before Pair Coeffs
ERROR: Pair_coeff command before simulation box is defined

Step 2 — Energy Minimisation

Why minimisation is not trivial for this system

Packmol produces geometrically valid but not energy-minimised coordinates. Bad contacts (atom separations < σ) generate forces of thousands of kcal/mol/Å that cause the integrator to diverge on Step 1 of dynamics.

The TIP4P Hydrogen Problem

Attempt v1 — SHAKE during minimise (k_SHAKE = 1,987,206 kcal/mol/Å²):

LAMMPS cannot run iterative SHAKE without a timestep. It substitutes each constraint with a stiff harmonic spring. With 14,410 O–H bonds:

\[F_\text{norm,floor} \approx \sqrt{14410} \times k_\text{SHAKE} \times \delta \approx 2100\ \text{kcal/mol/Å}\]

SD cannot descend through this floor. ftol = 1e-4 is unreachable.

Attempt v2 — No SHAKE, k = 0 bonds:

Hw has zero LJ (ε = 0, σ = 0). With k = 0 bonds and special_bonds zeroing 1–2 Coulomb, Hw feels only inter-molecular Coulomb attraction to neighbouring Ow — with no repulsive barrier. Result: topology collapse.

Step 0:    PotEng = -59,933    Fnorm = 3,862
Step 100:  PotEng = -775,396   Fnorm = 1,184,000,000   ← Hw fell into Ow
ERROR: Out of range atoms - cannot compute PPPM

Solution v5 — Temporary harmonic springs:

Use soft O–H springs during minimisation only (2000× softer than SHAKE). Spring forces → 0 as bonds relax to equilibrium. Fnorm converges cleanly. Reset to k = 0 before write_data.

# During minimisation ONLY:
bond_coeff   4   1000.0  0.8724    # k=1000, prevents Hw drift
angle_coeff  4    100.0  103.600   # k=100,  prevents H-O-H deformation

# After minimisation — reset before write_data:
bond_coeff   4   0.0  0.8724
angle_coeff  4   0.0  103.600

Why k = 1000 is the right value:

Coulomb force on Hw from neighbouring Ow at 2 Å ≈ 77 kcal/mol/Å. Spring force at 0.1 Å stretch = 1000 × 0.1 = 100 kcal/mol/Å > 77. ✓

The Complete Minimisation Script

# ==============================================================
# LAMMPS ENERGY MINIMIZATION v5
# System: Pd3 + OPC water, 21618 atoms
# ==============================================================
units           real
atom_style      full
boundary        p p p

# Styles BEFORE read_data
bond_style      harmonic
angle_style     harmonic
pair_style      lj/cut/tip4p/long 4 5 4 4 0.1594 12.0
kspace_style    pppm/tip4p 1.0e-5

read_data       data.lmp extra/special/per/atom 3

# pair_coeff AFTER read_data
pair_coeff      1 1   3.58509  2.48000    # Pd1-Pd1
pair_coeff      2 2   3.58509  2.48000    # Pd2-Pd2
pair_coeff      3 3   3.58509  2.48000    # Pd3-Pd3
pair_coeff      4 4   0.21280  3.16655    # Ow-Ow
pair_coeff      5 5   0.00000  0.00000    # Hw-Hw
pair_coeff      1 4   0.87352  2.82328    # Pd1-Ow
pair_coeff      2 4   0.87352  2.82328    # Pd2-Ow
pair_coeff      3 4   0.87352  2.82328    # Pd3-Ow
pair_coeff      1 5   0.00000  0.00000    # Pd1-Hw
pair_coeff      2 5   0.00000  0.00000    # Pd2-Hw
pair_coeff      3 5   0.00000  0.00000    # Pd3-Hw
pair_coeff      4 5   0.00000  0.00000    # Ow-Hw (explicit zero)
pair_modify     mix arithmetic
pair_modify     tail yes

set type 4 charge -1.3582

# Temporary springs (see rationale above)
bond_coeff      1   0.0    2.9440
bond_coeff      2   0.0    2.7476
bond_coeff      3   0.0    2.7476
bond_coeff      4   1000.0 0.8724    # TEMPORARY
angle_coeff     1   0.0    57.610
angle_coeff     2   0.0    57.610
angle_coeff     3   0.0    64.790
angle_coeff     4   100.0  103.600   # TEMPORARY

group           PD      type 1 2 3
group           WATER   type 4 5

neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes
comm_modify     cutoff 15.1

# Pd frozen via setforce (fix rigid/small conflicts with minimiser)
fix             FREEZE_PD   PD      setforce 0.0 0.0 0.0

thermo          100
thermo_style    custom step pe ke etotal press vol fnorm
thermo_modify   norm no flush yes

dump            dMin all dcd 500 minimize.dcd

# Stage 1: triage — clear repulsive contacts
min_style       sd
min_modify      dmax 0.01
minimize        0.0  50.0  1000  10000

# Stage 2: bulk relaxation
min_style       sd
min_modify      dmax 0.10
minimize        0.0  1.0  20000  200000

# Stage 3: tight convergence
min_style       sd
min_modify      dmax 0.10
minimize        0.0  1.0e-4  50000  500000

# Reset springs to zero before writing
bond_coeff      4   0.0   0.8724
angle_coeff     4   0.0   103.600

write_data      minimized.data pair ij
write_restart   minimized.restart

Convergence Achieved

Metric	Step 0	Final (~71,000 steps)
PotEng	−59,933 kcal/mol	−106,300 kcal/mol
Fnorm	3,859 kcal/mol/Å	~4 kcal/mol/Å
Pressure	40,509 atm	~−1,876 atm
Max force	73.3 kcal/mol/Å	~1.5 kcal/mol/Å

The residual pressure of −1,876 atm is corrected by the NPT barostat within the first ~10 ps of equilibration.

Minimisation Health Check (no VMD needed)

  
# 1. Did it converge on force tolerance (not max iterations)?
grep "Stopping criterion" log.lammps

# 2. Final force magnitude
grep "Force two-norm\|Force max component" log.lammps

# 3. Charge neutrality
awk '/^Atoms/{f=1;next} /^[A-Z]/{if(f)exit} f&&NF==9{s+=$4}
     END{printf "Net charge: %.4f e\n",s}' minimized.data

# 4. No NaN anywhere
grep -c "nan\|NaN" log.lammps | awk '{print ($1==0)?"NaN: PASS":"NaN: FAIL"}'

# 5. One-liner full report
echo "=== Health Check ===" && \
grep "Stopping criterion" log.lammps && \
grep "Force two-norm" log.lammps && \
awk '/^[[:space:]]*[0-9]/{pe=$2} END{print "Final PotEng:",pe}' log.lammps && \
awk '/^Atoms/{f=1;next}/^[A-Z]/{if(f)exit}f&&NF==9{s+=$4}
     END{printf "Net charge: %.4f e\n",s}' minimized.data

Step 3 — NPT Equilibration

Thermostat/Barostat Choice

The thesis states “Berendsen barostat + Nosé–Hoover thermostat + Langevin barostat” — three mutually contradictory methods. The correct implementation for a true isothermal-isobaric ensemble is fix npt (Martyna–Tobias–Klein Nosé–Hoover chain for both T and P).

SHAKE Cluster Building — Critical Syntax

# WRONG — builds 7205 frozen angles, 0 size-3 clusters:
fix SHAKE_WAT all shake 1.0e-4 20 0 b 4 a 4

# WRONG — group name is case-sensitive, 'water' ≠ 'WATER':
fix SHAKE_WAT water shake 1.0e-4 20 0 m 15.999 1.008

# CORRECT — mass-based, correct group name:
fix SHAKE_WAT WATER shake 1.0e-4 20 0 m 15.999 1.008

Expected log output when correct:

7205 = # of size 3 clusters
   0 = # of frozen angles

Two-Phase Equilibration Script

# ==============================================================
# LAMMPS NPT EQUILIBRATION
# From: minimized.restart
# T = 298.15 K, P = 1 atm, dt = 1 fs
# Phase 1: NVT 200 ps  (temperature equilibration from 0 K)
# Phase 2: NPT 1 ns    (pressure/density equilibration)
# ==============================================================
units           real
atom_style      full
boundary        p p p
bond_style      harmonic
angle_style     harmonic
pair_style      lj/cut/tip4p/long 4 5 4 4 0.1594 12.0
kspace_style    pppm/tip4p 1.0e-5

read_restart    minimized.restart

pair_coeff      1 1   3.58509  2.48000
pair_coeff      2 2   3.58509  2.48000
pair_coeff      3 3   3.58509  2.48000
pair_coeff      4 4   0.21280  3.16655
pair_coeff      5 5   0.00000  0.00000
pair_coeff      1 4   0.87352  2.82328
pair_coeff      2 4   0.87352  2.82328
pair_coeff      3 4   0.87352  2.82328
pair_coeff      1 5   0.00000  0.00000
pair_coeff      2 5   0.00000  0.00000
pair_coeff      3 5   0.00000  0.00000
pair_coeff      4 5   0.00000  0.00000
pair_modify     mix arithmetic
pair_modify     tail yes

bond_coeff      4   0.0   0.8724
angle_coeff     4   0.0  103.600

group           PD      type 1 2 3
group           WATER   type 4 5

neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes
comm_modify     cutoff 15.1

# SHAKE: true constraints now (timestep exists, no spring substitution)
fix   SHAKE_WAT WATER shake 1.0e-4 20 0 m 15.999 1.008
# Rigid/small: Pd cluster, performs own time integration
fix   RIGID_PD  PD    rigid/small molecule

# Velocities: minimisation leaves all at zero
velocity  all create 298.15 4928459 rot yes dist gaussian

# Phase 1: NVT
fix   NVT_WAT WATER nvt temp 298.15 298.15 100.0
timestep  1.0
thermo    1000
thermo_style  custom step temp press pe ke etotal vol density
thermo_modify norm no flush yes
dump  NPT_DCD all dcd 1000 npt_equil.dcd
run   200000
unfix NVT_WAT

# Phase 2: NPT
fix   NPT_WAT WATER npt temp 298.15 298.15 100.0 iso 1.0 1.0 1000.0
restart 100000 restart.npt.a restart.npt.b
run   1000000

write_restart   npt_equil.restart
write_data      npt_equil.data pair ij

HPC Submission — PARAM Shakti (IIT Kharagpur)

Partition Rules

Partition	Cores	Max walltime	Node policy
`shared`	1–36	3 days	Node sharing
`medium`	1–10 nodes × 40 cores	3 days	Exclusive
`large`	1–10 nodes × 40 cores	7 days	Exclusive, low priority
`gpu`	GPU nodes	—	Exclusive

Critical constraint: All jobs must run from /scratch/$USER. The /home directory is for source code and configuration only.

Key SBATCH Rules

# medium/large: always set ntasks-per-node=40 explicitly
#SBATCH --ntasks-per-node=40

Available LAMMPS Modules

  
module avail lammps
# Useful builds:
# apps/lammps/29Aug2024/oneapi-2021  (default, FFTW3, pure MPI)
# apps/lammps/29Aug2024/intel        (GCC 13 build, pure MPI)

Important finding: None of the available LAMMPS builds on PARAM Shakti include the OPENMP or INTEL acceleration packages. All runs are pure MPI. Check with:

  
lmp -h 2>&1 | grep -E "OPENMP|INTEL|USER-OMP"
# Empty output = packages not compiled in

Benchmark Results — 21618 Atoms, Pure MPI

All timings from Performance: line in LAMMPS log. Module: apps/lammps/29Aug2024/oneapi-2021 (FFTW3).

Cores	Nodes	Partition	ns/day	Notes
8	shared	shared	1.44	—
16	shared	shared	1.46	Near-zero scaling — socket boundary
36	shared	shared	5.80	—
40	1	medium	7.08	Full node + InfiniBand
80	2	medium	pending	—
160	4	medium	pending	—
280	7	medium	pending	—
400	10	medium	pending	—

Key observation: 8→16 cores on shared gave almost zero speedup (1.44 → 1.46 ns/day). Shared partition MPI tasks land on different sockets with no InfiniBand — communication overhead kills scaling. Always use medium partition for production.

Timing breakdown at 40 cores:

Pair    : 64%   ← LJ + real-space Coulomb
Kspace  : 18%   ← PPPM FFT (FFTW3 vs KISS: 2x faster)
Neigh   : 12%   ← neighbour list rebuilds (1604/1716 dangerous)
Comm    :  3%

Why dangerous builds are high: Neighbour list rebuild triggered when any atom moves > skin/2 = 1.0 Å. During NVT startup from 0 K, atoms are fast. Fix for equilibrated production:

neighbor        2.5 bin          # larger skin = fewer rebuilds
neigh_modify    delay 0 every 2 check yes  # check every 2 steps not 1

Expected gain: 8–12% improvement in production runs.

Lessons Learned — Ordered by Severity

Critical (will crash the run)

Parser order is strict. pair_coeff must come after read_data. pair_style, bond_style, angle_style must come before read_data.
Hw has zero LJ — never remove its bond restraint. Without a spring or SHAKE, Hw atoms collapse into neighbouring Ow Coulomb wells. No LJ repulsion stops them.
SHAKE is impossible without a timestep. During minimisation, LAMMPS substitutes SHAKE with springs of k = 1,987,206 kcal/mol/Å². This creates an Fnorm floor of ~2100 that SD cannot descend through. Use temporary harmonic bonds (k = 1000) instead.
fix rigid/small conflicts with the minimiser. Use fix setforce 0.0 0.0 0.0 on the cluster during minimisation. Switch to rigid/small for dynamics.
LAMMPS group names are case-sensitive. WATER ≠ water. A wrong group name gives ERROR: Could not find fix group ID.
Ow charge must be set manually. fftool writes q(Ow) = 0.0. Apply set type 4 charge -1.3582 after read_data or the system carries +9785 e net charge and PPPM diverges.

Moderate (wrong physics, may not crash)

pair_modify mix arithmetic does not override explicit pair_coeff. Explicit lines always win. Document cross-term mixing rule carefully.
Ow–Hw cross term is not zero by default. Arithmetic mixing gives ε(Ow-Hw) = 0.10640 kcal/mol. Set pair_coeff 4 5 0.0 0.0 explicitly.
Berendsen barostat ≠ NpT ensemble. Berendsen does not sample the correct isothermal-isobaric distribution. Use fix npt (Nosé–Hoover) for proper statistical mechanics.
etol = 0 in minimise. Non-zero energy tolerance causes false early exit when ΔE per step < etol even though forces are still large (caused an 8-step exit with Fnorm = 2076 in early runs).

Minor (warnings, not errors)

comm_modify cutoff 15.1 silences two TIP4P comm-cutoff warnings. The auto-corrected value is 15.0318 Å. Pre-declaring 15.1 Å prevents both warnings from firing.
Generated N of 10 mixed pair_coeff — always set all 15 pair coefficients (5 self + 10 cross) explicitly. Removes all auto-generation ambiguity.
Non-ASCII characters in print strings cause WARNING: Detected non-ASCII characters in input. Replace em-dashes (—) with -- and ± with +/-.
SHAKE b 4 a 4 builds frozen angles, not size-3 clusters. Use mass-based syntax m 15.999 1.008 for TIP4P-family water.

References

Heinz, H., Vaia, R. A., Farmer, B. L., & Naik, R. R. (2008). Accurate Simulation of Surfaces and Interfaces of FCC Metals Using 12-6 and 9-6 Lennard-Jones Potentials. J. Phys. Chem. C, 112(44), 17281–17290.
Izadi, S., Anandakrishnan, R., & Onufriev, A. V. (2014). Building Water Models: A Different Approach. J. Phys. Chem. Lett., 5(21), 3863–3871.
Thompson, A. P., et al. (2022). LAMMPS — A flexible simulation tool for particle-based materials modeling. Comput. Phys. Commun., 271, 108171.
Martyna, G. J., Tobias, D. J., & Klein, M. L. (1994). Constant pressure molecular dynamics algorithms. J. Chem. Phys., 101(5), 4177–4189.
LAMMPS Documentation (29 Aug 2024). https://docs.lammps.org
PARAM Shakti HPC Documentation, IIT Kharagpur. https://hpc.iitkgp.ac.in/HPCF/paramShakti