Crossovers are now representated

segmentmodel.py

@@ -47,12 +47,12 @@ class ConnectableElement():
     def __init__(self, connections=[]):
         self.connections = connections
 
-    def get_connections_and_locations(self, type_=None):
+    def get_connections_and_locations(self, type_=None, exclude=[]):
         """ Returns a list with each entry of the form:
             connection, location_in_self, location_in_other """
         ret = []
         for c in self.connections:
-            if type_ is None or c.type_ == type_:
+            if type_ is None or c.type_ == type_ and type_ not in exclude:
                 if   c.A.container is self:
                     ret.append( [c, c.A, c.B] )
                 elif c.B.container is self:
@@ -367,6 +367,44 @@ class DoubleStrandedSegment(Segment):
             end3 = end3.end3
         self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
 
+    def crossover(self, nt, other, other_nt, strands=None, other_strands=None):
+        """ Add a crossover between two helices """
+        ## Validate other, nt, other_nt
+        ##   TODO
+
+        ## Fix input
+        if strands is None and other_strands is None:
+            ## Regular double crossover
+            strands = ['fwd','rev']
+            other_strands = list(strands[::-1])
+        elif strands is None:
+            strands = []
+            for d in other_strands:
+                if d == 'fwd': strands.append('rev')
+                elif d == 'rev': strands.append('fwd')
+        elif other_strands == None:
+            other_strands = []
+            for d in strands:
+                if d == 'fwd': other_strands.append('rev')
+                elif d == 'rev': other_strands.append('fwd')
+
+        if isinstance(strands,str):
+            strands = [strands]
+        if isinstance(other_strands,str):
+            other_strands = [other_strands]
+
+        ## Create locations, connections and add to segments
+        c = nt/self.num_nts
+        print(self.name,nt,c)
+        assert(c >= 0 and c <= 1)
+        loc = Location( self, address=c, type_ = strands )
+
+        c = other_nt/other.num_nts
+        print(other.name,other_nt,c)
+        assert(c >= 0 and c <= 1)
+        other_loc = Location( other, address=c, type_ = other_strands )
+        self._connect(other, Connection( loc, other_loc, type_="crossover" )) # TODO: use multiple crossover types
+
     ## Real work
     def _connect_ends(self, end1, end2, type_, force_connection):
         ## TODO remove self?
@@ -520,12 +558,12 @@ class SegmentModel(ArbdModel):
         self.useNonbondedScheme( nbDnaScheme )
 
 
-    def get_connections(self,type_=None):
+    def get_connections(self,type_=None,exclude=[]):
         """ Find all connections in model, without double-counting """
         added=set()
         ret=[]
         for s in self.segments:
-            items = [e for e in s.get_connections_and_locations(type_) if e[0] not in added]
+            items = [e for e in s.get_connections_and_locations(type_,exclude=exclude) if e[0] not in added]
             added.update([e[0] for e in items])
             ret.extend( items )
         return ret
@@ -732,6 +770,19 @@ class SegmentModel(ArbdModel):
                 A.particle = B.particle = b
                 ...
 
+        for c,A,B in self.get_connections(exclude="intrahelical"):
+            s1,s2 = [l.container for l in (A,B)]
+            if A.particle is not None and B.particle is not None:
+                continue
+            assert( A.particle is None )
+            assert( B.particle is None )
+
+            ## TODO: offload the work here to s1
+            a1,a2 = [l.address   for l in (A,B)]
+            ## TODO: if existing particle of same type is very nearby, use that instead
+            A.particle = s1._generate_one_bead(a1,0)
+            B.particle = s2._generate_one_bead(a2,0)
+
         """ Some tests """
         for c,A,B in self.get_connections("intrahelical"):
             for l in (A,B):
@@ -787,7 +838,7 @@ class SegmentModel(ArbdModel):
 
         """ Add intrahelical bond potentials """
         if self.DEBUG: print("Adding intrahelical bond potentials")
-        dists = dict()          # built for later use
+        dists = dict()          # intrahelical distances built for later use
         intra_beads = self._get_intrahelical_beads() 
         if self.DEBUG: print("  Adding %d bonds" % len(intra_beads))
         for b1,b2 in intra_beads:
@@ -836,7 +887,7 @@ class SegmentModel(ArbdModel):
                     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                
+                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2
 
             angle = self.get_angle_potential(k,180)
             parent.add_angle( b1, b2, b3, angle )
@@ -908,6 +959,69 @@ class SegmentModel(ArbdModel):
             pot = self.get_bond_potential(4,18.5)
             self.add_bond(b1,b2, pot)
 
+        for c,A,B in self.get_connections("crossover"):
+            b1,b2 = [loc.particle for loc in (c.A,c.B)]
+            pot = self.get_bond_potential(4,18.5)
+            self.add_bond(b1,b2, pot)
+
+            if not local_twist:
+                ## TODO
+                ...
+            else:
+                ## TODO: fwd/rev
+                k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
+                if 'orientation_bead' in b1.__dict__:
+                    t0 = 90 + 60
+                    # if f1: t0 -= 120
+                    o = b1.orientation_bead
+                    pot = self.get_angle_potential(k,t0)
+                    self.add_angle( o,b1,b2, pot )
+                else:
+                    t0 = 90 + 60
+                    # if f2: t0 -= 120
+                    o = b2.orientation_bead
+                    pot = self.get_angle_potential(k,t0)
+                    self.add_angle( b1,b2,o, pot )
+
+
+                u1,u2 = [b.intrahelical_neighbors[-1] for b in (b1,b2)] # TODO: fix this
+                d1,d2 = [b.intrahelical_neighbors[0] for b in (b1,b2)] # TODO: fix this
+
+                ## TODO make k depend on sep
+                k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
+                t0 = 90
+                pot = self.get_dihedral_potential(k,t0)
+                if 'orientation_bead' in b1.__dict__:
+                    o1 = b1.orientation_bead
+                    if u2 is not None:
+                        self.add_dihedral( o1,b1,b2,u2, pot )
+                    elif d2 is not None:
+                        self.add_dihedral( o2,b2,b1,d2, pot )
+                if 'orientation_bead' in b2.__dict__:
+                    o2 = b2.orientation_bead
+                    if u1 is not None:
+                        self.add_dihedral( o2,b2,b1,u1, pot )
+                    elif d1 is not None:
+                        self.add_dihedral( o2,b2,b1,d1, pot )
+
+                if 'orientation_bead' in b1.__dict__ and 'orientation_bead' in b2.__dict__:
+                    if u1 is not None and u2 is not None:
+                        t0 = 0
+                        pot = self.get_dihedral_potential(k,t0)
+                        self.add_dihedral( u1,b1,b2,u2,  pot )
+                    elif d1 is not None and d2 is not None:
+                        t0 = 0
+                        pot = self.get_dihedral_potential(k,t0)
+                        self.add_dihedral( d1,b1,b2,d2, pot )
+                    elif d1 is not None and u2 is not None:
+                        t0 = 180
+                        pot = self.get_dihedral_potential(k,t0)
+                        self.add_dihedral( d1,b1,b2,u2, pot )
+                    elif u1 is not None and d2 is not None:
+                        t0 = 180
+                        pot = self.get_dihedral_potential(k,t0)
+                        self.add_dihedral( u1,b1,b2,d2, pot )
+
 
     # def get_bead(self, location):
     #     if type(location.container) is not list: