How to evaluate force-field energies on a frame#

Goal#

Compute bonded + non-bonded force-field energies (and optionally per-atom forces) for any coordinate snapshot, without going through an MD engine. Useful for ranking docking poses, validating a build, computing per-frame interaction energies, or scripting custom energy-based analyses.

Minimal example#

from moleculekit.molecule import Molecule
from ffevaluation.ffevaluate import FFEvaluate, loadParameters

mol = Molecule("./build/structure.prmtop")
mol.read("./build/structure.pdb")

prm = loadParameters("./build/structure.prmtop")          # parmed ParameterSet
ffev = FFEvaluate(mol, prm)
energies, forces, per_atom = ffev.calculate(mol.coords)

print("total energy:", energies.sum())                    # kcal/mol

FFEvaluate.calculate(coords) returns three arrays:

  • energies - (6, n_frames) per-term system totals in the order (bond, vdw, elec, angle, dihedral, improper).

  • forces - (n_atoms, 3, n_frames) per-atom force vector in kcal/mol/Å.

  • atmnrg - (n_atoms, 6, n_frames) per-atom energy contributions in the same term order.

FFEvaluate.calculateEnergies(coords) is a thin wrapper that calls calculate and discards both forces and atmnrg, returning only the energies (as a dict when formatted=True).

Both calculate and calculateEnergies accept an optional box=mol.box argument ((3, n_frames)). Pass it for periodic systems so non-bonded distances get wrapped through the box; omit it for vacuum / non-periodic snapshots.

Parameters that matter#

Parameter

What it does

mol

Molecule whose topology defines the bond / angle / dihedral graph. Frames in mol.coords define how many coordinate snapshots get evaluated.

prm

A parmed ParameterSet - load from prmtop with loadParameters("...") or build from a CHARMM .prm via loadParameters(prm_path).

betweensets

A tuple of two atom-selection strings - restricts non-bonded energy / force calc to interactions between the two sets. Computes only LJ + electrostatics, no intramolecular bonded terms.

cutoff

Non-bonded cutoff in Å. 0 (default) means no cutoff - exact electrostatics + LJ. Set to e.g. 9.0 to match a simulation’s cutoff.

rfa

Use with cutoff to enable the reaction-field approximation for electrostatics beyond the cutoff.

solventDielectric

Solvent dielectric used when rfa=True. Default 78.5 (water).

Common variations#

Per-frame energies on a trajectory#

mol.read("./output.xtc")                                  # multi-frame
all_energies, all_forces, _ = ffev.calculate(mol.coords)
# all_energies shape: (n_terms, n_frames)

Interaction energy between two selections#

ffev = FFEvaluate(mol, prm, betweensets=("resname LIG", "protein"))
e, _, _ = ffev.calculate(mol.coords)
print("protein-ligand interaction energy:", e.sum(axis=0))

Restricting to a pair of sets is much cheaper than the full system - no bonded terms are computed and the non-bonded loop is restricted to inter-set pairs.

Use a CHARMM parameter file instead of prmtop#

prm = loadParameters("./params.prm")                      # CHARMM .prm
ffev = FFEvaluate(mol, prm)

loadParameters dispatches by extension - only .prmtop, .prm, and .frcmod are recognised. Any other extension raises RuntimeError. For CHARMM systems built from a PSF, point loadParameters at the matching .prm; for AMBER pass the .prmtop directly.

Decomposed energy report#

e_by_term = ffev.calculateEnergies(mol.coords, formatted=True)
# Dict keyed by term name. Keys: bond, angle, dihedral, improper, vdw, elec, total.
for k, v in e_by_term.items():
    print(f"{k:>10}: {v:.3f} kcal/mol")

calculateEnergies(..., formatted=True) is the readable variant of calculate(...).

Gotchas#

  • FFEvaluate expects the same atom order in mol.coords and in the parameter set. If your mol has been reordered / filtered after the parmtop was loaded, the bond graph won’t match. Always load both from the same source.

  • With betweensets, the result excludes intra-set bonded and intra-set non-bonded energies. The reported number is the cross-set interaction only.

  • Default cutoff=0 is exact but O(N²) - on >100k-atom systems this is slow. Setting cutoff > 0 only short-circuits the per-pair maths inside the loop; the pair iteration itself is still O(N²) (there’s no neighbour list). Match your simulation’s cutoff (typically 9-12 Å) for production-scale evaluations, but don’t expect linear-scaling.

  • The output unit is kcal/mol regardless of the parameter file format.

See also#