diff --git a/mrdna/coords.py b/mrdna/coords.py
deleted file mode 100644
index 41628e83501c17098094fc12b7d4d4e6594ebcb6..0000000000000000000000000000000000000000
--- a/mrdna/coords.py
+++ /dev/null
@@ -1,204 +0,0 @@
-import numpy as np
-from scipy.optimize import newton
-
-def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
- ## Going through many iterations wasn't really needed
- tol = 1
- count = 0
-
- R = np.eye(3)
- comB = np.zeros([3,])
- cNext = coordsB
-
- while tol > 1e-6:
- q,cB,comA = _minimizeRmsd(cNext,coordsA, weights)
- R = R.dot(quaternion_to_matrix(q))
- assert( np.all(np.isreal( R )) )
-
- comB += cB
- cLast = cNext
- cNext = (coordsB-comB).dot(R)
-
- tol = np.sum(((cNext-cLast)**2)[:]) / np.max(np.shape(coordsB))
- if count > maxIter:
- Exception("Exceeded maxIter (%d)" % maxIter)
- count += 1
-
- print("%d iterations",count)
- return R, comB, comA
-
-
-def minimizeRmsd(coordsB, coordsA, weights=None):
- q,comA,comB = _minimizeRmsd(coordsB, coordsA, weights)
- assert( np.all(np.isreal( q )) )
- return quaternion_to_matrix(q),comA,comB
-
-
-## http://onlinelibrary.wiley.com/doi/10.1002/jcc.21439/full
-def _minimizeRmsd(coordsB, coordsA, weights=None):
- A = coordsA
- B = coordsB
-
- shapeA,shapeB = [np.shape(X) for X in (A,B)]
- for s in (shapeA,shapeB): assert( len(s) == 2 )
-
- A,B = [X.T if s[1] > s[0] else X for X,s in zip([A,B],(shapeA,shapeB))] # TODO: print warning
-
- shapeA,shapeB = [np.shape(X) for X in (A,B)]
- assert( shapeA == shapeB )
- for X,s in zip((A,B),(shapeA,shapeB)):
- assert( s[1] == 3 and s[0] >= s[1] )
-
- # if weights is None: weights = np.ones(len(A))
- if weights is None:
- comA,comB = [np.mean( X, axis=0 ) for X in (A,B)]
- else:
- assert( len(weights[:]) == len(B) )
- W = np.diag(weights)
- comA,comB = [np.sum( W.dot(X), axis=0 ) / np.sum(W) for X in (A,B)]
-
- A = np.array( A-comA )
- B = np.array( B-comB )
-
- if weights is None:
- s = A.T.dot(B)
- else:
- s = A.T.dot(W.dot(B))
-
- sxx,sxy,sxz = s[0,:]
- syx,syy,syz = s[1,:]
- szx,szy,szz = s[2,:]
-
- K = [[ sxx+syy+szz, syz-szy, szx-sxz, sxy-syx],
- [syz-szy, sxx-syy-szz, sxy+syx, sxz+szx],
- [szx-sxz, sxy+syx, -sxx+syy-szz, syz+szy],
- [sxy-syx, sxz+szx, syz+szy, -sxx-syy+szz]]
- K = np.array(K)
-
- # GA = np.trace( A.T.dot(W.dot(A)) )
- # GB = np.trace( B.T.dot(W.dot(B)) )
-
- ## Finding GA/GB can be done more quickly
- # I = np.eye(4)
- # x0 = (GA+GB)*0.5
- # vals = newtoon(lambda x: np.det(K-x*I), x0 = x0)
-
- vals, vecs = np.linalg.eig(K)
- i = np.argmax(vals)
- q = vecs[:,i]
-
- # RMSD = np.sqrt( (GA+GB-2*vals[i]) / len(A) )
- # print("CHECK:", K.dot(q)-vals[i]*q )
- return q, comB, comA
-
-def quaternion_to_matrix(q):
- assert(len(q) == 4)
-
- ## It looks like the wikipedia article I used employed a less common convention for q (see below
- ## http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_.E2.86.94_quaternion
- # q1,q2,q3,q4 = q
- # R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q4), 2*(q1*q3 + q2*q4)],
- # [ 2*(q1*q2 + q3*q4), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q4)],
- # [ 2*(q1*q3 - q2*q4), 2*(q1*q4 + q2*q3), 1-2*(q2*q2 + q1*q1)]]
-
- q = q / np.linalg.norm(q)
- q0,q1,q2,q3 = q
- R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q0), 2*(q1*q3 + q2*q0)],
- [ 2*(q1*q2 + q3*q0), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q0)],
- [ 2*(q1*q3 - q2*q0), 2*(q1*q0 + q2*q3), 1-2*(q2*q2 + q1*q1)]]
-
- return np.array(R)
-
-def quaternion_from_matrix( R ):
- e1 = R[0]
- e2 = R[1]
- e3 = R[2]
-
- # d1 = 0.5 * np.sqrt( 1+R[0,0]+R[1,1]+R[2,2] )
- # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
- # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
- # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
-
- d1 = 1+R[0,0]+R[1,1]+R[2,2]
- d2 = 1+R[0,0]-R[1,1]-R[2,2]
- d3 = 1-R[0,0]+R[1,1]-R[2,2]
- d4 = 1-R[0,0]-R[1,1]+R[2,2]
-
- maxD = max((d1,d2,d3,d4))
- d = 0.5 / np.sqrt(maxD)
-
- if d1 == maxD:
- return np.array(( 1.0/(4*d),
- d * (R[2,1]-R[1,2]),
- d * (R[0,2]-R[2,0]),
- d * (R[1,0]-R[0,1]) ))
- elif d2 == maxD:
- return np.array(( d * (R[2,1]-R[1,2]),
- 1.0/(4*d),
- d * (R[0,1]+R[1,0]),
- d * (R[0,2]+R[2,0]) ))
- elif d3 == maxD:
- return np.array(( d * (R[0,2]-R[2,0]),
- d * (R[0,1]+R[1,0]),
- 1.0/(4*d),
- d * (R[1,2]+R[2,1]) ))
- elif d4 == maxD:
- return np.array(( d * (R[1,0]-R[0,1]),
- d * (R[0,2]+R[2,0]),
- d * (R[1,2]+R[2,1]),
- 1.0/(4*d) ))
-
-def rotationAboutAxis(axis,angle, normalizeAxis=True):
- if normalizeAxis: axis = axis / np.linalg.norm(axis)
- angle = angle * 0.5 * np.pi/180
- cos = np.cos( angle )
- sin = np.sin( angle )
- q = [cos] + [sin*x for x in axis]
- return quaternion_to_matrix(q)
-
-def readArbdCoords(fname):
- coords = []
- with open(fname) as fh:
- for line in fh:
- coords.append([float(x) for x in line.split()[1:]])
- return np.array(coords)
-
-def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5):
- import MDAnalysis as mda
-
- usel = mda.Universe(psf, dcd)
- sel = usel.select_atoms("name D*")
-
- # r0 = ref.xyz[0,ids,:]
- 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 = np.mean(pos, axis=0)
- return pos
-
-def unit_quat_conversions():
- for axis in [[0,0,1],[1,1,1],[1,0,0],[-1,-2,0]]:
- for angle in np.linspace(-180,180,10):
- R = rotationAboutAxis(axis, angle)
- R2 = quaternion_to_matrix( quaternion_from_matrix( R ) )
- if not np.all( np.abs(R-R2) < 0.01 ):
- import pdb
- pdb.set_trace()
- quaternion_to_matrix( quaternion_from_matrix( R ) )
-
-
-if __name__ == "__main__":
- unit_quat_conversions()
-
diff --git a/mrdna/model/CanonicalNucleotideAtoms.py b/mrdna/model/CanonicalNucleotideAtoms.py
index fa0a2e98ca6bbe1ca19c1f84acd3858216058e52..7b2dcdce5f172e80aad27eeaa926ed9e17bf51ff 100644
--- a/mrdna/model/CanonicalNucleotideAtoms.py
+++ b/mrdna/model/CanonicalNucleotideAtoms.py
@@ -1,8 +1,8 @@
import json
import numpy as np
import copy
-from ..coords import rotationAboutAxis
-from .arbdmodel import Group, PointParticle, ParticleType
+from ..arbdmodel.coords import rotationAboutAxis
+from ..arbdmodel import Group, PointParticle, ParticleType
from .. import get_resource_path
seqComplement = dict(A='T',G='C')
diff --git a/mrdna/model/nbPot.py b/mrdna/model/nbPot.py
index c26c3b07f116bc67121642c0a7dbf5c1d28df374..2329bc31d322c3a1135a7a8b368146dcda99e899 100644
--- a/mrdna/model/nbPot.py
+++ b/mrdna/model/nbPot.py
@@ -3,9 +3,9 @@ import scipy.optimize as opt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter as savgol
-from mrdna.model.nonbonded import NonbondedScheme
-from mrdna import get_resource_path
-from mrdna.config import CACHE_DIR
+from ..arbdmodel.nonbonded import NonbondedScheme
+from .. import get_resource_path
+from ..config import CACHE_DIR
from scipy.optimize import curve_fit
from pathlib import Path
diff --git a/mrdna/readers/cadnano_segments.py b/mrdna/readers/cadnano_segments.py
index 304171843f0f975cdec8e492c8ac67b3652696a8..f5986074e8c633860c1dd6f041c26c41747c9246 100644
--- a/mrdna/readers/cadnano_segments.py
+++ b/mrdna/readers/cadnano_segments.py
@@ -5,7 +5,7 @@ import os,sys
from glob import glob
import re
-from ..coords import readArbdCoords, readAvgArbdCoords, rotationAboutAxis
+from ..arbdmodel.coords import readArbdCoords, readAvgArbdCoords, rotationAboutAxis
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from ..model.dna_sequence import m13 as m13seq
diff --git a/mrdna/readers/polygon_mesh.py b/mrdna/readers/polygon_mesh.py
index ec41a6774f6f5ef80373235acc999002bd4bd103..15ba8a26aa81b3430e85090de24425416a963b5e 100644
--- a/mrdna/readers/polygon_mesh.py
+++ b/mrdna/readers/polygon_mesh.py
@@ -1,7 +1,7 @@
import numpy as np
import sys
import re
-from ..coords import rotationAboutAxis, minimizeRmsd
+from ..arbdmodel.coords import rotationAboutAxis, minimizeRmsd
from ..version import maintainer
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
@@ -310,8 +310,8 @@ def write_pdb_psf(bases, prefix):
""" Function for debugging .ma parser """
- from ..model.arbdmodel import ArbdModel,ParticleType, PointParticle, Group
- from ..model.nonbonded import HarmonicBond
+ from ..arbdmodel import ArbdModel,ParticleType, PointParticle, Group
+ from ..arbdmodel.nonbonded import HarmonicBond
types = {c:ParticleType(c.upper(), mass=1,radius=1) for c in "bfuxyzs"}
bond=HarmonicBond(k=0,r0=1)
diff --git a/mrdna/readers/segmentmodel_from_lists.py b/mrdna/readers/segmentmodel_from_lists.py
index 277573d351e2eb327f79fd1f3413e1e5a8ed6a3c..475988982f5368eca84b9b03e63ca34e5ef8a2c3 100644
--- a/mrdna/readers/segmentmodel_from_lists.py
+++ b/mrdna/readers/segmentmodel_from_lists.py
@@ -6,7 +6,7 @@ import os,sys
import scipy
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
-from ..coords import quaternion_from_matrix
+from ..arbdmodel.coords import quaternion_from_matrix
from .. import get_resource_path
ref_stack_position = np.array((-2.41851735, -0.259761333, 3.39999978))
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 6f5efa809c6b7b118f87eb0ffbbe8732823200d0..a08670e324f1a3c8fc59c472da56deed41d64d21 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -2,9 +2,9 @@ import pdb
from pathlib import Path
import numpy as np
import random
-from .model.arbdmodel import PointParticle, ParticleType, Group, ArbdModel
-from .coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
-from .model.nonbonded import *
+from .arbdmodel import PointParticle, ParticleType, Group, ArbdModel
+from .arbdmodel.coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
+from .arbdmodel.nonbonded import *
from copy import copy, deepcopy
from .model.nbPot import nbDnaScheme
@@ -1545,14 +1545,16 @@ class SegmentModel(ArbdModel):
timestep=50e-6, cutoff=50,
decompPeriod=10000, pairlistDistance=None,
nonbondedResolution=0, DEBUG=0,
+ integrator='Brown',
debye_length = None,
extra_bd_file_lines = "",
):
self.DEBUG = DEBUG
if DEBUG > 0: print("Building ARBD Model")
ArbdModel.__init__(self,segments,
- dimensions, temperature, timestep, cutoff,
- decompPeriod, pairlistDistance=None,
+ dimensions, temperature,
+ timestep, integrator,
+ cutoff, decompPeriod, pairlistDistance=None,
nonbondedResolution = 0,
extra_bd_file_lines = extra_bd_file_lines
)
diff --git a/mrdna/simulate.py b/mrdna/simulate.py
index 348c452e901d0155354857195b25340f8dc40d1a..f0528d2b2c3fcaef07e7d243b08fa9e63be9e564 100644
--- a/mrdna/simulate.py
+++ b/mrdna/simulate.py
@@ -1,6 +1,6 @@
import os
import tempfile
-from .coords import readArbdCoords, readAvgArbdCoords
+from .arbdmodel.coords import readArbdCoords, readAvgArbdCoords
import shutil
from . import get_resource_path