Fixed spring constant for dihedral that spans crossovers, keeping helices...

Fixed spring constant for dihedral that spans crossovers, keeping helices parallel (note this differs from prototyped web service where 90 degree angle potentials were used)

Changed dna bead type clustering to have higher resolution

Streamlined some small parts of the code

segmentmodel.py

@@ -1654,10 +1654,9 @@ class SegmentModel(ArbdModel):
         data = np.array([b.num_nt for b in beads])[:,np.newaxis]
         order = int(2-np.log10(2*max_basepairs_per_bead)//1)
         try:
-            clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/80, criterion="distance")
+            clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/500, criterion="distance")
             cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
         except:
-            ## TODO debug
             clusters = data
             cluster_size = np.arange(len(data))+1
         for b,c in zip(beads,clusters):
@@ -1667,15 +1666,30 @@ class SegmentModel(ArbdModel):
         data = np.array([b.num_nt for b in beads])[:,np.newaxis]
         order = int(2-np.log10(max_nucleotides_per_bead)//1)
         try:
-            clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/40, criterion="distance")
+            clusters = hcluster.fclusterdata(data, float(max_nucleotides_per_bead)/500, criterion="distance")
             cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
         except:
-            ## TODO debug
             clusters = data
             cluster_size = np.arange(len(data))+1
         for b,c in zip(beads,clusters):
             _assign_bead_type(b, cluster_size[c-1], decimals=order)
 
+
+        # for bead in [b for s in segments for b in s]:
+        #     num_nt0 = bead.num_nt
+        #     # bead.num_nt = np.around( np.float32(num_nt), decimals=decimals )
+        #     key = (bead.type_.name[0].upper(), bead.num_nt)
+        #     if key in beadtype_s:
+        #         bead.type_ = beadtype_s[key]
+        #     else:
+        #         t = deepcopy(bead.type_)
+        #         t.__dict__["nts"] = bead.num_nt*2 if t.name[0].upper() in ("D","O") else bead.num_nt
+        #         # t.name = t.name + "%03d" % (t.nts*10**decimals)
+        #         t.name = t.name + "%.16f" % (t.nts)
+        #         print( "{} --> {} ({})".format(num_nt0, bead.num_nt, t.name) )
+        #         beadtype_s[key] = bead.type_ = t
+
+
         """ Update bead indices """
         self._countParticleTypes() # probably not needed here
         self._updateParticleOrder()
@@ -1853,30 +1867,28 @@ class SegmentModel(ArbdModel):
 
 
             k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+            a = None
             if u1 is not None and u2 is not None:
                 t0 = 0
-                k = k_fn( 0.5*(dists[b1][u1]+dists[b2][u2]) )
-                pot = self.get_dihedral_potential(k,t0)
-                self.add_dihedral( u1,b1,b2,u2,  pot )
+                a,b,c,d = (u1,b1,b2,u2)
             elif d1 is not None and d2 is not None:
                 t0 = 0
-                k = k_fn( 0.5*(dists[b1][d1]+dists[b2][d2]) )
-                pot = self.get_dihedral_potential(k,t0)
-                self.add_dihedral( d1,b1,b2,d2, pot )
+                a,b,c,d = (d1,b1,b2,d2 )
             elif d1 is not None and u2 is not None:
                 t0 = 180
-                k = k_fn( 0.5*(dists[b1][d1]+dists[b2][u2]) )
-                pot = self.get_dihedral_potential(k,t0)
-                self.add_dihedral( d1,b1,b2,u2, pot )
+                a,b,c,d = (d1,b1,b2,u2)
             elif u1 is not None and d2 is not None:
                 t0 = 180
-                k = k_fn( 0.5*(dists[b1][u1]+dists[b2][d2]) )
+                a,b,c,d = (u1,b1,b2,d2)
+
+            if a is not None:
+                k = k_fn( dists[b][a]+dists[c][d] )
                 pot = self.get_dihedral_potential(k,t0)
-                self.add_dihedral( u1,b1,b2,d2, pot )
+                self.add_dihedral( a,b,c,d,  pot )
 
 
             if local_twist:
-                k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
+                k = k_fn(1)
                 pot = self.get_angle_potential(k,120)
                 if 'orientation_bead' in b1.__dict__:
                     o = b1.orientation_bead
@@ -1933,11 +1945,13 @@ class SegmentModel(ArbdModel):
 
                 n1,n2,n3,n4 = (B1.particle, A1.particle, A2.particle, B2.particle)
 
-                ## <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
+                """
+                <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