sbmOpenMM Documentation

Introduction

Structure Based Models (SBMs) are representations of protein systems based on simplifications made over classical Molecular Dynamics (MD) force fields. Their are based on the energy landscape theory of protein folding and the principle of minimal frustration. The models maintain protein structures by focusing on chemical contacts formed at the native protein configuration, ignoring other non-native contacts. This allows for simpler force field definitions which capture essential protein dynamics at a much lower computational expense than traditional MD simulations.

sbmOpenMM is a Python library that offers flexibility to set up SBMs using the MD framework of OpenMM toolkit. It automates the creation of openmm.system classes that contain the necessary force field parameters to run molecular dynamics simulations using a protein structure and a contact map as the only necessary inputs.

sbmOpenMM is divided in three main classes:

1. geometry

2. models

3. system

The first class, geometry, contains methods to calculate the geometrical parameters from the input structures. These parameters are used to define the input conformation as the global minimum configuration in the potential energy function. The second class, models, allows to easily set up predefined SBM models, that encompass coarse grained, all atom and multi basin potentials. The third class, system, is the main class that holds all the methods to define, modify and create SBMs to be simulated with OpenMM.

The library is open-source and offers flexibility to create custom SBMs or to modify the predefined topology based models included in it.

SBM Models

The models class of sbmOpenMM contains three methods for automatic setting up predefined SBM potentials. It works by initializing a system class with the necessary force field parameters, derived from the input files, to set up one of the possible models which are detailed next:

Coarse grained, alpha-carbon (CA), model

The coarse grained method represents the protein system as beads centered at the alpha carbons of each residue in the protein. It uses harmonic potentials to hold the covalent connectivity and geometry of the beads. Torsional geometries are modeled with a periodic torsion potential. Native contacts are represented through the use of Lennard-Jones potentials that allow to form and break non-bonded interactions, permitting complete and local unfolding of the structures.

To create a CA model, call:

sbmOpenMM.models.getCAModel(pdb_file, contacts_file)

Here, pdb_file is the path to the PDB format structure of the protein and contacts_file is the path to the contact file containing only the CA atoms of the system. This last file should be numbered considering the CA atoms consecutively.

The force field equations are:

\[H_A = \sum_{bonds}V_{bond}+\sum_{angles}V_{angle}+\sum_{torsions}V_{torsion}+\sum_{contacts}V_{LJ_{12-10}}+\sum_{non-contacts}V_{LJ_{12}}\]

\[V_{bond} = \frac{k_b}{2}(r-r_0)^2\]

\[V_{angle} = \frac{k_a}{2}(\theta-\theta_0)^2\]

\[V_{torsion} = k_t(1-cos(\phi-\phi_0))+\frac{1}{2}(1-cos(3(\phi-\phi_0))))\]

\[V_{LJ_{12-10}} = \epsilon_{c}(5(\frac{\sigma_{ij}}{r})^{12}-6(\frac{\sigma_{ij}}{r})^{10})\]

\[V_{LJ_{12}} = \epsilon_{nc}(\frac{\sigma_{nc}}{r})^{12}\]

Here the default values are \(k_b=20000\ kJ/(mol \cdot nm^2)\), \(k_a=40\ kJ/(mol \cdot rad^2)\), \(k_t=1.0\ kJ/mol\), \(\epsilon_{c}=1.0\ kJ/mol\), \(\epsilon_{nc}=1.0\ kJ/mol\) and \(\sigma_{nc}=0.4\ nm\). The geometric parameters are set to the calculated structural values in the input structure, with \(r_0\) the equilibrium bond distance in nanometers, \(\theta_0\) the equilibrium angle length in radians, \(\phi_0\) the equilibrium torsional angle in radians and \(\sigma_{ij}\) the equilibrium contact distance in nanometers. The variable \(r\) represents, accordingly, the current bond or (non)contact distance in nanometers, \(\theta\) the current angle length in radians and \(\phi\) the current torsional angle in radians.

It is possible to use a \(V_{LJ_{12-10-6}}\) potential for the native contact interactions, defined as:

\[V_{LJ_{12-10-6}} = \epsilon_{c}(13(\frac{\sigma_{ij}}{r})^{12}-18(\frac{\sigma_{ij}}{r})^{10}+4(\frac{\sigma_{ij}}{r})^{6})\]

This potential gives a small energy barrier for contact formation/breaking that emulates a “desolvation effect”. To use this potential as the native contact energy function, instead of the \(V_{LJ_{12-10}}\) potential, give the option contact_force =’12-10-6’ to the sbmOpenMM.models.getCAModel() method.

Note that even if the units for the force constants are given in real physical units (e.g. \(kJ/mol\)), this is just to match the variables used by OpenMM. The models are not parametrized to equate this real physical values and comparison with experiments will require further adjustment to the energy unit system.

All-heavy-atoms (AA) model

The all-atom model represents the protein system with all its heavy atoms (i.e. excluding hydrogens). It uses harmonic potentials to hold the covalent connectivity, geometry and chirality of the protein residues. Periodic torsional potentials are used to maintain dihedral geometries of backbones and side chains. Native contacts are represented through the use of Lennard-Jones potentials that allow to form and break non-bonded interactions, permitting complete and local unfolding of the structures.

The method to create an AA model is:

sbmOpenMM.models.getAllAtomModel(pdb_file, contacts_file)

Here, pdb_file is the path to the PDB format structure of the protein and contacts_file is the path to the contact file containing only the non-hydorgen atoms of the protein system.

The force field equations are:

\[H_A = \sum_{bonds}V_{bond}+\sum_{angles}V_{angle}+\sum_{torsions}V_{torsion}+\sum_{impropers}V_{improper}+\sum_{planars}V_{planar}+\sum_{contacts}V_{LJ_{12-10}}+\sum_{non-contacts}V_{LJ_{12}}\]

\[V_{bond} = \frac{k_b}{2}(r-r_0)^2\]

\[V_{angle} = \frac{k_a}{2}(\theta-\theta_0)^2\]

\[V_{torsion} = k_t(1-cos(\phi-\phi_0))+\frac{1}{2}(1-cos(3(\phi-\phi_0))))\]

\[V_{improper} = \frac{k_i}{2}(\chi-\chi_{0})^2\]

\[V_{planar} = \frac{k_p}{2}(\chi-\chi_{0})^2\]

\[V_{LJ_{12-10}} = \epsilon_{c}(5(\frac{\sigma_{ij}}{r})^{12}-6(\frac{\sigma_{ij}}{r})^{10})\]

\[V_{LJ_{12}} = \epsilon_{nc}(\frac{\sigma_{nc}}{r})^{12}\]

Here the default values are \(k_b=10000\ kJ/(mol \cdot nm^2)\), \(k_a=80\ kJ/(mol \cdot rad^2)\), \(k_i=10.0\ kJ/(mol \cdot rad^2)\), \(k_p=20.0\ kJ/(mol \cdot rad^2)\), \(\epsilon_{nc}=0.1\ kJ/mol\) and \(\sigma_{nc}=0.25\ nm\). The values of the torsional \(k_t\) and native energy constant \(\epsilon_{c}\) are assigned by the following equations:

\[k_t=N_{atoms}/3N_{torsions}\ (kJ/mol)\]

\[k_c=2N_{atoms}/3N_{contacts}\ (kJ/mol)\]

Where \(N_{atoms}\) is the total number of atoms in the system, \(N_{torsions}\) is the total number of proper torsions assigned by the forcefield and \(N_{contacts}\) is the number of native contacts in the contact file definition. Additionally, the torsional energy constant \(k_t\) is further divided by classifying the torsions into backbone and sidechain groups. The assignment is carried out as:

\[k_{t}^{bb}=2k_t/3\]

\[k_{t}^{sc}=k_t/3\]

Here, \(k_{t}^{bb}\) and \(k_{t}^{sc}\) are the torsional energy constant for backbone and sidechain torsion groups, respectively. This grouping of torsions into backbone and side chains is the default behaviour of the sbmOpenMM.models.getAllAtomModel() method. It can be disabled by given the option group_by_bb_and_sc=False.

The geometric parameters are set to the calculated structural values in the input structure, with \(r_0\) the equilibrium bond distance in nanometers, \(\theta_0\) the equilibrium angle length in radians, \(\phi_0\) the equilibrium torsional angle in radians, \(\chi_0\) the equilibrium improper or planar equilibrium angle in radians and \(\sigma_{ij}\) the equilibrium contact distance in nanometers. The variable \(r\) represents, accordingly, the current bond or (non)contact distance in nanometers, \(\theta\) the current angle length in radians, \(\phi\) the current proper torsional angle in radians and \(\chi\) the equilibrium improper or planar torsional angles in radians.

Note that even if the units for the force constants are given in real physical units (e.g. \(kJ/mol\)), this is just to match the variables used by OpenMM. The models are not parametrized to equate this real physical values and comparison with experiments will require further adjustment to the energy unit system.

Multi basin model

The multi basin model automates the creation of a dual basin native contact potential. It receives as input two sbmOpenMM system classes, either CA or AA models, containing two different definitions of native contacts. One of the configurations is defined as the main model and the other is considered as the alternate model. All forcefield and topology parameters, different than the native contacts, are passed from the main configuration into the multi basin model. Then, the contacts are compared between the input configurations to define the sets of common and unique contacts. Common contacts with equilibrium length distances that differ more than a threshold are defined as dual basin and are assigned a special non-bonded Gaussian potential. The rest of the contacts are treated as single minima and are modeled with a Lennard-Jones (default) or a single basin Gaussian potential.

The multi basin Gaussian potential is defined as:

\[V_{Multi-basin} = \epsilon_{C}((1+(\frac{r_{ex}}{r})^{12})\prod_{minima}G(r,r_{0}^{\alpha})-1)\]

\[G(r,r_{0}^{\alpha}) = 1-exp(\frac{-(r-r_{0}^{\alpha})^2}{2\sigma^2})\]

\[\sigma^{2} = \frac{(r_{0}^{\alpha})^2}{50ln(2)}\]

Here, \(\epsilon_{C}\) is the native contact energy constant inherited from the main configuration, \(r_{ex}\) is the contact excluded volume radius, \(r_{0}^{\alpha}\) is the equilibrium distance for the \(alpha\)-th configuration and \(r\) is the current contact distance. \(\sigma\) is a parameter that modulates the well amplitude of the \(V_{Multi-basin}\) energy function. The single and double basin gaussian potential are distinguished by the number of \(r_{0}^{\alpha}\) parameters given.

The Lennard Jones contact potential is inherited accordingly from the CA or AA models used to build the multi basin SBM.

The method to create a multi basin model is:

sbmOpenMM.models.getMultiBasinModel(main_model, alternate_configuration=alternate_model)

Here, main_model and alternate_model are initialized sbmOpenMM system classes containing full force field parameter definitions.

Core classes

sbmOpenMM.geometry

core package of the sbmOpenMM package that contains the main sbmOpenMM classes.

The sbmOpenMM.core package contains the three sbmOpenMM main classes:

1. geometry

2. models

3. system

The first class, geometry, contains methods to calculate the geometrical parameters from the input structures. These parameters are used to define the input conformation as the global minimum configuration in the potential energy function. The second class, models, allows to easily set up predefined SBM models, that encompass coarse grained, all atom and multi basin potentials. The third class, system, is the main class that holds all the methods to define, modify and create SBMs to be simulated with OpenMM.

class geometry

A class to hold methods for calculating geometrical values given sets of atom coordinates.

angle(coord1, coord2, coord3)

Calculate the angle length between three (x,y,z) quantity coordinates.

Parameters

coord1simtk.unit.quantity.Quantity array

Vector for the first coordinate.

coord2simtk.unit.quantity.Quantity array

Vector for the second coordinate.

coord3simtk.unit.quantity.Quantity array

Vector for the third coordinate.

Returns

simtk.unit.quantity.Quantity

Quantity (value and unit) of the angle length in radians.

bond(coord1, coord2)

Calculate the distance length between two (x,y,z) quantity coordinates.

Parameters

coord1simtk.unit.quantity.Quantity array

Vector for the first coordinate.

coord2simtk.unit.quantity.Quantity array

Vector for the second coordinate.

Returns

simtk.unit.quantity.Quantity

Quantity (value and unit) of the distance length in nanometers.

position2Array(position, output_unit)

Converts an OpenMM position object quantity into a numpy array.

Parameters

positionsimtk.unit.quantity.Quantity

Array containing quantity objects [e.g. (x,y,z) array returned from positions].

output_unitsimtk.unit.unit.Unit

Unit in which to return the items of the array.

Returns

numpy.ndarray

A numpy array containing the quantity values converted to floats.

torsion(coord1, coord2, coord3, coord4)

Calculate the torsion angle length between four (x,y,z) quantity coordinates.

Parameters

coord1simtk.unit.quantity.Quantity array

Vector for the first coordinate.

coord2simtk.unit.quantity.Quantity array

Vector for the second coordinate.

coord3simtk.unit.quantity.Quantity array

Vector for the third coordinate.

coord4simtk.unit.quantity.Quantity array

Vector for the fourth coordinate.

Returns

simtk.unit.quantity.Quantity

Quantity (value and unit) of the torsion length in radians.

sbmOpenMM.models

core package of the sbmOpenMM package that contains the main sbmOpenMM classes.

The sbmOpenMM.core package contains the three sbmOpenMM main classes:

1. geometry

2. models

3. system

The first class, geometry, contains methods to calculate the geometrical parameters from the input structures. These parameters are used to define the input conformation as the global minimum configuration in the potential energy function. The second class, models, allows to easily set up predefined SBM models, that encompass coarse grained, all atom and multi basin potentials. The third class, system, is the main class that holds all the methods to define, modify and create SBMs to be simulated with OpenMM.

class models

A class to hold functions for the automated generation of default SBM models.

Methods

getAllAtomModel(structure_file, contact_file, kwarg)**	Creates an all atom SBM system class object with default initialized parameters.
getCAModel(structure_file, contact_file, kwarg)**	Creates an alpha-carbon only sbmOpenMM system class object with default initialized parameters.
getMultiBasinModel(main_model, alternate_configuration, kwarg)**	Creates a multi basin model from two sbmOpenMM.system class instances.

getAllAtomModel(structure_file, contact_file, default_parameters=True, default_forces=True, group_by_bb_and_sc=True, create_system=True, minimise=False, masses_per_element=False, radii_per_atom_type=False, ff_radii='amber', forcefield_file=None)

Initialises a default full-heavy-atom sbmOpenMM system class from a structure and a contact file defining the native contacts for the model. The system creation steps are:

1. Add the geometrical parameters for the model.

2. Add the default force field parameters for the model.

3. Create the default force objects.

4. Create the OpenMM system class.

The method can be used to generate an initialized sbmOpenMM system class, that only contains the geometrical parameters, by passing the option default_parameters as False. This is useful to store the geometrical values of bonds, angles, dihedrals, etc. in order to add custom parameters and forces.

The method can also be created without the initialisation of the forces classes by setting default_forces to False. This allows to load the default forcefield parameters and to modified them before creating the OpenMM force objects.

The method can also be stopped before creating the OpenMM system class using create_system=False.

Finally, a forcefield file can be given in order to read the forcefield parameters from it. You can give None to the contact file argument if you provide a forcefield file. Contacts definitions will be extracted from this file even if you give a path to a contact file.

Parameters

structure_filestring

Path to the input structure file.

contact_filestring

Path to the input native contact file. The file can be an output from the SMOG program (4 column contact file) or a two column file defining the atoms to be paired.

default_parametersboolean (True)

Whether to add default SBM All Atom forcefield parameters to the model.

default_forcesboolean (True)

Whether to initilize default SBM All Atom force objects. Set to False if the parameters will be different from the default ones.

group_by_bb_and_scboolean (True)

Wether to classify the torsions into backbone and side-chain to partition the torsional energy for each torsion in the system.

create_systemboolean (True)

If True the function will call the createSystemObject() method to create an OpenMM system object. If modifications to the default forcefield are necessary this option should be given False.

minimiseboolean (False)

Whether to minimise the system (with default options) if large forces are found.

masses_per_elementboolean (False)

Assign mass of the atoms according to their element.

radii_per_atom_typeboolean (False)

Assing radii per oplsaa atom type.

ff_radiistr (amber)

Which forcefield parameters to use to assign atom radii (amber or opslaa). Use together with radii_per_atom_type=True option.

forcefield_filestring

Path to the input forcefield file.

Returns

sbmsbmOpenMM.system

Initializes a sbmOpenMM.system class with default options for defining an All Atom SBM force field.

getCAModel(structure_file, contact_file, default_parameters=True, default_forces=True, create_system=True, contact_force='12-10', torsion_energy=1.0, contact_energy=1.0, minimise=False, residue_masses=False, residue_radii=False, forcefield_file=None)

Initialises a coarse-grained, carbon alpha (CA), sbmOpenMM system class from a structure and a contact file defining the native contacts for the coarse grained model.

The system creation steps are:

1. Add the geometrical parameters for the model.

2. Add the default force field parameters for the model.

3. Create the default force objects.

4. Create the OpenMM system class.

The method can be used to generate an initialized sbmOpenMM system class, that only contains the geometrical parameters, by passing the option default_parameters as False. This is useful to store the geometrical values of bonds, angles, dihedrals, etc. in order to add custom parameters and forces.

The method can also be created without the initialisation of the forces classes by setting default_forces to False. This allows to load the default forcefield parameters and to modified them before creating the OpenMM force objects.

The method can also be stopped before creating the OpenMM system class using create_system=False.

Finally, a forcefield file can be given in order to read the forcefield parameters from it. You can give None to the contact file argument if you provide a forcefield file. Contacts definitions will be extracted from this file even if you give a path to a contact file.

Parameters

structure_filestring

Path to the input structure file.

contact_filestring

Path to the input native contact file. The file can be an output from the SMOG program (4 column contact file) or a two column file defining the atoms to be paired.

default_parametersboolean (True)

Whether to add default SBM CA forcefield parameters to the model.

default_forcesboolean (True)

Whether to initilize default SBM CA force objects. Set to False if the parameters will be different from the default ones.

create_systemboolean (True)

If True the function will call the createSystemObject() method to create an OpenMM system object. If modifications to the default forcefield are necessary this option should be given False.

residue_massesboolean (False)

Set each alpha carbon atom mass to its average aminoacid residue mass.

residue_radiiboolean (False)

Set each alpha carbon atom radius to its statistical aminoacid residue radius.

forcefield_filestring

Path to the input forcefield file.

Returns

sbmsbmOpenMM.system

Initialized sbmOpenMM.system class with default options for defining a coarse-grained CA SBM force field.

getMultiBasinModel(main_model, alternate_configuration, double_minima_threshold=0.05, excluded_volume_radius=None, use_lennard_jones=True, create_system=True)

Generate a multibasin model from two sbmOpenMM initialized system classes. It currently supports only two configurations.

The two configurations (main and alternate) are compared at the level of native contacts to define common and unique contacts. If any common contact is defined to have equilibrium distances significantly different between the configurations, a multi basin Gaussian potential is added to explicitely consider both equilibrium distances in the non-bonded energy function. Unique contacts from the alternate configuration are added to the main configuration. All other bonded parameters are maintained from the main configuration.

Optionally, an excluded volume term can be given, with the excluded_volume_radius option, to control separately the sphere radius from the equilibrium contact distances.

The set of single-basin (or unique) contacts are added, by default, as Lennard-Jones potentials (use_lennard_jones=True) or can be added as Gaussian terms (use_lennard_jones=False) with separate control of the excluded volume term.

Finally, the method can be stopped before creating the OpenMM system class using create_system=False.

Parameters

main_modelsbmOpenMM.system

Configuration upon which the multi basin forces should be added.

alternate_configurationsbmOpenMM.system

Second configuration to define as a new minima in the native contact force object.

double_minima_thresholdfloat

The minimum equilibrium distance difference to consider a native contact having two different equilibrium distances. If shared contacts between the configurations have equilibrium distances that do not differ more than this value, only the main configuration equilibrium distance will be kept.

excluded_volume_radiusfloat

Radius of the excluded volume term in the Gaussian native-contact energy function.

use_lennard_jonesboolean (True)

Whether to use or not Lennard-Jones potential for the single-basin native contact set. If False a Gaussian function with separate control of the excluded volume term.

create_systemboolean (True)

If True the function will call the createSystemObject() method to create an OpenMM system object.

Returns

sbmsbmOpenMM.system

Initialized sbmOpenMM.system class with added multi basin contact potential.

Attributes

common_contactsset

Set containing the common contacts between configurations.

dual_basin_contactsset

Set containing the subset of common contacts that have a dual basin potential

unique_contactsdict

Dictionary containing the sets of unique contacts for each configuration. The keys are integers, 0 for the main configuration, 1 for the alternate configuration.

sbmOpenMM.system

core package of the sbmOpenMM package that contains the main sbmOpenMM classes.

The sbmOpenMM.core package contains the three sbmOpenMM main classes:

1. geometry

2. models

3. system

The first class, geometry, contains methods to calculate the geometrical parameters from the input structures. These parameters are used to define the input conformation as the global minimum configuration in the potential energy function. The second class, models, allows to easily set up predefined SBM models, that encompass coarse grained, all atom and multi basin potentials. The third class, system, is the main class that holds all the methods to define, modify and create SBMs to be simulated with OpenMM.

class system(structure_path, particles_mass=1.0)

A class containing methods and parameters for generating Structure Based Models (SBM) systems to be simulated using the OpenMM interface. It offers flexibility to create default and custom SBM systems and to easily modify their parameters.

Attributes

structure_pathstring

Path to the pdb or cif input file

structureopenmm.app.pdbfile.PDBFile or openmm.app.pdbxfile.PDBxFile

Object that holds the information of OpenMM PDB or CIF parsing methods.

topologyopenmm.app.topology.Topology

OpenMM topology of the model.

positionsunit.quantity.Quantity

Atomic positions of the model.

particles_massfloat or list

Mass of each particle. If float then uniform masses are given to all particles. If list per-particle masses are assigned.

model_typestring

String representing the model type: All-atom (AA), alpha-carbon (CA) or multi-basin variants (AA-MB, CA-MB).

atomslist

A list of the current atoms in the model. The items are simtk.openmm.app.topology.Atom initialised classes.

n_atomsint

Total numer of atoms in the model.

bondscollections.OrderedDict

A dict that uses bonds (2-tuple of simtk.openmm.app.topology.Atom objects) present in the model as keys and their forcefield properties as values.

bonds_indexeslist

A list containing the zero-based indexes of the atoms defining the bonds in the model.

n_bondsint

Total numer of bonds in the model.

anglescollections.OrderedDict

A dict that uses angles (3-tuple of simtk.openmm.app.topology.Atom objects) present in the model as keys and their forcefield properties as values.

angles_indexeslist

A list containing the zero-based indexes of the atoms defining the angles in the model.

n_anglesint

Total numer of angles in the model.

torsionscollections.OrderedDict

A dict that uses proper torsions (4-tuple of simtk.openmm.app.topology.Atom objects) present in the model as keys and their forcefield properties as values.

torions_indexeslist

A list containing the zero-based indexes of the atoms defining the torsions in the model.

n_torsionsint

Total numer of proper torsions in the model.

improperscollections.OrderedDict

A dict that uses improper torsions (4-tuple of simtk.openmm.app.topology.Atom objects) present in the model as keys and their forcefield properties as values.

impropers_indexeslist

A list containing the zero-based indexes of the atoms defining the imporpers in the model.

n_impropersint

Total numer of improper torsions in the model.

planarscollections.OrderedDict

A dict that uses planar torsions (4-tuple of simtk.openmm.app.topology.Atom objects) present in the model as keys and their forcefield properties as values.

planars_indexeslist

A list containing the zero-based indexes of the atoms defining the planars in the model.

n_planarsint

Total numer of planar torsions in the model.

contactscollections.OrderedDict

A dict that uses native contacts (2-tuple of simtk.openmm.app.topology.Atom objects present in the model as keys and their forcefield properties as values.

contacts_indexeslist

A list containing the zero-based indexes of the atoms defining the contacts in the model.

n_contactsint

Total numer of native contacts in the model.

torsions_groupdict

A dict that uses proper torsions two central atoms (2-tuple of simtk.openmm.app.topology.Atom objects) present in the model’s topology as keys and the number of torsions (int) that share these same middle bond atoms as values.

torsions_typedict

A dict that uses proper torsions (4-tuple of simtk.openmm.app.topology.Atom objects) present in the model as keys and, a string representing wheter the torsion is classified as ‘backbone’ or ‘sidechain’ as values.

energy_constantdict

A dict that holds the value for the different energy terms parameters used by different SBM models and forces.

harmonicBondForceopenmm.HarmonicBondForce

Stores the OpenMM HarmonicBondForce initialised-class. Implements an harmonic bond potential between pairs of particles, that depends quadratically on their distance.

harmonicAngleForceopenmm.HarmonicAngleForce

Stores the OpenMM HarmonicAngleForce initialised-class. Implements an harmonic angle potential between trios of particles, that depends quadratically on their angle length.

periodicTorsionForceopenmm.CustomTorsionForce

Stores the OpenMM CustomTorsionForce initialised-class. Implements a force potential that varies periodically with the value of the proper torsion angle.

generalPeriodicTorsionForceopenmm.CustomTorsionForce

Stores the OpenMM CustomTorsionForce initialised-class. Implements a general force potential that varies periodically with the value of the proper torsion angle.

harmonicImproperForceopenmm.CustomTorsionForce

Stores the OpenMM CustomTorsionForce initialised-class. Implements a force potential that varies quadratically with the value of the improper torsion angle.

harmonicPlanarForceopenmm.CustomTorsionForce

Stores the OpenMM CustomTorsionForce initialised-class. Implements a force potential that varies quadratically with the value of the planar torsion angle.

lj12_6contactForceopenmm.CustomBondForce

Stores the OpenMM CustomBondForce initialised-class to be applied to non-bonded interactions between native contact pairs. Implements a lennard-jones potential with exponents 12 and 6 for the repulsive and attractive componenets, respectively.

lj12_10contactForceopenmm.CustomBondForce

Stores the OpenMM CustomBondForce initialised-class to be applied to non-bonded interactions between native contact pairs. Implements a lennard-jones potential with exponents 12 and 10 for the repulsive and attractive components, respectively.

lj12_10_6contactForceopenmm.CustomBondForce

Stores the OpenMM CustomBondForce initialised-class to be applied to non-bonded interactions between native contact pairs. Implements a lennard-jones potential with exponents 12 and 10 for the repulsive and attractive components, respectively, and an additional 6-exponent term to model a “desolvation penalty” for forming/breaking contacts.

singleGaussianContactForceopenmm.CustomNonbondedForce

Stores the OpenMM CustomBondForce initialised-class to be applied to non-bonded interactions between native contact pairs with single minimum. Implements a mixed lennard-jones (repulsive) and gaussian potential (attractive) with separate control of equilibrium distance and excluded volume.

doubleGaussianContactForceopenmm.CustomNonbondedForce

Stores the OpenMM CustomBondForce initialised-class to be applied to non-bonded interactions between native contact pairs with two minima. Implements a mixed lennard-jones (repulsive) and gaussian potential (attractive) with separate control of the equilibrium distances and excluded volume.

ljRepulsionForceopenmm.CustomNonbondedForce

Stores the OpenMM CustomBondForce initialised-class to be applied to the non-bonded interactions between non-native contact pairs. Implements only a repulsive lennard-jones potential with exponent 12.

forceGroupscollections.OrderedDict

A dict that uses force names as keys and their corresponding force as values.

systemopenmm.System

Stores the OpenMM System initialised class. It stores all the forcefield information for the SBM model.

rf_epsilonfloat

Epsilon parameter used in the repulsion force object.

rf_sigmafloat

Sigma parameter used in the repulsion force object.

rf_cutofffloat

Cutoff value used for repulsion force interactions.

exclusionslist

List of added exclusions for the repuslion non-bonded term.

Methods

removeHydrogens()	Remove hydrogens from the input structure by using a regexpression pattern. Used specially for creating all atom (AA) models.
getCAlphaOnly()	Filter in only alpha carbon atoms from the input structure and updates the topology object to add new bonds between them. Used specially for creating alpha-carbon (CA) corse-grained models.
getAtoms()	Reads atoms from topology, adds them to the main class and sorts them into a dictionary to store their forcefield properties.
getBonds()	Reads bonds from topology, adds them to the main class and sorts them into a dictionary to store their forcefield properties.
getAngles()	Reads bonds from topology and creates a list of all possible angles, adds them to the main class and sorts them into a dictionary to store their forcefield properties.
getProperTorsions()	Using the created angles, usually by the getAngles() method, creates a list of all possible proper torsion dihedral angles, filtering out torsions based on residue-specific rules (only for the all-atom model). The torsions are then added to the main class and sorted into a dictionary to store their forcefield properties.
getImpropers()	Create improper torsions based on backbone and sidechain residue-specific rules (all-atom model only), adds them to the main class and sorts them into a dictionary to store their forcefield properties.
getPlanars()	Create planar torsions based on backbone and sidechain residue-specific rules (all-atom model only), adds them to the main class and sorts them into a dictionary to store their forcefield properties.
readContactFile()	Reads a file containing native contact information and adds it into the main class. The file can be smog-style (4 columns) or given as 2 columns. The format will be automatically detected.
setBondParameters()	Change the forcefield parameters for bonded terms.
setAngleParameters()	Change the forcefield parameters for angle terms.
setProperTorsionParameters()	Change the forcefield parameters for proper torsion terms.
setImproperParameters()	Change the forcefield parameters for improper torsion terms.
setPlanarParameters()	Change the forcefield parameters for planar torsion terms.
setNativeContactParameters()	Change the forcefield parameters for native contact terms.
setParticlesMasses()	Change the mass parameter for each atom in the system.
setParticlesRadii()	Change the excluded volume radius parameter for each atom in the system.
addHarmonicBondForces()	Creates an harmonic bonded force term for each bond in the main class using their defined forcefield parameters.
addHarmonicAngleForces()	Creates an harmonic angle force term for each angle in the main class using their defined forcefield parameters.
addPeriodicTorsionForces()	Creates a periodic torsion force term for each proper torsion in the main class using their defined forcefield parameters.
addGenericPeriodicTorsionForces()	Creates a periodic torsion force term for each proper torsion in the main class using their defined forcefield parameters.
addHarmonicImproperForces()	Creates an harmonic torsion force term for each improper torsion in the main class using their defined forcefield parameters. Used specially for simulating All Atom systems.
addHarmonicPlanarForces()	Creates an harmonic torsion force term for each planar torsion in the main class using their defined forcefield parameters. Used specially for simulating All Atom systems.
addLJ12_6ContactForces()	Creates a 12-6 Lennard-Jones bond potential for each native contact in the main class using their defined forcefield parameters. Used specially for simulating All Atom systems.
addLJ12_10ContactForces()	Creates a 12-10 Lennard-Jones bond potential for each native contact in the main class using their defined forcefield parameters. Used specially for simulating coarse-grained alpha-carbon systems.
addGaussianContactForces()	Creates a gaussian single and double basin bond potential for each native contact in the main class using their defined forcefield parameters. The contacts are recognized according two the number of parameters given as values in the attribute system. 3-parameters for single-basin potential and 4-paramters for dual-basin potentials.
addLJRepulsionForces()	Creates a repulsive-only 12 Lennard-Jones non-bonded potential specifying a exclusion list for bond, angle, torsion, and native contact terms.
groupTorsionsbyBBAndSC()	Groups proper torsions by backbone and sidechain torsions. Used exclusively for simulating all-atom SBM.
getAATorsionParameters()	Generates default periodic torsion forcefield parameters, for proper torsions, using pre-defined assignment schemes. Used exclusively for simulating all-atom SBM.
getAANativeContactParameters()	Generates default bonded contact forcefield parameters, for native contacts, using pre-defined assignment schemes. Used exclusively for simulating all-atom SBM.
createSystemObject()	Creates OpenMM system object adding particles, masses and forces. It also groups the added forces into Force-Groups for the sbmReporter class.
addParticles()	Add particles to the system OpenMM class instance.
addSystemForces()	Add forces to the system OpenMM class instance. It also save names for the added forces to include them in the reporter class.
dumpStructure()	Writes a structure file of the system in its current state.
dumpForceFieldData()	Writes to a file the parameters of the SBM forcefield.
loadForcefieldFromFile()	Loads forcefield parameters from a sbmOpenMM force field file written with the dumpForceFieldData() method.
setCAMassPerResidueType()	Sets alpha carbon atoms to their average residue mass. Used specially for modifying alpha-carbon (CA) corse-grained models.
setCARadiusPerResidueType()	Sets alpha carbon atoms to their average residue mass. Used specially for modifying alpha-carbon (CA) corse-grained models.
setAAMassPerAtomType()	Sets every atom mass to its element mass. Used specially for all-atom models.
setAARadiusPerAtomType()	Sets every atom radius to its element mass. Used specially for all-atom models.

addGaussianContactForces(self)

Creates an openmm.CustomBondForce() object with the equilibrium distances and parameters setted up in the “contacts” attribute. Two forces are created according to the number of parameters given:

If three parameters are given the custom bond force is initilized with the formula:

energy = epsilon*((1+(r0/r)^12)*(1-exp(-25*log(2)*(r-r0)^2/(r0*r0)))-1)

The force object is stored at the “singleGaussianContactForce” attribute.

The force parameters must be contained in self.contacts as follows: self.contacts is a dictionary:

• The keys are tuples for tww atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> rex (float)

  • second -> r0 (quantity)

  • third -> epsilon (float)

If four parameters are given the custom bond force is initilized with the formula:

energy = epsilon*((1+(r0/r)^12)*(1-exp(-25*log(2)*(r-r0)^2/(r0*r0)))*

(1-exp(-25*log(2)*(r-r1)^2/(r1*r1)))-1)

The force object is stored at the “doubleGaussianContactForce” attribute.

The force parameters must be contained in self.contacts as follows: self.contacts is a dictionary:

• The keys are tuples for two atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> rex (float)

  • second -> r0 (quantity)

  • second -> r1 (quantity)

  • third -> epsilon (float)

Parameters

None

Returns

None

addGeneralPeriodicTorsionForces(self)

Creates an openmm.CustomTorsionForce() object with the torsions and parameters setted up in the “torsions” attribute. The custom torsion force is initilized with the formula:

energy = k*(1+cos(n*(theta-theta0)))

The force object is stored at the “generalPeriodicTorsionForce” attribute.

The force parameters must be contained in self.torsions as follows: self.torsions is a dictionary:

• The keys are tuples for four atoms in self.topology.atoms attribute.

• The values are a list with the parameters as items:

• Each item is a tuple containing:

  • first -> theta0 (quantity)

  • second -> k (float)

  • third -> n (int)

Parameters

None

Returns

None

addHarmonicAngleForces(self)

Creates an openmm.HarmonicAngleForce() object with the angles and parameters setted up in the “angles” attribute. The force object is stored at the “harmonicAngleForce” attribute.

The force parameters must be contained in self.angles as follows: self.angles is dictionary:

• The keys are 3-tuples for three atoms in self.topology.atoms attribute.

• The values are a 2-tuple of parameters in the following order:

  • first -> angle0 (quantity)

  • second -> k (float)

Parameters

None

Returns

None

addHarmonicBondForces(self)

Creates an openmm.HarmonicBondForce() object with the bonds and parameters setted up in the “bonds” dictionary attribute. The force object is stored at the “harmonicBondForce” attribute.

The force parameters must be contained in self.bonds as follows: self.bonds is a dictionary:

• The keys are 2-tuples for two atom items in self.topology.atoms attribute.

• The values are a 2-tuple of parameters in the following order:

  • first -> bond0 (quantity)

  • second -> k (float)

Parameters

None

Returns

None

addHarmonicImproperForces(self)

Creates an openmm.CustomTorsionForce() object with the torsions and parameters setted up in the “impropers” attribute. The custom torsion force is initilized with the formula:

energy = 0.5*k*dtheta_torus^2 dtheta_torus = dtheta - floor(dtheta/(2*pi)+0.5)*(2*pi) dtheta = theta - theta0 pi = np.pi

The force object is stored at the “harmonicImproperForce” attribute.

The force parameters must be contained in self.torsions as follows: self.torsions is a dictionary:

• The keys are tuples for four atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> theta0 (quantity)

  • second -> k (float)

Parameters

None

Returns

None

addHarmonicPlanarForces(self)

Creates an openmm.CustomTorsionForce() object with the torsions and parameters setted up in the “planars” attribute. The custom torsion force is initilized with the formula:

energy = 0.5*k*dtheta_torus^2 dtheta_torus = dtheta - floor(dtheta/(2*pi)+0.5)*(2*pi) dtheta = theta - theta0 pi = np.pi

The force object is stored at the “harmonicPlanarForce” attribute.

The force parameters must be contained in self.torsions as follows: self.torsions is a dictionary:

• The keys are tuples for four atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> theta0 (quantity)

  • second -> k (float)

Parameters

None

Returns

None

addLJ12_10ContactForces(self)

Creates an openmm.CustomBondForce() object with the bonds and parameters setted up in the “contacts” attribute. The custom bond force is initilized with the formula:

energy = epsilon*(5*(sigma/r)^12-6*(sigma/r)^10)

The force object is stored at the “lj12_10contactForce” attribute.

The force parameters must be contained in self.contacts as follows: self.contacts is a dictionary:

• The keys are tuples for two atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> sigma (quantity)

  • second -> epsilon (float)

Parameters

None

Returns

None

addLJ12_10_6ContactForces(self)

Creates an openmm.CustomBondForce() object with the bonds and parameters setted up in the “contacts” attribute. The custom bond force is initilized with the formula:

energy = epsilon*(13*(sigma/r)^12-18*(sigma/r)^10+4*(sigma/r)^6)

The force object is stored at the “lj12_10_6contactForce” attribute.

The force parameters must be contained in self.contacts as follows: self.contacts is an ordered dictionary:

• The keys are tuples for two atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> sigma (quantity)

  • second -> epsilon (float)

Parameters

None

Returns

None

addLJ12_6ContactForces(self)

Creates an openmm.CustomBondForce() object with the bonds and parameters setted up in the “contacts” attribute. The custom bond force is initilized with the formula:

energy = epsilon*((sigma/r)^12-2*(sigma/r)^6)

The force object is stored at the “lj12_6contactForce” attribute.

The force parameters must be contained in self.contacts as follows: self.contacts is a dictionary:

• The keys are tuples for two atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> sigma (quantity)

  • second -> epsilon (float)

Parameters

None

Returns

None

addLJRepulsionForces(self, cutoff=None, bonded_exclusions_index=3, exclusions=None)

Creates an openmm.CustomNonbondedForce() object with the parameters sigma and epsilon given to this method. The custom non-bonded force is initilized with the formula:

energy = ‘epsilon*(sigma/r)^12; sigma=0.5*(sigma1+sigma2)’

The method adds exclusions for bonded atoms up until ‘bonded_exclusions_index’ bond orders (default 3) and also for all the pairs defined in the ‘contacts’ attribute. The force object is stored at the “ljRepulsionForce” attribute.

Parameters

epsilonfloat

Value of the epsilon constant in the energy function.

sigmafloat or list

Value of the sigma constant (in nm) in the energy function. If float the same sigma value is used for every particle. If list a unique parameter is given for each particle.

cutofffloat

The cutoff distance (in nm) being used for the nonbonded interactions.

bonded_exclusions_indexint

Specified number of bonds to consider when adding non-bonded exclusions.

Returns

None

addParticles(self)

Add a particle to the sbmOpenMM system for each atom in it. The mass of each particle is set up with the values in the ‘particles_mass’ attribute.

Parameters

None

Returns

None

addPeriodicTorsionForces(self)

Creates an openmm.CustomTorsionForce() object with the torsions and parameters setted up in the “torsions” attribute. The custom torsion force is initilized with the formula:

energy = k*(1-cos(theta-theta0)+0.5*(1-cos(3*(theta-theta0))))

The force object is stored at the “periodicTorsionForce” attribute.

The force parameters must be contained in self.torsions as follows: self.torsions is a dictionary:

• The keys are tuples for four atoms in self.topology.atoms attribute.

• The values are a tuple of parameters in the following order:

  • first -> theta0 (quantity)

  • second -> k (float)

Parameters

None

Returns

None

addSystemForces(self)

Adds generated forces to the sbmOpenMM system, also adding a force group to the ‘forceGroups’ attribute dictionary.

Parameters

None

Returns

None

checkBondDistances(self, threshold=0.25)

Searches for large bond distances for the atom pairs defined in the ‘bonds’ attribute. It raises an error when large bonds are found.

Parameters

thresholdfloat

Treshold to check for large bond distances.

Returns

None

checkLargeForces(self, threshold=1, minimise=False)

Prints the SBM system energies of the input configuration of the system. It optionally checks for large forces acting upon all particles in the SBM system and iteratively minimises the system configuration until no forces larger than a threshold are found.

Parameters

thresholdfloat

Treshold to check for large forces.

minimisefloat

Whether to iteratively minimise the system until all forces are lower or equal to the thresshold value.

Returns

None

createSystemObject(self, check_bond_distances=True, minimise=False, force_threshold=10, bond_threshold=0.25)

Creates an openmm.System() object using the force field parameters given to the SBM ‘system’ class. It adds particles, forces and creates a force group for each force object. Optionally the method can check for large bond distances (default) and minimise the atomic positions if large forces are found in any atom (default False).

Parameters

minimiseboolean

Wheter to minimise the system if large forces are found.

check_bond_distancesboolean

Wheter to check for large bond distances.

force_thresholdfloat

Treshold to check for large forces.

bond_thresholdfloat

Treshold to check for large bond distances.

Returns

None

dumpForceFieldData(self, output_file)

Writes a file containing the current forcefield parameters in the sbmOpenMM system.

Parameters

output_filestring

name of the output file.

Returns

None

dumpStructure(self, output_file)

Writes a PDB file containing the currently defined abmOpenMM system atoms.

Parameters

output_filestring

name of the PDB output file.

Returns

None

getAANativeContactParameters(self)

Calculates the contact force field parameters by distributing the total native contact energy equally among all the contacts.

The energy values are stored in the ‘energy_constant’ attribute, which contains the entry ‘C’ for the partitioned energy values.

Parameters

None

Returns

contact_parameterslist

List of the native-contact force constant parameters for each native contact in the SBM system.

getAATorsionParameters(self, group_by_bb_and_sc=True)

Calculates the torsion force field parameters by distributing the total torsional energy equally among all the torsions. Optionally, uses torsional classifications, in backbone or side-chain torsions, to group the energy partition. In the last case considers that BB torsion energies are twice the value of sidechain torsion energies.

The energy partioned values are stored in the ‘energy_constant’ attribute. It contains the entries ‘torsion’ for the equally partitioned torsion energies, or ‘BB’ and ‘SC’ for group-partioned torsion energies.

Parameters

group_by_bb_and_scboolean

Wheter to partion the torsional energy in backbone and side-chain groups.

Returns

torsion_parameterslist

List of the torsion force constant parameters for each torsion in the SBM system.

getAngles(self)

Adds angles to the sbmOpenMM system class based on the bondings in the topology object.

Parameters

None

Returns

None

getAtoms(self)

Adds atoms in the OpenMM topology instance to the sbmOpenMM system class.

Parameters

None

Returns

None

getBonds(self)

Adds bonds in the OpenMM topology instance to the sbmOpenMM system class.

Parameters

None

Returns

None

getCAlphaOnly(self)

Keeps in the SBM system only the alpha carbon atoms from the OpenMM topology.

Parameters

None

Returns

None

getImpropers(self, chains=None)

Adds improper torsions to the sbmOpenMM system class to maintain backbone chiralities.

Parameters

chainslist (None)

If this parameter is given, only add impropers if the chain id is in this list

Returns

None

getPlanars(self, chains=None)

Adds planar torsions to the sbmOpenMM system class to mantain side chain and backbone planar arrangements.

Parameters

chainslist (None)

If this parameter is given, only add planars if the chain id is in this list

Returns

None

getProperTorsions(self)

Adds proper torsions to the sbmOpenMM system class based on the bonded angles in it. Excludes special torsions for rings and backbones.

Parameters

None

Returns

None

groupTorsionsbyBBAndSC(self)

Classifies the torsions into backbone or sidechain, based on the atoms present in each torsion definition.

The classification is stored in the ‘torsions_type’ dictionary attribute.

Parameters

None

Returns

None

loadForcefieldFromFile(self, forcefield_file)

Loads force field parameters from a force field file written by the dumpForceFieldData() method into the sbmOpenMM system.

Parameters

forcefield_filestring

path to the force field file.

Returns

None

readContactFile(self, contact_file, shift=1)

Reads a file to add native contact information to the sbmOpenMM system class. The file format can be ‘smog’ like (i.e. 4 columns) or ‘2column’ like (i.e. one column for each atom). This is auto-dected by the method.

Parameters

None

Returns

None

removeHydrogens(self)

Removes all hydrogen atoms in the topology from the SBM system.

Parameters

None

Returns

None

setAAMassPerAtomType(self)

Sets the masses of the atoms to the mass its corresponding element.

Parameters

None

Returns

None

setAARadiusPerAtomType(self, ff_radii='amber')

Sets the excluded volume radii of the atoms to the sigma value of the oplsaa or amber (default) forcefield.

Parameters

None

Returns

None

setAngleParameters(self, angle_parameters)

Set the harmonic angle constant force parameters. The input can be a float, to set the same parameter for all force interactions, or a list, to define a unique parameter for each force interaction.

Parameters

angle_parametersfloat or list

Parameter(s) to set up for the harmonic angle forces.

Returns

None

setBondParameters(self, bond_parameters)

Set the harmonic bond constant force parameters. The input can be a float, to set the same parameter for all force interactions, or a list, to define a unique parameter for each force interaction.

Parameters

bond_parametersfloat or list

Parameter(s) to set up for the harmonic bond forces.

Returns

None

setCAMassPerResidueType(self)

Sets the masses of the alpha carbon atoms to the average mass of its amino acid residue.

Parameters

None

Returns

None

setCARadiusPerResidueType(self)

Sets the excluded volume radii of the alpha carbon atoms to characteristic radii of their corresponding amino acid residue.

Parameters

None

Returns

None

setImproperParameters(self, improper_parameters)

Set the harmonic torsion constant force parameters for impropers. The input can be a float, to set the same parameter for all force interactions, or a list, to define a unique paremater for each force interaction.

Parameters

improper_parametersfloat or list

Parameter(s) to set up for the harmonic torsion force parameters for impropers.

Returns

None

setNativeContactParameters(self, contact_parameters)

Set the native interactions constant force parameters for native contacts. The input can be a float, to set the same parameter for all force interactions, or a list, to define a unique parameter for each force interaction.

Parameters

contact_parametersfloat or list

Parameter(s) to set up for the native interaction force for native contacts.

Returns

None

setParticlesMasses(self, particles_mass)

Set the masses of the particles in the system. The input can be a float, to set the same mass for all particles, or a list, to define a unique mass for each particle.

Parameters

particles_massfloat or list

Mass(es) values to add for the particles in the sbmOpenMM system class.

Returns

None

setParticlesRadii(self, particles_radii)

Set the radii of the particles in the system. The input can be a float, to set the same radius for all particles, or a list, to define a unique radius for each particle.

Parameters

particles_radiifloat or list

Radii values to add for the particles in the sbmOpenMM system class.

Returns

None

setPlanarParameters(self, planar_parameters)

Set the harmonic torsion constant force parameters for planars. The input can be a float, to set the same parameter for all force interactions, or a list, to define a unique parameter for each force interaction.

Parameters

planar_parametersfloat or list

Parameter(s) to set up for the harmonic torsion force parameters for planars.

Returns

None

setProperTorsionParameters(self, torsion_parameters)

Set the periodic torsion constant force parameters. The input can be a float, to set the same parameter for all force interactions, or a list, to define a unique parameter for each force interaction.

Parameters

torsion_parametersfloat or list

Parameter(s) to set up for the periodic torsion forces.

Returns

None

Reporter class

sbmOpenMM.sbmReporter

reporter package of the sbmOpenMM package that contains the sbmReporter class.

The sbmOpenMM.reproter package contains the sbmReporter class.

sbmReporter is a special class of the OpenMM StateDataReporter class, that additionally accepts a sbmobject to print the SBM forcefield energies.

class sbmReporter(file, reportInterval, sbmObject=None, **kwargs)

A special case of the StateDataReporter class that outputs information about a simulation, such as energy and temperature, etc. to a file. This special reporter outputs the sbmOpenMM force group energies inside the sbmOpenMM system object.

It is used in the same way as the OpenMM StateDataReporter class, but it takes as additional input an instance of the sbmOpenMM object with the option ‘sbmObject’.