Improved atomic model backmapping for objects with long bonds and in general

atomicModel.py

@@ -531,20 +531,25 @@ class atomicModel():
             r0 = b.initialPosition
             r1 = simplerModelCoords[b.idx]
 
+            above = b.nodeAbove if b.nodeAbove is not None and \
+                    np.sum((b.nodeAbove.initialPosition-b.initialPosition)**2) < 10**2 else b
+            below = b.nodeBelow if b.nodeBelow is not None and \
+                    np.sum((b.nodeBelow.initialPosition-b.initialPosition)**2) < 10**2 else b
+
             if type(simplerModel).__name__ == "beadModelTwist" and \
                     b.orientationNode is not None:
-                above = b.nodeAbove if b.nodeAbove is not None else b
-                below = b.nodeBelow if b.nodeBelow is not None else b
                 o = b.orientationNode
                 assert(above != below)
 
-                x0 = o.initialPosition - r0
                 z0 = above.initialPosition - below.initialPosition
+                z1 = simplerModelCoords[above.idx] - simplerModelCoords[below.idx]
+                z0,z1 = [z/np.linalg.norm(z) for z in (z0,z1)]
 
-                r1 = simplerModelCoords[b.idx]
+                x0 = o.initialPosition - r0
                 x1 = simplerModelCoords[o.idx] - r1
-                z1 = simplerModelCoords[above.idx] - simplerModelCoords[below.idx]
-                x0,z0,x1,z1 = [x/np.linalg.norm(x) for x in (x0,z0,x1,z1)]
+                x0,x1 = [x-x.dot(z)*z for x,z in zip((x0,x1),(z0,z1))] # make x orthogonal to z
+                x0,x1 = [x/np.linalg.norm(x) for x in (x0,x1)]
+
                 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) )
 
@@ -555,8 +560,6 @@ class atomicModel():
                 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]
@@ -565,6 +568,7 @@ class atomicModel():
                     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)