Tutorials‎ > ‎

Force Fields

pDynamo employs standard biomolecular force fields, including AMBER, CHARMM and OPLS (see here for links). It can read files, prepared by external programs, that contain the force field definitions for a system, or, alternatively, set up the force field definitions for a system by itself.

External Files

pDynamo can read a number of third-party files that contain systems and their force field definitions, including those generated by the AMBER and CHARMM programs. A suitable set of commands for reading AMBER topology and coordinate files is:

from pBabel    import AmberTopologyFile_ToSystem, AmberCrdFile_ToCoordinates3
from pMolecule import NBModelFull

system = AmberTopologyFile_ToSystem ( "mySystem.top" )
system.coordinates3 = AmberCrdFile_ToCoordinates3 ( "mySystem.crd" )
system.DefineNBModel ( NBModelFull ( ) )
system.Summary ( )
system.Energy  ( )

The equivalent commands for CHARMM are a little more complicated as CHARMM requires a parameter file in addition to its PSF and coordinate files. Note that CHARMM PSF files should be in XPLOR format if a valid energy model is to be created for the system.

from pBabel    import CHARMMCRDFile_ToCoordinates3, CHARMMParameterFiles_ToParameters, CHARMMPSFFile_ToSystem
from pMolecule import NBModelFull

parameters          = CHARMMParameterFiles_ToParameters ( [ "myParameters.prm" ] )
system              = CHARMMPSFFile_ToSystem ( "mySystem.psfx", isXPLOR = True, parameters = parameters )
system.coordinates3 = CHARMMCRDFile_ToCoordinates3 ( "mySystem.chm" )
system.DefineNBModel ( NBModelFull ( ) )
system.Summary ( )
system.Energy  ( )

Force Field Definition Within pDynamo

A number of force field parameter sets are distributed with the pDynamo program, which are located in the parameters/forceFields directory of the distribution. The subdirectories within this directory correspond to different parameter sets so that, for example, the charmm/c36a2 directory contains the c36a2 revision of the CHARMM force field, whereas the opls/protein directory contains the OPLS-AA force field for proteins.

All force fields that pDynamo provides are assigned to a system in the same way. A two-step procedure is used, which consists of first determining the atom type of each atom in the system, followed by the definition of the individual force field terms and their parameters.

Atom typing itself can be done in two ways. The simplest, and the one that many biomolecular force fields employ, is to use the system's sequence information, notably an atom's label and the type of component, with its link and variant modifications, in which it occurs. This requires that the definitions of the atom types of the atoms in each of the components, links and variants in the system exist in the force field directory. A second, more general approach, that does not rely on sequence information, makes use of connectivity patterns, which are definitions that define a mapping between specific chemical motifs or substructures and force field atom types. Note that it is possible to use both methods to type the atoms of a system, in which case the more specific sequence-based method is applied first, followed up by the more general pattern-based method for the remaining untyped atoms.

The force field terms for the system are built once all the atoms have been typed. This is done by using the system's connectivity information to define all the bonds, angles, dihedrals and other terms that the force field requires, and then picking the parameters for each of the interactions based upon the types of the atoms that are involved.

It is not unusual when assigning a force field to a system to have untyped atoms or missing parameters for particular terms. In these cases the additional information that is required needs to be added to the parameter set by editing the appropriate definition files in the parameter set directory. It is beyond the scope of this tutorial to describe how to derive force field parameters and so it is assumed that the necessary definitions can be found.

All files in the parameter set are in YAML format. The atom types and their properties are in the file atomTypes.yaml, the definitions for atom typing that employ sequence information are in subdirectories with the names components, links and variants, whereas those that use connectivity patterns are in patterns.yaml. Likewise, there are files for each of the different types of force field terms that occur in the parameter set. Examples are harmonicBondParameters.yaml for harmonic bond terms, and lennardJonesParameters.yaml that specify the parameters for the Lennard-Jones interaction. Effort has been made to make the format of these files as easy to understand as possible and so a detailed description of the formats will not be given here. Nevertheless a few points are worth making:
  • The convention is used that the absence of a file indicates that that term is not to be calculated. Thus, for example, if harmonicBondParameters.yaml is absent, no harmonic bonds will be calculated and, similarly, if patterns.yaml is missing, no atom typing by pattern will be performed. Normally the only file that has to be present is atomTypes.yaml.
  • The atom typing files in the component, link and variant directories of the parameter set normally also specify the charges of the atoms, in addition to the atom/atom type mapping.
  • The specification of appropriate connectivity patterns for atom typing requires some care and there are several important points to remember: (i) patterns are applied to a system in the order in which they are specified in the file, which means that more specific patterns should occur before more general ones; (ii) hydrogens are normally omitted from MM patterns and are typed implicitly using the value of the hydrogenLabel keyword of the typed atom to which they are attached. Certain non-standard hydrogen types may require explicit treatment in a pattern, in which case they will be handled equivalently to the atoms of other elements; (iii) the charge of an atom type given in a pattern overrides that specified in the atom type definition and refers to the total charge of the atom and any implicit hydrogens it may have. The charge that will actually be assigned to the atom is the difference between the pattern charge and the sum of the charges on its implicit hydrogens — these being defined in the hydrogen's atom type definition; and (iv) defining bond patterns that match aromatic systems can be troublesome because the latter can occur in a number of different, albeit equivalent, canonical forms. It is best to use Undefined as the aromatic bond type in these cases.
  • pDynamo uses Angstroms, kJ mol-1 and radians as its default units of length, energy and angle, respectively, whereas many force fields use other units by default. To make life easier, the original, default values are kept in the YAML files, together with their units, and pDynamo makes the conversions that are necessary.