diff --git a/segmentmodel.py b/segmentmodel.py
index b00a29fb90caee5ec9efcc72f17b3b6af15aab02..1b7300ee7499b42c354f2ce69fa63a24c5f6d244 100644
--- a/segmentmodel.py
+++ b/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