Updated orientation of ssDNA in atomic model so that z-axis is aligned with...

Updated orientation of ssDNA in atomic model so that z-axis is aligned with vector connecting adjacent nts; untested

atomicModel.py

@@ -517,6 +517,9 @@ class atomicModel():
             helixCutoff = dsDnaHelixNeighborDist if b.type[0] == 'd' else ssDnaHelixNeighborDist
             allCutoff = dsDnaAllNeighborDist if b.type[0] == 'd' else ssDnaAllNeighborDist
 
+            r0 = b.initialPosition
+            r1 = simplerModelCoords[b.idx]
+
             if type(simplerModel).__name__ == "beadModelTwist" and \
                     b.orientationNode is not None:
                 above = b.nodeAbove if b.nodeAbove is not None else b
@@ -524,7 +527,6 @@ class atomicModel():
                 o = b.orientationNode
                 assert(above != below)
 
-                r0 = b.initialPosition
                 x0 = o.initialPosition - r0
                 z0 = above.initialPosition - below.initialPosition
 
@@ -534,7 +536,6 @@ class atomicModel():
                 x0,z0,x1,z1 = [x/np.linalg.norm(x) for x in (x0,z0,x1,z1)]
                 y0,y1 = [np.cross(z,x) for x,z in zip((x0,x1),(z0,z1))]
                 R,tmp1,tmp2 = minimizeRmsd( (x0,y0,z0), (x1,y1,z1) )
-                cOld,cNew = (r0,r1)
 
             else:
                 ids = simplerModel._getNeighborhoodIds(b, simplerModelCoords, helixCutoff, allCutoff)
@@ -542,6 +543,20 @@ class atomicModel():
                 posNew = np.array( [simplerModelCoords[i] for i in ids] )
                 R,cOld,cNew = minimizeRmsd( posOld, posNew )
 
+                if b.type[0] == 's':
+                    above = b.nodeAbove if b.nodeAbove is not None else b
+                    below = b.nodeBelow if b.nodeBelow is not None else b
+                    assert(above != below)
+                    z0 = above.initialPosition - below.initialPosition
+                    z1 = simplerModelCoords[above.idx] - simplerModelCoords[below.idx]
+
+                    zR = z0.dot(R)
+                    cross = np.cross(z1,zR) / (np.linalg.norm(z1)*np.linalg.norm(zR))
+                    angle = np.arcsin(np.linalg.norm(cross))*180/np.pi
+                    R2 = rotationAboutAxis(cross,angle)
+                    R = R.dot(R2)
+
+            cOld,cNew = (r0,r1)
             trans[b.idx] = (R,cOld,cNew)
 
         ## Optionally smooth orientations