af2rave.simulation module

Utilities for creating simulation box

class SimulationBox(filename: str, forcefield: ForceField = <openmm.app.forcefield.ForceField object>)[source]

Bases: object

Create a simulation box from a protein PDB box.

Parameters:
  • filename (str) – path to the pdb file

  • forcefield (OpenMM.app.ForceField) – forcefield to be used for adding hydrogens

add_disulfide_bond(resid: list[tuple[int, int]])[source]

By default the simulation box after create_box() does not have any disulfide bonds. This method will add back the bonds when necessary. This method currently only works with AMBER forcefields.

Parameters:

resid (list[tuple[int, int]]) – list of tuples of CYS residues that are bonded

create_box(**kwargs) tuple[list, Topology][source]

Generate the simulation box from a raw pdb file. Currently only soluble proteins are supported as we can only add water.

This function performs the following tasks: 1. use pdbfixer to add missing atoms, residues, and terminals 2. add hydrogen, at the given pH 3. solvate the system with water, neutralize the system by adding ions

Parameters:
  • pH (float) – float: pH of the system. Default is 7.0

  • padding (float) – padding around the protein. Default is 10. Unit: Angstrom.

  • water_model (str) – water model to be used. Default is ‘tip3p’

  • positiveIon (str) – positive ion used to neutralize the system. Default is ‘Na+’.

  • negativeIon (str) – negative ion used to neutralize the system. Default is ‘Cl-’

  • ionicStrength (float) – ionic strength of the system. Default is 0.0. Unit: molar

property map_atom_index

Map atom index from input PDB to the output PDB file. This is a callable property.

After adding hydrogen, the atom index will be changed. This function translates the old atom index to the new atom index.

Example:

box = SimulationBox(filename)
box.create_box()
box.map_atom_index(1)
Parameters:

index (AtomIndexLike: int or set[int] or list[int] or list[set[int]]) – atom index or a list of atom indices

save_pdb(filename: str)[source]

Write the simulation box to a pdb file.

Parameters:

filename (str) – path to the output pdb file

class TopologyMap(old_top: Topology | Topology | str, new_top: Topology | Topology | str, resid_offset=0)[source]

Bases: object

Create a mapping between the old and new topology.

Parameters:
  • old_top (openmm.app.Topology or mdtraj.Topology or str) – old topology, either a path or an OpenMM or MDTraj topology object

  • new_top (openmm.app.Topology or mdtraj.Topology or str) – new topology, either a path or an OpenMM or MDTraj topology object

  • resid_offset – Offset to be added to the residue id to align two topologies. Default is 0. For example, if a residue in the old topology has resid 1 and in the new topology has resid 101, then resid_offset = 100.

map_atom_index(index: int | set[int] | tuple[int] | list[int] | list[set[int]] | list[tuple[int]]) int | set[int] | tuple[int] | list[int] | list[set[int]] | list[tuple[int]][source]

Map atom index (0-based) from input PDB to the output PDB file.

After adding hydrogen, the atom index will be changed. This function translates the old atom index to the new atom index.

Parameters:

index (AtomIndexLike: int or set[int] or list[int] or list[set[int]]) – atom index or a list of atom indices

Tools for unbiased simulation

The UnbiasedSimulation class for running MD simulations.

class UnbiasedSimulation(pdb_file, **kwargs)[source]

Bases: object

UnbiasedSimulation class for running MD. An example use is:

import af2rave.simulation as af2sim
sim = af2sim.UnbiasedSimulation(<arguments>)
sim.run(steps)

These are the parameters the initialization function accepts:

Parameters:
  • pdb_file (str) – Path to OpenMM.app.pdbfile.PDBFile object

  • list_of_index (list[tuple[int, int]]) – List of indexes to calculate the CVs

  • forcefield (OpenMM.app.ForceField) – OpenMM.app.ForceField object. Default: Charmm36mFF

  • temp (float | openmm.unit.Quantity[unit=kelvin]) – Temperature of the system. Default: 310 K

  • pressure (float | openmm.unit.Quantity[unit=bar]) – Pressure of the system. Default: 1 Bar. Can be set to none to disable pressure coupling (NVT).

  • dt (float | openmm.unit.Quantity[unit=picoseconds]) – Time step of the integrator. Default: 0.002 ps

  • cutoff (float | openmm.unit.Quantity[unit=angstroms]) – Nonbonded cutoff. Default: 10.0 Angstroms

  • steps (int) – Simulation steps. Default: 50 million steps (100 ns)

  • cv_file (str) – File to write CVs to. Default: COLVAR.dat

  • cv_freq (int) – Frequency of CVs written to the file. Default: 50

  • xtc_file (str) – Trajectory file name. Default: traj.xtc

  • xtc_freq (int) – Frequency writing trajectory in steps. Default: 1000

  • xtc_reporter (mdtraj.reporters.XTCReporter) – XTCReporter object. This overrides traj_filename and traj_freq. Default: None

  • custom_forces – List of custom forces to be added to the system. Default: None

  • append (bool) – Appends to existing file. Default: False

  • progress_every (int) – Progress report to stdout frequency. Default: 1000 Can be set to None to suppress this.

property pos: list[Vec3]

Get the positions of the simulation.

Returns:

List of positions

Return type:

list[openmm.Vec3]

restart(steps: int = 50000000, restart_file: str = None) Simulation[source]

Run the simulation from given pdb file. Default: 50 million steps (100 ns).

Parameters:
  • steps (int) – Number of steps to run the simulation. Default: 50 million steps (100 ns)

  • restart_file (str) – Name of the checkpoint file. Default: {prefix}.chk

run(steps: int = 50000000) Simulation[source]

Run the simulation from given pdb file. Default: 50 million steps (100 ns).

Parameters:

steps (int) – Number of steps to run the simulation. Default: 50 million steps (100 ns)

save_checkpoint(filename: str = None) None[source]

Save the checkpoint of the simulation.

Parameters:

filename (str) – Name of the checkpoint file. Default: {prefix}_out.chk

save_pdb(filename: str = None) None[source]

Save the final state PDB of the simulation.

Parameters:

filename (str) – Name of the pdb file. Default: {prefix}_out.pdb

OpenMM reporters

This CVReporter generates a PLUMED-style COLVAR file of the features.

class CVReporter(file: str = 'COLVAR.dat', reportInterval: int = 100, list_of_indexes: list[tuple[int, int]] = None, append: bool = False)[source]

Bases: object

An OpenMM reporter that writes a PLUMED-style COLVAR file of the features. This reporter writes to the file every reportInterval steps. The first column is the number of steps instead of time. Distances are in the units of Angstorms.

This reporter calculates distances without a PBC. This is equivalent to NOPBC option in PLUMED. This should do the right thing most of the time.

Parameters:
  • file (str) – The name of the file to write the CVs to. Default: COLVAR.dat

  • reportInterval (int) – The interval at which to write the CVs. Default: 100

  • list_of_indexes (list[tuple[int, int]]) – The list of indexes to calculate the CVs. Default: None

  • append (bool) – Append to existing file

describeNextReport(simulation)[source]

OpenMM reporter method to describe the next report.

report(simulation, state)[source]

OpenMM reporter method to report the current state.

class MinimizationReporter(maxIter: int = 500)[source]

Bases: MinimizationReporter

report(iteration, x, grad, args)[source]

the report method is called every iteration of the minimization.

Args:
iteration (int): The index of the current iteration. This refers

to the current call to the L-BFGS optimizer. Each time the minimizer increases the restraint strength, the iteration index is reset to 0.

x (array-like): The current particle positions in flattened order:

the three coordinates of the first particle, then the three coordinates of the second particle, etc.

grad (array-like): The current gradient of the objective function

(potential energy plus restraint energy) with respect to the particle coordinates, in flattened order.

args (dict): Additional statistics about the current state of minimization.

In particular: “system energy”: the current potential energy of the system “restraint energy”: the energy of the harmonic restraints “restraint strength”: the force constant of the restraints (in kJ/mol/nm^2) “max constraint error”: the maximum relative error in the length of any constraint

Returns:

bool : Specify if minimization should be stopped.