Ctrl+K

HTMD FFEvaluate - Easy MM force-field evaluation#

HTMD provides the FFEvaluate module which allows the easy evaluation of molecular mechanics force-fields. This serves multiple purposes:

  • Application development: Many computational biology applications require the evaluation of energies between molecules. This allows faster development of applications based on those energies.

  • Educational: Most molecular dynamics codes are either closed source or written in a highly optimized manner making them hard to understand. The FFEvaluate code has been written in python and numba which makes it easy to read, improve and extend.

  • Debugging: FFEvaluate can provide a detailed breakdown of all energies and visualization of forces which helps debug problems in a MD simulation or parameterization.

We would like to acknowledge parmed, an awesome python library which FFEvaluate uses to parse forcefield parameters.

Quick example#

from ffevaluation.ffevaluate import FFEvaluate, loadParameters, viewForces
from moleculekit.molecule import Molecule
from ffevaluation.home import home
from os.path import join

datadir = home(dataDir='thrombin-ligand-amber')
# Load the parameters
prm = loadParameters(join(datadir, 'structure.prmtop'))
# Load the topology
mol = Molecule(join(datadir, 'structure.prmtop'))
# Create a FFEvaluate object
ffev = FFEvaluate(mol, prm)

# Load some coordinates
mol.read(join(datadir, 'structure.pdb'))
# Calculate energy and forces
energies, forces, atom_energies = ffev.calculate(mol.coords)
# Visualize the forces
mol.view(viewer="vmd")
viewForces(mol, forces * 0.1)
# Calculate interaction energy between the protein and the ligand
ffev = FFEvaluate(mol, prm, betweensets=('protein', 'resname MOL'))
energies = ffev.calculateEnergies(mol.coords)

Detailed explanation#

from ffevaluation.ffevaluate import FFEvaluate, loadParameters
from moleculekit.molecule import Molecule
from ffevaluation.home import home
from os.path import join

datadir = home(dataDir='thrombin-ligand-amber')

The loadParameters function serves as a utility function for loading AMBER and CHARMM forcefield parameters. It is simply a wrapper for the parmed loading functions as it can sometimes get a bit complex to obtain parameters from some files, but if you already know how to use parmed you can use it directly to obtain a ParameterSet object of your parameter files.

prm = loadParameters(join(datadir, 'structure.prmtop'))

2024-06-12 13:20:39,190 - root - WARNING - Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.

Now that we have loaded the parameters, we load the prmtop topology in a Molecule object to obtain the bond, angle, dihedral, improper and atomtype information.

mol = Molecule(join(datadir, 'structure.prmtop'))

The main class for the calculations is the FFEvaluate class. When instantiated, the class constructor converts and stores the parameters in a format friendly to the numba python calculations. This is only done once for each system & parameterization combination as it incurrs a small computational overhead. The energy calculations are then done with the class methods calculate and calculateEnergies which can be evaluated for multiple different configurations of the system without having to incurr the same overhead by re-converting the parameters at every call.

ffev = FFEvaluate(mol, prm)

Now let’s read in some coordinates into our Molecule to evaluate the energies and forces.

mol.read(join(datadir, 'structure.pdb'))
energies, forces, atom_energies = ffev.calculate(mol.coords)

As we can see, FFEvaluate returns:

  • The energy of the system broken down in bonded, angle, dihedral, improper, VdW and electrostatic terms in a (6, nframes) shaped array where the rows correspond to 0: bond 1: lennard-jones 2: electrostatic 3: angle 4: dihedral 5: improper. You can sum over axis 0 to get the total energy.

  • The forces of all atoms in a (natoms, 3, nframes) shaped array

  • The energies of each individual atom in a (natoms, 6, nframes) shaped array. These energies are not physically meaningful as they are approximate atom energies calculated as the sum of all non-bonded energies on that atom plus the fraction of the bonded energy terms (i.e. for bonds we divide equally the bond energy between the two atoms, for angles between all three atoms and for dihedrals/impropers between all 4 atoms). However they can be useful for identifying issues in a parameterization or high-energy areas such as collisions.

where natoms is the number of atoms in mol and nframe is the number of frames in the coordinates of mol (in this case only one).

Warning: The first call to calculate in a python session can take some time as numba will compile the functions. Subsequent calls in the same session will be much faster.

FFEvaluate helps us visualize the forces using VMD. This can be very helpful for debugging cases where the simulation crashes due to atoms flying away because of too large forces. In case the vectors of the forces show up too large in VMD we can simply scale the forces down as much as we want for visualization.

viewForces(mol, forces * 0.1)
# On systems with 1000s of atoms it can take a while to render all force vectors in VMD

To just obtain energies in a nicely formatted fashion, the class also provides the calculateEnergies utility method which performs the same calculations but returns only the energies in a python dictionary.

ffev.calculateEnergies(mol.coords)

{'angle': 724.149178661847,
 'bond': 313.3073023605492,
 'dihedral': 3217.575130422163,
 'elec': -9277.547950450653,
 'improper': 30.433080793505184,
 'vdw': -1312.6701018873812,
 'total': -6304.753360099971}

Lastly, we are often interested specifically in (non-bonded) interaction energies between two molecules or generally between two sets of atoms. FFEvaluate allows us to calculate these energies using the betweensets argument. This argument accepts a list or tuple of two atomselect strings which define the two atom sets.

ffev = FFEvaluate(mol, prm, betweensets=('protein', 'resname MOL'))
energies = ffev.calculateEnergies(mol.coords)

Disclaimers#

While FFEvaluate provides an easy and fast way of obtaining MM energies and forces, the user should be aware of some pitfalls.

  • The code is still in beta. There are some less-used or older force-field features which are not yet supported (i.e. Urey-Bradley terms, PRMTOP EXCLUDED_ATOMS_LIST)

  • It will never be as fast as professional MD codes. The purpose of FFEvaluate is to make an easily usable and understandable MM evaluation code. Therefore on systems with large amount of non-bonded interations such as solvated systems with many waters the code will take a long time to evaluate the energies and forces.

  • Speed improvements such as cell-lists for non-bonded interations and particle-mesh Ewald are not implemented (but could be in the future). FFEvaluate will calculate all non-bonded interactions explicitly.

  • Implicit solvent calculations are still missing but might be added in the future.

If you are interested in contributing code to the project, or want to report a bug, feel free to chat us up on the HTMD github issues!

On this page