polychrom.forces module

Detailed description of forces in polychrom

This module defines forces commonly used in polychrom. Most forces are implemented using custom forces in openmm. The force equations were generally derived such that the force and the first derivative both go to zero at the cutoff radius.

Parametrization of bond forces

Most of the bond forces are parametrized using two parametrs: bondLength and bondWiggleDistance. The parameter bondLength is length of the bond at rest, while bondWiggleDistance is the estension of the bond at which energy reaches 1kT.

Note that the actual standard deviation of the bond length is bondWiggleDistance/sqrt(2) for a harmonic bond force, and is bondWiggleDistance*sqrt(2) for constant force bonds, so if you are switching from harmonic bonds to constant force, you may choose to decrease the wiggleDistance by a factor of 2.

Note on energy equations

Energy equations are passed as strings to one of the OpenMM customXXXForce class (e.g. customNonbondedForce). Note two things. First, sub-equations are separated by semicolon, and are evaluated “bottom up”, last equation first. Second, equations seem much more scary than they actually are (see below).

All energy equations have to be continuous, and we strongly believe that the first derivative has to be continuous as well. As a result, all equations were carefully crafted to be smooth functions. This makes things more complicated. For example, a simple “k * abs(x-x0)” becomes “k * (sqrt((x-x0)^2 + a^2) - a)” where a is a small number (defined to be 0.01 for example).

All energy equations have to be calculatable in single precision. Any rounding error will throw you off. For example, you should never have sqrt(A - B) where A and B are expressions, and A >= B. Because by chance, due to rounding, you may and up with A slightly less than B, and you will receive NaN, and the whole simulation will blow up. Similarly, atan(very_large_number), while defined mathematically, could easily become NaN, because very_large_number may be larger than the largest allowable float.

Note that basically all nonbonded forces were written before OpenMM introduced a switching function http://docs.openmm.org/latest/api-python/generated/simtk.openmm.openmm.CustomNonbondedForce.html Therefore, we always manually sticth the value and the first derivative of the force to be 0 at the cutoff distance. For custom user-defined forces, it may be better to use switching function instead. This does not apply to custom external forces, there stitching is still necessary.

Force equations don’t have “if” statements, but it is possible to avoid them where they would be normally used. For example, “if a: b= b0 + c” can be replaced with “b = b0 + c * delta(a)”. Similarly “f(r) if r < r0; 0 otherwise” is just “f(r) * step(r0 - r)”. These examples appear frequently in the forces that we have. One of the finest examples of crafting complex forces with on-the-fly generation of force equation is in polychrom.forces.heteropolymer_SSW(). One of the best examples of optimizing complex forces using polynomials is in polychrom.forces.polynomial_repulsive().

polychrom.forces.angle_force(sim_object, triplets, k=1.5, theta_0=<Mock name='mock.pi' id='140368592802896'>, name='angle', override_checks=False)

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) – Stiffness of the bond. If list, then determines the stiffness of the i-th triplet Potential is k * alpha^2 * 0.5 * kT

• theta_0 (float or list of length N) – Equilibrium angle of the bond. By default it is np.pi.

• override_checks (bool) – If True then do not check that no bonds are repeated. False by default.

polychrom.forces.constant_force_bonds(sim_object, bonds, bondWiggleDistance=0.05, bondLength=1.0, quadraticPart=0.02, name='abs_bonds', override_checks=False)

Constant force bond force. Energy is roughly linear with estension after r=quadraticPart; before it is quadratic to make sure the force is differentiable.

Force is parametrized using the same approach as bond force: it reaches U=kT at extension = bondWiggleDistance

Note that, just as with bondForce, mean squared extension is actually larger than wiggleDistance by sqrt(2) factor.

Parameters

• bonds (iterable of (int, int)) – Pairs of particle indices to be connected with a bond.

• bondWiggleDistance (float) – Displacement at which bond energy equals 1 kT. Can be provided per-particle.

• bondLength (float) – The length of the bond. Can be provided per-particle.

• override_checks (bool) – If True then do not check that no bonds are repeated. False by default.

polychrom.forces.cylindrical_confinement(sim_object, r, bottom=None, k=0.1, top=9999, name='cylindrical_confinement')

As it says.

polychrom.forces.grosberg_angle(sim_object, triplets, k=1.5, name='grosberg_angle', override_checks=False)

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 N-long list of floats) – Synchronized with regular stiffness. Default value is very flexible, as in Grosberg paper. Default value maximizes entanglement length.

• override_checks (bool) – If True then do not check that no bonds are repeated. False by default.

polychrom.forces.grosberg_polymer_bonds(sim_object, bonds, k=30, name='grosberg_polymer', override_checks=False)

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.

Parameters

• k (float, optional) – Arbitrary parameter; default value as in Grosberg paper.

• override_checks (bool) – If True then do not check that no bonds are repeated. False by default.

polychrom.forces.grosberg_repulsive_force(sim_object, trunc=None, radiusMult=1.0, name='grosberg_repulsive', trunc_function='min(trunc1, trunc2)')

This is the fastest non-transparent repulsive force. (that preserves topology, doesn’t allow chain passing) 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.)

Parameters

• trunc (None, float or N-array of floats) – “transparency” values for each particular particle, which correspond to the truncation values in kT for the grosberg repulsion energy between a pair of such particles. Value of 1.5 yields frequent passing, 3 - average passing, 5 - rare passing.

• radiusMult (float (optional)) – Multiplier for the size of the force. To make scale the energy larger, set to be more than 1.

• trunc_function (str (optional)) – a formula to calculate the truncation between a pair of particles with transparencies trunc1 and trunc2 Default is min(trunc1, trunc2)

polychrom.forces.harmonic_bonds(sim_object, bonds, bondWiggleDistance=0.05, bondLength=1.0, name='harmonic_bonds', override_checks=False)

Adds harmonic bonds.

Bonds are parametrized in the following way.

• A length of a bond at rest is bondLength

• Bond energy equal to 1kT at bondWiggleDistance

Note that bondWiggleDistance is not the standard deviation of the bond extension: that is actually smaller by a factor of sqrt(2).

Parameters

• bonds (iterable of (int, int)) – Pairs of particle indices to be connected with a bond.

• bondWiggleDistance (float or iterable of float) – Distance at which bond energy equals kT. Can be provided per-particle. If 0 then set k=0.

• bondLength (float or iterable of float) – The length of the bond. Can be provided per-particle.

• override_checks (bool) – If True then do not check that no bonds are repeated. False by default.

polychrom.forces.heteropolymer_SSW(sim_object, interactionMatrix, monomerTypes, extraHardParticlesIdxs, repulsionEnergy=3.0, repulsionRadius=1.0, attractionEnergy=3.0, attractionRadius=1.5, selectiveRepulsionEnergy=20.0, selectiveAttractionEnergy=1.0, keepVanishingInteractions=False, name='heteropolymer_SSW')

A version of smooth square well potential that enables the simulation of heteropolymers. Every monomer is assigned a number determining its type, then one can specify additional attraction between the types with the interactionMatrix. Repulsion between all monomers is the same, except for extraHardParticles, which, if specified, have higher repulsion energy.

The overall potential is the same as in polychrom.forces.smooth_square_well()

Treatment of extraHard particles is the same as in polychrom.forces.selective_SSW()

This is an extension of SSW (smooth square well) force in which:

1. You can give monomerTypes (e.g. 0, 1, 2 for A, B, C) and interaction strengths between these types. The corresponding entry in interactionMatrix is multiplied by selectiveAttractionEnergy to give the actual additional depth of the potential well.

2. You can select a subset of particles and make them “extra hard”. See selective_SSW force for descrition.

Force summary

Potential is the same as smooth square well, with the following parameters for particles i and j:

• Attraction energy (i,j) = attractionEnergy + selectiveAttractionEnergy * interactionMatrix[i,j]

• Repulsion Energy (i,j) = repulsionEnergy + selectiveRepulsionEnergy; if (i) or (j) are extraHard

• Repulsion Energy (i,j) = repulsionEnergy; otherwise

param interactionMatrix

the EXTRA interaction strenghts between the different types. Only upper triangular values are used. See “Force summary” abovec

type interactionMatrix

np.array

param monomerTypes

the type of each monomer, starting at 0

type monomerTypes

list of int or np.array

param extraHardParticlesIdxs

the list of indices of the “extra hard” particles. The extra hard particles repel all other particles with extra selectiveRepulsionEnergy

type extraHardParticlesIdxs

list of int

param repulsionEnergy

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

type repulsionEnergy

float

param repulsionRadius

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

type repulsionRadius

float

param attractionEnergy

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

type attractionEnergy

float

param attractionRadius

the maximal range of the attractive part of the potential.

type attractionRadius

float

param selectiveRepulsionEnergy

the EXTRA repulsion energy applied to the “extra hard” particles

type selectiveRepulsionEnergy

float

param selectiveAttractionEnergy

the EXTRA attraction energy (prefactor for the interactionMatrix interactions)

type selectiveAttractionEnergy

float

param keepVanishingInteractions

a flag that determines whether the terms that have zero interaction are still added to the force. This can be useful when changing the force dynamically (i.e. switching interactions on at some point)

type keepVanishingInteractions

bool

polychrom.forces.polynomial_repulsive(sim_object, trunc=3.0, radiusMult=1.0, name='polynomial_repulsive')

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.

See the gist below with an example of the potential. https://gist.github.com/mimakaev/0327bf6ffe7057ee0e0625092ec8e318

Parameters

trunc (float) – the energy value around r=0

polychrom.forces.pull_force(sim_object, particles, force_vecs, name='Pull')

adds force pulling on each particle particles: list of particle indices force_vecs: list of forces [[f0x,f0y,f0z],[f1x,f1y,f1z], …] if there are fewer forces than particles forces are padded with forces[-1]

polychrom.forces.selective_SSW(sim_object, stickyParticlesIdxs, extraHardParticlesIdxs, repulsionEnergy=3.0, repulsionRadius=1.0, attractionEnergy=3.0, attractionRadius=1.5, selectiveRepulsionEnergy=20.0, selectiveAttractionEnergy=1.0, name='selective_SSW')

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.

This is a tunable version of SSW: a) You can specify the set of “sticky” particles. The sticky particles are attracted only to other sticky particles. b) You can smultaneously select a subset of particles and make them “extra hard”.

This force was used two-ways. First was to make a small subset of particles very sticky. In that case, it is advantageous to make the sticky particles and their neighbours “extra hard” and thus prevent the system from collapsing.

Another useage is to induce phase separation by making all B monomers sticky. In that case, extraHard particles may not be needed at all, because the system would not collapse on iteslf.

Parameters

• stickyParticlesIdxs (list of int) – the list of indices of the “sticky” particles. The sticky particles are attracted to each other with extra selectiveAttractionEnergy

• extraHardParticlesIdxs (list of int) – the list of indices of the “extra hard” particles. The extra hard particles repel all other particles with extra selectiveRepulsionEnergy

• 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

• attractionRadius (float) – the maximal range of the attractive part of the potential.

• selectiveRepulsionEnergy (float) – the EXTRA repulsion energy applied to the extra hard particles

• selectiveAttractionEnergy (float) – the EXTRA attraction energy applied to the sticky particles

polychrom.forces.smooth_square_well(sim_object, repulsionEnergy=3.0, repulsionRadius=1.0, attractionEnergy=0.5, attractionRadius=2.0, name='smooth_square_well')

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.

Parameters

• 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

• attractionRadius (float) – the radius of the attractive part of the potential. E(attractionRadius) = 0, E’(attractionRadius) = 0

polychrom.forces.spherical_confinement(sim_object, r='density', k=5.0, density=0.3, center=[0, 0, 0], invert=False, particles=None, name='spherical_confinement')

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

Parameters

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

• center ([float, float, float]) – The coordinates of the center of the sphere.

• invert (bool) – If True, particles are not confinded, but excluded from the sphere.

• particles (list of int) – The list of particles affected by the force. If None, apply the force to all particles.

polychrom.forces.spherical_well(sim_object, particles, r, center=[0, 0, 0], width=1, depth=1, name='spherical_well')

A spherical potential well, suited for example to simulate attraction to a lamina.

Parameters

• particles (list of int or np.array) – indices of particles that are attracted

• r (float) – Radius of the nucleus

• center (vector, optional) – center position of the sphere. This parameter is useful when confining chromosomes to their territory.

• width (float, optional) – Width of attractive well, nm.

• depth (float, optional) – Depth of attractive potential in kT NOTE: switched sign from openmm-polymer, because it was confusing. Now this parameter is really the depth of the well, i.e. positive = attractive, negative = repulsive

polychrom.forces.tether_particles(sim_object, particles, *, pbc=False, k=30, positions='current', name='Tethers')

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

Parameters

• particles (list of ints) – List of particles to be tethered (fixed in space). Negative values are allowed.

• pbc (Bool, optional) – If True, periodicdistance function is applied

• k (int, optional) – The steepness of the tethering potential. Values >30 will require decreasing potential, but will make tethering rock solid. Can be provided as a vector [kx, ky, kz].