Added models with twist

arbdmodel.py

@@ -765,8 +765,6 @@ component "data" value 3
 
         with open(self._dihedral_filename,'w') as fh:
             for b in self.get_dihedrals():
-                if type(b) is not str:
-                    b = b.filename()
                 item = tuple([p.idx for p in b[:-1]] + [str(b[-1])])
                 fh.write("DIHEDRAL %d %d %d %d %s\n" % item)
 

nonbonded.py

@@ -58,7 +58,7 @@ class TabulatedPotential(NonbondedScheme):
 
 ## Bonded potentials
 class HarmonicPotential():
-    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, prefix="potentials/"):
+    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
         self.k = k
         self.r0 = r0
         self.rRange = rRange
@@ -67,6 +67,7 @@ class HarmonicPotential():
         self.prefix = prefix
         self.periodic = False
         self.type_ = "None"
+        self.max_potential = max_potential
         self.kscale_ = None     # only used for 
 
     def filename(self):
@@ -98,14 +99,28 @@ class HarmonicPotential():
             f[f<-maxForce] = -maxForce
             u[0] = 0
             u[1:] = np.cumsum(f*np.diff(r))
+
+        if self.max_potential is not None:
+            f = np.diff(u)/np.diff(r)
+            ids = np.where( 0.5*(u[1:]+u[:-1]) > self.max_potential )[0]
+
+            w = np.sqrt(2*self.max_potential/self.k)
+            drAvg = 0.5*(np.abs(dr[ids]) + np.abs(dr[ids+1]))
+
+            f[ids] = f[ids] * np.exp(-(drAvg-w)/(w))
+            u[0] = 0
+            u[1:] = np.cumsum(f*np.diff(r))
+
+        u = u - np.min(u)
+
         np.savetxt( self.filename(), np.array([r, u]).T, fmt="%f" )
 
     def __hash__(self):
         assert(self.type_ != "None")
-        return hash((self.type_, self.k, self.r0, self.rRange, self.resolution, self.maxForce, self.prefix, self.periodic))
+        return hash((self.type_, self.k, self.r0, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
 
     def __eq__(self, other):
-        for a in ("type_", "k", "r0", "rRange", "resolution", "maxForce", "prefix", "periodic"):
+        for a in ("type_", "k", "r0", "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
             if self.__dict__[a] != other.__dict__[a]:
                 return False
         return True
@@ -116,20 +131,20 @@ class HarmonicPotential():
 #         self.kscale_ = 1.0
 
 class HarmonicBond(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, prefix="potentials/"):
-        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix)
+    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
+        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
         self.type_ = "bond"
         self.kscale_ = 1.0
 
 class HarmonicAngle(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, prefix="potentials/"):
-        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix)
+    def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
+        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
         self.type_ = "angle"
         self.kscale_ = (180.0/np.pi)**2
 
 class HarmonicDihedral(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(-180,180), resolution=0.1, maxForce=None, prefix="potentials/"):
-        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix)
+    def __init__(self, k, r0, rRange=(-180,180), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
+        HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
         self.periodic = True
         self.type_ = "dihedral"
         self.kscale_ = (180.0/np.pi)**2

segmentmodel.py

 import numpy as np
 from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
+from coords import rotationAboutAxis
 from nonbonded import *
 from copy import copy, deepcopy
 from nbPot import nbDnaScheme
 
+from scipy.special import erf
+import scipy.optimize as opt
+
 
 """
-# TODO:
+TODO:
+ - document
+ - 
 
 
 """
@@ -52,6 +58,13 @@ class Segment(ConnectableElement, Group):
                                   mass = 300,
                                   radius = 3,                 
                               )
+    orientation_particle = ParticleType("O",
+                                        diffusivity = 100,
+                                        mass = 300,
+                                        radius = 1,
+                                    )
+    # orientation_bond = HarmonicBond(10,2)
+    orientation_bond = HarmonicBond(30,1.5)
 
     ssDNA_particle = ParticleType("S",
                                   diffusivity = 43.5,
@@ -67,19 +80,24 @@ class Segment(ConnectableElement, Group):
         Group.__init__(self, name, children=[])
         ConnectableElement.__init__(self, connections=[])
 
+        self.start_orientation = None
+        self.twist_per_nt = 0
+
+        self.beads = [c for c in self.children] # self.beads will not contain orientation beads
+
         self._bead_model_generation = 0    # TODO: remove?
         self.segment_model = segment_model # TODO: remove?
 
-        # self.end5 = Location( self, address=0, type_= "end5" )
-        # self.end3 = Location( self, address=-1, type_ = "end3" )
-
         self.num_nts = num_nts
         if end_position is None:
             end_position = np.array((0,0,self.distance_per_nt*num_nts)) + start_position
         self.start_position = start_position
         self.end_position = end_position
 
-    def _generate_one_bead(self, pos, nts):
+    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead ):
+        raise NotImplementedError
+
+    def _generate_one_bead(self, pos, nts, orientation=None):
         raise NotImplementedError
 
     def _assign_particles_to_locations(self):
@@ -87,49 +105,52 @@ class Segment(ConnectableElement, Group):
 
     def get_all_consecutive_beads(self, number):
         assert(number >= 1)
-        ## Assume that consecutive beads in self.children are bonded
+        ## Assume that consecutive beads in self.beads are bonded
         ret = []
-        for i in range(len(self.children)-number+1):
-            tmp = [self.children[i+j] for j in range(0,number)]
+        for i in range(len(self.beads)-number+1):
+            tmp = [self.beads[i+j] for j in range(0,number)]
             ret.append( tmp )
         return ret
 
     def get_beads_before_bead(self, bead, number, inclusive=False):
-        ## Assume that consecutive beads in self.children are bonded
-        i = self.children.index(bead)
-        l = len(self.children)
+        ## Assume that consecutive beads in self.beads are bonded
+        i = self.beads.index(bead)
+        l = len(self.beads)
         if i-number < 0:
             raise Exception("Not enough beads after bead")
         
         start = 1
         if inclusive: start = 0
-        return [self.children[i-j] for j in range(start,number)]
+        return [self.beads[i-j] for j in range(start,number)]
 
     def get_beads_after_bead(self, bead, number, inclusive=False):
-        ## Assume that consecutive beads in self.children are bonded
-        i = self.children.index(bead)
-        l = len(self.children)
+        ## Assume that consecutive beads in self.beads are bonded
+        i = self.beads.index(bead)
+        l = len(self.beads)
         if i+number >= l:
             raise Exception("Not enough beads after bead")
         
         start = 1
         if inclusive: start = 0
-        return [self.children[i+i] for j in range(start,number)]
+        return [self.beads[i+i] for j in range(start,number)]
 
     # def get_bead_pairs_within(self, cutoff):
     #     for b1,b2 in self.get_all_consecutive_beads(self, number)
 
 
-    def _generate_beads(self, bead_model, max_nts_per_bead=4):
+    def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
         
         """ Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
 
+        ## TODO: decide whether to remove bead_model argument
+
         # self._bead_model_generation += 1
         # self._bead_model_max_nts_per_bead = max_nts_per_bead
 
         direction = self.end_position - self.start_position
-        num_beads = (self.num_nts // max_nts_per_bead) + 1
+        num_beads = self._get_num_beads( max_basepairs_per_bead, max_nucleotides_per_bead )
         nts_per_bead = float(self.num_nts)/num_beads
+        twist_per_bead = nts_per_bead * self.twist_per_nt
 
         last = None
 
@@ -140,9 +161,24 @@ class Segment(ConnectableElement, Group):
 
             s = i*float(nts_per_bead)/(self.num_nts) # contour
             pos = direction * s + self.start_position
+            if self.start_orientation is not None:
+                axis = self.start_orientation.dot( np.array((0,0,1)) )
+                orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=False )
+                ## TODO: ensure this is correct
+                # orientation = self.start_orientation.dot(orientation) # .dot( self.start_orientation )
+                orientation = orientation.dot( self.start_orientation )
+            else:
+                orientation = None
 
-            b = self._generate_one_bead(pos,nts)
+            b = self._generate_one_bead(pos,nts,orientation)
             self.children.append(b)
+            self.beads.append(b) # don't add orientation bead
+
+            if "orientation_bead" in b.__dict__:
+                o = b.orientation_bead
+                self.children.append(o)
+                self.add_bond(b,o, Segment.orientation_bond)
+                
             # if last is not None:
             #     self.add_bond( i=last, j=b, bond="ssdna" )
             # last = b
@@ -163,15 +199,26 @@ class DoubleStrandedSegment(Segment):
     def __init__(self, name, num_nts, start_position = np.array((0,0,0)),
                  end_position = None, 
                  segment_model = None,
-                 twist = None,
+                 local_twist = False,
+                 num_turns = None,
                  start_orientation = None):
-
-        self.distance_per_nt = 5
+        
+        self.helical_rise = 10.44
+        self.distance_per_nt = 3.4
         Segment.__init__(self, name, num_nts, 
                          start_position,
                          end_position, 
                          segment_model)
 
+        self.local_twist = local_twist
+        if num_turns is None:
+            num_turns = float(num_nts) / self.helical_rise
+        self.twist_per_nt = float(360 * num_turns) / num_nts
+
+        if start_orientation is None:
+            start_orientation = np.array(((1,0,0),(0,1,0),(0,0,1)))
+        self.start_orientation = start_orientation
+
         self.nicks = []
 
         self.start5 = Location( self, address=0, type_= "end5" )
@@ -210,13 +257,28 @@ class DoubleStrandedSegment(Segment):
 
         end1.container._connect( end2.container, Connection( end1, end2, type_="intrahelical" ) )
 
-    def _generate_one_bead(self, pos, nts):
-        return PointParticle( Segment.dsDNA_particle, pos, nts,
-                              num_nts=nts, parent=self )
+    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
+        return (self.num_nts // max_basepairs_per_bead) + 1
+
+    def _generate_one_bead(self, pos, nts, orientation=None):
+        if self.local_twist:
+            assert(orientation is not None)
+            # opos = pos + np.array((2,0,0)).dot(orientation)
+            opos = pos + orientation.dot( np.array((Segment.orientation_bond.r0,0,0)) )
+            o = PointParticle( Segment.orientation_particle, opos, nts,
+                               num_nts=nts, parent=self )
+            bead = PointParticle( Segment.dsDNA_particle, pos, nts,
+                                  num_nts=nts, parent=self, orientation_bead=o )
+
+        else:
+            bead = PointParticle( Segment.dsDNA_particle, pos, nts,
+                                  num_nts=nts, parent=self )
+        return bead
+
 
     def _assign_particles_to_locations(self):
-        self.start3.particle =  self.start5.particle = self.children[0]
-        self.end3.particle =  self.end5.particle = self.children[-1]
+        self.start3.particle =  self.start5.particle = self.beads[0]
+        self.end3.particle   =  self.end5.particle   = self.beads[-1]
 
     def _generate_atomic(self, atomic_model):
         ...
@@ -277,7 +339,10 @@ class SingleStrandedSegment(Segment):
 
         self._connect( other.container, Connection( my_end, other, type_="intrahelical" ) )
 
-    def _generate_one_bead(self, pos, nts):
+    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
+        return (self.num_nts // max_nucleotides_per_bead) + 1
+
+    def _generate_one_bead(self, pos, nts, orientation=None):
         return PointParticle( Segment.ssDNA_particle, pos, nts,
                               num_nts=nts, parent=self )
 
@@ -290,9 +355,9 @@ class SingleStrandedSegment(Segment):
     
 
 class SegmentModel(ArbdModel):
-    def __init__(self, segments = [], 
-                 max_basepairs_per_bead = 7,
-                 max_nucleotides_per_bead = 4,
+    def __init__(self, segments=[], local_twist=True,
+                 max_basepairs_per_bead=7,
+                 max_nucleotides_per_bead=4,
                  dimensions=(1000,1000,1000), temperature=291,
                  timestep=50e-6, cutoff=50, 
                  decompPeriod=10000, pairlistDistance=None, 
@@ -304,13 +369,13 @@ class SegmentModel(ArbdModel):
                            decompPeriod, pairlistDistance=None,
                            nonbondedResolution=0)
 
-        self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
-        self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
+        # self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
+        # self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
         self.children = self.segments = segments
 
         self._bonded_potential = dict() # cache bonded potentials
 
-        self._generate_bead_model(segments, max_nucleotides_per_bead, max_nucleotides_per_bead)
+        self._generate_bead_model(segments, max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
     
     def _get_intrahelical_beads(self):
         ret = []
@@ -393,15 +458,16 @@ class SegmentModel(ArbdModel):
 
     #     ret
 
-    def _get_potential(self, type_, kSpring, d):
+    def _get_potential(self, type_, kSpring, d, max_potential = None):
         key = (type_,kSpring,d)
         if key not in self._bonded_potential:
             if type_ == "bond":
-                self._bonded_potential[key] = HarmonicBond(kSpring,d)
+                self._bonded_potential[key] = HarmonicBond(kSpring,d, max_potential=max_potential)
             elif type_ == "angle":
-                self._bonded_potential[key] = HarmonicAngle(kSpring,d) # , resolution = 1, maxForce=0.1)
+                self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
+                # , resolution = 1, maxForce=0.1)
             elif type_ == "dihedral":
-                self._bonded_potential[key] = HarmonicDihedral(kSpring,d)
+                self._bonded_potential[key] = HarmonicDihedral(kSpring,d, max_potential=max_potential)
             else:
                 raise Exception("Unhandled potential type '%s'" % type_)
         return self._bonded_potential[key]
@@ -409,18 +475,51 @@ class SegmentModel(ArbdModel):
         return self._get_potential("bond", kSpring, d)
     def get_angle_potential(self, kSpring, d):
         return self._get_potential("angle", kSpring, d)
-    def get_dihedral_potential(self, kSpring, d):
-        return self._get_potential("dihedral", kSpring, d)
+    def get_dihedral_potential(self, kSpring, d, max_potential=None):
+        while d > 180: d-=360
+        while d < -180: d+=360
+        return self._get_potential("dihedral", kSpring, d, max_potential)
 
 
+    def _getParent(self, *beads ):
+        ## TODO: test
+        if np.all( [b1.parent == b2.parent 
+                    for b1,b2 in zip(beads[::2],beads[1::2])] ):
+            return beads[0].parent
+        else:
+            return self
+
+    def _get_twist_spring_constant(self, sep):
+        """ sep in nt """
+        kT = 0.58622522         # kcal/mol
+        twist_persistence_length = 90  # set semi-arbitrarily as there is a large spread in literature
+        ## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
+        ##   Assume A is small
+        ## int[B_] :=  Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
+        ##             Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]
+
+        ## Actually, without assumptions I get fitFun below
+        ## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
+        ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
+        Lp = twist_persistence_length/0.34 
+
+        fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - np.exp(-sep/Lp)
+        k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
+        return k[0][0] * 2*kT*0.00030461742
+        
+
     def _generate_bead_model(self, segments,
                              max_basepairs_per_bead = 7,
-                             max_nucleotides_per_bead = 4):
+                             max_nucleotides_per_bead = 4,
+                             local_twist=False
+    ):
 
 
         """ Generate beads """
         for s in segments:
-            s._generate_beads( max_nucleotides_per_bead )
+            if local_twist:
+                s.local_twist = True
+            s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
 
         """ Combine beads at junctions as needed """
         for s in segments:
@@ -441,9 +540,11 @@ class SegmentModel(ArbdModel):
                         t.__dict__["nts"] = b.num_nts*2
                     elif key[0] == "S":
                         t.__dict__["nts"] = b.num_nts
+                    elif key[0] == "O":
+                        t.__dict__["nts"] = b.num_nts
                     else:
                         raise Exception("TODO")
-                    print(t.nts)
+                    # print(t.nts)
                     t.name = t.name + "%03d" % (100*t.nts)
                     beadtype_s[key] = b.type_ = t
 
@@ -451,13 +552,12 @@ class SegmentModel(ArbdModel):
         """ Add intrahelical bond potentials """
         dists = dict()          # built for later use
         for b1,b2 in self._get_intrahelical_beads():
+            parent = self._getParent(b1,b2)
             if b1.parent == b2.parent:
                 sep = 0.5*(b1.num_nts+b2.num_nts)
-                parent = b1.parent
             else:
                 sep = 1
-                parent = self
-            
+
             conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
             if b1.type_.name[0] == "D" and b1.type_.name[0] == "D":
                 elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
@@ -493,10 +593,7 @@ class SegmentModel(ArbdModel):
             else:
                 sep += 1
                 
-            if b1.parent == b2.parent and b2.parent == b3.parent:
-                parent = b1.parent
-            else:
-                parent = self
+            parent = self._getParent(b1,b2,b3)
 
             kT = 0.58622522         # kcal/mol
             if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
@@ -504,6 +601,8 @@ class SegmentModel(ArbdModel):
                 ## <cos(q)> ~ 1 - 3/4A                                                                                            
                 ## where A = k_spring / (2 kT)                                                                                    
                 k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+                if local_twist:
+                    k *= 0.5    # halve because orientation beads have similar springs
             else:
                 ## TODO: get correct number from ssDNA model
                 k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2                
@@ -528,11 +627,48 @@ class SegmentModel(ArbdModel):
         for b1 in beads:
             exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
         for b1,b2 in exclusions:
-            if b1.parent == b2.parent:
-                parent = b1.parent
-            else:
-                parent = self
+            parent = self._getParent(b1,b2)
             parent.add_exclusion( b1, b2 )
+
+        """ Twist potentials """
+        if local_twist:
+
+            ## TODO: decide whether to add bond here
+            # """Add bonds between orientation bead and parent"""
+            # for s in self.segments:
+            #     for b,o in zip(s.children[::2],s.children[1::2]):
+            #         s.add_bond(
+                    
+
+            for b1 in beads:
+                if "orientation_bead" not in b1.__dict__: continue
+                for b2,sep in dists[b1]:
+                    if "orientation_bead" not in b2.__dict__: continue
+
+                    p1,p2 = [b.parent for b in (b1,b2)]
+                    o1,o2 = [b.orientation_bead for b in (b1,b2)]
+
+                    parent = self._getParent( b1, b2 )
+                    
+
+                    """ Add heuristic 90 degree potential to keep orientation bead orthogonal """
+                    k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+                    pot = self.get_angle_potential(k,90)
+                    parent.add_angle(o1,b1,b2, pot)
+                    parent.add_angle(b1,b2,o2, pot)
+                                        
+                    ## TODO: improve this
+                    twist_per_nt = 0.5 * (p1.twist_per_nt + p2.twist_per_nt)
+                    angle = sep*twist_per_nt
+                    if angle > 360 or angle < -360:
+                        raise Exception("The twist between beads is too large")
+                        
+                    k = self._get_twist_spring_constant(sep)
+                    pot = self.get_dihedral_potential(k,angle,max_potential=1)
+                    parent.add_dihedral(o1,b1,b2,o2, pot)
+
+                
+                
         
         """ Add connection potentials """
         ## TODO

segmentmodel_example.py

@@ -14,7 +14,7 @@ PATH=/home/cmaffeo2/anaconda3/bin:$PATH; source activate cadnano
 
 if __name__ == "__main__":
 
-    seg1 = DoubleStrandedSegment("strand", num_nts = 46)
+    seg1 = DoubleStrandedSegment("strand", num_nts = 46, local_twist=True)
 
     seg2 = SingleStrandedSegment("strand", 
                                  start_position = seg1.end_position + np.array((5,0,5)),
@@ -34,7 +34,10 @@ if __name__ == "__main__":
     seg1.connect_end5(seg3)
     seg1.connect_start3(seg3)
 
-    model = SegmentModel( [seg1, seg2, seg3],
+    model = SegmentModel( [seg1, seg2, seg3], 
+                          local_twist = True,
+                          max_basepairs_per_bead = 5,
+                          max_nucleotides_per_bead = 5,
                           dimensions=(5000,5000,5000),
                       )