MM Models

pDynamo employs standard biomolecular force fields, namely Amber, CHARMM and OPLS (see here for links), implemented in the MMModelAMBER, MMModelCHARMM and MMModelOPLS classes, respectively. To define these models, 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.

Files from third party MM programs

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

from pBabel import ImportCoordinates3 , \

ImportSystem

system = ImportSystem ( "mySystem.top" )

system.coordinates3 = ImportCoordinates3 ( "mySystem.crd" )

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 CHARMMParameterFileReader , \

ImportCoordinates3 , \

ImportSystem

parameters = CHARMMParameterFileReader.PathsToParameters ( [ "myParameters.prm" ] )

system = ImportSystem ( "mySystem.psfx" ,

isXPLOR = True ,

parameters = parameters )

system.coordinates3 = ImportCoordinates3 ( "mySystem.chm" )

Reading Gromacs files is even more complicated as Gromacs topology files have not yet been incorporated into pDynamo's import and export mechanism, and so explicit file readers must be used. The commands are:

from pBabel import GromacsDefinitionsFileReader , \

GromacsParameterFileReader , \

ImportCoordinates3

parameters = GromacsParameterFileReader.PathToParameters ( "mySystem.top" )

system = GromacsDefinitionsFileReader.PathToSystem ( "mySystem.top" ,

parameters = parameters )

system.coordinates3 = ImportCoordinates3 ( "mySystem.gro" )

The reading of Amber and CHARMM force field files will return systems with instances of the MMModelAMBER and MMModelCHARMM classes, respectively. By contrast, the MM models of the systems returned by reading Gromacs files will vary, as these files can contain definitions for all of the models that pDynamo implements.

Force field definition within pDynamo

It is often convenient to use external files to define a system's MM model. For example, it is common to use a third party program to set up a simulation system and perform long MM molecular dynamics simulations, before switching to pDynamo to carry out a QC/MM study.

pDynamo can, of course, also be employed directly. To do this, a number of force field parameter sets are included with the 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, which lists each atom's label, the type of component in which it occurs, together with any link and variant modifications. This procedure 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.