From 76f781fd80512d0c004c74238858bf8cbf8e393d Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 18 Dec 2018 13:48:36 -0600
Subject: [PATCH] Fixed code adding dihedral angles between adjacent crossovers

---
 mrdna/segmentmodel.py | 16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 7782a80..a6d5632 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -2317,9 +2317,6 @@ class SegmentModel(ArbdModel):
         for crossovers in self.get_consecutive_crossovers():
             if local_twist: break
             ## filter crossovers
-            new_cl = []
-            lastParticle = None
-            crossovers = new_cl
             for i in range(len(crossovers)-2):
                 c1,A1,B1,dir1 = crossovers[i]
                 c2,A2,B2,dir2 = crossovers[i+1]
@@ -2333,10 +2330,6 @@ class SegmentModel(ArbdModel):
 
                 """
                 <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
@@ -2348,7 +2341,7 @@ class SegmentModel(ArbdModel):
                     k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
                     return k[0][0] * 2*kT*0.00030461742
 
-                k = get_spring( max(sep,1.0) )
+                k = get_spring( max(sep,2) )
                 t0 = sep*s1.twist_per_nt # TODO weighted avg between s1 and s2
                 
                 # pdb.set_trace()
@@ -2363,9 +2356,10 @@ class SegmentModel(ArbdModel):
 
                 # if n2.idx == 0:
                 #     print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
-                if sep == 0:
-                    pot = self.get_angle_potential(k,t0)
-                    self.add_angle( n1,n2,n4, pot )
+                if sep == 0 and n1 is not n4:
+                    # pot = self.get_angle_potential(k,t0)
+                    # self.add_angle( n1,n2,n4, pot )
+                    pass
                 else:
                     pot = self.get_dihedral_potential(k,t0)
                     self.add_dihedral( n1,n2,n3,n4, pot )
-- 
GitLab