Removed mdtraj dependency

dnarbd/coords.py

@@ -164,63 +164,29 @@ def readArbdCoords(fname):
     return np.array(coords)
 
 def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5):
-    import mdtraj as md
-    ## align trajectory to pdb
-    ref = md.load(pdb, top=pdb)
-    sel = md.load(dcd, top=pdb)
-    # return ref.xyz[0,:,:]*10
+    import MDAnalysis as mda
 
-    # ref = sel[-1]
-    ids = ref.topology.select("name =~ 'D.*'")
-    assert(len(ids) > 3)
+    usel = mda.Universe(psf, dcd)
+    sel = usel.select_atoms("name D*")
 
     # r0 = ref.xyz[0,ids,:]
-    r0 = sel.xyz[-1,ids,:]
-    t = -1                      # in case dcd_frames < 3
-    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:
+    ts = usel.trajectory[-1]
+    r0 = sel.positions
+    pos = []
+    for t in range(ts.frame-1,-1,-1):
+        usel.trajectory[t]
+        R,comA,comB = minimizeRmsd(sel.positions,r0)
+        r = np.array( [(r-comA).dot(R)+comB for r in sel.positions] )
+        rmsd = np.mean( (r-r0)**2 )
+        r = np.array( [(r-comA).dot(R)+comB for r in usel.atoms.positions] )
+        pos.append( r )
+        if rmsd > 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 readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5):
-#     import MDAnalysis as md
-#     ## align trajectory to pdb
-#     from pdb import set_trace
-#     set_trace()
-#     uref = md.Universe(psf, pdb)
-#     usel = md.Universe(psf, dcd)
-
-#     ref = uref.select_atoms("name 'd*'")
-#     sel = usel.select_atoms("name 'd*'")
-
-#     # r0 = ref.xyz[0,ids,:]
-#     ts = usel.trajectory[-1]
-#     ts.frame
-#     r0 = sel.positions
-#     for t in range(ts.frame-1,-1,-1):
-#         print(t)
-#         usel.trajectory[t]
-#         R,comA,comB = minimizeRmsd(sel.positions,r0)
-#         sel.positions = np.array( [(r-comA).dot(R)+comB for r in sel.positions] )
-#         rmsd = np.mean( (sel.positions-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
+    return pos
 
 def unit_quat_conversions():
     for axis in [[0,0,1],[1,1,1],[1,0,0],[-1,-2,0]]:

dnarbd/segmentmodel.py

@@ -744,7 +744,7 @@ class Segment(ConnectableElement, Group):
             b = self._generate_one_bead( self.nt_pos_to_contour(0), 0)
             existing_beads = [b] + existing_beads
 
-        if existing_beads[-1].get_nt_position(self)-(self.num_nt-1) < -0.5:
+        if existing_beads[-1].get_nt_position(self)-(self.num_nt-1) < -0.5 or len(existing_beads)==1:
             b = self._generate_one_bead( self.nt_pos_to_contour(self.num_nt-1), 0)
             existing_beads.append(b)
         assert(len(existing_beads) > 1)

dnarbd/simulate.py

@@ -76,7 +76,7 @@ def multiresolution_simulation( model, output_name,
         model._clear_beads()
         model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=False )
         model.simulate( outputPrefix = output_prefix, numSteps=fine_steps, **simargs )
-        coordinates = readAvgArbdCoords('%s.psf' % output_name,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix )
+        coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix )
 
 
         """ Atomic simulation """