diff --git a/cadnano_segments.py b/cadnano_segments.py
index d0adeea2b822735fc997b5081a0e0a1e74fce940..b8149c2ac45f13dddbade565d95154332d1b6cd9 100644
--- a/cadnano_segments.py
+++ b/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()
diff --git a/run.py b/run.py
index 71f2f4f118a0a15fe6a182b9a9c3022ce51dd7ab..fd24ebc674b58885ca43087ae4d23536765312c0 100644
--- a/run.py
+++ b/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
diff --git a/simulate.py b/simulate.py
new file mode 100644
index 0000000000000000000000000000000000000000..933db19e0cd13cdc92e25c5841df6862a9e01a3a
--- /dev/null
+++ b/simulate.py
@@ -0,0 +1,89 @@
+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