Additional scripts to analyze the data
This is a wrapper above a GPUassisted molecular dynamics package Openmm.
You can find extensive description of openmm classes here: https://simtk.org/api_docs/openmm/api10/annotated.html
The main library is located in openmmlib.Simulation
A file knotAnalysis.py contains code needed to calculate Alexander’s polynomials, a measure of knot complexity. It also uses OpenMM to help unwrap polymer and reduce time needed to calculate Alexander’s polynomial
A file polymerScalings.py has some utilities to calculate Rg(s) and Pc(s) for polymer conformations in a fast and efficient way.
A file contactmaps.py has code to quickly calculate contacts within a polymer structure, and organize them as contactmaps. It is used by polymerScalings.
A file pymol_show has utilities to visualize polymer conformations using pymol
Polymer configuration is represented as a Nx3 numpy array of coordinates. Start/end of chains/rings are not directly specified in the file, and have to be added through method setLayout
Input file may have the simplistic format, described in txtToJoblib.py (first line with number of particles, then N lines with three floats corresponding to x,y,z coordinates each). Input file can also be any of the output files.
Output file format is a dictionary, saved with joblib.dump. Nx3 data array is stored under the key “data”. The rest of the dictionary consists of metadata, describing details of the simulation. This metadata is for informative purpose only, and is ignored by the code.
All forces of the system are added by the code to the self.ForceDict dictionary. After the force is added, the class (or user) can get it out of the self.ForceDict and modify the force. Once the simulation is started, all forces are automatically applied and cannot be modified.
Two types of bond forces are harmonic bond force and FENEtype bond as described in Grosberg papers. Individual bonds can be added using addBond, while polymer bonds can be added using addHarmonicPolymerBonds, etc.
Grosberg repulsive force  a faster and better implementation of repulsive force; and LennardJones Force, that can be attractive and allows to specify extra attraction/repulsion between any pairs of particles.
There are also polynomial repulsive and attractive forces which are faster due to a shorter cutoff radius.
Stiffness force can be harmonic, or the “special” Grosberg force, kept for compatibility with the systems used in Grosberg forces
External forces include spherical confinement, cylindrical confinement, attraction to “lamina” surface of a cylinder, gravity, etc.
A sample line of information printed on each step looks like this:
minim bl=5 (i) pos[1]=[99.2 52.1 52.4] dr=0.54 kin=3.65 pot=3.44 Rg=107.312 SPS=144:
Let’s go over it step by step.
minim  simulation name (sim default, minim  energy minimization. Other name can be provided in self.name).
bl=5  name of a current block
(i) indicate that velocity reinitialization was done at this step. You will simultaneously see that Ek is more than 2.4
pos[1] is a position of a first monomer
dr is sqrt(mean square displacement) of monomers, i.e. how much did a monomer shift on average.
kin=3.65 pot=3.44  energies: kinetic, potential
Rg=107.312  current radius of gyration (size of the system)
SPS  steps per second
Functions depend on each other, and have to be applied in certain groups
load — Mandatory
setup — Mandatory
setLayout — Mandatory, after self.load()
saveFolder — Optional (default is folder with the code)
self.add___Force() — Use any set of forces
Then go all the addBond, and other modifications and tweaks of forces.
Before running actual simulation, it is advised to resolve all possible conflict by doing energyMinimization
doBlock — the actual simulation
save — saves conformation
Select GPU (“0” or “1”)  setup
Select chain/ring  setLayout
Select timestep or collision rate  Simulation
Base class for openmm simulations
Methods
Sets up the important lowlevel parameters of the platform. Mandatory to run.
Parameters:  platform : string, optional
PBC : bool, optional
PBCbox : (float,float,float), optional
GPU : “0” or “1”, optional
integrator : “langevin”, “variableLangevin”, “verlet”, “variableVerlet”,
verbose : bool, optional
errorTol : float, optional


sets the folder where to save data.
Parameters:  folder : string


Sets configuration of the chains in the system. This information is later used by the chainforming methods, e.g. addHarmonicPolymerBonds() and addStiffness(). This method supersedes the less general getLayout().
Parameters:  chains: list of tuples :


loads data from file. Accepts text files, joblib files or pure data as Nx3 or 3xN array
If h5dictKey is specified, uses new h5dictbased I/O
Parameters:  filename : joblib file, or text file name, or Nx3 or 3xN numpy array
center : bool, optional
h5dictKey : int or str, optional
masses : array


Saves conformation plus some metadata. Metadata is not interpreted by this library, and is for your reference
If data is saved to the .vtf format, the same filename should be specified each time. .vtf format is then viewable by VMD.
Initializes an HDF5based storage for multiple conformation.
When the storage is initialized, save() by default saves to the storage.
If “r+” mode is chosen, simulation is automaticaly set up to continue from existing conformation. In that case, last step is automatically determined, and data from secondtolast file is loaded.
Note
Continue mode does not copy forces and layouts!
Parameters:  filename : str
mode : :


Sets particle positions
Parameters:  data : Nx3 arrayline


Sets up domains for the simulation. Also, pickles domain vector to “domains.dat”.
Parameters:  domains : boolean array or None
filename : str or None


Adds bond between two particles, allows to specify parameters
Parameters:  i,j : int
bondWiggleDistance : float
bondType : “Harmonic” or “Grosberg”
verbose : bool


Adds harmonic bonds connecting polymer chains wiggleDist controls the distance at which energy of the bond equals kT
Adds FENE bonds according to HalversonGrosberg paper. (Halverson, Jonathan D., et al. “Molecular dynamics simulation study of
nonconcatenated ring polymers in a melt. I. Statics.” The Journal of chemical physics 134 (2011): 204904.)
This method has a repulsive potential buildin, so that Grosberg bonds could be used with truncated potentials. Is of no use unless you really need to simulate Grosbergtype system.
Parameters:  k : float, optional


Adds harmonic angle bonds. k specifies energy in kT at one radian If k is an array, it has to be of the length N. Xth value then specifies stiffness of the angle centered at monomer number X. Values for ends of the chain will be simply ignored.
Parameters:  k : float or list of length N


Adds stiffness according to the Grosberg paper. (Halverson, Jonathan D., et al. “Molecular dynamics simulation study of
nonconcatenated ring polymers in a melt. I. Statics.” The Journal of chemical physics 134 (2011): 204904.)
Parameters are synchronized with normal stiffness
If k is an array, it has to be of the length N. Xth value then specifies stiffness of the angle centered at monomer number X. Values for ends of the chain will be simply ignored.
Parameters:  k : float or Nlong list of floats


Adds a special force which could be use for very efficient resolution of crossings Use this force if your monomers are all “on top of each other” E.g. if you start your simulations with fractional brownyan motion with h < 0.4
This is the fastest nontransparent repulsive force. Done according to the paper: (Halverson, Jonathan D., et al. “Molecular dynamics simulation study of
nonconcatenated ring polymers in a melt. I. Statics.” The Journal of chemical physics 134 (2011): 204904.)
This is a simple polynomial repulsive potential. It has the value of trunc at zero, stays flat until 0.60.7 and then drops to zero together with its first derivative at r=1.0.
Parameters:  trunc : float


This is a simple and fast polynomial force that looks like a smoothed version of the squarewell potential. The energy equals repulsionEnergy around r=0, stays flat until 0.60.7, then drops to zero together with its first derivative at r=1.0. After that it drop down to attractionEnergy and gets back to zero at r=`attractionRadius`.
The energy function is based on polynomials of 12th power. Both the function and its first derivative is continuous everywhere within its domain and they both get to zero at the boundary.
Parameters:  repulsionEnergy: float :
repulsionRadius: float :
attractionEnergy: float :
attractionEnergy: float :


This is almost the same potential as in addSmoothSquareWellTailedForce. The only difference is that the attractive part of the potential flattens out to the value of tailEnergy at r=`attractionRadius` and then goes quadratically to zero at tailRadius. Please, refer to the documentation for addSmoothSquareWellForce for the details of the repulsive and attractive parts of the potential.
Parameters:  kwargs: :
tailEnergy: :
tailRadius: :


Adds a lennardjones force, that allows for mutual attraction. This is the slowest force out of all repulsive.
Note
This is the only force that allows for socalled “exceptions’. Exceptions allow you to change parameters of the force for a specific pair of particles. This can be used to create shortrange attraction between pairs of particles. See manual for Openmm.NonbondedForce.addException.
Parameters:  cutoff : float, optional
domains : bool, optional
epsilonRep : float, optional
epsilonAttr : float, optional
blindFraction : float, 0<x<1
sigmaRep, sigmaAttr: float, optional :


A softened version of lennardJones force. Now we’re moving to polynomial forces, so go there instead.
used to exclude a bunch of particles from calculation of nonbonded force
Parameters:  particles : list


Adds attractive shortrange interaction of strength epsilon between particles i,j and a few neighboring particles requires :py:func:’LennardJones Force<Simulation.addLennardJonesForce>’
Parameters:  i,j : int
epsilon : float
sigma : float, optional
length : int, optional, default = 3


Constrain particles to be within a sphere. With no parameters creates sphere with density .3
Parameters:  r : float or “density”, optional
k : float, optional
density : float, optional, <1


Excludes particles from a sphere of radius r at certain position.
Attracts one domain to the lamina. Infers radius from spherical confinement, that has to be initialized already.
Parameters:  width : float, optional
depth : float, optional
r : float, optional


tethers particles in the ‘particles’ array. Increase k to tether them stronger, but watch the system!
Parameters:  particles : list of ints
k : int, optional


adds force pulling downwards in z direction When using cutoff, acts only when z>cutoff
Initializes particles velocities
Parameters:  mult : float, optional


Sends particle coordinates to OpenMM system. If system has exploded, this is
used in the code to reset coordinates.
Reinitializes the OpenMM context object. This should be called if lowlevel parameters, such as forces, have changed.
Parameters:  mult : float, optional


A wrapper to the buildin OpenMM Local Energy Minimization
Runs system at smaller timestep and higher collision rate to resolve possible conflicts.
Now we’re moving towards local energy minimization, this is here for backwards compatibility.
performs one block of simulations, doing steps timesteps, or steps_per_block if not specified.
Parameters:  steps : int or None
increment : bool, optional
num : int or None, optional


A subclass used to do simulations for the Metaphase project
Methods
Attracts a subset of particles to the core, repells the rest from the core
Limits position of a set of particles in z coordinate
Parameters:  particles : list
zCoordinates : list, or tuple of length 2
k : float, optional
useOtherAxis : “x”,”y” or “z”, optional


contains some experimental features
Methods
quickly loads a set of repulsive chains, possibly adds spherical confinement