System

A class containing methods and parameters for generating TOPO coarse-grained systems to be simulated with OpenMM. It is typically constructed via topo.models.buildCoarseGrainModel(), which sets bonds, angles, torsions, Yukawa electrostatics, and structure-based non-bonded forces.

See also

The TOPO model: theory and force field for the physics of each force these methods add, and Using TOPO from Python for how to build and inspect a system from Python.

class topo.core.system(structure_path: str, model: str = 'topo')[source]

A class for generating coarse-grained (CG) systems that can be simulated using the OpenMM interface. It allows for the creation of both default and custom CG systems and easy modification of their parameters.

Parameters:

  • structure_path (str) – The path to the input PDB or CIF file. Required.

  • model (str, optional, default='topo') – The model to be used. Currently only ‘topo’ (topology-based for folded proteins) is supported.

Variables:

  • structure (mm.app.pdbfile.PDBFile or mm.app.pdbxfile.PDBxFile) – An object that holds the information of the parsed PDB or CIF file.

  • topology (mm.app.topology.Topology) – The topology of the model, as generated by OpenMM.

  • positions (unit.quantity.Quantity) – The atomic positions of the model.

  • particles_mass (float or list) – The mass of each particle. If a float is provided, all particles will have the same mass. If a list is provided, per-particle masses will be assigned.

  • particles_charge (list) – The charge of each particle.

  • particle_rmin_2 (float) – Per-particle collision radius Rmin/2 (nm) – the bead’s excluded-volume (vdW) radius. Feeds only dumpForceFieldData(); the contact force itself uses the per-pair rmin_matrix.

  • atoms (list) – A list of the current atoms in the model. The items are mm.app.topology.atoms initialised classes.

  • n_atoms (int) – The total number of atoms in the model.

  • bonds (collections.OrderedDict) – A dictionary that uses bonds (2-tuple of mm.app.topology.bonds objects) present in the model as keys and their forcefield properties as values.

  • bonds_indexes (list) – A list containing the zero-based indexes of the atoms defining the bonds in the model.

  • n_bonds (int) – The total number of bonds in the model.

  • bonded_exclusions_index (int) – The exclusion rule for nonbonded force. =2 for topo model (excludes 1-2 and 1-3 bonded neighbours, preserving native contacts at sequence separation abs(i - j) >= 3)

  • harmonicBondForce (mm.HarmonicBondForce) – The mm.HarmonicBondForce object that implements a harmonic bond potential between pairs of particles, that depends quadratically on their distance.

  • n_angles (int) – The total number of angles in the model.

  • gaussianAngleForce (mm.CustomAngleForce) – The mm.CustomAngleForce object that implements a Gaussian angle bond potential between pairs of three particles.

  • n_torsions (int) – Total number of torsion angles in the model.

  • gaussianTorsionForce (mm.CustomTorsionForce) – Stores the OpenMM CustomTorsionForce initialised-class. Implements a Gaussian torsion angle bond potential between pairs of four particles.

  • yukawaForce (mm.CustomNonbondedForce) – Stores the OpenMM CustomNonbondedForce initialized-class. Implements the Debye-Huckle potential.

  • custom_non_bonded_force (mm.CustomNonbondedForce) – Stores the OpenMM CustomNonbondedForce initialized-class. Implements the structure-based (native + non-native) contact potential.

  • forceGroups (collections.OrderedDict) – A dict that uses force names as keys and their corresponding force as values.

  • system (mm.System) – Stores the OpenMM System initialised class. It stores all the forcefield information for the topo model.

loadForcefieldFromFile()

Loads forcefield parameters from a force field file written with the dumpForceFieldData() method.

__init__(structure_path: str, model: str = 'topo')[source]

Initialises the TOPO OpenMM system class.

Parameters:

  • structure_path (string [requires]) – Name of the input PDB or CIF file

  • model (str [optional, default='topo']) – Model name. Currently only ‘topo’ is supported.

Return type:

None

getCAlphaOnly() → None[source]

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) coarse-grained models.

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

Return type:

None

Notes

TODO: check that residue indices are consecutive, e.g.:

np_indices = np.array(atom_indices)
if len(np_indices) > 1:
    if (np_indices[1:] - np_indices[:-1]) SOME_CONDITION_TO_CHECK_FOR_CONSECUTIVE:
        warnings.warn('atom indices are not monotonically increasing')
    if len(np.unique(np_indices)) < len(np_indices):
        warnings.warn('atom_indices are not unique')
getAtoms()[source]

Reads atoms from topology, adds them to the main class and sorts them into a dictionary to store their forcefield properties.

After getCAlphaOnly, C-alpha atoms are stored on self.topology only. We need to add them to atoms attribute and system also. Adds atoms in the OpenMM topology instance to the TOPO system class.

Return type:

None

getBonds(except_chains=None)[source]

Reads bonds from topology, adds them to the main class and sorts them into a dictionary to store their forcefield properties.

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

Parameters:

except_chains (String [optional])

Return type:

None

getAngles()[source]

Reads bonds from topology, adds them to the main class and sorts them into a dictionary to store their forcefield properties.

Angles are built from bond (bondedTo instance). This function modify self.angles, self.angles_indexes and self.n_angles

getTorsions()[source]

Reads bonds from topology, adds them to the main class and sorts them into a dictionary to store their forcefield properties.

Torsion Angles are built from angles and bondedTo instance. This function modify self.torsions, self.torsions_indexes and self.n_torsions

setBondForceConstants() → None[source]

Change the forcefield parameters for bonded terms.

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.

Return type:

None

setParticlesMass(particles_mass)[source]

Change the mass parameter for each atom in the system.

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_mass (float or list) – Mass(es) values to add for the particles in the TOPO system class.

Return type:

None

setParticlesRadii(particles_radii)[source]

Set each particle’s collision radius Rmin/2 (stored in particle_rmin_2).

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. These per-particle Rmin/2 values feed only dumpForceFieldData(); the contact force uses the per-pair rmin_matrix.

Parameters:

particles_radii (float or list) – Per-particle collision radius Rmin/2 (nm) for the TOPO system class.

Return type:

None

setParticlesCharge(particles_charge)[source]

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

Parameters:

particles_charge (float or list) – Charge values to add for the particles in the TOPO system class.

Return type:

None

addHarmonicBondForces() → None[source]

Creates a harmonic bonded force term for each bond in the main class using their defined forcefield parameters.

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

openMM uses harmonic bond that has energy term of form:

\[E= \frac{1}{2}k(r-r_0)^2\]

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) (measured in unit of kj/mol/nm^2)

Return type:

None

addGaussianAngleForces() → None[source]

Add Gaussian functional form of angle. Note that in openMM log is neutral logarithm.

Angle potential takes a bimodal (double-Gaussian, log-sum-exp) form:

\[U_\mathrm{angle}(\theta) = -\frac{1}{\gamma} \ln\left[ e^{-\gamma\left[k_\alpha(\theta-\theta_\alpha)^2 + \epsilon_\alpha\right]} + e^{-\gamma k_\beta(\theta-\theta_\beta)^2} \right]\]

with basins at \(\theta_\alpha = 91.7^\circ\) and \(\theta_\beta = 130.0^\circ\). See the model-theory documentation for all constants.

addPeriodicTorsionForce() → None[source]

Add the sequence-dependent periodic torsion potential.

Each backbone dihedral uses a periodic torsion with four periodicities:

\[U_\mathrm{torsion}(\varphi) = \sum_{n=1}^{4} k_{D,n}\left[1 + \cos(n\varphi - \delta_n)\right]\]

The force constants \(k_{D,n}\) and phases \(\delta_n\) depend on the two central residues of the dihedral and are read from topo/parameters/data/dihedral_params.csv (scaled by a 0.756 calibration factor; see topo.parameters.dihedral.load_dihedral_params()).

addYukawaForces(use_pbc: bool) → None[source]

Creates a nonbonded force term for electrostatic interaction DH potential.

Creates an mm.CustomNonbondedForce() object implementing the screened-Coulomb (Debye-Huckel / Yukawa) potential:

\[U^\mathrm{el}_{ij}(r) = f\,\frac{q_i q_j}{\epsilon_r\, r}\, e^{-r/l_D}\]

where \(f = 1/4\pi\epsilon_0 = 138.935458\ \mathrm{kJ\,nm\,mol^{-1}\,e^{-2}}\) is the Coulomb constant in MD units, \(\epsilon_r = 78.5\) is the relative dielectric of water, and \(l_D = 1.0\ \mathrm{nm}\) is the Debye screening length (about 100 mM monovalent salt). A 2.0 nm cutoff with a switching function at 1.8 nm is used.

The force object is stored at the yukawaForce attribute.

Parameters:

use_pbc ((bool) whether use PBC, cutoff periodic boundary condition)

Return type:

None

addCustomNonBondedForce(rmin_matrix, energy_matrix, use_pbc)[source]

Add the structure-based (native + non-native) contact non-bonded force.

Creates an mm.CustomNonbondedForce() implementing a 12-10-6 (Go-type) pairwise potential, with per-pair well position \(R_{ij}\) and well depth \(\varepsilon_{ij}\) supplied as tabulated functions of the two particles’ ids:

\[U^\mathrm{nb}_{ij}(r) = \varepsilon_{ij}\left[ 13\left(\frac{R_{ij}}{r}\right)^{12} - 18\left(\frac{R_{ij}}{r}\right)^{10} + 4\left(\frac{R_{ij}}{r}\right)^{6} \right]\]

The well minimum is at \(r = R_{ij}\), where \(U^\mathrm{nb}_{ij} = -\varepsilon_{ij}\), so \(R_{ij}\) is the preferred distance and \(\varepsilon_{ij}\) is the well depth. Both matrices are produced by topo.utils.nonbonded.build_nonbonded_interaction():

  • Native contacts (residue pairs in contact in the input structure): \(R_{ij}\) is the native Ca-Ca distance and \(\varepsilon_{ij}\) is the sum of hydrogen-bond, backbone-sidechain and (domain-scaled) sidechain-sidechain energies.

  • Non-native pairs: a soft excluded-volume repulsion – a negligible well depth with \(R_{ij} = R_{\mathrm{min}/2,i} + R_{\mathrm{min}/2,j}\) (sum rule).

A 2.0 nm cutoff with a switching function at 1.8 nm is used, and pairs two or fewer bonds apart (1-2 and 1-3 neighbours) are excluded. The force object is stored at the custom_non_bonded_force attribute.

Parameters:

  • rmin_matrix (numpy.ndarray) – Symmetric (n_atoms x n_atoms) matrix of well positions \(R_{ij}\) (Rmin) in nm. Looked up per pair via R_table(id1, id2).

  • energy_matrix (numpy.ndarray) – Symmetric (n_atoms x n_atoms) matrix of well depths \(\varepsilon_{ij}\) in kJ/mol. Looked up via eps_table(id1, id2).

  • use_pbc (bool) – If True, use CutoffPeriodic; otherwise CutoffNonPeriodic.

Return type:

None

createSystemObject(check_bond_distances: bool = True, minimize: bool = False, check_large_forces: bool = True, force_threshold: float = 10.0, bond_threshold: float = 0.5) → None[source]

Creates OpenMM system object adding particles, masses and forces. It also groups the added forces into Force-Groups for the topoReporter class.

Creates an mm.System() object using the force field parameters given to the ‘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 minimize the atomic positions if large forces are found in any atom (default False).

Parameters:

  • minimize (boolean (False)) – Whether to minimize the system if large forces are found.

  • check_bond_distances (boolean (True)) – Whether to check for large bond distances.

  • check_large_forces (boolean (False)) – Whether to print force summary of force groups

  • force_threshold (float (10.0)) – Threshold to check for large forces.

  • bond_threshold (float (0.5)) – Threshold to check for large bond distances.

Return type:

None

checkBondDistances(threshold: float = 0.5) → None[source]

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

The comparison is strict (>): a bond exactly at threshold is allowed, so the 0.5 nm default (the nucleic equilibrium bond length, bond_length_nucleic) does not spuriously reject a valid RNA/nucleic structure; only bonds strictly longer than threshold raise.

Parameters:

threshold ((float, default=0.5 nm)) – Threshold to check for large bond distances.

Return type:

None

checkLargeForces(minimize: bool = False, threshold: float = 10) → None[source]

Prints the topo system energies of the input configuration of the system and checks for large forces acting upon its particles.

With minimize=True the system is iteratively minimized until no force larger than threshold remains (updating self.positions). With minimize=False the structure is left untouched, but if any force still exceeds threshold a RuntimeWarning is emitted so an unstable / clashing input is not passed over silently.

Parameters:

  • threshold ((float, default=10)) – Threshold to check for large forces (kJ/mol/nm).

  • minimize ((bool, default= False)) – Whether to iteratively minimize the system until all forces are lower or equal to the threshold value. When False, forces above threshold only trigger a warning (no minimization).

Return type:

None

addParticles() → None[source]

Add particles to the system OpenMM class instance.

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

addSystemForces() → None[source]

Add forces to the system OpenMM class instance. It also save names for the added forces to include them in the reporter class.

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

dumpStructure(output_file: str) → None[source]

Writes a structure file of the system in its current state.

Writes a PDB file containing the currently defined CG system atoms and its positions.

Parameters:

output_file (string) – name of the PDB output file.

Return type:

None

dumpTopology(output_file: str) → None[source]

Writes a topology file of the system in PSF format, this is used for visualization and post-analysis.

Writes a file containing the current topology in the TOPO system. This file contains topology of system, used in visualization and analysis.

Here, we used parmed to load openMM topology and openMM system to create Structure object in parmed. Because parmed doesn’t automatically recognize charge, mass of atoms by their name. We need to set charge, mass back to residues properties.

Parameters:

output_file (string [requires]) – name of the output PSF file.

Return type:

None

dumpForceFieldData(output_file: str) → None[source]

Writes to a file the parameters of the forcefield.

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

Parameters:

output_file (string [requires]) – name of the output file.

Return type:

None

setCAMassPerResidueType()[source]

Sets alpha carbon atoms to their average residue mass. Used specially for modifying alpha-carbon (CA) coarse-grained models.

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

Return type:

None

setCARadiusPerResidueType() → None[source]

Set each CA bead’s excluded-volume radius to the fixed per-AA Rmin_2 from model_parameters[self.model].

Warning

These are the fixed per-type radii used for rigid ribosome scenery beads. A mobile protein chain (nascent chain, folded-protein sim) must not use them – its Rmin/2 is the per-residue, structure-derived Karanicolas-Brooks value (build_nonbonded_interaction(return_rmin_2=True)), set via setParticlesRadii(). This method is therefore no longer called by the model builders; it is retained only for scenery/diagnostic use.

Return type:

None

setCAChargePerResidueType() → None[source]

Sets the charge of the alpha carbon atoms to characteristic charge of their corresponding amino acid residue.

Return type:

None

static _setParameters(term, parameters)[source]

General function to set up or change force field parameters. protected method, can be called only inside class system.

Parameters:

  • term (dict) – Dictionary object containing the set of degrees of freedom (DOF) to set up attributes to (e.g. bonds attribute)

  • parameters (integer or float or list) – Value(s) for the specific forcefield parameters. If integer or float, sets up the same value for all the DOF in terms. If a list is given, sets a unique parameter for each DOF.

Return type:

None