OPLS: verifying the xnn implementation against OpenMM and the 1996 paper#
The xnn OPLS model is a clean-room implementation of the published
functional form (Jorgensen, Maxwell & Tirado-Rives, J. Am. Chem. Soc.
118, 11225, 1996, eqs 1-4). This notebook verifies it two independent
ways:
numerical parity with OpenMM — the same parameter library and topology (atom types assigned from the parameter file’s SMARTS templates, bonds perceived from the coordinates) are assembled into an OpenMM
System(an entirely independent energy/force engine) and evaluated on randomized conformations of butane, ethanol and ethylene. Energies and forces must agree to the precision of the unit constants (~1e-7 kJ/mol);the paper’s own numbers — relaxed torsional scans reproduce the OPLS-AA column of Table 1, and the Fourier/Ryckaert-Bellemans torsion conversion reproduces the dual-form rows of Table 2 of the L-OPLS paper (Siu, Pluhackova & Böckmann, JCTC 8, 1459, 2012).
Unit conventions are part of what is verified: xnn stores OPLS
parameters in kcal/mol with the thermochemical calorie (4.184 kJ exactly)
and uses the CODATA Coulomb constant, matching the GROMACS/OpenMM
ecosystem in which OPLS parameters are distributed.
0. Setup#
import warnings
warnings.filterwarnings("ignore")
import math
import numpy as np
import torch
torch.set_default_dtype(torch.float64)
torch.manual_seed(0)
rng = np.random.default_rng(1)
import openmm as mm
import openmm.unit as u
import xnn
from xnn.common.data import structure_to_graph
from xnn.common.models import ForceStressOutput
from xnn.ffnn.models import OPLS, builtin_library, rb_to_fourier, fourier_to_rb
from xnn.ffnn.models.oplslib import (resolve_angle_type, resolve_bond_type,
resolve_dihedral_type,
resolve_improper_type)
EV_TO_KJ = 96.48533212331
print("xnn:", xnn.__version__, "| torch:", torch.__version__,
"| openmm:", mm.version.version)
xnn: 0.1.0 | torch: 2.5.1+cu121 | openmm: 8.6.0.dev-c6173db
1. An OpenMM twin of the same library + topology#
The translation is direct: harmonic bonds/angles carry a factor 2 (OPLS
writes k (x-x0)^2, OpenMM k/2 (x-x0)^2), the Fourier torsion becomes
four phased periodic torsions, the V2 improper a phase-180 periodic
torsion, and 1,4 pairs become explicit exceptions. The one subtlety is the
Lennard-Jones combination rule: OpenMM’s NonbondedForce hard-codes
Lorentz-Berthelot mixing, so OPLS geometric mixing needs a
CustomNonbondedForce (charges stay on the NonbondedForce).
def openmm_twin(lib, top, impropers=(), fudge_lj=0.5, fudge_qq=0.5):
"""Assemble an OpenMM Context evaluating the same OPLS model.
``impropers`` are the (i, j, center, l) quadruples the xnn model placed
(``model.impropers``); their parameters resolve by class pattern.
"""
at = lib.atom_types
cls = [at[n]["cls"] for n in top.types]
cls_oop = [lib.cls(n, "oop") for n in top.types]
system = mm.System()
for name in top.types:
system.addParticle(at[name]["mass"])
nb = mm.NonbondedForce() # Coulomb + 1,4 exceptions
nb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
lj = mm.CustomNonbondedForce( # OPLS geometric mixing
"4*eps*((sig/r)^12-(sig/r)^6); sig=sqrt(sig1*sig2);"
" eps=sqrt(eps1*eps2)")
lj.addPerParticleParameter("sig")
lj.addPerParticleParameter("eps")
for name in top.types:
nb.addParticle(at[name]["charge"], 0.1, 0.0)
lj.addParticle([max(at[name]["sigma"], 1e-6) * 0.1,
at[name]["epsilon"] * 4.184])
for i, j in top.exclusions + top.pairs14:
nb.addException(i, j, 0.0, 0.1, 0.0)
lj.addExclusion(i, j)
for i, j in top.pairs14:
ti, tj = at[top.types[i]], at[top.types[j]]
nb.addException(i, j, fudge_qq * ti["charge"] * tj["charge"],
max(math.sqrt(ti["sigma"] * tj["sigma"]), 1e-6) * 0.1,
fudge_lj * math.sqrt(ti["epsilon"] * tj["epsilon"])
* 4.184, replace=True)
nb.setForceGroup(0); lj.setForceGroup(4)
system.addForce(nb); system.addForce(lj)
bond = mm.HarmonicBondForce(); bond.setForceGroup(1)
for i, j in top.bonds:
bt = lib.bond_types[resolve_bond_type(lib.bond_types, cls[i], cls[j])]
bond.addBond(i, j, bt["r0"] * 0.1, 2 * bt["k"] * 4.184 * 100)
system.addForce(bond)
ang = mm.HarmonicAngleForce(); ang.setForceGroup(2)
for i, j, k in top.angles:
a = lib.angle_types[resolve_angle_type(lib.angle_types,
cls[i], cls[j], cls[k])]
ang.addAngle(i, j, k, math.radians(a["theta0"]), 2 * a["k"] * 4.184)
system.addForce(ang)
tors = mm.PeriodicTorsionForce(); tors.setForceGroup(3)
const = 0.0
for i, j, k, l in top.dihedrals:
key = resolve_dihedral_type(lib.dihedral_types,
cls[i], cls[j], cls[k], cls[l])
v = lib.dihedral_types[key]["v"]
const += v[0] * 4.184 # the constant V0 offset
for n, (vn, ph) in enumerate(zip(v[1:], (0.0, math.pi, 0.0, math.pi)),
start=1):
if vn:
tors.addTorsion(i, j, k, l, n, ph, 0.5 * vn * 4.184)
for i, j, k, l in impropers:
key = resolve_improper_type(lib.improper_types, cls_oop[i], cls_oop[j],
cls_oop[k], cls_oop[l])
v2 = lib.improper_types[key]["v2"]
tors.addTorsion(i, j, k, l, 2, math.pi, 0.5 * v2 * 4.184)
system.addForce(tors)
ctx = mm.Context(system, mm.VerletIntegrator(1e-3),
mm.Platform.getPlatformByName("Reference"))
return ctx, const
def openmm_eval(ctx, const, pos, groups=None):
"""Energy (kJ/mol) and forces (kJ/mol/nm) from the OpenMM twin."""
ctx.setPositions(np.asarray(pos) * 0.1)
kw = {"groups": groups} if groups is not None else {}
st = ctx.getState(getEnergy=True, getForces=True, **kw)
e = st.getPotentialEnergy().value_in_unit(u.kilojoule_per_mole)
f = st.getForces(asNumpy=True).value_in_unit(
u.kilojoule_per_mole / u.nanometer)
return e + (const if groups is None or 3 in groups else 0.0), f
2. Energy and force parity on randomized conformations#
Three gas-phase molecules cover every term of the force field: butane
(bonds, angles, Fourier torsions, 1,4 scaling), ethanol (heteroatoms,
alcohol torsions, a zero-LJ hydroxyl hydrogen) and ethylene (V2
impropers, wildcard X-CM-CM-X torsion). Ten heavily jittered
conformations each.
def molecule_set():
butane = (np.array([[0.0, 0.0, 0.0], [1.53, 0.0, 0.0], [2.05, 1.44, 0.0],
[3.58, 1.44, 0.0],
[-0.4, -0.5, 0.9], [-0.4, -0.5, -0.9], [-0.4, 1.0, 0.0],
[1.93, -0.52, 0.88], [1.93, -0.52, -0.88],
[1.65, 1.96, -0.88], [1.65, 1.96, 0.88],
[3.98, 0.44, 0.0], [3.98, 1.96, 0.88],
[3.98, 1.96, -0.88]]),
[6] * 4 + [1] * 10)
ethanol = (np.array([[0.0, 0.0, 0.0], [1.512, 0.0, 0.0], [2.0, 1.32, 0.0],
[-0.39, -0.51, 0.89], [-0.39, -0.51, -0.89],
[-0.39, 1.02, 0.0], [1.90, -0.52, 0.88],
[1.90, -0.52, -0.88], [2.60, 1.30, 0.7]]),
[6, 6, 8, 1, 1, 1, 1, 1, 1])
ethylene = (np.array([[0.0, 0.0, 0.0], [1.34, 0.0, 0.0],
[-0.54, 0.94, 0.0], [-0.54, -0.94, 0.0],
[1.88, 0.94, 0.0], [1.88, -0.94, 0.0]]),
[6, 6, 1, 1, 1, 1])
return {"butane": butane, "ethanol": ethanol, "ethylene": ethylene}
lib = builtin_library("oplsaa")
print(f"{'molecule':<10} {'max |dE| (kJ/mol)':>18} {'max |dF| (kJ/mol/nm)':>22}")
for name, (pos0, z) in molecule_set().items():
# atom types from the library's SMARTS templates, bonds perceived from
# the coordinates, impropers placed at the trigonal centers
opls = OPLS.from_atoms((pos0, z), lib, cutoff=100.0)
top = opls.topology
model = ForceStressOutput(opls)
ctx, const = openmm_twin(lib, top, opls.impropers)
de = df = 0.0
for _ in range(10):
pos = pos0 + 0.1 * rng.standard_normal(pos0.shape)
out = model(structure_to_graph(
{"pos": torch.tensor(pos), "atomic_numbers": torch.tensor(z)},
cutoff=100.0))
e_o, f_o = openmm_eval(ctx, const, pos)
de = max(de, abs(float(out["energy"]) * EV_TO_KJ - e_o))
df = max(df, np.abs(out["forces"].detach().numpy() * EV_TO_KJ * 10
- f_o).max())
print(f"{name:<10} {de:>18.2e} {df:>22.2e}")
assert de < 1e-6 and df < 1e-5
molecule max |dE| (kJ/mol) max |dF| (kJ/mol/nm)
butane 9.67e-08 2.89e-07
ethanol 7.29e-08 1.35e-06
ethylene 1.70e-07 3.53e-07
Agreement at ~1e-7 kJ/mol in the energy and ~1e-6 kJ/mol/nm in every force
component — the residual is the rounding of OpenMM’s Coulomb constant
(138.935456), i.e. the two engines agree to the precision at which the
physical constants themselves are written down.
3. Term-by-term decomposition#
xnn returns the energy decomposition directly; OpenMM force groups give
the same split (1,4 interactions live in the NonbondedForce exceptions,
so they are compared together with the Coulomb group).
pos0, z = molecule_set()["butane"]
model = OPLS.from_atoms((pos0, z), lib, cutoff=100.0)
top = model.topology
ctx, const = openmm_twin(lib, top, model.impropers)
pos = pos0 + 0.1 * rng.standard_normal(pos0.shape)
out = model(structure_to_graph(
{"pos": torch.tensor(pos), "atomic_numbers": torch.tensor(z)},
cutoff=100.0))
rows = [
("bonds", float(out["e_bond"]) * EV_TO_KJ, {1}),
("angles", float(out["e_angle"]) * EV_TO_KJ, {2}),
("torsions (+impropers)",
(float(out["e_torsion"]) + float(out["e_improper"])) * EV_TO_KJ, {3}),
("Lennard-Jones (direct)", float(out["e_lj"]) * EV_TO_KJ, {4}),
("Coulomb + all 1,4",
(float(out["e_coulomb"]) + float(out["e_lj14"])
+ float(out["e_coulomb14"])) * EV_TO_KJ, {0}),
]
print(f"{'term':<24} {'xnn (kJ/mol)':>14} {'OpenMM':>12} {'diff':>10}")
for label, e_x, grp in rows:
e_o, _ = openmm_eval(ctx, const, pos, groups=grp)
print(f"{label:<24} {e_x:>14.6f} {e_o:>12.6f} {e_x - e_o:>10.2e}")
assert abs(e_x - e_o) < 1e-6
term xnn (kJ/mol) OpenMM diff
bonds 694.914753 694.914753 -2.59e-10
angles 86.283288 86.283288 -5.34e-11
torsions (+impropers) 22.001321 22.001321 -2.24e-12
Lennard-Jones (direct) 0.316461 0.316461 -3.07e-12
Coulomb + all 1,4 11.832974 11.832974 -9.39e-08
4. Anchors from the papers themselves#
Relaxed ethane barrier (Table 1 of Jorgensen 1996; the oplsaa-1996
built-in carries the paper’s original alkane torsions):
from ase import Atoms
from ase.constraints import FixInternals
from ase.optimize import BFGS
from xnn.common.deploy import XNNCalculator
d, dh, ang = 1.529, 1.09, math.radians(110.7)
pos = [[0.0, 0.0, 0.0], [d, 0.0, 0.0]]
for base, sign, off in ((0.0, -1.0, 60.0), (d, 1.0, 0.0)):
for k in range(3):
phi = math.radians(off + 120.0 * k)
pos.append([base - sign * dh * math.cos(math.pi - ang),
dh * math.sin(math.pi - ang) * math.cos(phi),
dh * math.sin(math.pi - ang) * math.sin(phi)])
z6 = [6, 6] + [1] * 6
model = OPLS.from_atoms((pos, z6), "oplsaa-1996", cutoff=30.0)
energies = {}
for target in (60.0, 0.0):
at = Atoms(numbers=z6, positions=pos)
at.calc = XNNCalculator(ForceStressOutput(model), cutoff=model.cutoff)
at.set_dihedral(2, 0, 1, 5, target, indices=[5, 6, 7])
at.set_constraint(FixInternals(dihedrals_deg=[[target, [2, 0, 1, 5]]]))
BFGS(at, logfile=None).run(fmax=1e-5, steps=500)
energies[target] = at.get_potential_energy() / 4.3364104242e-2
barrier = energies[0.0] - energies[60.0]
print(f"relaxed ethane rotational barrier: {barrier:.2f} kcal/mol "
f"(Table 1: 3.01)")
assert abs(barrier - 3.01) < 0.02
relaxed ethane rotational barrier: 3.01 kcal/mol (Table 1: 3.01)
The torsion-form conversion against Table 2 of the L-OPLS paper, which
lists the same hexane CT-CT-CT-CT torsion in both Ryckaert-Bellemans and
Fourier form (kJ/mol) — including the constant V0 that makes the two
match exactly:
rb = [0.518787, -0.230192, 0.896807, -1.49134, 0.0, 0.0]
fourier = [-0.305938, 2.697394, -0.896807, 0.74567, 0.0]
print("rb_to_fourier:", [round(x, 6) for x in rb_to_fourier(rb)])
print("Siu Table 2: ", fourier)
assert max(abs(a - b) for a, b in zip(rb_to_fourier(rb), fourier)) < 1e-6
assert max(abs(a - b) for a, b in zip(fourier_to_rb(fourier), rb)) < 1e-6
print("both directions agree to 1e-6 kJ/mol")
rb_to_fourier: [-0.305938, 2.697394, -0.896807, 0.74567, -0.0]
Siu Table 2: [-0.305938, 2.697394, -0.896807, 0.74567, 0.0]
both directions agree to 1e-6 kJ/mol
Summary#
check |
result |
|---|---|
energy parity vs OpenMM (butane / ethanol / ethylene, 10 random conformations each) |
max deviation ~1e-7 kJ/mol |
force parity vs OpenMM |
max deviation ~1e-6 kJ/mol/nm |
term-by-term decomposition vs OpenMM force groups |
~1e-7 kJ/mol per term |
relaxed ethane barrier vs Jorgensen 1996 Table 1 |
3.01 vs 3.01 kcal/mol |
Fourier / Ryckaert-Bellemans conversion vs Siu 2012 Table 2 |
exact to the table’s digits |
Together with examples/ffnn/opls/opls_conformational_energetics.ipynb
(the full Table 1 reproduction) this establishes that the xnn OPLS
implementation is numerically equivalent to the reference MD engines and
faithful to the published parameterizations.