Commit 2175845a authored by cmaffeo2's avatar cmaffeo2
Browse files

tmp commit of older work; changed how duplicate bonds are handled

parent 81546af9
...@@ -57,12 +57,13 @@ class Transformable(): ...@@ -57,12 +57,13 @@ class Transformable():
return obj return obj
class Parent(): class Parent():
def __init__(self, children=[]): def __init__(self, children=[], remove_duplicate_bonded_terms=False):
self.children = [] self.children = []
for x in children: for x in children:
self.add(x) self.add(x)
# self.children = children # self.children = children
self.remove_duplicate_bonded_terms = remove_duplicate_bonded_terms
self.bonds = [] self.bonds = []
self.angles = [] self.angles = []
self.dihedrals = [] self.dihedrals = []
...@@ -130,32 +131,64 @@ class Parent(): ...@@ -130,32 +131,64 @@ class Parent():
ret = self.bonds ret = self.bonds
for c in self.children: for c in self.children:
if isinstance(c,Parent): ret.extend( c.get_bonds() ) if isinstance(c,Parent): ret.extend( c.get_bonds() )
if self.remove_duplicate_bonded_terms:
return list(set(ret))
else:
return ret return ret
def get_angles(self): def get_angles(self):
ret = self.angles ret = self.angles
for c in self.children: for c in self.children:
if isinstance(c,Parent): ret.extend( c.get_angles() ) if isinstance(c,Parent): ret.extend( c.get_angles() )
if self.remove_duplicate_bonded_terms:
return list(set(ret))
else:
return ret return ret
def get_dihedrals(self): def get_dihedrals(self):
ret = self.dihedrals ret = self.dihedrals
for c in self.children: for c in self.children:
if isinstance(c,Parent): ret.extend( c.get_dihedrals() ) if isinstance(c,Parent): ret.extend( c.get_dihedrals() )
if self.remove_duplicate_bonded_terms:
return list(set(ret))
else:
return ret return ret
def get_impropers(self): def get_impropers(self):
ret = self.impropers ret = self.impropers
for c in self.children: for c in self.children:
if isinstance(c,Parent): ret.extend( c.get_impropers() ) if isinstance(c,Parent): ret.extend( c.get_impropers() )
if self.remove_duplicate_bonded_terms:
return list(set(ret))
else:
return ret return ret
def get_exclusions(self): def get_exclusions(self):
ret = self.exclusions ret = self.exclusions
for c in self.children: for c in self.children:
if isinstance(c,Parent): ret.extend( c.get_exclusions() ) if isinstance(c,Parent): ret.extend( c.get_exclusions() )
if self.remove_duplicate_bonded_terms:
return list(set(ret))
else:
return ret return ret
## Removed because prohibitively slow
# def remove_duplicate_terms(self):
# for key in "bonds angles dihedrals impropers exclusions".split():
# self.remove_duplicate_item(key)
# def remove_duplicate_item(self, dict_key, existing=None):
# if existing is None: existing = []
# ret = [i for i in list(set(self.__dict__[dict_key])) if i not in existing]
# self.__dict__[dict_key] = ret
# existing.extend(ret)
# for c in self.children:
# if isinstance(c,Parent):
# ret = ret + c.remove_duplicate_item(dict_key, existing)
# return ret
def __iter__(self): def __iter__(self):
## TODO: decide if this is the nicest way to do it! ## TODO: decide if this is the nicest way to do it!
"""Depth-first iteration through tree""" """Depth-first iteration through tree"""
...@@ -354,10 +387,11 @@ class Group(Transformable, Parent, Child): ...@@ -354,10 +387,11 @@ class Group(Transformable, Parent, Child):
def __init__(self, name=None, children = [], parent=None, def __init__(self, name=None, children = [], parent=None,
position = np.array((0,0,0)), position = np.array((0,0,0)),
orientation = np.array(((1,0,0),(0,1,0),(0,0,1))), orientation = np.array(((1,0,0),(0,1,0),(0,0,1))),
remove_duplicate_bonded_terms = False,
): ):
Transformable.__init__(self, position, orientation) Transformable.__init__(self, position, orientation)
Child.__init__(self, parent) # Initialize Child first Child.__init__(self, parent) # Initialize Child first
Parent.__init__(self, children) Parent.__init__(self, children, remove_duplicate_bonded_terms)
self.name = name self.name = name
self.isClone = False self.isClone = False
...@@ -393,9 +427,9 @@ class Group(Transformable, Parent, Child): ...@@ -393,9 +427,9 @@ class Group(Transformable, Parent, Child):
class PdbModel(Transformable, Parent): class PdbModel(Transformable, Parent):
def __init__(self, children=[], dimensions=None): def __init__(self, children=[], dimensions=None, remove_duplicate_bonded_terms=False):
Transformable.__init__(self,(0,0,0)) Transformable.__init__(self,(0,0,0))
Parent.__init__(self, children) Parent.__init__(self, children, remove_duplicate_bonded_terms)
self.dimensions = dimensions self.dimensions = dimensions
self.particles = [p for p in self] self.particles = [p for p in self]
self.cacheInvalid = True self.cacheInvalid = True
...@@ -554,8 +588,8 @@ class PdbModel(Transformable, Parent): ...@@ -554,8 +588,8 @@ class PdbModel(Transformable, Parent):
class ArbdModel(PdbModel): class ArbdModel(PdbModel):
def __init__(self, children, dimensions=(1000,1000,1000), temperature=291, timestep=50e-6, cutoff=50, decompPeriod=10000, pairlistDistance=None, nonbondedResolution=0.1): def __init__(self, children, dimensions=(1000,1000,1000), temperature=291, timestep=50e-6, cutoff=50, decompPeriod=10000, pairlistDistance=None, nonbondedResolution=0.1, remove_duplicate_bonded_terms=True):
PdbModel.__init__(self, children, dimensions) PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms)
self.temperature = temperature self.temperature = temperature
self.timestep = timestep self.timestep = timestep
......
...@@ -531,8 +531,8 @@ def run_simulation_protocol( output_name, job_id, json_data, ...@@ -531,8 +531,8 @@ def run_simulation_protocol( output_name, job_id, json_data,
remove_long_bonds=False, remove_long_bonds=False,
gpu = 0, gpu = 0,
directory=None, directory=None,
coarse_steps = 1e5+1, coarse_steps = 1e7+1,
fine_steps = 1e5+1 fine_steps = 1e7+1
): ):
ret = None ret = None
......
...@@ -19,8 +19,8 @@ if __name__ == '__main__': ...@@ -19,8 +19,8 @@ if __name__ == '__main__':
print("WARNING: skipping unreadable json file {}".format(f)) print("WARNING: skipping unreadable json file {}".format(f))
continue continue
run_simulation_protocol( out, "job-"+out+"-", data, gpu=1 )
# 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("WARNING: failed to simulate {}".format(f))
...@@ -735,7 +735,7 @@ class Segment(ConnectableElement, Group): ...@@ -735,7 +735,7 @@ class Segment(ConnectableElement, Group):
## Stop if we found the 3prime end ## Stop if we found the 3prime end
if l.on_fwd_strand == is_fwd and l.type_ == "3prime": if l.on_fwd_strand == is_fwd and l.type_ == "3prime":
print(" found end at",l) # print(" found end at",l)
return pos, None, None, None, None return pos, None, None, None, None
## Check location connections ## Check location connections
...@@ -745,15 +745,15 @@ class Segment(ConnectableElement, Group): ...@@ -745,15 +745,15 @@ class Segment(ConnectableElement, Group):
## Found a location on the same strand? ## Found a location on the same strand?
if l.on_fwd_strand == is_fwd: if l.on_fwd_strand == is_fwd:
print(" passing through",l) # print(" passing through",l)
print("from {}, connection {} to {}".format(nt_pos,l,B)) # print("from {}, connection {} to {}".format(nt_pos,l,B))
Bpos = B.get_nt_pos() Bpos = B.get_nt_pos()
return pos, B.container, Bpos, B.on_fwd_strand, 0.5 return pos, B.container, Bpos, B.on_fwd_strand, 0.5
## Stop at other strand crossovers so basepairs line up ## Stop at other strand crossovers so basepairs line up
elif c.type_ == "crossover": elif c.type_ == "crossover":
if nt_pos == pos: continue if nt_pos == pos: continue
print(" pausing at",l) # print(" pausing at",l)
return pos, l.container, pos+(2*is_fwd-1), is_fwd, 0 return pos, l.container, pos+(2*is_fwd-1), is_fwd, 0
import pdb import pdb
...@@ -805,14 +805,15 @@ class Segment(ConnectableElement, Group): ...@@ -805,14 +805,15 @@ class Segment(ConnectableElement, Group):
self.beads = [] self.beads = []
if True: if True:
print("WARNING: DEBUG") ## TODO: remove this if duplicates are never found
# print("Searching for duplicate particles...")
## Remove duplicates, preserving order ## Remove duplicates, preserving order
tmp = [] tmp = []
for c in new_children: for c in new_children:
if c not in tmp: if c not in tmp:
tmp.append(c) tmp.append(c)
else: else:
print(" duplicate particle found!") print(" DUPLICATE PARTICLE FOUND!")
new_children = tmp new_children = tmp
for b in new_children: for b in new_children:
...@@ -878,7 +879,7 @@ class Segment(ConnectableElement, Group): ...@@ -878,7 +879,7 @@ class Segment(ConnectableElement, Group):
# import pdb # import pdb
# pdb.set_trace() # pdb.set_trace()
print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2])) # print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
e_ds = eb2.get_contour_position(self) - eb1.get_contour_position(self) e_ds = eb2.get_contour_position(self) - eb1.get_contour_position(self)
num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead ) num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead )
ds = e_ds / (num_beads+1) ds = e_ds / (num_beads+1)
...@@ -1636,7 +1637,7 @@ class SegmentModel(ArbdModel): ...@@ -1636,7 +1637,7 @@ class SegmentModel(ArbdModel):
# if s1.name in ("49-3","49-4") and s2.name in ("49-3","49-4"): # if s1.name in ("49-3","49-4") and s2.name in ("49-3","49-4"):
# if s1.name in ("51-2","51-3") and s2.name in ("51-2","51-3"): # if s1.name in ("51-2","51-3") and s2.name in ("51-2","51-3"):
# pdb.set_trace() # 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
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
...@@ -2040,6 +2041,8 @@ class SegmentModel(ArbdModel): ...@@ -2040,6 +2041,8 @@ class SegmentModel(ArbdModel):
pot = self.get_dihedral_potential(k,t0) pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( n1,n2,n3,n4, pot ) self.add_dihedral( n1,n2,n3,n4, pot )
# ## remove duplicate potentials; ## TODO ensure that they aren't added twice in the first place?
# self.remove_duplicate_terms()
def walk_through_helices(segment, direction=1, processed_segments=None): def walk_through_helices(segment, direction=1, processed_segments=None):
""" """
...@@ -2182,7 +2185,7 @@ class SegmentModel(ArbdModel): ...@@ -2182,7 +2185,7 @@ class SegmentModel(ArbdModel):
if locs is None: continue if locs is None: continue
# for pos, is_fwd in locs: # for pos, is_fwd in locs:
for l in locs: for l in locs:
print("Tracing",l) # print("Tracing",l)
# TODOTODO # TODOTODO
pos = seg.contour_to_nt_pos(l.address, round_nt=True) pos = seg.contour_to_nt_pos(l.address, round_nt=True)
is_fwd = l.on_fwd_strand is_fwd = l.on_fwd_strand
......
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