diff --git a/cadnano_segments.py b/cadnano_segments.py index 32c2104ee2ddecd36b9c7d4002eb4d0746c6dd4c..626ceea25462861fcf827261fc693adb4c7a418c 100644 --- a/cadnano_segments.py +++ b/cadnano_segments.py @@ -526,15 +526,15 @@ def run_simulation_protocol( output_name, job_id, json_data, sequence=None, remove_long_bonds=False, gpu = 0, - directory=None + directory=None, + coarse_steps = 1e6, + fine_steps = 1e6 ): - coarseSteps = 1e5+1 - fineSteps = 1e5+1 - ret = None d_orig = os.getcwd() try: + directory = "/dev/shm/origami-job-robot_v1.9_bent--0f4axf4k" if directory is None: import tempfile directory = tempfile.mkdtemp(prefix='origami-%s-' % job_id, dir='/dev/shm/') @@ -583,18 +583,18 @@ def run_simulation_protocol( output_name, job_id, json_data, 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*coarseSteps, **simargs ) + 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*coarseSteps, **simargs ) + 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=coarseSteps, **simargs ) + # model.simulate( outputPrefix = output_prefix, numSteps=coarse_steps, **simargs ) coordinates = readArbdCoords('%s.0.restart' % full_output_prefix) @@ -605,7 +605,7 @@ def run_simulation_protocol( output_name, job_id, json_data, 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=fineSteps, **simargs ) + # 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) @@ -615,7 +615,7 @@ def run_simulation_protocol( output_name, job_id, json_data, 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=fineSteps, **simargs ) + 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 ) @@ -625,7 +625,7 @@ def run_simulation_protocol( output_name, job_id, json_data, model._update_segment_positions(coordinates) model._clear_beads() model._generate_atomic_model(scale=0.25) - model.atomic_simulate( outputPrefix = full_output_prefix ) + model.atomic_simulate( outputPrefix = output_prefix ) ret = directory diff --git a/run.py b/run.py index 03b1bb6fd10aa5f2fe09775db60825bcfcafe1e8..499f49c5a4bf0be1c78ce66af4602363f1e387b6 100644 --- a/run.py +++ b/run.py @@ -18,7 +18,9 @@ if __name__ == '__main__': except: print("WARNING: skipping unreadable json file {}".format(f)) continue - try: - run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 ) - except: - print("ERROR: failed to simulate {}".format(f)) + + run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 ) + # try: + # run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 ) + # except: + # print("ERROR: failed to simulate {}".format(f)) diff --git a/segmentmodel.py b/segmentmodel.py index 5875ecbb02ec049acbe45bdac3ccbfa4ea40b7f0..8ba44617991778594b81138b8b2e678b8d461fa6 100644 --- a/segmentmodel.py +++ b/segmentmodel.py @@ -668,6 +668,9 @@ class Segment(ConnectableElement, Group): # conn_locs = self.get_contour_sorted_connections_and_locations() # locs = [A for c,A,B in conn_locs] # existing_beads = [l.particle for l in locs if l.particle is not None] + # if self.name == "51-2": + # pdb.set_trace() + existing_beads = {l.particle for l in self.locations if l.particle is not None} existing_beads = sorted( list(existing_beads), key=lambda b: b.get_contour_position(self) ) @@ -722,6 +725,8 @@ class Segment(ConnectableElement, Group): last = eb1 for j in range(num_beads): s = ds*(j+1) + s0 + # if self.name in ("51-2","51-3"): + # print(" adding bead at {}".format(s)) b = self._generate_one_bead(s,nts) last.make_intrahelical_neighbor(b) @@ -732,6 +737,8 @@ class Segment(ConnectableElement, Group): if eb2.parent == self: tmp_children.append(eb2) + # if self.name in ("51-2","51-3"): + # pdb.set_trace() self._rebuild_children(tmp_children) def _regenerate_beads(self, max_nts_per_bead=4, ): @@ -1292,7 +1299,11 @@ class SegmentModel(ArbdModel): beads = list(beads) ## Skip beads that are None (some locations were not assigned a particle to avoid double-counting) - beads = [b for b in beads if b is not None] + # beads = [b for b in beads if b is not None] + beads = list(filter(lambda x: x is not None, beads)) + if isinstance(self, DoubleStrandedSegment): + beads = list(filter(lambda x: x.type_.name[0] == "D", beads)) + contours = [] for b in beads: try: @@ -1316,10 +1327,31 @@ class SegmentModel(ArbdModel): """ Get positions """ positions = bead_coordinates[ids,:].T - try: - tck, u = interpolate.splprep( positions, u=contours, s=0, k=3 ) - except: - tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 ) + ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) + tck = ret[0][0] + + ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) + tck2 = ret[0][0] + + assert( tck[0][0] == tck2[0][0] ) + + # if len(beads) < 8: + # ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) + # tck = ret[0][0] + # if ret[2] > 0: + # pdb.set_trace() + + # else: + # try: + # ret = interpolate.splprep( positions, u=contours, s=0, k=3, full_output=1 ) + # tck = ret[0][0] + # if ret[2] > 0: + # pdb.set_trace() + # except: + # ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) + # tck = ret[0][0] + # if ret[2] > 0: + # pdb.set_trace() s.position_spline_params = tck @@ -1404,6 +1436,8 @@ class SegmentModel(ArbdModel): ## Loop through all connections, generating beads at appropriate locations for c,A,B in self.get_connections("intrahelical"): s1,s2 = [l.container for l in (A,B)] + # if s1.name in ("51-2","51-3") and s2.name in ("51-2","51-3"): + # pdb.set_trace() print("Working on {}".format(c)) ## TODO be more elegant! # if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue @@ -1498,7 +1532,7 @@ class SegmentModel(ArbdModel): b1,b2 = [l.particle for l in (A,B)] if b1 is b2: ## already handled by Segment._generate_beads - continue + pass else: for b in (b1,b2): assert( b is not None ) b1.make_intrahelical_neighbor(b2)