From 09b759a351fa5930f5bae87efe63e234ccdce53e Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Wed, 5 Sep 2018 18:21:18 +0000
Subject: [PATCH] Removed mdtraj dependency
---
dnarbd/coords.py | 64 ++++++++++--------------------------------
dnarbd/segmentmodel.py | 2 +-
dnarbd/simulate.py | 2 +-
3 files changed, 17 insertions(+), 51 deletions(-)
diff --git a/dnarbd/coords.py b/dnarbd/coords.py
index 70b6216..b74906d 100644
--- a/dnarbd/coords.py
+++ b/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]]:
diff --git a/dnarbd/segmentmodel.py b/dnarbd/segmentmodel.py
index 314bcab..eabccb5 100644
--- a/dnarbd/segmentmodel.py
+++ b/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)
diff --git a/dnarbd/simulate.py b/dnarbd/simulate.py
index 011fb7a..8d9cc11 100644
--- a/dnarbd/simulate.py
+++ b/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 """
--
GitLab