Cleaned up a few bugs

segmentmodel.py

@@ -685,7 +685,8 @@ class DoubleStrandedSegment(Segment):
                  segment_model = None,
                  local_twist = False,
                  num_turns = None,
-                 start_orientation = None):
+                 start_orientation = None,
+                 twist_persistence_length = 90 ):
         
         self.helical_rise = 10.44
         self.distance_per_nt = 3.4
@@ -702,6 +703,7 @@ class DoubleStrandedSegment(Segment):
         if start_orientation is None:
             start_orientation = np.eye(3) # np.array(((1,0,0),(0,1,0),(0,0,1)))
         self.start_orientation = start_orientation
+        self.twist_persistence_length = twist_persistence_length
 
         self.nicks = []
 
@@ -1572,11 +1574,11 @@ class SegmentModel(ArbdModel):
                     new_cl.append(cl)
                     lastParticle = A.particle
             crossovers = new_cl
-            for i in range(len(crossovers)-1):
+            for i in range(len(crossovers)-2):
                 c1,A1,B1,dir1 = crossovers[i]
-                c2,A2,B2,dir2 = crossovers[j]
+                c2,A2,B2,dir2 = crossovers[i+1]
                 s1,s2 = [l.container for l in (A1,A2)]
-                sep = A1.get_nt_position(s1) - A2.get_nt_position(s2)
+                sep = A1.particle.get_nt_position(s1) - A2.particle.get_nt_position(s2)
                 sep = np.abs(sep)
 
                 n1,n2,n3,n4 = (B1.particle, A1.particle, A2.particle, B2.particle)
@@ -1589,28 +1591,28 @@ class SegmentModel(ArbdModel):
 
                 ## 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 = s1.twistPersistenceLength/0.34  # set semi-arbitrarily as there is a large spread in literature
+                Lp = s1.twist_persistence_length/0.34  # set semi-arbitrarily as there is a large spread in literature
 
-                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) - exp(-sep/Lp)
-                k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
+                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) )
                 k = k[0][0] * 2*kT*0.00030461742
                         
                 t0 = sep*(360.0/s1.twist_per_nt)
                 
                 # pdb.set_trace()
-                if c1.on_fwd_strand: t0 -= 120
+                if A1.on_fwd_strand: t0 -= 120
                 if dir1 != dir2:
-                    c2_on_fwd = not c2.on_fwd_strand
+                    A2_on_fwd = not A2.on_fwd_strand
                 else:
-                    c2_on_fwd = c2.on_fwd_strand
-                if c2_on_fwd: t0 += 120
+                    A2_on_fwd = A2.on_fwd_strand
+                if A2_on_fwd: t0 += 120
    
                 # t0 = (t0 % 360
 
                 # if n2.idx == 0:
                 #     print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
                 pot = self.get_dihedral_potential(k,t0)
-                dihedrals.add( n1,n2,n3,n4, pot )
+                self.add_dihedral( n1,n2,n3,n4, pot )
                 
 
     def walk_through_helices(segment, direction=1, processed_segments=None):