Home Building pd3.ff — A Force Field File for the Pd₃ Cluster in fftool
Post
Cancel

Building pd3.ff — A Force Field File for the Pd₃ Cluster in fftool

This post documents every decision made in constructing pd3.ff, the fftool-compatible force field file for a zero-valent palladium trimer (Pd₃) cluster. Every number in the file is traced to a clearly identified source. The goal is a file that can be dropped directly into an fftool workflow alongside opc.ff for Pd₃/water simulations.


What is pd3.ff?

fftool reads a plain-text .ff file that defines 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 and the relevant excerpt is:

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

That single line fixes two things that are easy to get wrong: the column is sigma (the zero-crossing distance, not the potential minimum), and the energy unit is kJ/mol throughout the file — for epsilon, bond force constants, and angle force constants alike.


The Final File

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# INTERFACE Force Field for Pd3 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 (AMBER/CHARMM/CVFF/OPLS-AA form):
#   r0   = 2.819 A   (r at potential minimum)
#   eps0 = 6.15 kcal/mol (well depth)
# Conversion to fftool form:
#   sigma   = r0 / 2^(1/6) = 2.5114 A
#   epsilon = eps0 * 4.184 = 25.7316 kJ/mol
#
# CLUSTER GEOMETRY
# Locally optimised via DFT (QM) calculation on the bare Pd3 cluster.
# D3h equilateral triangle: r(Pd-Pd) = 2.6760 A, all angles = 60.0 deg

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

BONDS
# i    j    pot   re/A    kr/kJmol-1A-2
Pd1  Pd2  cons  2.6760  0.0000
Pd1  Pd3  cons  2.6760  0.0000
Pd2  Pd3  cons  2.6760  0.0000

ANGLES
# i    j    k    pot   th/deg   ka/kJmol-1rad-2
Pd1  Pd2  Pd3  harm  60.0000  0.0000
Pd2  Pd1  Pd3  harm  60.0000  0.0000
Pd1  Pd3  Pd2  harm  60.0000  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 with no conversion needed.

1.2 Charge

The Pd₃ cluster is zero-valent (Pd⁰). No ionic character and no partial charge model is applied: q = 0.0000 e for all three 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 paper presents the INTERFACE force field, which 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}\]

Why INTERFACE FF and not UFF? The Universal Force Field (UFF, Rappé 1992) assigns generic LJ parameters estimated from ionisation potentials and polarisabilities. The INTERFACE FF was explicitly fitted to condensed-phase properties of bulk Pd metal — lattice parameter, cohesive energy, surface energy — making it far more appropriate for a zero-valent metallic cluster in a solvated MD simulation.

The critical unit conversion — why sigma ≠ r0:

Heinz Table 1 uses the AMBER/CHARMM functional form:

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

where $r_0$ is the distance at the potential minimum.

fftool uses the standard 4-epsilon form (confirmed by the fftool README):

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

where $\sigma$ is the distance where E = 0 (the zero-crossing). Equating the two minima gives:

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

Converting the well depth to kJ/mol (fftool’s required energy unit):

\[\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. This makes these values the most directly applicable to this simulation — they were obtained on the actual system being studied, not transferred from a different chemical context.

The QM-optimised geometry for Pd₃ is a D₃h equilateral triangle:

\[r(\text{Pd–Pd}) = 2.6760\ \text{Å} \quad \text{(all three bonds equal)}\]

The 60.0° angles follow analytically from the three equal bond lengths — this is the defining property of an equilateral triangle and requires no separate optimisation.

This value is consistent with the range reported in the DFT literature for bare Pd₃ clusters (2.53–2.68 Å across different functionals and basis sets), and falls close to the B3LYP/LANL2DZ value of 2.57 Å reported by Moc, Musaev & Morokuma (2000).

Why pot = cons with kr = 0.0000?

The Pd₃ cluster is treated as a rigid body 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 — but applies no spring force. The force constant kr = 0.0000 kJ/mol/Ų is correct and intentional: the integrator enforces the geometry at every step, not a spring.


Section 3 — ANGLES

All three interior angles of the equilateral triangle are exactly 60.0° by definition of D₃h symmetry. All three i-j-k triplets are listed explicitly so fftool can write a complete angle topology:

1
2
3
Pd1  Pd2  Pd3  harm  60.0000  0.0000
Pd2  Pd1  Pd3  harm  60.0000  0.0000
Pd1  Pd3  Pd2  harm  60.0000  0.0000

The angle force constant ka = 0.0000 kJ/mol/rad² is correct for the same reason as the bond force constants: fix rigid/small makes all spring constants physically irrelevant.


Unit System Summary

Quantityfftool unitSource unitConversion
sigmaÅÅ$\sigma = r_0 / 2^{1/6}$
epsilonkJ/molkcal/mol$\times\ 4.184$
bond reÅÅnone (QM direct)
bond krkJ/mol/Ų0 (rigid)
angle thdegreesdegreesnone
angle kakJ/mol/rad²0 (rigid)

Parameter Table

ParameterValueUnitSource
Mass (Pd)106.42AMUIUPAC 2021
Charge0.0000ePd⁰ neutral cluster
sigma2.5114ÅHeinz 2008 Table 1, converted
epsilon25.7316kJ/molHeinz 2008 Table 1, converted
r(Pd–Pd)2.6760ÅLocal DFT optimisation of bare Pd₃ cluster
∠Pd–Pd–Pd60.0000degD₃h geometry, follows from equal bonds

References

  1. fftool — Padua, A. fftool: Tool to build force field input files for molecular simulation. Zenodo, 2015. DOI: 10.5281/zenodo.18618 · GitHub

  2. 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

  3. Pd₃ cluster geometry (literature context) — Moc, J.; Musaev, D. G.; Morokuma, K. Adsorption of Multiple H₂ Molecules on Pd₃ and Pd₄ Clusters. J. Phys. Chem. A 2000, 104, 11606–11614. DOI: 10.1021/jp0022104 · ACS page

This post is licensed under CC BY 4.0 by the author.