Moved run_simulate_protocols outside of cadnano_segments since more modular;...

Moved run_simulate_protocols outside of cadnano_segments since more modular; increased density of coarsest DNA model

cadnano_segments.py

@@ -449,37 +449,29 @@ class cadnano_part(SegmentModel):
         for hid,segs in self.helices.items():
             occ = self.strand_occupancies[hid]
             for i in range(len(segs)-1):
-                    seg1,seg2 = [segs[j] for j in (i,i+1)]
-                    r1,r2 = [self.helix_ranges[hid][j] for j in (i,i+1)]
-                    if r1[1]+1 == r2[0]:
-                        ## TODO: handle nicks that are at intrahelical connections(?)
-                        zmid1 = int(0.5*(r1[0]+r1[1]))
-                        zmid2 = int(0.5*(r2[0]+r2[1]))
-                        # if seg1.name == "1-1" or seg2.name == "1-1":
-                        #     # if seg1.name == "19-3" or seg2.name == "19-3":
-                        #     import pdb
-                        #     pdb.set_trace()
-
-                        ## TODO: validate
-                        if zmid1 in occ[0] and zmid2 in occ[0]:
-                            seg1.connect_end3(seg2.start5)
-
-                        if zmid1 in occ[1] and zmid2 in occ[1]:
-                            if zmid1 in occ[0]:
-                                end = seg1.end5
-                            else:
-                                end = seg1.start5
-                            if zmid2 in occ[0]:
-                                seg2.connect_start3(end)
-                            else:
-                                seg2.connect_end3(end)
+                seg1,seg2 = [segs[j] for j in (i,i+1)]
+                r1,r2 = [self.helix_ranges[hid][j] for j in (i,i+1)]
+                if r1[1]+1 == r2[0]:
+                    ## TODO: handle nicks that are at intrahelical connections(?)
+                    zmid1 = int(0.5*(r1[0]+r1[1]))
+                    zmid2 = int(0.5*(r2[0]+r2[1]))
+
+                    ## TODO: validate
+                    if zmid1 in occ[0] and zmid2 in occ[0]:
+                        seg1.connect_end3(seg2.start5)
+
+                    if zmid1 in occ[1] and zmid2 in occ[1]:
+                        if zmid1 in occ[0]:
+                            end = seg1.end5
+                        else:
+                            end = seg1.start5
+                        if zmid2 in occ[0]:
+                            seg2.connect_start3(end)
+                        else:
+                            seg2.connect_end3(end)
                         
     def _add_crossovers(self):
         for h1,f1,z1,h2,f2,z2 in self.xover_list:
-            # if (h1 == 52 or h2 == 52) and z1 == 221:
-            if (h1 == 15 and z1 == 230) or h2 == 15 and z2 == 231:
-                import pdb
-                pdb.set_trace()
             seg1, nt1 = self._get_segment_nucleotide(h1,z1,not f1)
             seg2, nt2 = self._get_segment_nucleotide(h2,z2,f2)
             ## TODO: use different types of crossovers
@@ -548,8 +540,8 @@ def read_model(json_data, sequence=None):
     """ Read in data """
     part = decode_cadnano_part(json_data)
     model = cadnano_part(part,
-                         max_basepairs_per_bead = 7,
-                         max_nucleotides_per_bead = 4,
+                         max_basepairs_per_bead = 5,
+                         max_nucleotides_per_bead = 5,
                          local_twist=False)
     model._generate_strands() # TODO: move into model creation
 
@@ -576,120 +568,6 @@ def read_model(json_data, sequence=None):
 
     return model
 
-
-def run_simulation_protocol( output_name, job_id, json_data,
-                             sequence=None,
-                             remove_long_bonds=False,
-                             gpu = 0,
-                             directory=None,
-                             coarse_steps = 1e7+1,
-                             fine_steps = 1e7+1
-                         ):
-
-    ret = None
-    d_orig = os.getcwd()
-    try:
-        if directory is None:
-            import tempfile
-            directory = tempfile.mkdtemp(prefix='origami-%s-' % job_id, dir='/dev/shm/')
-        elif not os.path.exists(directory):
-            os.makedirs(directory)
-        os.chdir(directory)
-
-        output_directory = "output"
-        
-        """ Read in data """
-        part = decode_cadnano_part(json_data)
-        model = cadnano_part(part,
-                             max_basepairs_per_bead = 7, 
-                             max_nucleotides_per_bead = 4,
-                             local_twist=False)
-        model._generate_strands() # TODO: move into model creation
-
-        # TODO        
-        # try:
-        #     model.set_cadnano_sequence()
-        # finally:
-        #     ...
-        #     if sequence is not None and len() :
-        #         model.strands[0].set_sequence(seq)
-        
-        if sequence is None or len(sequence) == 0:
-            ## default m13mp18
-            sequence = list(m13seq)
-            try:
-                model.strands[0].set_sequence(sequence)
-            except:
-                ...
-        else:
-            model.strands[0].set_sequence(sequence)
-
-        for s in model.segments:
-            s.randomize_unset_sequence()
-
-        # TODO: add the output directory in to the readArbdCoords functions or make it an attribute of the model object 
-
-        """ Coarse simulation """
-        # simargs = dict(timestep=200e-6, outputPeriod=1e5, gpu=gpu )
-        simargs = dict(timestep=200e-6, outputPeriod=1e5, gpu=gpu, arbd=arbd )
-        
-        if remove_long_bonds:
-            output_prefix = "%s-0" % output_name
-            full_output_prefix = "%s/%s" % (output_directory,output_prefix)
-            ## TODO Remove long bonds
-            model.simulate( outputPrefix = output_prefix, numSteps=0.05*coarse_steps, **simargs )
-            coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
-            output_prefix = "%s-1" % output_name
-            model._update_segment_positions(coordinates)
-            # model._clear_beads()
-            # # model._generate_bead_model( 7, 4, False )
-            # model._generate_bead_model( 7, 4, False )
-            model.simulate( outputPrefix = output_prefix, numSteps=0.95*coarse_steps, **simargs )
-        else:
-            output_prefix = "%s-1" % output_name
-            full_output_prefix = "%s/%s" % (output_directory,output_prefix)
-            model.simulate( outputPrefix = output_prefix, numSteps=coarse_steps, **simargs )
-        coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
-
-
-        """ Fine simulation """ 
-        output_prefix = "%s-2" % output_name
-        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
-        simargs['timestep'] = 50e-6
-        model._update_segment_positions(coordinates)
-        model._clear_beads()
-        model._generate_bead_model( 1.03, 1.03, local_twist=True, escapable_twist=True )
-        model.simulate( outputPrefix = output_prefix, numSteps=fine_steps, **simargs )
-        coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
-
-
-        """ Freeze twist """
-        output_prefix = "%s-3" % output_name
-        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
-        model._update_segment_positions(coordinates)
-        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 )
-
-
-        """ Atomic simulation """
-        output_prefix = "%s-4" % output_name
-        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
-        model._update_segment_positions(coordinates)
-        model._clear_beads()
-        model._generate_atomic_model(scale=0.25)
-        model.atomic_simulate( outputPrefix = output_prefix )
-
-        ret = directory
-        
-    except:
-        raise
-    finally:
-        os.chdir(d_orig)
-        
-    return ret
-
 # pynvml.nvmlInit()
 # gpus = range(pynvml.nvmlDeviceGetCount())
 # pynvml.nvmlShutdown()

run.py

@@ -2,7 +2,8 @@
 import sys
 from glob import glob
 import re
-from cadnano_segments import read_json_file, run_simulation_protocol
+from cadnano_segments import read_json_file, read_model
+from simulate import run_simulation_protocol
 
 if __name__ == '__main__':
     if len(sys.argv) > 1:
@@ -15,12 +16,18 @@ if __name__ == '__main__':
         out = re.match('json/(.*).json',f).group(1)
         try:
             data = read_json_file(f)
+            model = read_model(data)
         except:
             print("WARNING: skipping unreadable json file {}".format(f))
             continue
             
-        # run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
-        try:
-            run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
-        except:
-            print("WARNING: failed to simulate {}".format(f))
+        gpu = 0
+        directory = "run.out.7bd92d9.twist-angles"
+        run_simulation_protocol( out, "job-"+out+"-", model, gpu=gpu,
+                                 directory = directory)
+        # try:
+        #     run_simulation_protocol( out, "job-"+out+"-", model, gpu=gpu,
+        #                              directory = directory)
+        # except:
+        #     print("WARNING: failed to simulate json file {}".format(f))
+        #     continue

simulate.py

+import os
+from coords import readArbdCoords, readAvgArbdCoords
+
+arbd="/home/cmaffeo2/development/cuda/arbd.dbg/src/arbd" # reduced the mem footprint cause vmd
+namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multicore-CUDA/namd2"
+
+def run_simulation_protocol( output_name, job_id, model,
+                             remove_long_bonds=False,
+                             gpu = 0,
+                             directory=None,
+                             coarse_steps = 5e7+1,
+                             fine_steps = 5e7+1
+                         ):
+
+    ret = None
+    d_orig = os.getcwd()
+    try:
+        if directory is None:
+            import tempfile
+            directory = tempfile.mkdtemp(prefix='origami-%s-' % job_id, dir='/dev/shm/')
+        elif not os.path.exists(directory):
+            os.makedirs(directory)
+        os.chdir(directory)
+
+        output_directory = "output"
+        
+        # TODO: add the output directory in to the readArbdCoords functions or make it an attribute of the model object 
+
+        """ Coarse simulation """
+        # simargs = dict(timestep=200e-6, outputPeriod=1e5, gpu=gpu )
+        simargs = dict(timestep=200e-6, outputPeriod=1e5, gpu=gpu, arbd=arbd )
+
+        model._clear_beads()
+        model._generate_bead_model( 5, 5, local_twist=False )
+
+        if remove_long_bonds:
+            output_prefix = "%s-0" % output_name
+            full_output_prefix = "%s/%s" % (output_directory,output_prefix)
+            ## TODO Remove long bonds
+            model.simulate( outputPrefix = output_prefix, numSteps=0.05*coarse_steps, **simargs )
+            coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
+            output_prefix = "%s-1" % output_name
+            model._update_segment_positions(coordinates)
+            model.simulate( outputPrefix = output_prefix, numSteps=0.95*coarse_steps, **simargs )
+        else:
+            output_prefix = "%s-1" % output_name
+            full_output_prefix = "%s/%s" % (output_directory,output_prefix)
+            model.simulate( outputPrefix = output_prefix, numSteps=coarse_steps, **simargs )
+        coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
+
+
+        """ Fine simulation """ 
+        output_prefix = "%s-2" % output_name
+        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
+        simargs['timestep'] = 50e-6
+        model._update_segment_positions(coordinates)
+        model._clear_beads()
+        model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
+        model.simulate( outputPrefix = output_prefix, numSteps=fine_steps, **simargs )
+        coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
+
+
+        """ Freeze twist """
+        output_prefix = "%s-3" % output_name
+        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
+        model._update_segment_positions(coordinates)
+        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 )
+
+
+        """ Atomic simulation """
+        output_prefix = "%s-4" % output_name
+        full_output_prefix = "%s/%s" % (output_directory,output_prefix)
+        model._update_segment_positions(coordinates)
+        model._clear_beads()
+        # model._generate_atomic_model(scale=0.25)
+        model._generate_atomic_model(scale=1.0)
+        model.atomic_simulate( outputPrefix = output_prefix )
+
+        ret = directory
+        
+    except:
+        raise
+    finally:
+        os.chdir(d_orig)
+        
+    return ret