Commit b670c7f0 authored by cmaffeo2's avatar cmaffeo2
Browse files

Conditionally update coordinates using linear interpolation when there are not many beads

parent 9063e58f
...@@ -526,15 +526,15 @@ def run_simulation_protocol( output_name, job_id, json_data, ...@@ -526,15 +526,15 @@ def run_simulation_protocol( output_name, job_id, json_data,
sequence=None, sequence=None,
remove_long_bonds=False, remove_long_bonds=False,
gpu = 0, gpu = 0,
directory=None directory=None,
coarse_steps = 1e6,
fine_steps = 1e6
): ):
coarseSteps = 1e5+1
fineSteps = 1e5+1
ret = None ret = None
d_orig = os.getcwd() d_orig = os.getcwd()
try: try:
directory = "/dev/shm/origami-job-robot_v1.9_bent--0f4axf4k"
if directory is None: if directory is None:
import tempfile import tempfile
directory = tempfile.mkdtemp(prefix='origami-%s-' % job_id, dir='/dev/shm/') 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, ...@@ -583,18 +583,18 @@ def run_simulation_protocol( output_name, job_id, json_data,
output_prefix = "%s-0" % output_name output_prefix = "%s-0" % output_name
full_output_prefix = "%s/%s" % (output_directory,output_prefix) full_output_prefix = "%s/%s" % (output_directory,output_prefix)
## TODO Remove long bonds ## 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) coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
output_prefix = "%s-1" % output_name output_prefix = "%s-1" % output_name
model._update_segment_positions(coordinates) model._update_segment_positions(coordinates)
# model._clear_beads() # model._clear_beads()
# # model._generate_bead_model( 7, 4, False ) # # model._generate_bead_model( 7, 4, False )
# 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: else:
output_prefix = "%s-1" % output_name output_prefix = "%s-1" % output_name
full_output_prefix = "%s/%s" % (output_directory,output_prefix) 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) coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
...@@ -605,7 +605,7 @@ def run_simulation_protocol( output_name, job_id, json_data, ...@@ -605,7 +605,7 @@ def run_simulation_protocol( output_name, job_id, json_data,
model._update_segment_positions(coordinates) model._update_segment_positions(coordinates)
model._clear_beads() model._clear_beads()
model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=True ) 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) 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, ...@@ -615,7 +615,7 @@ def run_simulation_protocol( output_name, job_id, json_data,
model._update_segment_positions(coordinates) model._update_segment_positions(coordinates)
model._clear_beads() model._clear_beads()
model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=False ) 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 ) 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, ...@@ -625,7 +625,7 @@ def run_simulation_protocol( output_name, job_id, json_data,
model._update_segment_positions(coordinates) model._update_segment_positions(coordinates)
model._clear_beads() model._clear_beads()
model._generate_atomic_model(scale=0.25) model._generate_atomic_model(scale=0.25)
model.atomic_simulate( outputPrefix = full_output_prefix ) model.atomic_simulate( outputPrefix = output_prefix )
ret = directory ret = directory
......
...@@ -18,7 +18,9 @@ if __name__ == '__main__': ...@@ -18,7 +18,9 @@ if __name__ == '__main__':
except: except:
print("WARNING: skipping unreadable json file {}".format(f)) print("WARNING: skipping unreadable json file {}".format(f))
continue continue
try:
run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 ) run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
except: # try:
print("ERROR: failed to simulate {}".format(f)) # run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
# except:
# print("ERROR: failed to simulate {}".format(f))
...@@ -668,6 +668,9 @@ class Segment(ConnectableElement, Group): ...@@ -668,6 +668,9 @@ class Segment(ConnectableElement, Group):
# conn_locs = self.get_contour_sorted_connections_and_locations() # conn_locs = self.get_contour_sorted_connections_and_locations()
# locs = [A for c,A,B in conn_locs] # locs = [A for c,A,B in conn_locs]
# existing_beads = [l.particle for l in locs if l.particle is not None] # 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 = {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) ) existing_beads = sorted( list(existing_beads), key=lambda b: b.get_contour_position(self) )
...@@ -722,6 +725,8 @@ class Segment(ConnectableElement, Group): ...@@ -722,6 +725,8 @@ class Segment(ConnectableElement, Group):
last = eb1 last = eb1
for j in range(num_beads): for j in range(num_beads):
s = ds*(j+1) + s0 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) b = self._generate_one_bead(s,nts)
last.make_intrahelical_neighbor(b) last.make_intrahelical_neighbor(b)
...@@ -732,6 +737,8 @@ class Segment(ConnectableElement, Group): ...@@ -732,6 +737,8 @@ class Segment(ConnectableElement, Group):
if eb2.parent == self: if eb2.parent == self:
tmp_children.append(eb2) tmp_children.append(eb2)
# if self.name in ("51-2","51-3"):
# pdb.set_trace()
self._rebuild_children(tmp_children) self._rebuild_children(tmp_children)
def _regenerate_beads(self, max_nts_per_bead=4, ): def _regenerate_beads(self, max_nts_per_bead=4, ):
...@@ -1292,7 +1299,11 @@ class SegmentModel(ArbdModel): ...@@ -1292,7 +1299,11 @@ class SegmentModel(ArbdModel):
beads = list(beads) beads = list(beads)
## Skip beads that are None (some locations were not assigned a particle to avoid double-counting) ## 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 = [] contours = []
for b in beads: for b in beads:
try: try:
...@@ -1316,10 +1327,31 @@ class SegmentModel(ArbdModel): ...@@ -1316,10 +1327,31 @@ class SegmentModel(ArbdModel):
""" Get positions """ """ Get positions """
positions = bead_coordinates[ids,:].T positions = bead_coordinates[ids,:].T
try: ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 )
tck, u = interpolate.splprep( positions, u=contours, s=0, k=3 ) tck = ret[0][0]
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 )
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 s.position_spline_params = tck
...@@ -1404,6 +1436,8 @@ class SegmentModel(ArbdModel): ...@@ -1404,6 +1436,8 @@ class SegmentModel(ArbdModel):
## Loop through all connections, generating beads at appropriate locations ## Loop through all connections, generating beads at appropriate locations
for c,A,B in self.get_connections("intrahelical"): for c,A,B in self.get_connections("intrahelical"):
s1,s2 = [l.container for l in (A,B)] 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)) print("Working on {}".format(c))
## TODO be more elegant! ## TODO be more elegant!
# if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue # if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue
...@@ -1498,7 +1532,7 @@ class SegmentModel(ArbdModel): ...@@ -1498,7 +1532,7 @@ class SegmentModel(ArbdModel):
b1,b2 = [l.particle for l in (A,B)] b1,b2 = [l.particle for l in (A,B)]
if b1 is b2: if b1 is b2:
## already handled by Segment._generate_beads ## already handled by Segment._generate_beads
continue pass
else: else:
for b in (b1,b2): assert( b is not None ) for b in (b1,b2): assert( b is not None )
b1.make_intrahelical_neighbor(b2) b1.make_intrahelical_neighbor(b2)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment