polychrom.simulation module¶
Creating a simulation: Simulation class¶
Both initialization and running the simulation is done by interacting with an instance
of polychrom.simulation.Simulation class.
Overall parameters¶
Overall technical parameters of a simulation are generally initialized in the constructor of the
Simulation class. polychrom.simulation.Simulation.__init__() . This includes
Techcnical parameters not affecting the output of simulations
Platform (cuda (usually), opencl, or CPU (slow))
GPU index
reporter (where to save results): see :py:mod`polychrom.hdf5_reporter`
Parameters affecting the simulation
number of particles
integrator (we usually use variable Langevin) + error tolerance of integrator
collision rate
Whether to use periodic boundary conditions (PBC)
timestep (if using non-variable integrator)
Parameters that are changed rarely, but may be useful
particle mass, temperature and length scale
kinetic energy at which to raise an error
OpenMM precision
Rounding before saving (default is to 0.01)
Starting conformation is loaded using polychrom.simulation.Simulation.set_data() method.
Many tools for creating starting conformations are in polychrom.starting_conformations
Adding forces¶
Forces define the main aspects of a given simulation. Polymer connectivity, confinement, crosslinks, tethering monomers, etc. are all defined as different forces acting on the particles.
Typicall used forces are listed in polychrom.forces module. Forces out of there can be added using
polychrom.simulation.Simulation.add_force() method.
Forces and their parameters are an essential part of nearly any polymer simulations. Some forces have just a few paramters (e.g. spherical confinement just needs a radius), while other forces may have lots of parameters and can define complex structures. For example, harmonidBondForce with a specially-created bond list was used to create a backbone-plectoneme conformation in Caulobacter simulations (Le et al, Science 2013). Same harmonic bonds that change over time are used to simulate loop extrusion as in (Fudenberg, 2016).
Some forces need to be added together. Those include forces defining polymer connectivity. Those forces are combined
into forcekits. Forcekits are defined in polychrom.forcekits module. The only example
of a forcekit for now is defining polymer connectivity using bonds, polymer stiffness, and inter-monomer interaction
(“nonbonded force”).
Some forces were written for openmm-polymer library and were not fully ported/tested into the polychrom library.
Those forces reside in polychrom.legacy.forces module. Some of them can be used as is, and some of them
would need to be copied to your code and potentially conformed to the new style of defining forces. This includes
accepting simulation object as a parameter, and having a .name attribute.
Defining your own forces¶
Each force in polychrom.forces is a simple function that wraps creation of an openmm force object. Users
can create new forces in the script defining their simulation and add them using add_force method. Good examples of
forces are in polychrom.forces - all but harmonic bond force use custom forces, and provide explanations of
why particular energy function was chosen. Description of the module polychrom.forces has some important
information about adding new forces.
Running a simulation¶
To run a simulation, you call polychrom.simulation.Simulation.doBlock() method in a loop.
Unless specified otherwise, this would save a conformation into a defined reporter. Terminating a
simulation is not necessary; however, terminating a reporter using reporter.dump_data() is needed for
the hdf5 reporter. This all can be viewed in the example script.
- class polychrom.simulation.Simulation(**kwargs)[source]¶
Bases:
objectThis is a base class for creating a Simulation and interacting with it. All general simulation parameters are defined in the constructor. Forces are defined in
polychrom.forcesmodule, and are added usingpolychrom.simulation.Simulation.add_force()method.- __init__(**kwargs)[source]¶
All numbers here are floats. Units specified in a parameter.
- Parameters
N (int) – number of particles
error_tol (float, optional) – Error tolerance parameter for variableLangevin integrator Values of around 0.01 are reasonable for a “nice” simulation (i.e. simulation with soft forces etc). Simulations with strong forces may need 0.001 or less OpenMM manual recommends 0.001, but our forces tend to be “softer” than theirs
timestep (number) – timestep in femtoseconds. Mandatory for non-variable integrators. Ignored for variableLangevin integrator. Value of 70-80 are appropriate
collision_rate (number) –
collision rate in inverse picoseconds. values of 0.01 or 0.05 are often used. Consult with lab members on values.
In brief, equilibrium simulations likely do not care about the exact dynamics you’re using, and therefore can be simulated in a “ballistic” dynamics with col_rate of around 0.001-0.01.
Dynamical simulations and active simulations may be more sensitive to col_rate, though this is still under discussion/investigation.
Johannes converged on using 0.1 for loop extrusion simulations, just to be safe.
PBCbox ((float,float,float) or False; default:False) – Controls periodic boundary conditions If PBCbox is False, do not use periodic boundary conditions If intending to use PBC, then set PBCbox to (x,y,z) where x,y,z are dimensions of the bounding box for PBC
GPU (GPU index as a string ("0" for first, "1" for second etc.)) – Machines with 1 GPU automatically select their GPU.
integrator ("langevin", "variableLangevin", "verlet", "variableVerlet",) – “brownian”, or tuple containing Integrator from Openmm class reference and string defining integrator type. For user-defined integrators, specify type “brownian” to avoid checking if kinetic energy exceeds max_Ek.
mass (number or np.array) – Particle mass (default 100 amu)
temperature (simtk.units.quantity(units.kelvin), optional) – Temperature of the simulation. Devault value is 300 K.
verbose (bool, optional) – If True, prints a lot of stuff in the command line.
length_scale (float, optional) – The geometric scaling factor of the system. By default, length_scale=1.0 and harmonic bonds and repulsive forces have the scale of 1 nm.
max_Ek (float, optional) – raise error if kinetic energy in (kT/particle) exceeds this value
platform (string, optional) – Platform to use: CUDA (preferred fast GPU platform) OpenCL (maybe slower GPU platofrm, does not need CUDA installed) CPU (medium speed parallelized CPU platform) reference (slow CPU platform for debug)
verbose – Shout out loud about every change.
precision (str, optional (not recommended to change)) – single is the default now, because mixed is much slower on 3080 and other new GPUs If you are using double precision, it will be slower by a factor of 10 or so.
save_decimals (int or False, optional) – Round to this number of decimals before saving.
Falseis no rounding. Default is 2. It gives maximum error of 0.005, which is nearly always harmless but saves up to 40% of storage space (0.6 of the original) Using one decimal is safe most of the time, and reduces storage to 40% of int32. NOTE that using periodic boundary conditions will make storage advantage less.reporters (list, optional) – List of reporters to use in the simulation.
- dist(i, j)[source]¶
Calculates distance between particles i and j.
Added for convenience, and not for production code. Not for use in large for-loops.
- do_block(steps=None, check_functions=[], get_velocities=False, save=True, save_extras={})[source]¶
performs one block of simulations, doing steps timesteps, or steps_per_block if not specified.
- Parameters
steps (int or None) – Number of timesteps to perform.
increment (bool, optional) – If true, will not increment self.block and self.steps counters
check_functions (list of functions, optional) – List of functions to call every block. coordinates are passed to a function. If a function returns False, simulation is stopped.
get_velocities (bool, default=False) – If True, will return velocities
save (bool, defualt=True) – If True, save results of this block
save_extras (dict) – A dict of (key, value) with additional info to save
- init_positions()[source]¶
Sends particle coordinates to OpenMM system. If system has exploded, this is
used in the code to reset coordinates.
- init_velocities(temperature='current')[source]¶
Initializes particles velocities
- Parameters
temperature (temperature to set velocities (default: temerature of the simulation)) –
- initialize(**kwargs)[source]¶
Initialize, particles, velocities for the first time. Only need to use this function if your system has no forces (free Brownian particles). Otherwise _apply_force() will execute these lines to add particles to the system, initialize their positions/velocities, initialize the context.
- local_energy_minimization(tolerance=0.3, maxIterations=0, random_offset=0.02)[source]¶
A wrapper to the build-in OpenMM Local Energy Minimization
See caveat below
- Parameters
tolerance (float) –
It is something like a value of force below which the minimizer is trying to minimize energy to. see openmm documentation for description
Value of 0.3 seems to be fine for most normal forces.
maxIterations (int) –
Maximum # of iterations for minimization to do. default: 0 means there is no limit
This is relevant especially if your simulation does not have a well-defined energy minimum (e.g. you want to simulate a collapse of a chain in some potential). In that case, if you don’t limit energy minimization, it will attempt to do a whole simulation for you. In that case, setting a limit to the # of iterations will just stop energy minimization manually when it reaches this # of iterations.
random_offset (float) –
A random offset to introduce after energy minimization. Should ideally make your forces have realistic values.
For example, if your stiffest force is polymer bond force with “wiggle_dist” of 0.05, setting this to 0.02 will make separation between monomers realistic, and therefore will make force values realistic.
See why do we need it in the caveat below.
Notes
If using variable langevin integrator after minimization, a big error may happen in the first timestep. The reason is that enregy minimization makes all the forces basically 0. Variable langevin integrator measures the forces and assumes that they are all small - so it makes the timestep very large, and at the first timestep it overshoots completely and energy goes up a lot.
The workaround for now is to randomize positions after energy minimization
- reinitialize()[source]¶
Reinitializes the OpenMM context object. This should be called if low-level parameters, such as parameters of forces, have changed
- set_data(data, center=False, random_offset=1e-05, report=True)[source]¶
Sets particle positions
- Parameters
data (Nx3 array-like) – Array of positions
center (bool or "zero", optional) – Move center of mass to zero before starting the simulation if center == “zero”, then center the data such as all positions are positive and start at zero
random_offset (float or None) – add random offset to each particle Recommended for integer starting conformations and in general
report (bool, optional) – If set to False, will not report this action to reporters.