diff --git a/polychrom/extra_openmm_forces.py b/polychrom/extra_forces.py similarity index 100% rename from polychrom/extra_openmm_forces.py rename to polychrom/extra_forces.py diff --git a/polychrom/forcekits.py b/polychrom/forcekits.py new file mode 100644 index 0000000..7349a15 --- /dev/null +++ b/polychrom/forcekits.py @@ -0,0 +1,70 @@ +from . import openmm_forces + +def polymerChains( + sim_object, + chains=[(0, None, False)] + + bondForceFunc=openmm_forces.harmonicBonds, + bondForceKwargs={'wiggleDist'=0.05, + 'bondLength'=1.0}, + + angleForceFunc=openmm_forces.harmonicBonds, + angleForceKwargs={'k'=0.05}, + + nonbondedForceFunc=openmm_forces.polynomialRepulsiveForce, + nonbondedForceKwargs={'trunc'=3.0, + 'radiusMult'=1.}, + + exceptBonds=True, +): + """Adds harmonic bonds connecting polymer chains + + Parameters + ---------- + 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. + + exceptBonds : bool + If True then do not calculate non-bonded forces between the + particles connected by a bond. True by default. + """ + + bonds = [] + for start, end, is_ring in chains: + end = sim_object.N if (end is None) else end + + bonds += [(j, j+1) for j in range(start, end - 1)] + triples += [(j - 1, j, j + 1) for j in range(start + 1, end - 1)] + if is_ring: + bonds.append((start, end-1)) + triplets.append((int(end - 2), int(end - 1), int(start))) + triplets.append((int(end - 1), int(start), int(start + 1))) + + bondForceFunc(sim_object, bonds, **bondForceKwargs) + + if angleForceFunc is not None: + angleForceFunc(sim_object, triplets, **angleForceKwargs) + + if nonbondedForceFunc is not None: + nb_force = nonbondedForceFunc(sim_object, **nonbondedForceKwargs) + + if exceptBonds: + exc = list(set([tuple(i) for i in np.sort(np.array(bonds), axis=1)])) + if hasattr(nb_force, "addException"): + print('Add exceptions for {0} force'.format(i)) + for pair in exc: + nb_force.addException(int(pair[0]), int(pair[1]), 0, 0, 0, True) + + # The built-in LJ nonbonded force uses "exclusions" instead of "exceptions" + elif hasattr(nb_force, "addExclusion"): + print('Add exclusions for {0} force'.format(i)) + for pair in exc: + nb_force.addExclusion(int(pair[0]), int(pair[1])) + + print("Number of exceptions:", len(bonds)) + + diff --git a/polychrom/openmm_forces.py b/polychrom/forces.py similarity index 69% rename from polychrom/openmm_forces.py rename to polychrom/forces.py index 0831bc3..d26dbdc 100644 --- a/polychrom/openmm_forces.py +++ b/polychrom/forces.py @@ -1,149 +1,121 @@ import simtk.openmm as openmm import simtk.unit as units +import itertools + nm = units.meter * 1e-9 fs = units.second * 1e-15 ps = units.second * 1e-12 import numpy as np -def _initHarmonicBondForce(sim_object): - "Internal, inits harmonic forse for polymer and non-polymer bonds" - if "HarmonicBondForce" not in list(sim_object.forceDict.keys()): - sim_object.forceDict["HarmonicBondForce"] = openmm.HarmonicBondForce() - sim_object.bondType = "Harmonic" - - -def _initAbsBondForce(sim_object): - "inits abs(x) FENE bond force" - if "AbsBondForce" not in list(sim_object.forceDict.keys()): - force = "(1. / ABSwiggle) * ABSunivK * "\ - "(sqrt((r-ABSr0 * ABSconlen)* "\ - " (r - ABSr0 * ABSconlen) + ABSa * ABSa) - ABSa)" - - bondforceAbs = openmm.CustomBondForce(force) - bondforceAbs.addPerBondParameter("ABSwiggle") - bondforceAbs.addPerBondParameter("ABSr0") - bondforceAbs.addGlobalParameter("ABSunivK", sim_object.kT / sim_object.conlen) - bondforceAbs.addGlobalParameter("ABSa", 0.02 * sim_object.conlen) - bondforceAbs.addGlobalParameter("ABSconlen", sim_object.conlen) - sim_object.forceDict["AbsBondForce"] = bondforceAbs - - - -def addBond(sim_object, - i, j, # particles connected by bond - bondWiggleDistance=0.2, - # Flexibility of the bond, - # measured in distance at which energy equals kT - distance=None, # Equilibrium length of the bond, default = sim_object.length_scale - bondType=None, # Harmonic, Grosberg, ABS - verbose=None): # Set this to False or True to override sim_object.verbose for this function - # and don't want to contaminate output by 10000 messages - """Adds bond between two particles, allows to specify parameters - - Parameters - ---------- - - i,j : int - Particle numbers - bondWiggleDistance : float - Average displacement from the equilibrium bond distance +def _to_array_1d(scalar_or_array, arrlen, dtype=float): + if not hasattr(bondLength, "__iter__"): + outarr = np.full(arrlen, scalar_or_array, dtype) + else: + outarr = np.asarray(scalar_or_array, dtype=dtype) + + if len(outarr) != arrlen + raise ValueError('The length of the array differs from the expected one!') + + return outarr - bondType : "Harmonic" or "abs" - 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 +def harmonicBonds(sim_object, + bonds, + bondWiggleDistance=0.05, + bondLength=1.0): + """Adds harmonic bonds + Parameters + ---------- + + bonds : iterable of (int, int) + Pairs of particle indices to be connected with a bond. + bondWiggleDistance : float or iterable of float + Average displacement from the equilibrium bond distance. + Can be provided per-particle. + bondLength : float or iterable of float + The length of the bond. + Can be provided per-particle. """ + + bondForce = openmm.HarmonicBondForce() - if not hasattr(sim_object, "bondLengths"): - sim_object.bondLengths = [] - - if verbose is None: - verbose = sim_object.verbose - if (i >= sim_object.N) or (j >= sim_object.N): - raise ValueError("\nCannot add bond with monomers %d,%d that"\ - "are beyound the polymer length %d" % (i, j, sim_object.N)) - - bondSize = float(bondWiggleDistance) - - if distance is None: - distance = sim_object.length_scale - else: - distance = sim_object.length_scale * distance - distance = float(distance) - - if not hasattr(sim_object, "kbondScalingFactor"): # caching kbondScalingFactor - performance improvement... - t = (2 * sim_object.kT / (sim_object.conlen) ** 2) / (units.kilojoule_per_mole / nm ** 2) - print(t) - sim_object.kbondScalingFactor = float((2 * sim_object.kT / (sim_object.conlen) ** 2) / (units.kilojoule_per_mole / nm ** 2)) - - kbondScalingFactor = sim_object.kbondScalingFactor #... not to calculate it eevry time we add bond - # this will be an integer, so we don't have to deal with slow simtk.units every time we add a bond - - if bondType is None: - bondType = sim_object.bondType - - if bondType.lower() == "harmonic": - _initHarmonicBondForce(sim_object) - kbond = kbondScalingFactor / (bondSize ** 2) # using kbondScalingFactor because force accepts parameters with units - sim_object.forceDict["HarmonicBondForce"].addBond(int(i), int(j), float(distance), float(kbond)) - sim_object.bondLengths.append([int(i), int(j), float(distance), float(bondSize)]) - elif bondType.lower() == "abs": - _initAbsBondForce(sim_object) - sim_object.forceDict["AbsBondForce"].addBond(int(i), int( - j), [float(bondWiggleDistance), float(distance)]) # force is initialized to accept floats already - sim_object.bondLengths.append([int(i), int(j), float(distance), float(bondSize)]) - else: - sim_object._exitProgram("Bond type not known") - if verbose == True: - print("%s bond added between %d,%d, wiggle %lf dist %lf" % ( - bondType, i, j, float(bondWiggleDistance), float(distance))) - - - -def harmonicPolymerBonds(sim_object, - wiggleDist=0.05, - bondLength=1.0, - exceptBonds=True): - """Adds harmonic bonds connecting polymer chains + bondLength = _to_array_1d(bondLength, len(bonds)) * sim_object.length_scale + bondWiggleDistance = _to_array_1d(bondWiggleDistance, len(bonds)) * sim_object.length_scale + + # using kbondScalingFactor because force accepts parameters with units + kbond = sim_object.kbondScalingFactor / (bondWiggleDistance ** 2) + + for bond_idx, (i, j) in enumerate(bonds): + if (i >= sim_object.N) or (j >= sim_object.N): + raise ValueError("\nCannot add bond with monomers %d,%d that"\ + "are beyound the polymer length %d" % (i, j, sim_object.N)) + + bondForce.addBond(int(i), + int(j), + float(bondLength[bond_idx]), + float(kbond[bond_idx])) + + sim_object.forceDict["HarmonicBondForce"] = bondForce + + return bondForce + + +def FENEBonds(sim_object, + bonds, + bondWiggleDistance=0.05, + bondLength=1.0): + """Adds harmonic bonds Parameters ---------- - - wiggleDist : float - Average displacement from the equilibrium bond distance + + bonds : iterable of (int, int) + Pairs of particle indices to be connected with a bond. + bondWiggleDistance : float + Average displacement from the equilibrium bond distance. + Can be provided per-particle. bondLength : float - The length of the bond - exceptBonds : bool - If True then do not calculate non-bonded forces between the - particles connected by a bond. True by default. + The length of the bond. + Can be provided per-particle. """ + + forceExpr = "(1. / ABSwiggle) * ABSunivK * "\ + "(sqrt((r-ABSr0 * ABSconlen)* "\ + " (r - ABSr0 * ABSconlen) + ABSa * ABSa) - ABSa)" + bondforceAbs = openmm.CustomBondForce(force) + + bondforceAbs.addPerBondParameter("ABSwiggle") + bondforceAbs.addPerBondParameter("ABSr0") + bondforceAbs.addGlobalParameter("ABSunivK", sim_object.kT / sim_object.conlen) + bondforceAbs.addGlobalParameter("ABSa", 0.02 * sim_object.conlen) + bondforceAbs.addGlobalParameter("ABSconlen", sim_object.conlen) + + bondLength = _to_array_1d(bondLength, len(bonds)) * sim_object.length_scale + bondWiggleDistance = _to_array_1d(bondWiggleDistance, len(bonds)) * sim_object.length_scale + + for bond_idx, (i, j) in enumerate(bonds): + if (i >= sim_object.N) or (j >= sim_object.N): + raise ValueError("\nCannot add bond with monomers %d,%d that"\ + "are beyound the polymer length %d" % (i, j, sim_object.N)) + + bondforceAbs.addBond(int(i), + int(j), + [float(bondWiggleDistance[bond_idx]), + float(bondLength[bond_idx])]) + + sim_object.forceDict["AbsBondForce"] = bondforceAbs + + return bondforceAbs - for start, end, isRing in sim_object.chains: - for j in range(start, end - 1): - addBond(sim_object, j, j + 1, wiggleDist, - distance=bondLength, - bondType="Harmonic", verbose=False) - if exceptBonds: - sim_object.bondsForException.append((j, j + 1)) - - if isRing: - addBond(sim_object, start, end - 1, wiggleDist, - distance=bondLength, bondType="Harmonic") - if exceptBonds: - sim_object.bondsForException.append((start, end - 1)) - if sim_object.verbose == True: - print("ring bond added", start, end - 1) - - -def angleForce(sim_object, k=1.5): +def angleForce( + sim_object, + triplets, + k=1.5): """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 @@ -155,29 +127,30 @@ def angleForce(sim_object, k=1.5): k : float or list of length N Stiffness of the bond. - If list, then determines stiffness of the bond at monomer i. + If list, then determines the stiffness of the i-th triplet Potential is k * alpha^2 * 0.5 * kT """ - try: - k[0] - except: - k = np.zeros(sim_object.N, float) + k + + k = _to_array_1d(k, len(triples)) + stiffForce = openmm.CustomAngleForce( "kT*angK * (theta - 3.141592) * (theta - 3.141592) * (0.5)") - sim_object.forceDict["AngleForce"] = stiffForce - for start, end, isRing in sim_object.chains: - for j in range(start + 1, end - 1): - stiffForce.addAngle(j - 1, j, j + 1, [float(k[j])]) - if isRing: - stiffForce.addAngle(int(end - 2),int( end - 1), int(start), [float(k[end - 1]) ]) - stiffForce.addAngle(int(end - 1),int( start), int(start + 1), [float(k[start]) ]) - + stiffForce.addGlobalParameter("kT", sim_object.kT) stiffForce.addPerAngleParameter("angK") + + for triplet_idx, (p1, p2, p3) in enumerate(triplets): + stiffForce.addAngle(p1, p2, p3, [k[triplet_idx]]) + + sim_object.forceDict["AngleForce"] = stiffForce + + return stiffForce - -def polynomialRepulsiveForce(sim_object, trunc=3.0, radiusMult=1.): +def polynomialRepulsiveForce( + sim_object, + trunc=3.0, + radiusMult=1.): """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. @@ -197,9 +170,7 @@ def polynomialRepulsiveForce(sim_object, trunc=3.0, radiusMult=1.): "rsc4 = rsc2 * rsc2;" "rsc2 = rsc * rsc;" "rsc = r / REPsigma * REPrmin;") - sim_object.forceDict["Nonbonded"] = openmm.CustomNonbondedForce( - repul_energy) - repulforceGr = sim_object.forceDict["Nonbonded"] + repulforceGr = openmm.CustomNonbondedForce(repul_energy) repulforceGr.addGlobalParameter('REPe', trunc * sim_object.kT) repulforceGr.addGlobalParameter('REPsigma', radius) @@ -213,6 +184,11 @@ def polynomialRepulsiveForce(sim_object, trunc=3.0, radiusMult=1.): repulforceGr.addParticle(()) repulforceGr.setCutoffDistance(nbCutOffDist) + + sim_object.forceDict["Nonbonded"] = repulforceGr + + return repulforceGr + def smoothSquareWellForce(sim_object, repulsionEnergy=3.0, repulsionRadius=1., @@ -263,9 +239,8 @@ def smoothSquareWellForce(sim_object, "rshft = (r - REPsigma - ATTRdelta) / ATTRdelta * rmin12" ) - sim_object.forceDict["Nonbonded"] = openmm.CustomNonbondedForce( - energy) - repulforceGr = sim_object.forceDict["Nonbonded"] + + repulforceGr = openmm.CustomNonbondedForce( energy) repulforceGr.addGlobalParameter('REPe', repulsionEnergy * sim_object.kT) repulforceGr.addGlobalParameter('REPsigma', repulsionRadius * sim_object.conlen) @@ -281,7 +256,12 @@ def smoothSquareWellForce(sim_object, repulforceGr.addParticle(()) repulforceGr.setCutoffDistance(nbCutOffDist) + + sim_object.forceDict["Nonbonded"] = repulforceGr + + return repulforceGr + def selectiveSSWForce(sim_object, stickyParticlesIdxs, extraHardParticlesIdxs, @@ -395,7 +375,8 @@ def selectiveSSWForce(sim_object, float(i in extraHardParticlesIdxs))) sim_object.forceDict["Nonbonded"] = repulforceGr - + + return repulforceGr def cylindricalConfinement(sim_object, r, bottom=None, k=0.1, top=9999): @@ -417,7 +398,6 @@ def cylindricalConfinement(sim_object, r, bottom=None, k=0.1, top=9999): "step(r-CYLaa) * CYLkb * (sqrt((r-CYLaa)*(r-CYLaa) + CYLt*CYLt) - CYLt);" "r = sqrt(x^2 + y^2 + CYLtt^2)") - sim_object.forceDict["CylindricalConfinement"] = extforce2 for i in range(sim_object.N): extforce2.addParticle(i, []) extforce2.addGlobalParameter("CYLkb", k * sim_object.kT / nm) @@ -429,7 +409,12 @@ def cylindricalConfinement(sim_object, r, bottom=None, k=0.1, top=9999): extforce2.addGlobalParameter("CYLaa", (r - 1. / k) * nm) extforce2.addGlobalParameter("CYLt", (1. / (10 * k)) * nm) extforce2.addGlobalParameter("CYLtt", 0.01 * nm) + + sim_object.forceDict["CylindricalConfinement"] = extforce2 + + return extforce2 + def sphericalConfinement(sim_object, r="density", # radius... by default uses certain density k=5., # How steep the walls are @@ -454,14 +439,12 @@ def sphericalConfinement(sim_object, spherForce = openmm.CustomExternalForce( "step(r-SPHaa) * SPHkb * (sqrt((r-SPHaa)*(r-SPHaa) + SPHt*SPHt) - SPHt) " ";r = sqrt(x^2 + y^2 + z^2 + SPHtt^2)") - sim_object.forceDict["SphericalConfinement"] = spherForce for i in range(sim_object.N): spherForce.addParticle(i, []) if r == "density": r = (3 * sim_object.N / (4 * 3.141592 * density)) ** (1 / 3.) - sim_object.sphericalConfinementRadius = r if sim_object.verbose == True: print("Spherical confinement with radius = %lf" % r) # assigning parameters of the force @@ -469,9 +452,12 @@ def sphericalConfinement(sim_object, spherForce.addGlobalParameter("SPHaa", (r - 1. / k) * nm) spherForce.addGlobalParameter("SPHt", (1. / k) * nm / 10.) spherForce.addGlobalParameter("SPHtt", 0.01 * nm) - return r - - + + ## TODO: move 'r' elsewhere?.. + sim_object.sphericalConfinementRadius = r + sim_object.forceDict["SphericalConfinement"] = spherForce + + return spherForce def tetherParticles(sim_object, particles, k=30, positions="current"): @@ -488,12 +474,9 @@ def tetherParticles(sim_object, particles, k=30, positions="current"): Values >30 will require decreasing potential, but will make tethering rock solid. """ - if "Tethering Force" not in sim_object.forceDict: - tetherForce = openmm.CustomExternalForce( - " TETHkb * ((x - TETHx0)^2 + (y - TETHy0)^2 + (z - TETHz0)^2)") - sim_object.forceDict["Tethering Force"] = tetherForce - else: - tetherForce = sim_object.forceDict["Tethering Force"] + + tetherForce = openmm.CustomExternalForce( + "TETHkb * ((x - TETHx0)^2 + (y - TETHy0)^2 + (z - TETHz0)^2)") # assigning parameters of the force tetherForce.addGlobalParameter("TETHkb", k * sim_object.kT / nm) @@ -511,6 +494,11 @@ def tetherParticles(sim_object, particles, k=30, positions="current"): if sim_object.verbose == True: print("particle %d tethered! " % i) + sim_object.forceDict["Tethering Force"] = tetherForce + + return tetherForce + + def pullForce(sim_object, particles, forces): """ adds force pulling on each particle @@ -518,7 +506,6 @@ def pullForce(sim_object, particles, forces): forces: list of forces [[f0x,f0y,f0z],[f1x,f1y,f1z], ...] if there are fewer forces than particles forces are padded with forces[-1] """ - import itertools pullForce = openmm.CustomExternalForce( "PULLx * x + PULLy * y + PULLz * z") pullForce.addPerParticleParameter("PULLx") @@ -528,3 +515,5 @@ def pullForce(sim_object, particles, forces): force = [float(i) * (sim_object.kT / sim_object.conlen) for i in force] pullForce.addParticle(num, force) sim_object.forceDict["PullForce"] = pullForce + + return pullForce diff --git a/polychrom/openmm_simulator.py b/polychrom/simulaton.py similarity index 91% rename from polychrom/openmm_simulator.py rename to polychrom/simulaton.py index 63127bb..81818bf 100755 --- a/polychrom/openmm_simulator.py +++ b/polychrom/simulaton.py @@ -17,12 +17,15 @@ fs = units.second * 1e-15 ps = units.second * 1e-12 + class integrationFailError(Exception): pass + class eKExceedsError(Exception): pass + class Simulation(): """Base class for openmm simulations @@ -184,10 +187,12 @@ def __init__(self, **kwargs): self.kB = units.BOLTZMANN_CONSTANT_kB * \ units.AVOGADRO_CONSTANT_NA # Boltzmann constant self.kT = self.kB * self.temperature # thermal energy + # All masses are the same, # unless individual mass multipliers are specified in self.load() - self.bondsForException = [] self.conlen = 1. * nm * self.length_scale + + self.kbondScalingFactor = float((2 * self.kT / (self.conlen) ** 2) / (units.kilojoule_per_mole / nm ** 2)) self.system = openmm.System() self.PBC = kwargs["PBC"] @@ -201,10 +206,6 @@ def __init__(self, **kwargs): self.forceDict = {} # Dictionary to store forces - - - - def saveFolder(self, folder): """ sets the folder where to save data. @@ -219,46 +220,13 @@ def saveFolder(self, folder): os.mkdir(folder) self.folder = folder + def _exitProgram(self, line): print(line) print("--------------> Bye <---------------") exit() - def setChains(self, chains=[(0, None, 0)]): - """ - 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(). - - Parameters - ---------- - - 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. - """ - - if not hasattr(self, "N"): - raise ValueError("Load the chain first, or provide chain length") - - self.chains = [i for i in chains] # copy - for i in range(len(self.chains)): - start, end, isRing = self.chains[i] - end = self.N if (end is None) else end - self.chains[i] = (start, end, isRing) - - def getChains(self): - "returns configuration of chains" - return self.chains - - - - def save(self, filename=None, mode="joblib"): """Saves conformation plus some metadata. Metadata is not interpreted by this library, and is for your reference @@ -291,6 +259,7 @@ def getData(self): "Returns an Nx3 array of positions" return np.asarray(self.data / nm, dtype=np.float32) + def getScaledData(self): """Returns data, scaled back to PBC box """ if self.PBC != True: @@ -302,6 +271,7 @@ def getScaledData(self): assert toRet.min() >= 0 return toRet + def setData(self, data, center=False, random_offset = 1e-5): """Sets particle positions @@ -344,10 +314,6 @@ def setData(self, data, center=False, random_offset = 1e-5): self.data = units.Quantity(data, nm) - - if not hasattr(self, "chains"): - self.setChains() - if hasattr(self, "context"): self.initPositions() @@ -363,6 +329,7 @@ def RG(self): data = data - np.mean(data, axis=0)[None,:] return np.sqrt(np.sum(np.var(np.array(data), 0))) + def dist(self, i, j): """ Calculates distance between particles i and j @@ -371,6 +338,7 @@ def dist(self, i, j): dif = data[i] - data[j] return np.sqrt(sum(dif ** 2)) + def _applyForces(self): """Adds all particles to the system. Then applies all the forces in the forcedict. @@ -384,30 +352,17 @@ def _applyForces(self): for mass in self.masses: self.system.addParticle(mass) - print("Number of exceptions:", len(self.bondsForException)) - - if len(self.bondsForException) > 0: - exc = list(set([tuple(i) for i in np.sort(np.array(self.bondsForException), axis=1)])) - for i in list(self.forceDict.keys()): # Adding exceptions + for i in list(self.forceDict.keys()): # Adding forces force = self.forceDict[i] - if hasattr(force, "addException"): - print('Add exceptions for {0} force'.format(i)) - for pair in exc: - force.addException(int(pair[0]), - int(pair[1]), 0, 0, 0, True) - elif hasattr(force, "addExclusion"): - print('Add exclusions for {0} force'.format(i)) - for pair in exc: - # force.addExclusion(*pair) - force.addExclusion(int(pair[0]), int(pair[1])) if hasattr(force, "CutoffNonPeriodic") and hasattr(force, "CutoffPeriodic"): if self.PBC: force.setNonbondedMethod(force.CutoffPeriodic) - print("Using periodic boundary conditions!!!!") + print("Using periodic boundary conditions!!!") else: force.setNonbondedMethod(force.CutoffNonPeriodic) + print("adding force ", i, self.system.addForce(self.forceDict[i])) self.context = openmm.Context(self.system, self.integrator, self.platform, self.properties) @@ -415,6 +370,7 @@ def _applyForces(self): self.initVelocities() self.forcesApplied = True + def initVelocities(self, temperature="current"): """Initializes particles velocities @@ -448,6 +404,7 @@ def initPositions(self): eP = self.context.getState(getEnergy=True).getPotentialEnergy() / self.N / self.kT print("Particles loaded. Potential energy is %lf" % eP) + def reinitialize(self): """Reinitializes the OpenMM context object. This should be called if low-level parameters, @@ -595,7 +552,7 @@ def printStats(self): minmedmax = lambda x: (x.min(), np.median(x), x.mean(), x.max()) - print("\n Statistics: number of particles: %d, number of chains: %d\n" % (self.N, len(self.chains))) + print("\n Statistics: number of particles: %d\n" % (self.N, )) print("Statistics for particle position") print(" mean position is: ", np.mean( pos, axis=0), " Rg = ", self.RG()) @@ -622,10 +579,10 @@ def printStats(self): print() print("Statistics for the system:") print(" Forces are: ", list(self.forceDict.keys())) - print(" Number of exceptions: ", len(self.bondsForException)) print() print("Potential Energy Ep = ", eP / self.N / self.kT) + def show(self, shifts=[0., 0.2, 0.4, 0.6, 0.8], scale="auto"): """shows system in rasmol by drawing spheres draws 4 spheres in between any two points (5 * N spheres total) diff --git a/tests/bondLength.py b/tests/bondLength.py index ba6c96d..8ee5219 100755 --- a/tests/bondLength.py +++ b/tests/bondLength.py @@ -2,7 +2,7 @@ # Code written by: Maksim Imakaev (imakaev@mit.edu) import matplotlib.pyplot as plt import numpy as np -from openmmlib import Simulation +from polychrom import simulation, forces, forcekits, starting_conformation import polymerutils import os @@ -12,24 +12,30 @@ def exampleOpenmm(): Please follow comments along the text for explanations. """ - # ----------- Initializing general simulation parameters--------- - - # Initialization for simulations with constant environment, using default integrator (langevin) - # Fine-tune timestep and thermostat parameters so that your simulation does not blow up, - # But is going as fast as possible. You might need to increase timestep, but don't let - # your kinetic energy be above 1.6 in the "steady" regime - - # a = Simulation(timestep=80, thermostat=0.005) - # a.setup(platform="cuda", verbose=True) # Switch to platform="OpenCL" if you don't have cuda - - # Alternative initialization for dynamic simulations with strong forces, - # which would automatically adjusts timestep - # This is relevant, for example, for simulations of polymer collapse - # If simulation blows up, decrease errorTol by a factor of two and try again - a = Simulation(thermostat=0.02) # timestep not necessary for variableLangevin - a.setup(platform="cuda", integrator="variableLangevin", errorTol=0.06, verbose=True) + sim = simulation.Simulation( + platform="cuda", # Switch to platform="OpenCL" if you don't have cuda + + # + # Fine-tune timestep and thermostat parameters so that your simulation does not blow up, + # But is going as fast as possible. You might need to increase timestep, but don't let + # your kinetic energy be above 1.6 in the "steady" regime + # integrator="langevin", + # timestep=80, + # thermostat=0.005 + + # Alternative initialization for dynamic simulations with strong forces, + # which would automatically adjusts timestep + # This is relevant, for example, for simulations of polymer collapse + # If simulation blows up, decrease errorTol by a factor of two and try again + # a variableLangevin thermostat automatically calculates the timestep + # for a given error level + integrator="variableLangevin", + errorTol=0.06, + thermostat=0.02, + + verbose=True) - a.saveFolder("trajectory") # folder where to save trajectory + sim.saveFolder("trajectory") # folder where to save trajectory # ------- Creation of the initial conformation----------- @@ -37,66 +43,68 @@ def exampleOpenmm(): # polymer = polymerutils.load("globule") # loads compact polymer conformation of the length 6000 - # polymer = polymerutils.grow_rw(8000, 50, method="standard") + # polymer = starting_conformations.grow_cubic(8000, 50, method="standard") # grows a compact polymer ring of a length 8000 in a 50x50x50 box - # polymer = polymerutils.create_spiral(r1=4, r2=10, N=8000) + # polymer = starting_conformations.create_spiral(r1=4, r2=10, N=8000) # Creates a compact polymer arranged in a cylinder of radius 10, 8000 monomers long - polymer = polymerutils.create_random_walk(1, 8000) + polymer = starting_conformations.create_random_walk(1, 8000) # Creates an extended "random walk" conformation of length 8000 - a.load(polymer, center=True) # loads a polymer, puts a center of mass at zero + sim.setData(polymer, center=True) # loads a polymer, puts a center of mass at zero - # -----------Initialize conformation of the chains-------- - # By default the library assumes you have one polymer chain - # If you want to make it a ring, or more than one chain, use self.setChains - # self.setChains([(0,50,1),(50,None,0)]) will set a 50-monomer ring and a chain from monomer 50 to the end # -----------Adding forces --------------- -# a.addSphericalConfinement(density=0.85, k=1) + forcekits.polymerChains( + sim, + # By default the library assumes you have one polymer chain + # If you want to make it a ring, or more than one chain, provide a chains parameter + # chains = [(0, 50, True), (50, None, False)], # set a 50-monomer ring and a chain from monomer 50 to the end + + #bondForceFunc=forces.FENEBonds, # uncomment to bind particles in the chains with + # a constant force (i.e. use a linear potential instead of harmonic) + bondForceKwags={'wiggleDist':0.05}, # Bond distance will fluctuate +- 0.05 on average + + #angleForceFunc=None, # uncomment to disable stiffness + angleForceKwargs={'k':4}, + # K is more or less arbitrary, k=4 corresponds to presistence length of 4, + # k=1.5 is recommended to make polymer realistically flexible; k=8 is very stiff + + #nonbondedForceFunc=None, # uncomment to disable particle repulsion + nonbondedForceKwargs={'trunc':1.5}, + ) + +# forces.sphericalConfinement(sim, density=0.85, k=1) # Specifying density is more intuitive than radius # k is the slope of confinement potential, measured in kT/mon # set k=5 for harsh confinement # and k = 0.2 or less for collapse simulation - a.addHarmonicPolymerBonds(wiggleDist=0.05) - # Bond distance will fluctuate +- 0.05 on average - - #a.addGrosbergRepulsiveForce(trunc=50) - # this will resolve chain crossings and will not let chain cross anymore - - # a.addGrosbergRepulsiveForce(trunc=5) - # this will let chains cross sometimes - - #a.addStiffness(k=4) - # K is more or less arbitrary, k=4 corresponds to presistence length of 4, - # k=1.5 is recommended to make polymer realistically flexible; k=8 is very stiff # If your simulation does not start, consider using energy minimization below - # a.localEnergyMinimization() + # sim.localEnergyMinimization() # A very efficient algorithm to reach local energy minimum # Use it to minimize energy if you're doing diffusion simulations # If you're simulating dynamics of collapse or expansion, please do not use it - #a.energyMinimization(stepsPerIteration=10) # increase to 100 for larger or more complex systems + #sim.energyMinimization(stepsPerIteration=10) # increase to 100 for larger or more complex systems # An algorithm to start a simulation # Which works only with langevin integrator (but will not throw an error otherwise) # Decreases a timestep, and then increases it slowly back to normal # -----------Running a simulation --------- - a.save() # save original conformation + sim.save() # save original conformation allDists = [] for _ in xrange(10): # Do 10 blocks - a.doBlock(2000) # Of 2000 timesteps each - data = a.getData() + sim.doBlock(2000) # Of 2000 timesteps each + data = sim.getData() dists = np.sqrt(np.sum(np.diff(data, axis=0)**2, axis=1)) allDists.append(dists) - allDists = np.concatenate(allDists) plt.hist(allDists,100) plt.title("deviation should be around 0.05") diff --git a/tests/bondLength.py~ b/tests/bondLength.py~ deleted file mode 100755 index 00b8735..0000000 --- a/tests/bondLength.py~ +++ /dev/null @@ -1,108 +0,0 @@ -# (c) 2013 Massachusetts Institute of Technology. All Rights Reserved -# Code written by: Maksim Imakaev (imakaev@mit.edu) -import matplotlib.pyplot as plt -import numpy as np -from openmmlib import Simulation -import polymerutils -import os - -def exampleOpenmm(): - """ - An example script which generates an extended polymer, and lets it collapse to a sphere. - Please follow comments along the text for explanations. - """ - - # ----------- Initializing general simulation parameters--------- - - # Initialization for simulations with constant environment, using default integrator (langevin) - # Fine-tune timestep and thermostat parameters so that your simulation does not blow up, - # But is going as fast as possible. You might need to increase timestep, but don't let - # your kinetic energy be above 1.6 in the "steady" regime - - # a = Simulation(timestep=80, thermostat=0.005) - # a.setup(platform="cuda", verbose=True) # Switch to platform="OpenCL" if you don't have cuda - - # Alternative initialization for dynamic simulations with strong forces, - # which would automatically adjusts timestep - # This is relevant, for example, for simulations of polymer collapse - # If simulation blows up, decrease errorTol by a factor of two and try again - a = Simulation(thermostat=0.02) # timestep not necessary for variableLangevin - a.setup(platform="cuda", integrator="variableLangevin", errorTol=0.06, verbose=True) - - a.saveFolder("trajectory") # folder where to save trajectory - - - # ------- Creation of the initial conformation----------- - - # polymer = polymerutils.load("globule") - # loads compact polymer conformation of the length 6000 - - # polymer = polymerutils.grow_rw(8000, 50, method="standard") - # grows a compact polymer ring of a length 8000 in a 50x50x50 box - - # polymer = polymerutils.create_spiral(r1=4, r2=10, N=8000) - # Creates a compact polymer arranged in a cylinder of radius 10, 8000 monomers long - - polymer = polymerutils.create_random_walk(1, 8000) - # Creates an extended "random walk" conformation of length 8000 - - a.load(polymer, center=True) # loads a polymer, puts a center of mass at zero - - # -----------Initialize conformation of the chains-------- - # By default the library assumes you have one polymer chain - # If you want to make it a ring, or more than one chain, use self.setChains - # self.setChains([(0,50,1),(50,None,0)]) will set a 50-monomer ring and a chain from monomer 50 to the end - - - # -----------Adding forces --------------- -# a.addSphericalConfinement(density=0.85, k=1) - # Specifying density is more intuitive than radius - # k is the slope of confinement potential, measured in kT/mon - # set k=5 for harsh confinement - # and k = 0.2 or less for collapse simulation - - a.addHarmonicPolymerBonds(wiggleDist=0.05) - # Bond distance will fluctuate +- 0.05 on average - - #a.addGrosbergRepulsiveForce(trunc=50) - # this will resolve chain crossings and will not let chain cross anymore - - # a.addGrosbergRepulsiveForce(trunc=5) - # this will let chains cross sometimes - - #a.addStiffness(k=4) - # K is more or less arbitrary, k=4 corresponds to presistence length of 4, - # k=1.5 is recommended to make polymer realistically flexible; k=8 is very stiff - - # If your simulation does not start, consider using energy minimization below - - # a.localEnergyMinimization() - # A very efficient algorithm to reach local energy minimum - # Use it to minimize energy if you're doing diffusion simulations - # If you're simulating dynamics of collapse or expansion, please do not use it - - #a.energyMinimization(stepsPerIteration=10) # increase to 100 for larger or more complex systems - # An algorithm to start a simulation - # Which works only with langevin integrator (but will not throw an error otherwise) - # Decreases a timestep, and then increases it slowly back to normal - - # -----------Running a simulation --------- - - a.save() # save original conformation - allDists = [] - for _ in xrange(10): # Do 10 blocks - a.doBlock(2000) # Of 2000 timesteps each - data = a.getData() - dists = np.sqrt(np.sum(np.diff(data, axis=0)**2, axis=1)) - allDists.append(dists) - - - allDists = np.concatenate(allDists) - plt.hist(allDists,100) - - plt.show() - - - -exampleOpenmm() -exit()