Pd3 Cluster in OPC Water: Complete MD Simulation Guide — LAMMPS and GROMACS

What this post covers: A complete, production-ready walkthrough of setting up a Pd3 palladium nanocluster solvation simulation — from raw force field files through 100 ns production runs — in both LAMMPS and GROMACS. Every error encountered, fixed, and explained. Written as a living debug log so others skip weeks of trial and error.

Bottom line up front: GROMACS is 14× faster than LAMMPS for this system on CPU-only hardware. Both engines produce physically equivalent results. Choose LAMMPS if you need custom force field flexibility; choose GROMACS if you need throughput.

---

1. System Description

A single Pd3 palladium nanocluster solvated in OPC water inside a 60 × 60 × 60 Å periodic cubic box, replicating the methodology from:

“Interaction of Palladium Clusters in Different Solvents” (thesis)

Property	Value
Pd3 cluster atoms	3 (rigid asymmetric triangle)
Water molecules	7,205 OPC (4-site TIP4P-family)
Total atoms (LAMMPS)	21,618 (M-site implicit)
Total atoms (GROMACS)	28,823 (M-site explicit)
Box	60 × 60 × 60 Å, periodic in xyz
Temperature	298.15 K
Pressure	1 atm / 1 bar
Production target	100 ns NPT

---

2. Force Field Parameters

2.1 Pd3 Cluster — INTERFACE Force Field

Source: Heinz et al., J. Phys. Chem. C 2008, 112(44), 17281–17290.

Three distinct atom labels (Pd1, Pd2, Pd3) are required by fftool to map the triangle topology. All three carry identical physical parameters.

# pd3.ff -- fftool format
ATOMS
# type   element  mass      charge  pot  sigma/A   epsilon/[kcal/mol]
Pd1      Pd       106.4200  0.0000  lj   2.48000   3.58509
Pd2      Pd       106.4200  0.0000  lj   2.48000   3.58509
Pd3      Pd       106.4200  0.0000  lj   2.48000   3.58509

BONDS
# i    j    pot   re/A     ka
Pd1   Pd2   cons  2.9440   0.0
Pd1   Pd3   cons  2.7476   0.0
Pd2   Pd3   cons  2.7476   0.0

ANGLES
# i    j    k    pot   th/deg  ka
Pd1   Pd2   Pd3  harm  57.61   0.0
Pd2   Pd1   Pd3  harm  57.61   0.0
Pd1   Pd3   Pd2  harm  64.79   0.0

Unit conversion — epsilon:

\[\varepsilon = \frac{15.0\ \text{kJ/mol}}{4.184} = 3.58509\ \text{kcal/mol}\]

Geometry check: 57.61 + 57.61 + 64.79 = 180.01° ✓

2.2 OPC Water — Optimal Point Charge Model

Source: Izadi et al., J. Phys. Chem. Lett. 2014, 5(21), 3863–3871.

OPC is a 4-site rigid water model extending TIP4P. The virtual M-site carries the negative charge at distance d_OM = 0.1594 Å from Ow along the H-O-H bisector.

# opc.ff -- fftool format
ATOMS
# type   element  mass      charge   pot  sigma/A   epsilon/[kcal/mol]
Ow       O        15.9990   0.0000   lj   3.16655   0.21280
Hw       H        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

Unit conversion — epsilon:

\[\varepsilon = \frac{0.890359\ \text{kJ/mol}}{4.184} = 0.21280\ \text{kcal/mol}\]

Charge: q(Ow) = -2 × 0.6791 = -1.3582 e

Critical note: fftool writes q(Ow) = 0.0. The negative charge must be applied manually after reading the data file.

2.3 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 zero (Hw has no LJ site). The Ow–Hw cross term must also be explicitly set to zero — arithmetic mixing auto-generates ε = 0.10640 kcal/mol which is physically wrong for TIP4P-family models.

---

3. System Preparation with fftool + Packmol

3.1 Input File Formats

fftool accepts .zmat (Z-matrix / internal coordinate) format rather than plain .xyz. The Z-matrix encodes connectivity and atom ordering that fftool needs to write the LAMMPS topology.

pd3.zmat     -- Pd3 cluster Z-matrix (defines Pd1, Pd2, Pd3 connectivity)
opc.zmat     -- OPC water Z-matrix (defines Ow, Hw bond/angle topology)

3.2 Step-by-Step fftool + Packmol Workflow

Step 1 — Generate Packmol input and force field files

fftool 1 pd3.zmat 7205 opc.zmat -b 60.0 -u r -m a

Flag	Meaning
1 pd3.zmat	1 copy of the Pd3 cluster
7205 opc.zmat	7205 OPC water molecules
-b 60.0	Cubic box of side 60.0 Å
-u r	Real units (kcal/mol, Å). Without this, fftool defaults to reduced units and miscalculates LJ parameters
-m a	Arithmetic Lorentz-Berthelot mixing (σ_ij = (σ_i+σ_j)/2). Without this, fftool silently falls back to geometric mixing, breaking INTERFACE FF cross-term compatibility

Outputs: pack.inp (Packmol input), field.lmp (force field reference)

Step 2 — Pack molecules with Packmol

packmol < pack.inp

Using the custom pack.inp below — do not use the auto-generated one if you need precise cluster centering:

# pack.inp -- Pd3 cluster centred in OPC water box
tolerance 2.0
filetype xyz
output simbox.xyz

# Structure 1: Pd3 cluster fixed at box centre
structure pd3_pack.xyz
  number 1
  center
  fixed 30.0000 30.0000 30.0000 0.0 0.0 0.0
end structure

# Structure 2: 7205 OPC water molecules
# 1.5 A PBC padding: fills 1.5 to 58.5 A in each direction
# Exclusion sphere of 5.5 A around cluster centre prevents overlap
structure opc_pack.xyz
  number 7205
  inside box 1.5000 1.5000 1.5000 58.5000 58.5000 58.5000
  outside sphere 30.0000 30.0000 30.0000 5.5
end structure

Why these Packmol parameters:

• tolerance 2.0 — minimum 2.0 Å between any two atoms from different molecules (avoids bad LJ contacts)
• fixed 30.0 30.0 30.0 0.0 0.0 0.0 — pins the Pd3 centre-of-mass to the exact box centre with zero rotation
• inside box 1.5 1.5 1.5 58.5 58.5 58.5 — the 1.5 Å inset from each face creates a periodic image buffer; water molecules placed within 1.5 Å of the box face would have their images overlapping on the other side
• outside sphere 30.0 30.0 30.0 5.5 — 5.5 Å exclusion radius around the cluster; prevents water placement inside the first Pd solvation shell during packing

Output: simbox.xyz — all 21,618 atom positions

Step 3 — Generate LAMMPS data file

fftool 1 pd3.zmat 7205 opc.zmat -b 60.0 -u r -m a -l

The -l flag tells fftool to read simbox.xyz and write:

Output	Contents
data.lmp	Full LAMMPS data file (atoms, bonds, angles, coefficients)
in.lmp	Suggested LAMMPS input script (use as reference only)

All four flags must be identical between Step 1 and Step 3. The -m a flag in particular must be present in both calls — Step 1 uses it to compute cross-term mixing rules into field.lmp, and Step 3 uses it to write consistent pair_coeff lines into data.lmp. Omitting -m a in Step 3 silently generates geometric cross-terms that conflict with the arithmetic terms used in Step 1.

3.2 LAMMPS Data File Structure

LAMMPS data file

21618 atoms
5 atom types        # Pd1, Pd2, Pd3, Ow, Hw
14413 bonds
4 bond types        # Pd1-Pd2, Pd1-Pd3, Pd2-Pd3, Ow-Hw
7208 angles
4 angle types       # 3 Pd triangle + H-O-H

0 60 xlo xhi
0 60 ylo yhi
0 60 zlo zhi

Masses
1 106.42    # Pd1
2 106.42    # Pd2
3 106.42    # Pd3
4 15.999    # Ow
5 1.008     # Hw

Bond Coeffs # harmonic
1 0 2.944    # Pd1-Pd2 (constrained, k=0)
2 0 2.7476   # Pd1-Pd3
3 0 2.7476   # Pd2-Pd3
4 0 0.8724   # Ow-Hw

Angle Coeffs # harmonic
1 0 57.61    # Pd triangle angles
2 0 57.61
3 0 64.79
4 0 103.6    # H-O-H

---

4. LAMMPS Simulation

4.1 The Mandatory Parser Order

LAMMPS is a single-pass parser. Every style must be registered before the parser hits its coefficient section inside the data file. Violating this order produces hard errors.

bond_style / angle_style   BEFORE read_data  (Bond/Angle Coeffs)
pair_style                 BEFORE read_data  (Pair Coeffs)
kspace_style               BEFORE read_data  (safe)
read_data                  defines box + atom types
pair_coeff                 AFTER  read_data  (needs atom types)
set type 4 charge          AFTER  read_data

Violations produce:

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

4.2 Why TIP4P Hydrogen Atoms Need Special Treatment

Hw has zero LJ (ε = 0, σ = 0). With k = 0 bonds and special_bonds zeroing the 1-2 intramolecular Coulomb:

• Hw feels only intermolecular Coulomb attraction to neighbouring Ow
• No LJ repulsion, no bond spring = 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: Temporary harmonic O-H springs during minimisation only.

\[k_\text{spring} = 1000\ \text{kcal/mol/Å}^2\]

This overcomes the intermolecular Coulomb attraction (F_\text{Coulomb} \approx 77\ \text{kcal/mol/Å} at 2 Å) while being 2000× softer than SHAKE springs (k_\text{SHAKE} = 1{,}987{,}206). Springs are reset to k = 0 before write_data.

4.3 Complete LAMMPS 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 (mandatory)
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

# Charge correction (fftool writes q(Ow)=0)
set type 4 charge -1.3582

# Temporary springs: prevent Hw collapse, allow Fnorm convergence
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 k=1000

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 k=100

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

# fix setforce: freeze Pd cluster during minimisation
# (fix rigid/small conflicts with the 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

Run:

mpirun -np 8 lmp -in in.minimization -log minimize.log

Convergence result:

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

4.4 LAMMPS NPT Equilibration Script

# ==============================================================
# LAMMPS NPT EQUILIBRATION
# Phase 1: NVT 200 ps (heat from 0 K)
# Phase 2: NPT 1 ns  (equilibrate density)
# ==============================================================

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 (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

# Minimisation left all velocities = 0; re-initialise
velocity  all create 298.15 4928459 rot yes dist gaussian

# Phase 1: NVT 200 ps
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 1 ns
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

Run:

mpirun -np 40 lmp -in in.npt_equil.lammps -log npt_equil.log

4.5 SHAKE Syntax — Critical Gotchas

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

# WRONG 2 -- group name is case-sensitive:
fix SHAKE water shake 1.0e-4 20 0 m 15.999 1.008
# ERROR: Could not find fix group ID water

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

Expected log when correct:

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

---

5. GROMACS Simulation

5.1 Converting LAMMPS Data to GROMACS Format

GROMACS needs .gro (coordinates) and .top (topology). The conversion requires three steps because GROMACS makes the OPC M-site an explicit atom while LAMMPS handles it implicitly.

# convert_lmp_to_gmx.py -- reads minimized.data, writes system files
# Usage: python3 convert_lmp_to_gmx.py minimized.data

import sys
import numpy as np

def parse_lammps_data(fname):
    atoms, bonds, angles = [], [], []
    masses, box = {}, {}
    section = None
    with open(fname) as f:
        for line in f:
            line = line.split('#')[0].strip()
            if not line:
                continue
            if 'xlo xhi' in line:
                lo,hi = float(line.split()[0]), float(line.split()[1])
                box['x'] = hi - lo
            elif 'ylo yhi' in line:
                lo,hi = float(line.split()[0]), float(line.split()[1])
                box['y'] = hi - lo
            elif 'zlo zhi' in line:
                lo,hi = float(line.split()[0]), float(line.split()[1])
                box['z'] = hi - lo
            if line in ('Masses','Atoms','Bonds','Angles',
                        'Velocities','Bond Coeffs','Angle Coeffs'):
                section = line; continue
            if section == 'Masses':
                p = line.split()
                if len(p) >= 2:
                    try: masses[int(p[0])] = float(p[1])
                    except: pass
            elif section == 'Atoms':
                p = line.split()
                if len(p) >= 7:
                    atoms.append({'id':int(p[0]),'mol_id':int(p[1]),
                        'type':int(p[2]),'charge':float(p[3]),
                        'x':float(p[4]),'y':float(p[5]),'z':float(p[6])})
    atoms.sort(key=lambda a: a['id'])
    return atoms, box

def assign_names(atoms):
    TYPE_NAME = {1:'Pd1',2:'Pd2',3:'Pd3',4:'OW',5:'HW'}
    TYPE_RES  = {1:'PD3',2:'PD3',3:'PD3',4:'SOL',5:'SOL'}
    for a in atoms:
        a['name']    = TYPE_NAME.get(a['type'],f"X{a['type']}")
        a['resname'] = TYPE_RES.get(a['type'],'UNK')
    return atoms

def write_gro(atoms, box, outfile='system.gro'):
    ANG2NM = 0.1
    resmap, resnum, prev = {}, 0, None
    for a in atoms:
        if a['mol_id'] != prev:
            resnum += 1; resmap[a['mol_id']] = resnum; prev = a['mol_id']
    with open(outfile,'w') as f:
        f.write('Pd3 + OPC water minimized from LAMMPS\n')
        f.write(f'{len(atoms):5d}\n')
        for a in atoms:
            rn = resmap[a['mol_id']] % 100000
            an = a['id'] % 100000
            f.write(f"{rn:5d}{a['resname']:<5s}{a['name']:>5s}{an:5d}"
                    f"{a['x']*ANG2NM:8.3f}{a['y']*ANG2NM:8.3f}"
                    f"{a['z']*ANG2NM:8.3f}\n")
        f.write(f"  {box['x']*ANG2NM:.5f}"
                f"  {box['y']*ANG2NM:.5f}"
                f"  {box['z']*ANG2NM:.5f}\n")

atoms, box = parse_lammps_data(sys.argv[1])
atoms = assign_names(atoms)
write_gro(atoms, box)

5.2 Adding the OPC M-site (MW virtual atom)

GROMACS needs explicit MW coordinates:

# add_mw_sites.py
# Usage: python3 add_mw_sites.py system.gro system_opc4.gro
import sys, numpy as np
A_VSITE = 0.14806  # d_OM=0.1594A derived weight

def parse_gro(f):
    lines = open(f).readlines()
    return lines[0], int(lines[1]), lines[2:2+int(lines[1])], lines[2+int(lines[1])]

def parse_line(l):
    return (int(l[0:5]), l[5:10].strip(), l[10:15].strip(),
            int(l[15:20]), float(l[20:28]), float(l[28:36]), float(l[36:44]))

def fmt_line(rn,resname,atom,an,x,y,z):
    return f"{rn%100000:5d}{resname:<5s}{atom:>5s}{an%100000:5d}{x:8.3f}{y:8.3f}{z:8.3f}\n"

title, natoms, atomlines, boxline = parse_gro(sys.argv[1])
parsed = [parse_line(l) for l in atomlines]
out, new_an = [], 0
i = 0
while i < len(parsed):
    rn,resname,atom,an,x,y,z = parsed[i]
    if resname == 'SOL' and atom == 'OW':
        ow = np.array([x,y,z]); i+=1
        _,_,_,_,hx,hy,hz = parsed[i]; hw1 = np.array([hx,hy,hz]); i+=1
        _,_,_,_,hx,hy,hz = parsed[i]; hw2 = np.array([hx,hy,hz])
        mw = ow + A_VSITE*(hw1-ow) + A_VSITE*(hw2-ow)
        for name,pos in [('OW',ow),('HW1',hw1),('HW2',hw2),('MW',mw)]:
            new_an += 1
            out.append(fmt_line(rn,'SOL',name,new_an,*pos))
    else:
        new_an += 1
        out.append(fmt_line(rn,resname,atom,new_an,x,y,z))
    i += 1
with open(sys.argv[2],'w') as f:
    f.write(title); f.write(f'{len(out):5d}\n')
    f.writelines(out); f.write(boxline)
print(f"Input: {natoms} atoms  Output: {len(out)} atoms  MW added: {len(out)-natoms}")

5.3 GROMACS Topology Files

pd3.itp:

; pd3.itp -- Pd3 cluster (INTERFACE FF, Heinz 2008)
[ moleculetype ]
PD3     1

[ atoms ]
;nr  type  resnr  res   atom  cgnr  charge    mass
  1  Pd1    1    PD3   Pd1    1    0.0000  106.4200
  2  Pd2    1    PD3   Pd2    1    0.0000  106.4200
  3  Pd3    1    PD3   Pd3    1    0.0000  106.4200

[ constraints ]
;ai  aj  funct  b0(nm)
  1   2    1    0.29440    ; Pd1-Pd2: 2.9440 A
  1   3    1    0.27476    ; Pd1-Pd3: 2.7476 A
  2   3    1    0.27476    ; Pd2-Pd3: 2.7476 A

opc.itp:

; opc.itp -- OPC 4-site water (Izadi 2014)
[ moleculetype ]
SOL     2

[ atoms ]
;nr  type  resnr  res   atom  cgnr  charge    mass
  1  OW     1    SOL    OW    1    0.00000  15.9990
  2  HW     1    SOL   HW1    1    0.67910   1.0080
  3  HW     1    SOL   HW2    1    0.67910   1.0080
  4  MW     1    SOL    MW    1   -1.35820   0.0000

[ constraints ]
;ai  aj  funct  b0(nm)
  1   2    1   0.08724    ; O-H = 0.8724 A
  1   3    1   0.08724
  2   3    1   0.13712    ; H-H derived from geometry

[ virtual_sites3 ]
; MW placed at d_OM=0.1594 A along bisector
; weight a = 0.1594/(2*0.8724*cos(51.8 deg)) = 0.14806
  4   1   2   3    1   0.14806  0.14806

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

system.top:

[ defaults ]
; nbfunc  comb-rule  gen-pairs  fudgeLJ  fudgeQQ
    1         2         yes        1.0      1.0

[ atomtypes ]
; name  at_num  mass      charge  ptype  sigma(nm)   epsilon(kJ/mol)
Pd1     46    106.4200   0.0000    A    0.24800    15.00000
Pd2     46    106.4200   0.0000    A    0.24800    15.00000
Pd3     46    106.4200   0.0000    A    0.24800    15.00000
OW       8     15.9990  -1.3582    A    0.31666     0.89036
HW       1      1.0080   0.6791    A    0.00000     0.00000
MW       0      0.0000   0.0000    D    0.00000     0.00000

[ nonbond_params ]
; ai   aj  funct  sigma(nm)    epsilon(kJ/mol)
; Pd-Ow: sigma=(0.248+0.31666)/2=0.28233, eps=sqrt(15.0*0.89036)=3.6548
Pd1   OW    1    0.28233    3.65480
Pd2   OW    1    0.28233    3.65480
Pd3   OW    1    0.28233    3.65480
Pd1   HW    1    0.00000    0.00000
Pd2   HW    1    0.00000    0.00000
Pd3   HW    1    0.00000    0.00000

#include "pd3.itp"
#include "opc.itp"

[ system ]
Pd3 cluster in OPC water

[ molecules ]
PD3    1
SOL    7205

5.4 GROMACS MDP Files

minim.mdp — Energy Minimisation:

integrator      = steep
emtol           = 100.0        ; Fmax < 100 kJ/mol/nm = 1.44 kcal/mol/A
emstep          = 0.01
nsteps          = 50000
nstxout         = 500
nstlog          = 100
nstenergy       = 100
cutoff-scheme   = Verlet
nstlist         = 10
coulombtype     = PME
rcoulomb        = 1.2          ; 12 A (matches LAMMPS)
pme-order       = 4
fourierspacing  = 0.12
ewald-rtol      = 1e-5
vdwtype         = Cut-off
rvdw            = 1.2
DispCorr        = EnerPres     ; = pair_modify tail yes in LAMMPS
constraint-algorithm = LINCS
constraints          = none
gen-vel         = no
pbc             = xyz

nvt.mdp — NVT 200 ps:

integrator      = md
nsteps          = 100000       ; 100,000 x 2 fs = 200 ps
dt              = 0.002        ; 2 fs (vs 1 fs in LAMMPS)
nstxout-compressed = 5000
nstenergy       = 500
nstlog          = 500
cutoff-scheme   = Verlet
nstlist         = 10
coulombtype     = PME
rcoulomb        = 1.2
pme-order       = 4
fourierspacing  = 0.12
ewald-rtol      = 1e-5
vdwtype         = Cut-off
rvdw            = 1.2
DispCorr        = EnerPres
tcoupl          = V-rescale
tc-grps         = PD3 SOL
tau-t           = 0.1   0.1
ref-t           = 298.15 298.15
pcoupl          = no
constraint-algorithm = LINCS
constraints          = h-bonds
lincs-iter      = 1
lincs-order     = 4
gen-vel         = yes
gen-temp        = 298.15
gen-seed        = 4928459
pbc             = xyz

npt.mdp — NPT equilibration + production:

integrator      = md
nsteps          = 500000       ; 1 ns equil (change to 50000000 for 100 ns)
dt              = 0.002
nstxout-compressed = 5000
nstenergy       = 500
nstlog          = 500
cutoff-scheme   = Verlet
nstlist         = 10
coulombtype     = PME
rcoulomb        = 1.2
pme-order       = 4
fourierspacing  = 0.12
ewald-rtol      = 1e-5
vdwtype         = Cut-off
rvdw            = 1.2
DispCorr        = EnerPres
tcoupl          = Nose-Hoover
tc-grps         = PD3 SOL
tau-t           = 0.5   0.5
ref-t           = 298.15 298.15
pcoupl              = Parrinello-Rahman
pcoupltype          = isotropic
ref-p               = 1.0
tau-p               = 2.0
compressibility     = 4.5e-5
constraint-algorithm = LINCS
constraints          = h-bonds
lincs-iter      = 1
lincs-order     = 4
gen-vel         = no
pbc             = xyz

5.5 Complete GROMACS Workflow

# 0. Environment
module load apps/gromacs/2022/cpu
source /home/apps/gromacs/gromacs-2022.2/installCPUIMPI/bin/GMXRC.bash
GMX="gmx_CPUIMPI"

# 1. Convert LAMMPS data to GROMACS format
python3 convert_lmp_to_gmx.py minimized.data    # writes system.gro
python3 add_mw_sites.py system.gro system_opc4.gro  # adds MW virtual sites

# 2. Minimisation
${GMX} grompp -f minim.mdp -c system_opc4.gro -p system.top \
              -o minim.tpr -maxwarn 2
mpirun -np 8 ${GMX} mdrun -deffnm minim -ntomp 5 -pin on

# 3. NVT 200 ps
${GMX} grompp -f nvt.mdp -c minim.gro -p system.top \
              -o nvt.tpr -maxwarn 2
mpirun -np 8 ${GMX} mdrun -deffnm nvt -ntomp 5 -pin on

# 4. NPT equilibration 1 ns
${GMX} grompp -f npt.mdp -c nvt.gro -t nvt.cpt \
              -p system.top -o npt_equil.tpr -maxwarn 2
mpirun -np 8 ${GMX} mdrun -deffnm npt_equil -ntomp 5 -pin on

# 5. Check equilibration
${GMX} energy -f npt_equil.edr << 'EOF'
Temperature
Pressure
Density
0
EOF

# 6. Production 100 ns (change nsteps in npt.mdp to 50000000)
sed 's/nsteps.*=.*500000/nsteps          = 50000000/' npt.mdp > prod.mdp
${GMX} grompp -f prod.mdp -c npt_equil.gro -t npt_equil.cpt \
              -p system.top -o prod.tpr -maxwarn 2
mpirun -np 8 ${GMX} mdrun -deffnm prod -ntomp 5 -pin on -cpt 60

---

6. LAMMPS vs GROMACS — Detailed Comparison

6.1 Unit Systems

Quantity	LAMMPS real	GROMACS	Conversion
Energy	kcal/mol	kJ/mol	× 4.184
Length	Å	nm	÷ 10
Time	fs	ps	÷ 1000
Mass	AMU	AMU	1:1
Temperature	K	K	1:1
Pressure	atm	bar	× 1.01325
Force	kcal/mol/Å	kJ/mol/nm	× 4.184 × 10 = × 41.84
Velocity	Å/fs	nm/ps	× 0.1
Charge	e	e	1:1

Example conversions from this simulation:

LAMMPS Fmax    =  1.5 kcal/mol/A  x 41.84  = 62.8 kJ/mol/nm
GROMACS Fmax   = 60.5 kJ/mol/nm   / 41.84  =  1.4 kcal/mol/A
-> Same physical force OK

LAMMPS PotEng  = -106,300 kcal/mol x 4.184 = -444,759 kJ/mol
GROMACS Epot   = -432,261 kJ/mol
-> Difference = 12,498 kJ/mol from temporary springs in LAMMPS OK

6.2 Electrostatics

Feature	LAMMPS	GROMACS
Method	pppm/tip4p	PME
Grid	50×50×50 (KISS FFT / FFTW3)	52×52×52 (MKL)
Accuracy	1.0e-5	ewald-rtol = 1e-5
M-site handling	Implicit (computed per step)	Explicit virtual_sites3 atom
Real-space cutoff	12.0 Å	1.2 nm (= 12 Å)

6.3 Water Constraints

Feature	LAMMPS	GROMACS
Algorithm	SHAKE (iterative)	SETTLE (analytical exact)
During minimisation	Replaced by k=1,987,206 springs	LINCS/SETTLE work natively
Minimisation steps needed	~71,000	51
O-H bond	shake m 15.999 1.008	constraints = h-bonds
H-O-H angle	Enforced by SHAKE cluster	Enforced by SETTLE

6.4 Rigid Pd3 Cluster

Feature	LAMMPS	GROMACS
During minimisation	fix setforce 0.0 0.0 0.0	No special fix needed (LINCS)
During dynamics	fix rigid/small molecule	[constraints] in pd3.itp
Thermostat	Implicit (water collisions)	Separate tc-grps = PD3 SOL
Time integration	rigid/small does own NVE	Standard leap-frog via LINCS

6.5 Thermostat and Barostat

Feature	LAMMPS	GROMACS
NVT thermostat	fix nvt ... 100.0 (Nose-Hoover)	V-rescale (heating phase)
NPT thermostat	fix npt ... 100.0	Nose-Hoover, tau-t=0.5 ps
NPT barostat	fix npt iso 1.0 1.0 1000.0	Parrinello-Rahman, tau-p=2 ps
Tdamp	100 fs (LAMMPS)	500 fs (GROMACS, 0.5 ps)
Pdamp	1000 fs (LAMMPS)	2000 fs (GROMACS, 2.0 ps)
Ensemble	Correct NpT ✓	Correct NpT ✓

Note: GROMACS uses longer damping constants because the 2 fs timestep vs 1 fs in LAMMPS requires proportionally larger damping for stability. Both are physically equivalent in the long-time limit.

6.6 Timestep

	LAMMPS	GROMACS
Timestep	1 fs	2 fs
Why	SHAKE iterative (conservative)	LINCS/SETTLE analytical (safe at 2 fs)
Impact	Baseline throughput	2× throughput before any other change

6.7 Output Formats

Data	LAMMPS	GROMACS
Trajectory	.dcd (binary, ~4 bytes/atom/frame)	.xtc (compressed, ~2 bytes/atom/frame)
Energy	thermo in log	.edr binary (analysed with gmx energy)
Restart	.restart binary	.cpt checkpoint binary
Data/topology	.data	.gro + .top
Log	log.lammps	.log

---

7. Performance Benchmarks — PARAM Shakti (IIT Kharagpur)

System: 21,618 atoms (LAMMPS) / 28,823 atoms (GROMACS) Module: apps/lammps/29Aug2024/oneapi-2021 / apps/gromacs/2022/cpu Hardware: Intel Skylake, 40 cores/node, InfiniBand

7.1 LAMMPS Benchmarks (pure MPI, 1 fs timestep)

Cores	Nodes	Partition	ns/day	Scaling
8	shared	shared	1.44	baseline
16	shared	shared	1.46	~1.0× (zero scaling!)
36	shared	shared	5.80	4.0×
40	1	medium	7.08	4.9×

Key finding: 8→16 cores on the shared partition gave essentially zero speedup because LAMMPS MPI tasks landed on different NUMA nodes without InfiniBand. Always use the medium partition.

LAMMPS wall time breakdown at 40 cores:

Pair    : 64%  (LJ + real-space Coulomb)
Kspace  : 18%  (PPPM FFT -- FFTW3 is 2x faster than KISS FFT)
Neigh   : 12%  (neighbour list rebuilds)
Comm    :  3%

Note: OPENMP and INTEL packages are not compiled into any LAMMPS build on PARAM Shakti. Only pure MPI is available.

7.2 GROMACS Benchmark (MPI+OMP hybrid, 2 fs timestep)

Cores	Config	ns/day
40	8 MPI × 5 OMP	100.2

14× faster than LAMMPS on identical hardware.

7.3 Source of the 14× Speedup

Factor	Speedup
2 fs vs 1 fs timestep	2.0×
SETTLE vs SHAKE (no iterations)	~2.0×
Native OpenMP (LAMMPS had none)	~2.0×
PME vs PPPM (MKL-optimised FFT)	~1.5×
Combined	~12–14×

7.4 Practical Impact on 100 ns Simulation

Engine	ns/day	Days for 100 ns	Nodes needed	Queue strategy
LAMMPS (1 node)	7	14 days	4–5 chained jobs	medium, multiple
LAMMPS (10 nodes)	~40	2.5 days	10 nodes	medium, Priority wait
GROMACS (1 node)	100	1 day	1 node	medium, instant

---

8. HPC Setup — PARAM Shakti

8.1 Partition Rules

Partition	Cores	Walltime	Policy
shared	1–36	3 days	Node sharing
medium	1–10 nodes × 40	3 days	Exclusive
large	1–10 nodes × 40	7 days	Exclusive, low priority

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

8.2 LAMMPS Slurm Script (medium, 1 node)

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

module load apps/lammps/29Aug2024/oneapi-2021
cd /scratch/20ch91r15/simulations/pd3_h2o/02_npt_equil
mpirun -np 40 lmp -in in.npt_equil.lammps -log npt.log

8.3 GROMACS Slurm Script (medium, 1 node, 8×5 hybrid)

#!/bin/bash
#SBATCH --job-name=gmx_prod
#SBATCH --partition=medium
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=8
#SBATCH --cpus-per-task=5
#SBATCH --time=03-00:00:00
#SBATCH --output=gmx_%j.out

module load apps/gromacs/2022/cpu
source /home/apps/gromacs/gromacs-2022.2/installCPUIMPI/bin/GMXRC.bash

export OMP_NUM_THREADS=5
export OMP_PROC_BIND=close
export OMP_PLACES=cores

cd /scratch/20ch91r15/simulations/pd3_h2o/02_gmx_equil
mpirun -np 8 gmx_CPUIMPI mdrun -deffnm prod -ntomp 5 -pin on -cpt 60

Note: gmx_CPUIMPI is compiled with real Intel MPI, not thread-MPI. Use mpirun -np N for ranks. Do NOT use -ntmpi N — that flag only works with thread-MPI builds.

---

9. Verification — No VMD Required

9.1 LAMMPS Health Check

# Convergence criterion
grep "Stopping criterion" log.lammps
# Must show "force tolerance", NOT "max iterations"

# Force magnitudes
grep "Force two-norm\|Force max component" log.lammps

# 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

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

9.2 GROMACS Health Check

# Minimisation converged
grep "converged to Fmax" minim.log

# NPT equilibration properties
gmx_CPUIMPI energy -f npt_equil.edr << 'EOF'
Temperature
Pressure
Density
0
EOF
# Targets: T=298.15+/-5 K, P=1+/-200 bar, density~0.997 g/cm3

# Performance
grep "Performance" npt_equil.log

---

10. Lessons Learned

Critical (will crash the run)

1. LAMMPS parser order is strict. pair_coeff after read_data; pair_style, bond_style, angle_style before read_data.

2. Hw has zero LJ — never remove its O-H restraint. Without a spring or constraint, Hw collapses into neighbouring Ow.

3. SHAKE during minimisation creates an Fnorm floor of ~2,100 (k=1,987,206 springs × 14,410 bonds). Use k=1000 harmonic bonds instead.

4. fix rigid/small conflicts with the LAMMPS minimiser. Use fix setforce 0.0 0.0 0.0 on the cluster during minimisation.

5. LAMMPS group names are case-sensitive. WATER ≠ water.

6. Ow charge must be set manually. fftool writes q(Ow) = 0.0. Apply set type 4 charge -1.3582 after read_data.

7. GROMACS [ defaults ] must be the first section in system.top. Without it, [ atomtypes ] produces Invalid order for directive.

8. MW virtual site must be explicit in system.gro for GROMACS. LAMMPS handles it implicitly; GROMACS needs a 4th atom per water. Atom count mismatch: 21,618 (LAMMPS) vs 28,823 (GROMACS).

9. gmx_CPUIMPI uses real MPI, not thread-MPI. Use mpirun -np N for rank count. -ntmpi N throws Fatal error: Setting number of thread-MPI ranks... GROMACS compiled without thread-MPI.

10. etol = 0 in LAMMPS minimise. Non-zero etol triggers false early exit (caused 8-step exit with Fnorm=2,076 in early runs).

Moderate (wrong physics)

1. Ow–Hw LJ cross term is not zero by default. Arithmetic mixing generates ε = 0.10640 kcal/mol. Set pair_coeff 4 5 0.0 0.0.

2. Berendsen barostat ≠ correct NpT ensemble. Use fix npt (LAMMPS) or Parrinello-Rahman (GROMACS).

3. SHAKE b 4 a 4 builds frozen angles, not size-3 clusters. Use mass-based m 15.999 1.008 for TIP4P water.

4. comm_modify cutoff 15.1 silences paired TIP4P comm-cutoff warnings in LAMMPS.

---

11. References

1. Heinz, H., Vaia, R. A., Farmer, B. L., & Naik, R. R. (2008). Accurate Simulation of Surfaces and Interfaces of FCC Metals. J. Phys. Chem. C, 112(44), 17281–17290.

2. Izadi, S., Anandakrishnan, R., & Onufriev, A. V. (2014). Building Water Models: A Different Approach. J. Phys. Chem. Lett., 5(21), 3863–3871.

3. Thompson, A. P., et al. (2022). LAMMPS — A flexible simulation tool. Comput. Phys. Commun., 271, 108171.

4. Martyna, G. J., Tobias, D. J., & Klein, M. L. (1994). Constant pressure molecular dynamics algorithms. J. Chem. Phys., 101(5), 4177–4189.

5. LAMMPS Documentation (29 Aug 2024). https://docs.lammps.org

6. GROMACS Manual (2022.2). https://manual.gromacs.org

7. PARAM Shakti HPC, IIT Kharagpur. https://hpc.iitkgp.ac.in/HPCF/paramShakti