Additional scripts to analyze the data

Openmm-lib - a wrapper around Openmm to use with polymer simulations


This is a wrapper above a GPU-assisted molecular dynamics package Openmm.

You can find extensive description of openmm classes here:

The main library is located in openmmlib.Simulation

A file 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 has some utilities to calculate Rg(s) and Pc(s) for polymer conformations in a fast and efficient way.

A file 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

Input/Output file format

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 (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.

Implemented forces

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 FENE-type bond as described in Grosberg papers. Individual bonds can be added using addBond, while polymer bonds can be added using addHarmonicPolymerBonds, etc.

Nonbonded (inter-monomer) force can be of many types. Three Lennard-Jones-based forces:
simple repulsife force U = 1/r^12;

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.

Information, printed on current step

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

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

Frequently-used settings - where to specify them?

Select GPU (“0” or “1”) - setup

Select chain/ring - setLayout

Select timestep or collision rate - Simulation

class openmmlib.Simulation(timestep=80, thermostat=0.001, temperature=Quantity(value=300, unit=kelvin), verbose=False, velocityReinitialize=True, name='sim', length_scale=1.0, mass_scale=1.0)[source]

Base class for openmm simulations


setup(platform='CUDA', PBC=False, PBCbox=None, GPU='default', integrator='langevin', verbose=True, errorTol=None, precision='mixed')[source]

Sets up the important low-level parameters of the platform. Mandatory to run.


platform : string, optional

Platform to use

PBC : bool, optional

Use periodic boundary conditions, default:False

PBCbox : (float,float,float), optional

Define size of the bounding box for PBC

GPU : “0” or “1”, optional

Switch to another GPU. Mostly unnecessary. Machines with 1 GPU automatically select right GPU. Machines with 2 GPUs select GPU that is less used.

integrator : “langevin”, “variableLangevin”, “verlet”, “variableVerlet”,

“brownian”, optional Integrator to use (see Openmm class reference)

verbose : bool, optional

Shout out loud about every change.

errorTol : float, optional

Error tolerance parameter for variableLangevin integrator Values of 0.03-0.1 are reasonable for “nice” simulation Simulations with strong forces may need 0.01 or less


sets the folder where to save data.


folder : string

folder to save the data

setChains(chains=[(0, None, 0)])[source]

Sets configuration of the chains in the system. This information is later used by the chain-forming methods, e.g. addHarmonicPolymerBonds() and addStiffness(). This method supersedes the less general getLayout().


chains: list of tuples :

The list of chains in format [(start, end, isRing)]. The particle range should be semi-open, i.e. a chain (0,3,0) links the particles 0, 1 and 2. If bool(isRing) is True than the first and the last particles of the chain are linked into a ring. The default value links all particles of the system into one chain.


returns configuration of chains

setLayout(mode='chain', chains=None, Nchains=1)[source]

Deprecated method. Use setChains instead


returns configuration of chains

load(filename, center=False, h5dictKey=None, masses=None)[source]

loads data from file. Accepts text files, joblib files or pure data as Nx3 or 3xN array

If h5dictKey is specified, uses new h5dict-based I/O


filename : joblib file, or text file name, or Nx3 or 3xN numpy array

Input filename or array with data

center : bool, optional

Move center of mass to zero before starting the simulation

h5dictKey : int or str, optional

Indicates that you need to load from an h5dict

masses : array

Masses of each atom, measured in self.mass (default: 100 AMU, but could be modified by self.mass_scale)

save(filename=None, mode='auto')[source]

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.

mode : str
“h5dict” : use build-in h5dict storage “joblib” : use joblib storage “txt” : use text file storage “vtf” : append data to the .vtf trajectory file “auto” : use h5dict if initialized, joblib otherwise
filename : str or None
Filename not needed for h5dict storage (use initStorage command) for joblib and txt, if filename not provided, it is created automatically.
initStorage(filename, mode='w-')[source]

Initializes an HDF5-based 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 second-to-last file is loaded.


Continue mode does not copy forces and layouts!


filename : str

Filename of an h5dict storage file

mode : :

‘w’ - Create file, truncate if exists ‘w-‘ - Create file, fail if exists (default) ‘r+’ - Continue simulation, file must exist.


Returns an Nx3 array of positions


Returns data, scaled back to PBC box


Sets particle positions


data : Nx3 array-line

Array of positions with distance ~1 between connected atoms.


Runs automatically to offset data if it is integer based

Returns:Gyration ratius in units of length (bondlength). :
Returns:Distance to the furthest from the origin particle. :
dist(i, j)[source]

Calculates distance between particles i and j

useDomains(domains=None, filename=None)[source]

Sets up domains for the simulation. Also, pickles domain vector to “domains.dat”.


domains : boolean array or None

N-long array with domain vector

filename : str or None

Filename with pickled domain vector

addBond(i, j, bondWiggleDistance=0.2, distance=None, bondType=None, verbose=None)[source]

Adds bond between two particles, allows to specify parameters


i,j : int

Particle numbers

bondWiggleDistance : float

Average displacement from the equilibrium bond distance

bondType : “Harmonic” or “Grosberg”

Type of bond. Distance and bondWiggleDistance can be specified for harmonic bonds only

verbose : bool

Set this to False if you’re in verbose mode and don’t want to print “bond added” message


Adds harmonic bonds connecting polymer chains wiggleDist controls the distance at which energy of the bond equals kT


Adds FENE bonds according to Halverson-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.)

This method has a repulsive potential build-in, so that Grosberg bonds could be used with truncated potentials. Is of no use unless you really need to simulate Grosberg-type system.


k : float, optional

Arbitrary parameter; default value as in Grosberg paper.


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.


k : float or list of length N

Stiffness of the bond. If list, then determines stiffness of the bond at monomer i. Potential is k * alpha^2 * 0.5 * kT


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.


k : float or N-long list of floats

Synchronized with regular stiffness. Default value is very flexible, as in Grosberg paper. Default value maximizes entanglement length.


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

addGrosbergRepulsiveForce(trunc=None, radiusMult=1.0)[source]

This is the fastest non-transparent 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.)
trunc : None or float
truncation energy in kT, used for chain crossing. Value of 1.5 yields frequent passing, 3 - average passing, 5 - rare passing.
addPolynomialRepulsiveForce(trunc=3.0, radiusMult=1.0)[source]

This is a simple polynomial repulsive potential. It has the value of trunc at zero, stays flat until 0.6-0.7 and then drops to zero together with its first derivative at r=1.0.


trunc : float

the energy value around r=0

addSmoothSquareWellForce(repulsionEnergy=3.0, repulsionRadius=1.0, attractionEnergy=0.5, attractionRadius=2.0)[source]

This is a simple and fast polynomial force that looks like a smoothed version of the square-well potential. The energy equals repulsionEnergy around r=0, stays flat until 0.6-0.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.


repulsionEnergy: float :

the heigth of the repulsive part of the potential. E(0) = repulsionEnergy

repulsionRadius: float :

the radius of the repulsive part of the potential. E(repulsionRadius) = 0, E’(repulsionRadius) = 0

attractionEnergy: float :

the depth of the attractive part of the potential. E(repulsionRadius/2 + attractionRadius/2) = attractionEnergy

attractionEnergy: float :

the maximal range of the attractive part of the potential.

addSmoothSquareWellTailedForce(repulsionEnergy=3.0, repulsionRadius=1.0, attractionEnergy=0.5, attractionRadius=2.0, tailEnergy=0.1, tailRadius=3.0)[source]

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.


kwargs: :

same as in addSmoothSquareWellForce.

tailEnergy: :

the depth of the tail part of the potential.

tailRadius: :

the maximal range of the tail part of the potential.

addLennardJonesForce(cutoff=2.5, domains=False, epsilonRep=0.24, epsilonAttr=0.27, blindFraction=-1, sigmaRep=None, sigmaAttr=None)[source]

Adds a lennard-jones force, that allows for mutual attraction. This is the slowest force out of all repulsive.


This is the only force that allows for so-called “exceptions’. Exceptions allow you to change parameters of the force for a specific pair of particles. This can be used to create short-range attraction between pairs of particles. See manual for Openmm.NonbondedForce.addException.


cutoff : float, optional

Cutoff value. Default is good.

domains : bool, optional

Use domains, defined by :py:func:’setDomains <Simulation.setDomains>’

epsilonRep : float, optional

Epsilon (attraction strength) for LJ-force for all particles (except for domain) in kT

epsilonAttr : float, optional

Epsilon for attractive domain (if domains are used) in kT

blindFraction : float, 0<x<1

Fraction of particles that are “transparent” - used here instead of truncation

sigmaRep, sigmaAttr: float, optional :

Radius of particles in the LJ force. For advanced fine-tuning.

addSoftLennardJonesForce(epsilon=0.42, trunc=2, cutoff=2.5)[source]

A softened version of lennard-Jones force. Now we’re moving to polynomial forces, so go there instead.


used to exclude a bunch of particles from calculation of nonbonded force


particles : list

List of particles for whom to exclude nonbonded force.

addInteraction(i, j, epsilon, sigma=None, length=3)[source]

Adds attractive short-range interaction of strength epsilon between particles i,j and a few neighboring particles requires :py:func:’LennardJones Force<Simulation.addLennardJonesForce>’


i,j : int

Interacting particles

epsilon : float

LJ strength

sigma : float, optional

LJ length. If you increase it past 1.5, note the cutoff!

length : int, optional, default = 3

Number of particles around i,j that also attract each other

addCylindricalConfinement(r, bottom=None, k=0.1, top=9999)[source]

As it says.

addSphericalConfinement(r='density', k=5.0, density=0.3)[source]

Constrain particles to be within a sphere. With no parameters creates sphere with density .3


r : float or “density”, optional

Radius of confining sphere. If “density” requires density, or assumes density = .3

k : float, optional

Steepness of the confining potential, in kT/nm

density : float, optional, <1

Density for autodetection of confining radius. Density is calculated in particles per nm^3, i.e. at density 1 each sphere has a 1x1x1 cube.

excludeSphere(r=5, position=(0, 0, 0))[source]

Excludes particles from a sphere of radius r at certain position.

addLaminaAttraction(width=1, depth=1, r=None)[source]

Attracts one domain to the lamina. Infers radius from spherical confinement, that has to be initialized already.


width : float, optional

Width of attractive layer next to the lamina, nm.

depth : float, optional

Depth of attractive potential in kT

r : float, optional

Radius of an attractive cage. If not specified, inferred from previously defined spherical potential.

tetherParticles(particles, k=30, positions='current')[source]

tethers particles in the ‘particles’ array. Increase k to tether them stronger, but watch the system!


particles : list of ints

List of particles to be tethered (fixed in space)

k : int, optional

Steepness of the tethering potential. Values >30 will require decreasing potential, but will make tethering rock solid.

addGravity(k=0.1, cutoff=None)[source]

adds force pulling downwards in z direction When using cutoff, acts only when z>cutoff

addPullForce(particles, forces)[source]

adds force pulling on each particle


Initializes particles velocities


mult : float, optional

Multiply velosities by this. Is good for a cold/hot start.


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 low-level parameters, such as forces, have changed.


mult : float, optional

mult to be passed to

:py:func:’initVelocities <Simulation.initVelocities>’

localEnergyMinimization(tolerance=0.3, maxIterations=0)[source]

A wrapper to the build-in OpenMM Local Energy Minimization

energyMinimization(stepsPerIteration=100, maxIterations=1000, failNotConverged=True)[source]

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.

doBlock(steps=None, increment=True, num=None, reinitialize=True)[source]

performs one block of simulations, doing steps timesteps, or steps_per_block if not specified.


steps : int or None

Number of timesteps to perform. If not specified, tries to infer it from self.steps_per_block

increment : bool, optional

If true, will not increment self.steps counter

num : int or None, optional

If specified, will split the block in num subblocks. Default value is 10.


Prints detailed statistics of a system. Will be run every 50 steps

show(shifts=[0.0, 0.2, 0.4, 0.6, 0.8])[source]

shows system in rasmol by drawing spheres draws 4 spheres in between any two points (5 * N spheres total)

A subclass used to do simulations for the Metaphase project


addAttractionToTheCore(k, r0, coreParticles=[])[source]

Attracts a subset of particles to the core, repells the rest from the core

fixParticlesZCoordinate(particles, zCoordinates, k=0.3, useOtherAxis='z', mode='abs', gap=None)[source]

Limits position of a set of particles in z coordinate


particles : list

List of particles to be fixed.

zCoordinates : list, or tuple of length 2

If has length of particles, then should contain all Z coordinates If has length 2, then contains z coordinates of first and Nth particles, and the rest is approximated linearly.

k : float, optional

Strength of attraction, measured in kT/(bondlength)

useOtherAxis : “x”,”y” or “z”, optional

Apply the same in the other dimension

class openmmlib.ExperimentalSimulation(timestep=80, thermostat=0.001, temperature=Quantity(value=300, unit=kelvin), verbose=False, velocityReinitialize=True, name='sim', length_scale=1.0, mass_scale=1.0)[source]

contains some experimental features


quickLoad(data, mode='chain', Nchains=1, trunc=None, confinementDensity='NoConfinement')[source]

quickly loads a set of repulsive chains, possibly adds spherical confinement

createWalls(left=None, right=None, k=0.5)[source]

creates walls at x = left, x = right, x direction only

addSphericalWell(r=10, depth=1)[source]

pushes particles towards a boundary of a cylindrical well to create uniform well coverage

class openmmlib.YeastSimulation(timestep=80, thermostat=0.001, temperature=Quantity(value=300, unit=kelvin), verbose=False, velocityReinitialize=True, name='sim', length_scale=1.0, mass_scale=1.0)[source]

This class is maintained by Geoff to do simulations for the Yeast project


addNucleolus(k=1, r=None)[source]