HTMD provides the FFEvaluate module which allows the easy evaluation of molecular mechanics force-fields. This serves multiple purposes:
We would like to acknowledge parmed, an awesome python library which FFEvaluate uses to parse forcefield parameters.
from htmd.ffevaluation.ffevaluate import FFEvaluate, loadParameters, viewForces from htmd.molecule.molecule import Molecule from htmd.home import home from os.path import join datadir = join(home(dataDir='test-ffevaluate'), '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 htmd.ffevaluation.ffevaluate import FFEvaluate, loadParameters from htmd.molecule.molecule import Molecule from htmd.home import home from os.path import join datadir = join(home(dataDir='test-ffevaluate'), '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:
0: bond 1: lennard-jones 2: electrostatic 3: angle 4: dihedral 5: improper. You can sum over axis 0 to get the total energy.
(natoms, 3, nframes)shaped array
(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.
FFEvaluateis 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.
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!