Theory · 8 methods

Force fields

What CAPS computes, as it computes it: the equation, its symbols with CAPS's defaults, when to use the method, the source file, the tests that check it, and any departure from the cited method.

DL_POLY export (FIELD, CONFIG, CONTROL)

The typed system as DL_POLY 4 input in units of kcal/mol, written the way the DL_POLY force-field tools write it: harmonic bonds and angles with DL_POLY's ½k (twice the force field's k), class II bonds and angles as quartic terms (2K2, 3K3, 4K4), torsions as cos terms (OPLS triples as cos3) with the 1-4 electrostatic and van der Waals scale factors on exactly one term per 1-4 pair (a zero-amplitude term where the force field has no torsion), impropers in the rule's own atom order (explicit rules before wildcards; the outer atoms read by index, backwards, then otherwise), class II out-of-plane terms as harmonic inversions with the centre first, and every van der Waals type pair written out (lj, nm 9-6 for class II, buck). Consecutive identical molecules share one molecular type (nummols N); CONFIG holds the atoms about the cell centre; CONTROL is a generic NVT run with the same cut-off.

$$\begin{gathered}\text{harm:}\ \ \tfrac{1}{2}k(r - r_0)^2 \qquad \text{quar:}\ \ \tfrac{1}{2}k(r - r_0)^2 + \tfrac{1}{3}k'(r - r_0)^3 + \tfrac{1}{4}k''(r - r_0)^4 \\ \text{cos:}\ \ A\left[1 + \cos(m\phi - \delta)\right] \qquad \text{cos3:}\ \ \tfrac{1}{2}\left[A_1(1 + \cos\phi) + A_2(1 - \cos 2\phi) + A_3(1 + \cos 3\phi)\right]\end{gathered}$$

DL_POLY 4 user manual forms

Symbol	Meaning	In CAPS
units	kcal/mol, Å, degrees	kcal
1-4 scale	electrostatic, van der Waals	the force field's

When to use it

To run the same model in DL_POLY. For PCFF, CVFF, COMPASS and OPLS 2005 the FIELD file equals, term for term, the one the DL_POLY force-field tools write for the same structure (bench/ff/check_dlpoly.py: OPLS 2005 110 of 111 polymers, PCFF 104 of 105, CVFF 100 of 110 — the ten where the reference's own amide/sulfone torsions depend on the atom order —, COMPASS 7 of 9; chlorine keeps its standard atomic weight where the reference's templates give 35.065).

Source: core/src/dlpoly.cpp (write_dlpoly); impropers_dlpoly in core/src/ffdef.cpp
Tested by: Dlpoly.FieldConfigControl; bench/ff/check_dlpoly.py against the reference FIELD files
Departure from the reference: DL_POLY has no class II cross terms, separate 1-4 Lennard-Jones parameters or combined bending–torsion: those are left out and named in the notes

Universal force field (UFF)

Every element from H to Lr, typed from its bonds, hybridisation and oxidation state. Used for clean-up and for any structure the other force fields cannot type (metals, silicon, phosphorus, halogens).

$$\begin{gathered}E = \sum k(r - r_0)^2 + \sum K\left[C_0 + C_1\cos\theta + C_2\cos 2\theta\right] + \sum \tfrac{1}{2}V\left[1 - \cos(n\phi_0)\cos(n\phi)\right] \\ \quad + \sum D\left[\left(\frac{x}{r}\right)^{12} - 2\left(\frac{x}{r}\right)^{6}\right] + E_\text{inversion}\end{gathered}$$

Rappé et al. 1992

Symbol	Meaning	In CAPS
r₀	natural bond length r_i + r_j + r_BO − r_EN	from the atom types
k	bond constant 332.06 Z_i Z_j / r₀³	from the atom types
D, x	van der Waals well and distance	geometric means

When to use it

Geometry clean-up of any molecule; runs on structures outside GAFF's scope. UFF has no charges of its own: use QEq when electrostatics matter.

Source: core/src/uff.cpp · core/src/uff_params.inc
Tested by: Uff.TableCoversEveryElementToLawrencium · Uff.ForcesMatchFiniteDifferences · Uff.CleanUpGivesTextbookGeometry
Departure from the reference: parameter table typed in from the paper's Tables I and III and the QEq paper (Rh6+3 Z1 = 3.508 as printed); a wall below 30° on the periodic angle terms, not in the paper; LAMMPS export checked term by term

References

Rappé, A. K., Casewit, C. J., Colwell, K. S., Goddard, W. A., Skiff, W. M., "UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations", J. Am. Chem. Soc. 114, 10024–10035 (1992). doi:10.1021/ja00051a040
Rappé, A. K., Goddard, W. A., "Charge equilibration for molecular dynamics simulations", J. Phys. Chem. 95, 3358–3363 (1991). doi:10.1021/j100161a070

Charge equilibration (QEq)

Charges that equalise the electronegativity of every atom at fixed total charge, for any element. The energy below is minimised; the linear system is solved by conjugate gradients.

$$\begin{gathered}E(q) = \sum_i \chi_i q_i + \frac{1}{2}\sum_i J_i q_i^2 + \sum_{i<j} q_i q_j K(r_{ij}) \\ K = \frac{14.40}{\sqrt{r^2 + a_{ij}^2}}\ \text{eV}, \qquad a_{ij} = 7.20\left(\frac{1}{J_i} + \frac{1}{J_j}\right)\ \text{Å}\end{gathered}$$

Rappé & Goddard 1991, with the Ohno–Klopman kernel

Symbol	Meaning	In CAPS
χ, J	electronegativity and idempotential	UFF / QEq table
K	shielded Coulomb kernel	Ohno–Klopman, tapered to zero at the cut-off
r_c	cut-off	10 Å (half the cell at most)

When to use it

Charges for structures with metals, silicon or other elements without force-field charges, such as fillers and coupling agents. Like the original QEq it gives ionic crystals large charges (quartz O ≈ −1.6 e).

Source: core/src/qeq.cpp
Tested by: QEq.MoleculesAreNeutralWithSensibleSigns · QEq.PeriodicCrystals
Departure from the reference: the Ohno–Klopman kernel with a 7th-order taper replaces the paper's Slater-orbital Coulomb integrals; J of hydrogen does not depend on its charge

References

Rappé, A. K., Goddard, W. A., "Charge equilibration for molecular dynamics simulations", J. Phys. Chem. 95, 3358–3363 (1991). doi:10.1021/j100161a070

General AMBER force field (GAFF)

Atom types for organic molecules assigned by rules in antechamber's order (rings, aromaticity, conjugation), with AMBER functional forms: harmonic bonds and angles, cosine torsions, 12–6 Lennard-Jones, 1–4 interactions scaled by ½ (LJ) and 1/1.2 (Coulomb).

$$\begin{gathered}E = \sum k_b(r - r_0)^2 + \sum k_\theta(\theta - \theta_0)^2 + \sum \frac{V_n}{2}\left[1 + \cos(n\phi - \gamma)\right] \\ \quad + \sum 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right] + \sum \frac{q_i q_j}{4\pi\epsilon_0 r}\end{gathered}$$

Wang et al. 2004

Symbol	Meaning	In CAPS
types	c3, ca, hc, ha … assigned by rules	data/typing/gaff*
1–4 scaling	LJ and Coulomb	½ and 1/1.2

When to use it

Organic polymers and additives. CAPS's built-in default for pure hydrocarbons; the library also holds GAFF2 and the other published force fields, with CAPS typing rules.

Source: core/src/typing.cpp · core/src/ffdef.cpp · data/forcefields
Tested by: Field.TypesPolystyrene · tests/test_typing.cpp · bench/ff/check_data_lammps.py
Departure from the reference: none

References

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., Case, D. A., "Development and testing of a general amber force field", J. Comput. Chem. 25, 1157–1174 (2004). doi:10.1002/jcc.20035

Gasteiger–Marsili charges (PEOE)

Partial equalisation of orbital electronegativity: charge flows along each bond from the less to the more electronegative atom, in six iterations damped by half each time, so it stays local.

$$\chi_i = a_i + b_i q_i + c_i q_i^2, \qquad \Delta q_i^{(k)} = \left(\frac{1}{2}\right)^{k}\sum_j \frac{\chi_j - \chi_i}{\chi^+}$$

Gasteiger & Marsili 1980

Symbol	Meaning	In CAPS
a, b, c	per element and hybridisation	H, C, N, O (sp³, sp², sp), F, Cl, Br, I, sp³ S
χ⁺	χ of the cation (H: 20.02)	—
k	iteration	1 … 6

When to use it

Quick charges for organic molecules. Other elements are refused; use QEq for them.

Source: core/src/grow.cpp (gasteiger_ch)
Tested by: Gasteiger.HeteroatomsGetTheirOwnParameters
Departure from the reference: parameters of Table 1 plus the sp³ sulfur extension

References

Gasteiger, J., Marsili, M., "Iterative partial equalization of orbital electronegativity–-a rapid access to atomic charges", Tetrahedron 36, 3219–3228 (1980). doi:10.1016/0040-4020(80)80168-2

Water models and their massless sites

Three-site models (SPC, SPC/E, SPC/Fw, TIP3P and its variants) carry their charges on O and H. Four-site models (TIP4P, TIP4P-Ew, TIP4P/2005, TIP4P/Ice, OPC) move the negative charge to a massless site M on the H–O–H bisector. TIP5P puts it on two massless lone pairs L out of the molecular plane, on the side away from the hydrogens. The sites are placed from O and H at every evaluation and their forces handed back to them by the chain rule; for the out-of-plane lone pairs the cross product adds one more term to the virial.

$$\begin{gathered}M = (1-\alpha)\,O + \tfrac{\alpha}{2}(H_1 + H_2), \qquad \alpha = \frac{d_{OM}}{r_{OH}\cos(\theta/2)} \\ L_\pm = O + a\,(\mathbf{r}_1 + \mathbf{r}_2) \pm c\,(\mathbf{r}_1 \times \mathbf{r}_2), \qquad a = -\frac{d_{OL}\cos(\phi/2)}{2 r_{OH}\cos(\theta/2)}, \qquad c = \frac{d_{OL}\sin(\phi/2)}{r_{OH}^2 \sin\theta}\end{gathered}$$

r₁, r₂ the O–H vectors, θ the H–O–H angle, φ the L–O–L angle; for TIP5P (r_OH 0.9572 Å, θ 104.52°, d_OL 0.70 Å, φ 109.47°) a = −0.344908 and c = 0.64438 Å⁻¹, GROMACS's virtual_sites3 function 4 (c in nm⁻¹ = 6.4438); the cross part scales as λ² when the atoms are scaled by λ, hence the extra virial term

Symbol	Meaning	In CAPS
d_OM	O–M distance	Å
d_OL	O–lone pair distance	Å
φ	L–O–L angle	°
a, c	out-of-plane site constants	—, 1/Å

When to use it

Solvating rubbers, fillers or biomolecules, water uptake and sorption, hydrated interfaces. TIP5P reproduces water's density maximum near 4 °C; it exports to GROMACS only (LAMMPS has no five-site water). Four-site models export to LAMMPS with the tip4p pair styles (M left out, its charge on O) and to GROMACS with M as a virtual site.

Source: core/src/water.cpp (water_models, apply_water_model, water_forcefield), core/src/field.cpp (virtual sites)
Tested by: tests/test_solvate.cpp Water.ModelsGeometryAndSites (TIP4P/2005's M 0.1546 Å from O, LAMMPS tip4p styles); Water.Tip5pLonePairsForcesAndVirial (lone pairs at 0.70 Å and 109.47° away from the hydrogens, GROMACS's a and c, forces against finite differences, the virial against −dE/dλ: 0.98082 with the λ² term, 0.634 without); TIP5P in GROMACS on 512 waters: Coulomb 34.371 against CAPS's 34.382 kcal/mol (PME), Lennard-Jones with its tail 6162.68 against 6162.69
Departure from the reference: rigid models are run with SHAKE/RATTLE on CAPS's side (GROMACS uses SETTLE); bond and angle constants are given only for flexible runs

References

Berendsen, H. J. C., Grigera, J. R., Straatsma, T. P., "The missing term in effective pair potentials", J. Phys. Chem. 91, 6269–6271 (1987). doi:10.1021/j100308a038
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., Klein, M. L., "Comparison of simple potential functions for simulating liquid water", J. Chem. Phys. 79, 926–935 (1983). doi:10.1063/1.445869
Abascal, J. L. F., Vega, C., "A general purpose model for the condensed phases of water: TIP4P/2005", J. Chem. Phys. 123, 234505 (2005). doi:10.1063/1.2121687
Mahoney, M. W., Jorgensen, W. L., "A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions", J. Chem. Phys. 112, 8910–8922 (2000). doi:10.1063/1.481505

OPLS charges from bond increments

OPLS 2005 gives each atom a charge key (an OPLS number chosen by its chemical environment) and each bonded pair of keys an increment that moves charge from one atom to the other. Summing over the bonds gives charges that balance bond by bond, so every molecule keeps its formal charge. OPLS-AA 2024's fixed per-type charges sometimes do not add up where two groups meet (a substituted ring, a benzylic junction); CAPS then uses OPLS 2005's typing and increments for the charges, and OPLS-AA 2024's types for everything else, and says so.

$$q_i = \sum_{j \text{ bonded to } i}\delta_{ij}, \qquad \delta_{ji} = -\delta_{ij}, \qquad \delta \text{ from the OPLS 2005 table by the pair's keys (its first entry for a pair)}$$

OPLS 2005 as in Banks et al. 2005

Symbol	Meaning	In CAPS
key	the atom's OPLS number from its environment	charge rules in the OPLS 2005 typing file
δ_ij	bond charge increment	6179 pairs

When to use it

Automatic whenever OPLS-AA 2024's own charges miss the formal charge; by hand with Field › Use OPLS 2005's charges or charges "increments" in recipes and Python.

Source: core/src/ffdef.cpp (parameterize, companion_charges)
Tested by: Typing.Opls2005ChargesFromBondIncrements · Typing.CompanionChargesFromOpls2005
Departure from the reference: charges verified against the reference program's OPLS 2005 on 110 polymers; the combination with OPLS-AA 2024 parameters is CAPS's choice

References

Banks, J. L., Beard, H. S., Cao, Y., Cho, A. E., Damm, W., Farid, R., Felts, A. K., Halgren, T. A., Mainz, D. T., Maple, J. R., others, "Integrated Modeling Program, Applied Chemical Theory (IMPACT)", J. Comput. Chem. 26, 1752–1780 (2005). doi:10.1002/jcc.20292

Force fields by group and crystal potentials

Each group of molecules (a filler and its matrix, each component of a blend) is typed and parameterised by its own force field with its own charges and mixing rule; the Lennard-Jones pairs between groups follow a cross rule chosen for them, or pairs given explicitly, and every cross pair is written out in the engine files. A crystal group can instead take a literature many-body potential (Tersoff, Stillinger–Weber, Vashishta, EAM) that LAMMPS reads from its file. Units of the LAMMPS files: automatic (metal — eV, ps, bar — when a potential is read by LAMMPS in metal units only: AIREBO, AIREBO-M, REBO, MEAM; the published file then unchanged and every force-field energy parameter divided by 23.060549, LAMMPS's own factor; real — kcal/mol, fs, atm — otherwise, where LAMMPS converts Tersoff, SW, Vashishta, GW and EAM files itself), or real with AIREBO / REBO (a copy of the file CAPS converts: A, B and the LJ, torsion and Morse ε; the splines are dimensionless), or metal for any system. Both are exact: the physics does not depend on the unit system, only LAMMPS's electric-conversion constants differ in their last digits. MEAM maps each element to a library entry (by default the first of its atomic number) and keeps the parameter file's index order.

$$\begin{gathered}\epsilon_{ab} = \sqrt{\epsilon_a\epsilon_b}\ \text{ or }\ \frac{\epsilon_a + \epsilon_b}{2}, \qquad \sigma_{ab} = \frac{\sigma_a + \sigma_b}{2}\ \text{ or }\ \sqrt{\sigma_a\sigma_b} \\ \text{sixth power:}\ \ \sigma_{ab} = \left(\frac{\sigma_a^6 + \sigma_b^6}{2}\right)^{1/6}, \qquad \epsilon_{ab} = \frac{2\sqrt{\epsilon_a\epsilon_b}\,\sigma_a^3\sigma_b^3}{\sigma_a^6 + \sigma_b^6}\end{gathered}$$

between groups; inside a group its force field's own rule

Symbol	Meaning	In CAPS
a, b	atom types of two different groups	each group's own parameters
crystal site	one type per element, standard atomic mass, no charge	UFF's well depth D and minimum x across; 12-6 σ = x/2^(1/6), or 9-6 σ = x beside PCFF / COMPASS
potential file	LAMMPS's own file, converted from metal to real units by LAMMPS	the library (data/potentials) or the user's

When to use it

Composites and blends where one force field cannot describe every part: silica, SiC, graphene / CNT or a metal in a PCFF, COMPASS, OPLS or GAFF polymer; blend components parameterised by different force fields. What one simulation cannot hold is refused with the reason: different 1-4 scalings (unless the first group's is taken for all), 9-6 with 12-6 (unless 9-6 sites enter cross pairs as 12-6 with the same ε and minimum); MEAM in real units (LAMMPS reads it in metal units only). A crystal under a literature potential runs in LAMMPS; CAPS's own runs hold it still.

Source: core/src/ffmerge.cpp · core/src/manybody.cpp · capi caps_field_assign_groups · data/potentials
Tested by: FFMerge.SameForceFieldInTwoPartsIsExact · FFMerge.DifferentFamiliesAreRefusedOrMergedAsAsked · FFMerge.ManyBodyGroup · FFMerge.PotentialLibrary · FFMerge.AireboWritesMetalUnits · FFMerge.MeamMapsLibraryEntries
Departure from the reference: Si Tersoff slab under hexane: LAMMPS van der Waals beside Tersoff equals CAPS's for GAFF2 (366.6761), OPLS 2005 (296.6776), PCFF (133.1490 / 133.149) and COMPASS (84.3140 / 84.3139 kcal/mol); GAFF2 + OPLS polystyrene in GROMACS 573.546, LAMMPS 573.549 kcal/mol; (6,6) nanotube under AIREBO in hexane, metal units: van der Waals beside AIREBO GAFF2 3.49464, PCFF 16.67939, COMPASS 1.26178, OPLS 2005 −2.13761 kcal/mol in LAMMPS and CAPS; PS melt real vs metal export identical to 1e-9 in bonded and van der Waals terms; AIREBO, AIREBO-M, REBO in real units with CAPS's converted file: energy equal to the metal run ×23.060549 to 1e-13, forces to 2e-12 kcal/mol/Å, pressure to every digit; MEAM (SiC set) in the composite equals the slab alone (atom_style atomic) to 1e-10 eV, van der Waals beside it equals CAPS's (GAFF2 366.67611, PCFF 133.14904 kcal/mol)

References

Rappé, A. K., Casewit, C. J., Colwell, K. S., Goddard, W. A., Skiff, W. M., "UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations", J. Am. Chem. Soc. 114, 10024–10035 (1992). doi:10.1021/ja00051a040