diff --git a/bin/enrgmd b/bin/enrgmd
new file mode 100755
index 0000000000000000000000000000000000000000..6360704555092c60ec10f3e8526e5000cfc17b78
--- /dev/null
+++ b/bin/enrgmd
@@ -0,0 +1,80 @@
+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import argparse
+import re
+import pathlib, os
+import shutil
+from mrdna import get_resource_path
+
+parser = argparse.ArgumentParser(prog="mrdna",
+ description="""Easy multiresolution simulations of DNA nanotechnology objects using ARBD""")
+
+parser.add_argument('-o','--output-prefix', type=str, default=None,
+ help="Name for your job's output")
+parser.add_argument('-d','--directory', type=str, default=None,
+ help='Directory for simulation; does not need to exist yet')
+parser.add_argument('--backbone-scale', type=float, default=1.0,
+ help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
+
+parser.add_argument('input_file', type=str,
+ help="""Any of the following:
+ (1) a cadnano JSON file;
+ (2) a vHelix Maya (.ma) file;
+ (3) an atomic PDB file""")
+
+args = parser.parse_args()
+
+if __name__ == '__main__':
+
+ infile = pathlib.Path(args.input_file)
+ try:
+ prefix = infile.stem
+ extension = ".".join(infile.suffixes)
+ except:
+ raise Exception("Unrecognized input file '{}'".format(infile))
+
+ if extension == '.json':
+ from mrdna.readers import read_cadnano as read_model
+ elif extension == '.ma':
+ from mrdna.readers import read_vhelix as read_model
+ elif extension == '.pdb':
+ from mrdna.readers import read_atomic_pdb as read_model
+ else:
+ raise Exception("Unrecognized input file '{}'".format(infile))
+
+ model = read_model( str(infile) )
+
+ if args.output_prefix is not None:
+ prefix = args.output_prefix
+ directory = args.directory
+
+ d_orig = os.getcwd()
+ try:
+ if directory is None:
+ directory='.'
+ elif not os.path.exists(directory):
+ os.makedirs(directory)
+ os.chdir(directory)
+
+ print("""Generating your ENRG MD simulation files...
+If your simulation results in any publications, please cite the following:
+http://dx.doi.org/10.1093/nar/gkw155
+""")
+ model._clear_beads()
+ model._generate_atomic_model(scale=args.backbone_scale)
+ model.write_atomic_ENM( prefix )
+ model.atomic_simulate( output_name = prefix )
+ try:
+ shutil.copytree( get_resource_path("charmm36.nbfix"), "charmm36.nbfix" )
+ except FileExistsError:
+ ...
+ try:
+ os.makedirs('output')
+ except FileExistsError:
+ ...
+
+ except:
+ raise
+ finally:
+ os.chdir(d_orig)
diff --git a/mrdna/readers/segmentmodel_from_lists.py b/mrdna/readers/segmentmodel_from_lists.py
index c9e476bbd4cc42cad12ace3dd4c2ca673d7f51ac..a83645be41b99b2692cad1bae63218fe87e193ed 100644
--- a/mrdna/readers/segmentmodel_from_lists.py
+++ b/mrdna/readers/segmentmodel_from_lists.py
@@ -5,6 +5,7 @@ import numpy as np
import os,sys
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+from ..coords import quaternion_from_matrix
from .. import get_resource_path
def _three_prime_list_to_five_prime(three_prime):
@@ -81,26 +82,43 @@ def SegmentModelFromPdb(*args,**kwargs):
## Build map from residue index to helix index
-def set_splines(seg, coordinates, hid, hmap, hrank):
+def set_splines(seg, coordinates, hid, hmap, hrank, fwd, orientation=None):
maxrank = np.max( hrank[hmap==hid] )
if maxrank == 0:
ids = np.where((hmap == hid))[0]
pos = np.mean( [coordinates[r,:] for r in ids ], axis=0 )
coords = [pos,pos]
contours = [0,1]
+ if orientation is not None:
+ ids = np.where((hmap == hid) * fwd)[0]
+ assert( len(ids) == 1 )
+ q = quaternion_from_matrix( orientation[ids[0]] )
+ quats = [q, q]
+
else:
- coords,contours = [[],[]]
+ coords,contours,quats = [[],[],[]]
for rank in range(int(maxrank)+1):
ids = np.where((hmap == hid) * (hrank == rank))[0]
coords.append(np.mean( [coordinates[r,:] for r in ids ], axis=0 ))
contours.append( float(rank)/maxrank )
+ if orientation is not None:
+ ids = np.where((hmap == hid) * (hrank == rank) * fwd)[0]
+ assert(len(ids) == 1)
+ quats.append( quaternion_from_matrix( orientation[ids[0]] ) )
+ ## TODO: average quaterion with reverse direction
+
coords = np.array(coords)
seg.set_splines(contours,coords)
+ if orientation is not None:
+ quats = np.array(quats)
+ seg.set_orientation_splines(contours,quats)
+
seg.start_position = coords[0,:]
seg.end_position = coords[-1,:]
-def model_from_basepair_stack_3prime(coordinates, basepair, stack, three_prime, sequence=None):
+def model_from_basepair_stack_3prime(coordinates, basepair, stack, three_prime,
+ sequence=None, orientation=None):
"""
Creates a SegmentModel object from lists of each nucleotide's
basepair, its stack (on 3' side) and its 3'-connected nucleotide
@@ -180,7 +198,7 @@ def model_from_basepair_stack_3prime(coordinates, basepair, stack, three_prime,
for hid in double_stranded_helices:
seg = DoubleStrandedSegment(name=str(hid),
num_bp = np.sum(hmap==hid)//2)
- set_splines(seg, coordinates, hid, hmap, hrank)
+ set_splines(seg, coordinates, hid, hmap, hrank, fwd, orientation)
assert(hid == len(doubleSegments))
doubleSegments.append(seg)
@@ -190,7 +208,7 @@ def model_from_basepair_stack_3prime(coordinates, basepair, stack, three_prime,
for hid in single_stranded_helices:
seg = SingleStrandedSegment(name=str(hid),
num_nt = np.sum(hmap==hid))
- set_splines(seg, coordinates, hid, hmap, hrank)
+ set_splines(seg, coordinates, hid, hmap, hrank, fwd, orientation)
assert(hid == len(doubleSegments) + len(singleSegments))
singleSegments.append(seg)
diff --git a/mrdna/readers/segmentmodel_from_pdb.py b/mrdna/readers/segmentmodel_from_pdb.py
index 4bce635f313a232bcf7848de12cf0077dcf14d62..3cbb683de5a95be8f0fa76ea565e449b7cf12b76 100644
--- a/mrdna/readers/segmentmodel_from_pdb.py
+++ b/mrdna/readers/segmentmodel_from_pdb.py
@@ -12,6 +12,7 @@ import MDAnalysis as mda
import MDAnalysis.analysis.align as align
from MDAnalysis.analysis.distances import distance_array
+from ..model.CanonicalNucleotideAtoms import CanonicalNucleotideFactory
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from .. import get_resource_path
from .segmentmodel_from_lists import model_from_basepair_stack_3prime
@@ -284,12 +285,18 @@ def SegmentModelFromPdb(*args,**kwargs):
## Find sequence
seq = [resname_to_key[rn]
for rn in u.select_atoms("name C1'").resnames]
-
+
+ ## Get orientations
+ orientations = []
+ for R in transforms:
+ orientations.append( R.T.dot(CanonicalNucleotideFactory.DefaultOrientation) )
+
return model_from_basepair_stack_3prime(centers,
bps,
stacks,
three_prime,
seq,
+ orientations
)
if __name__ == "__main__":
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 47294aa71cb73192e795ba255988462320f69b07..16872d4a9448af9bbf32f97a4d38f042b74d6580 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -389,6 +389,10 @@ class Segment(ConnectableElement, Group):
tck, u = interpolate.splprep( coords.T, u=contours, s=0, k=1)
self.position_spline_params = tck
+ def set_orientation_splines(self, contours, quaternions):
+ tck, u = interpolate.splprep( quaternions.T, u=contours, s=0, k=1)
+ self.quaternion_spline_params = tck
+
def _set_splines_from_ends(self):
self.quaternion_spline_params = None
coords = np.array([self.start_position, self.end_position])
diff --git a/setup.py b/setup.py
index 4a43ce4dd38ff22e5f14b6b5b517a3d9e314aef7..e339ebf4072cdd6d140da637a1196f7b4ec65288 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,7 @@ setuptools.setup(
url="https://gitlab.engr.illinois.edu/tbgl/tools/mrdna",
packages=setuptools.find_packages(),
include_package_data=True,
- scripts=['bin/mrdna'],
+ scripts=['bin/mrdna','bin/enrgmd'],
install_requires=(
'numpy>=1.14',
'scipy>=1.1',