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.
from ffevaluation.ffevaluate import FFEvaluate, loadParameters, viewForces from moleculekit.molecule import Molecule from htmd.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 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)
from ffevaluation.ffevaluate import FFEvaluate, loadParameters from moleculekit.molecule import Molecule from ffevaluation.home import home from os.path import join datadir = home(dataDir='test-thrombin-ligand-amber')
loadParameters function serves as a utility function for loading
CHARMM forcefield parameters. It is simply a wrapper
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'))
Now that we have loaded the parameters, we load the
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
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
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.
natoms is the number of atoms in
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
calculateEnergies utility method which performs the
same calculations but returns only the energies in a python dictionary.
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)
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!