diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index be162ec19fc2ec59f2604731d15d316da9942b4b..7782a80b08d0932076baca65d8b69720c36c08b7 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -2319,11 +2319,6 @@ class SegmentModel(ArbdModel):
             ## filter crossovers
             new_cl = []
             lastParticle = None
-            for cl in crossovers:
-                c,A,B,d = cl
-                if A.particle is not lastParticle:
-                    new_cl.append(cl)
-                    lastParticle = A.particle
             crossovers = new_cl
             for i in range(len(crossovers)-2):
                 c1,A1,B1,dir1 = crossovers[i]
@@ -2332,6 +2327,8 @@ class SegmentModel(ArbdModel):
                 sep = A1.particle.get_nt_position(s1) - A2.particle.get_nt_position(s2)
                 sep = np.abs(sep)
 
+                assert(sep >= 0)
+
                 n1,n2,n3,n4 = (B1.particle, A1.particle, A2.particle, B2.particle)
 
                 """
@@ -2346,10 +2343,12 @@ class SegmentModel(ArbdModel):
                 ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
                 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) - np.exp(-sep/Lp)
-                k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
-                k = k[0][0] * 2*kT*0.00030461742
-                        
+                def get_spring(sep):
+                    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
+
+                k = get_spring( max(sep,1.0) )
                 t0 = sep*s1.twist_per_nt # TODO weighted avg between s1 and s2
                 
                 # pdb.set_trace()
@@ -2364,8 +2363,12 @@ class SegmentModel(ArbdModel):
 
                 # if n2.idx == 0:
                 #     print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
-                pot = self.get_dihedral_potential(k,t0)
-                self.add_dihedral( n1,n2,n3,n4, pot )
+                if sep == 0:
+                    pot = self.get_angle_potential(k,t0)
+                    self.add_angle( n1,n2,n4, pot )
+                else:
+                    pot = self.get_dihedral_potential(k,t0)
+                    self.add_dihedral( n1,n2,n3,n4, pot )
 
         # ## remove duplicate potentials; ## TODO ensure that they aren't added twice in the first place? 
         # self.remove_duplicate_terms()