Building pd4.ff — A Force Field File for the Pd₄ Cluster in fftool

This post documents every decision made in constructing pd4.ff, the fftool-compatible force field file for a zero-valent palladium tetramer (Pd₄) cluster. It is a companion to the pd3-ff post and follows the same structure. The LJ parameters are identical between the two files — the differences lie entirely in the geometry: Pd₄ is a tetrahedron, not a triangle, and that introduces one non-obvious angle convention that must be handled carefully. —

What is pd4.ff?

fftool reads a plain-text .ff file defining non-bonded (Lennard-Jones), bonded (stretch, bend), and charge parameters for each atom type in a molecule. The format is documented in the fftool README:

“potential parameters, namely Lennard-Jones sigma (angstrom) and epsilon (kJ mol⁻¹)”

The energy unit is kJ/mol throughout — for epsilon, bond force constants, and angle force constants alike.

The Final File

# INTERFACE Force Field for Pd4 cluster
#
# LJ PARAMETERS
# Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R.
# J. Phys. Chem. C 2008, 112, 17281-17290. DOI: 10.1021/jp801931d
# Table 1, Pd 12-6 LJ: r0 = 2.819 A, eps0 = 6.15 kcal/mol
# Conversion: sigma = 2.5114 A, epsilon = 25.7316 kJ/mol
#
# CLUSTER GEOMETRY
# Locally optimised via DFT (QM) calculation on the bare Pd4 cluster.
# Geometry: Td regular tetrahedron
#   Z-matrix bonds: 2.6759, 2.6754 (primary), 2.6751 A
#   Primary equilibrium bond distance: 2.6754 A
#   3D Pd-Pd-Pd angle: 109.4712 deg = arccos(-1/3)

ATOMS
# type  m/amu     q/e     pot  sigma/A  eps/kJmol-1
Pd1  Pd1  106.42  0.0000  lj   2.5114   25.7316
Pd2  Pd2  106.42  0.0000  lj   2.5114   25.7316
Pd3  Pd3  106.42  0.0000  lj   2.5114   25.7316
Pd4  Pd4  106.42  0.0000  lj   2.5114   25.7316

BONDS
# i    j    pot   re/A    kr/kJmol-1A-2
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

ANGLES
# i    j    k    pot   th/deg    ka/kJmol-1rad-2
Pd1  Pd2  Pd3  harm  109.4712  0.0000
Pd1  Pd2  Pd4  harm  109.4712  0.0000
Pd3  Pd2  Pd4  harm  109.4712  0.0000
Pd1  Pd3  Pd2  harm  109.4712  0.0000
Pd1  Pd3  Pd4  harm  109.4712  0.0000
Pd2  Pd3  Pd4  harm  109.4712  0.0000
Pd1  Pd4  Pd2  harm  109.4712  0.0000
Pd1  Pd4  Pd3  harm  109.4712  0.0000
Pd2  Pd4  Pd3  harm  109.4712  0.0000
Pd2  Pd1  Pd3  harm  109.4712  0.0000
Pd2  Pd1  Pd4  harm  109.4712  0.0000
Pd3  Pd1  Pd4  harm  109.4712  0.0000

Section 1 — ATOMS

1.1 Mass

Palladium has a standard atomic weight of 106.42 AMU (IUPAC 2021). Written directly into the mass column for all four atoms, no conversion needed.

1.2 Charge

The Pd₄ cluster is zero-valent (Pd⁰). No partial charge model is applied: q = 0.0000 e for all four atoms.

1.3 Lennard-Jones Parameters — Source and Conversion

Source: Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R. J. Phys. Chem. C 2008, 112, 17281–17290. DOI: 10.1021/jp801931d · ACS page

This is the same source used for pd3.ff. The INTERFACE force field fits 12-6 Lennard-Jones parameters for fcc metals against bulk lattice constants, surface energies, and elastic moduli. For palladium, Table 1, 12-6 LJ column:

\[r_0 = 2.819\ \text{Å}, \quad \varepsilon_0 = 6.15\ \text{kcal/mol}\]

The functional form mismatch — why sigma ≠ r0:

Heinz Table 1 uses the AMBER/CHARMM form where $r_0$ is the distance at the potential minimum:

\[E(r) = \varepsilon_0 \left[\left(\frac{r_0}{r}\right)^{12} - 2\left(\frac{r_0}{r}\right)^{6}\right]\]

fftool uses the standard 4-epsilon form where $\sigma$ is the zero-crossing distance (confirmed by the fftool README):

\[E(r) = 4\varepsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]\]

The minimum of the fftool form sits at $r_\text{min} = 2^{1/6}\sigma$. Equating the two minima:

\[\sigma = \frac{r_0}{2^{1/6}} = \frac{2.819}{1.12246} = 2.5114\ \text{Å}\]

Converting the well depth to kJ/mol:

\[\varepsilon = 6.15\ \text{kcal/mol} \times 4.184 = 25.7316\ \text{kJ/mol}\]

Cross-check against Table 1 columns A and B (form $E = A/r^{12} - B/r^6$):

\[A = \varepsilon_0 \cdot r_0^{12} = 6.15 \times 2.819^{12} = 1{,}548{,}874 \approx 1{,}549{,}000\ \checkmark\] \[B = 2\varepsilon_0 \cdot r_0^{6} = 2 \times 6.15 \times 2.819^{6} = 6{,}173\ \checkmark\]

Both match the table entries exactly, confirming the raw values were read correctly before conversion.

Section 2 — BONDS

Geometry: Locally Optimised QM Calculation

The bond lengths used in this file come from a locally performed DFT optimisation of the bare Pd₄ cluster. These values are the most directly applicable to this simulation — obtained on the actual system being studied rather than transferred from a different chemical context.

The QM-optimised geometry for Pd₄ is a Tᵈ regular tetrahedron. The Z-matrix used to build the structure encodes the four atoms via the following internal coordinates:

Z-matrix variable	Value	Description
r12	2.6759 Å	bond Pd1–Pd2
r13	2.6754 Å	bond Pd1–Pd3 (primary value)
r14	2.6754 Å	bond Pd1–Pd4
a132	60.0000°	Z-matrix construction angle
a142	59.9900°	Z-matrix construction angle
d1423	60.0000°	Z-matrix dihedral

The three distinct bond lengths span 2.6751–2.6759 Å, a spread of only 0.0008 Å — well within the numerical noise of any DFT optimisation. The primary equilibrium distance 2.6754 Å is used for all six edges of the tetrahedron in the BONDS section, which preserves the Tᵈ symmetry of the force field.

This value sits in excellent agreement with the DFT literature for bare Pd₄ clusters. German, Efremenko, and Sheintuch (J. Phys. Chem. A 2001, 105, 11312) report the Pd₄ C₃ᵥ triplet ground state (Table 1 of that paper) with apex-to-base bonds of 2.604 Å and base-to-base bonds of 2.714 Å at the B3LYP/LANL2DZ level. The present locally optimised value of 2.6754 Å falls between these two, consistent with a near-regular tetrahedron at a slightly different level of theory.

Why all six bonds use the same `re`

A regular tetrahedron has six equivalent edges. The sub-0.001 Å variation across the three Z-matrix bond lengths is an artefact of the Z-matrix construction order — each atom is placed relative to the previously placed ones, accumulating small numerical differences. These differences carry no physical meaning. Using the primary value 2.6754 Å for all six bonds enforces Tᵈ symmetry in the force field and avoids introducing spurious symmetry breaking into the rigid-body dynamics.

Why `pot = cons` with `kr = 0.0000`

The Pd₄ cluster is treated as a rigid body throughout the simulation using LAMMPS fix rigid/small. The cons potential type tells fftool to write a constrained bond entry in the topology — establishing the connectivity and the equilibrium distance — while applying no spring force. The force constant kr = 0.0000 kJ/mol/Å² is correct and intentional: the integrator enforces every Pd–Pd distance at every timestep, not a spring.

Section 3 — ANGLES

The Z-matrix angle vs the 3D bond angle — a critical distinction

This is the most important technical point specific to Pd₄ that does not arise for Pd₃.

The Z-matrix angles a132 = 60.00° and a142 = 59.99° are internal coordinate construction variables, not physical Pd–Pd–Pd bond angles. They describe the angle between bond vectors in the sequential Z-matrix build — they have no direct geometric meaning once the 3D Cartesian coordinates are generated.

The actual Pd–Pd–Pd angle in a regular tetrahedron is:

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

This is the angle that fftool measures when it reads the 3D structure from Packmol output and tries to assign angle terms. The fftool README states a tolerance of ±15° around the equilibrium angle in the .ff file. If 60° were written in the ANGLES section, fftool would look for angles near 60° in a structure that has angles of 109.47° — a 49.47° discrepancy — and would silently drop all 12 angle terms from the topology. This is why 109.4712° must appear in the ANGLES section regardless of what the Z-matrix contains.

Why 12 angle entries

A regular tetrahedron has 4 atoms. For each central atom j, there are $\binom{3}{2} = 3$ unique i-j-k triplets (choosing 2 from the 3 remaining atoms to flank j). With 4 central atoms: $4 \times 3 = 12$ unique angle entries. All are listed explicitly so fftool can write a complete topology.

Pd1  Pd2  Pd3  harm  109.4712  0.0000
Pd1  Pd2  Pd4  harm  109.4712  0.0000
Pd3  Pd2  Pd4  harm  109.4712  0.0000
  ... (9 more, one for each atom as central j)

The angle force constant ka = 0.0000 kJ/mol/rad² is correct for the same reason as the bond force constants: fix rigid/small holds the geometry rigidly at every timestep, making any spring constant physically irrelevant.

Unit System Summary

Quantity	fftool unit	Source unit	Conversion
sigma	Å	Å	$\sigma = r_0 / 2^{1/6}$
epsilon	kJ/mol	kcal/mol	$\times\ 4.184$
bond re	Å	Å	none (QM direct)
bond kr	kJ/mol/Å²	—	0 (rigid)
angle th	degrees	degrees	arccos(−1/3), not Z-matrix
angle ka	kJ/mol/rad²	—	0 (rigid)

Parameter Table

Parameter	Value	Unit	Source
Mass (Pd)	106.42	AMU	IUPAC 2021
Charge	0.0000	e	Pd⁰ neutral cluster
sigma	2.5114	Å	Heinz 2008 Table 1, converted
epsilon	25.7316	kJ/mol	Heinz 2008 Table 1, converted
r(Pd–Pd)	2.6754	Å	Local DFT optimisation of bare Pd₄ cluster (primary)
r(Pd–Pd) range	2.6751–2.6759	Å	Full Z-matrix spread (0.0008 Å)
∠Pd–Pd–Pd (ff)	109.4712	deg	arccos(−1/3), exact Tᵈ geometry
∠Pd–Pd–Pd (Z-matrix)	60.00 / 59.99	deg	Z-matrix construction variables only

References

fftool — Padua, A. fftool: Tool to build force field input files for molecular simulation. Zenodo, 2015. DOI: 10.5281/zenodo.18618 · GitHub
LJ parameters (INTERFACE FF) — Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R. Accurate Simulation of Surfaces and Interfaces of Face-Centered Cubic Metals Using 12-6 and 9-6 Lennard-Jones Potentials. J. Phys. Chem. C 2008, 112, 17281–17290. DOI: 10.1021/jp801931d · ACS page
Pd₄ cluster geometry (literature context) — German, E. D.; Efremenko, I.; Sheintuch, M. Hydrogen Interactions with a Pd₄ Cluster: Triplet and Singlet States and Transition Probability. J. Phys. Chem. A 2001, 105, 11312–11326. DOI: 10.1021/jp012796+ · ACS page