Updated run script to use smaller rmsd

run.py

@@ -11,6 +11,7 @@ from beadModelTwist import beadModelTwist
 from atomicModel import atomicModel
 
 import mdtraj as md
+from coords import minimizeRmsd # mdtraj.superpose not working
 
 import multiprocessing as mp
 
@@ -21,9 +22,9 @@ twistLp = 90
 
 # arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd"
 arbd="/home/cmaffeo2/development/cuda/arbd.dbg/src/arbd" # reduced the mem footprint cause vmd
-# namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multicore-CUDA/namd2"
+namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multicore-CUDA/namd2"
 # namd="/home/wilsja/namd/NAMD_CVS-2016-11-09_Linux-x86_64-multicore-CUDA/namd2"
-namd="/home/cmaffeo2/development/namd-bin/NAMD_2.12_Linux-x86_64-multicore-CUDA/namd2"
+# namd="/home/cmaffeo2/development/namd-bin/NAMD_2.12_Linux-x86_64-multicore-CUDA/namd2"
 
 def readSequenceFile(sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat'):
     seq = []
@@ -41,27 +42,30 @@ def readArbdCoords(fname):
             coords.append([float(x) for x in line.split()[1:]])
     return np.array(coords)
 
-def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=20):
+def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5):
     ## align trajectory to pdb
-    ref = md.load(pdb, top=psf)
-    sel = md.load(dcd, top=psf)
+    ref = md.load(pdb, top=pdb)
+    sel = md.load(dcd, top=pdb)
     # ref = sel[-1]
-    sel.superpose(ref)
-
-    ## Find where rmsd compared to last frame increases 
-    ref = sel[-1]
-    rmsds = md.rmsd(sel,ref)*10
-    t0 = np.where(rmsds > rmsdThreshold)[0]
-    if len(t0) == 0:
-        t0 = 0
-    else:
-        t0 = t0[-1]
-
-    pos = sel.xyz[t0:]
-    pos =  np.mean(pos, axis=0)
+    ids = ref.topology.select("name =~ 'd.*'")
+    assert(len(ids) > 3)
+
+    # r0 = ref.xyz[0,ids,:]
+    r0 = sel.xyz[-1,ids,:]
+    for t in range(len(sel.xyz)-2,-1,-1):
+        # print(t)
+        R,comA,comB = minimizeRmsd(sel.xyz[t,ids,:],r0)
+        sel.xyz[t,:,:] = np.array( [(r-comA).dot(R)+comB for r in sel.xyz[t]] )
+        rmsd = np.mean( (sel.xyz[t,ids,:]-r0)**2 )
+        if rmsd > (0.1*rmsdThreshold)**2:
+            break
+    t0=t+1
+    print( "Averaging coordinates in %s after frame %d" % (dcd, t0) )
+
+    pos = sel.xyz[t0:,:,:]
+    pos = np.mean(pos, axis=0)
     return 10*pos               # convert to Angstroms
 
-
 def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEulerZ=0,
                    sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat'):
     # ----------------- #