Added support for reading in orientation; added enrgmd script for output...

Added support for reading in orientation; added enrgmd script for output enrgmd files without running simulations

bin/enrgmd

+#! /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)

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)

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__":

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])

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',