diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index b6d4a6c2e56f3a7fd4ca80f5d7a06fd24573644e..28b96ddad7ce06f995921f8c30d09d0fd76325dc 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -1052,97 +1052,110 @@ class Segment(ConnectableElement, Group):
assert(len(old_beads) == len(self.beads))
- def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
+ def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead, one_bead_per_monomer=False):
""" Generate beads (positions, types, etc) and bonds, angles, dihedrals, exclusions """
## TODO: decide whether to remove bead_model argument
## (currently unused)
- ## First find points between-which beads must be generated
- # 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 == "S001":
- # pdb.set_trace()
+ tmp_children = []
- # pdb.set_trace()
- existing_beads0 = { (l.particle, l.particle.get_contour_position(self,l.address))
- for l in self.locations if l.particle is not None }
- existing_beads = sorted( list(existing_beads0), key=lambda bc: bc[1] )
+ if one_bead_per_monomer:
+ last = None
+ for i in range(self.num_nt):
+ s = self.nt_pos_to_contour(i)
+ b = self._generate_one_bead(s,1)
- # if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
- # pdb.set_trace()
- # return
-
- for b,c in existing_beads:
- assert(b.parent is not None)
-
- ## Add ends if they don't exist yet
- ## TODOTODO: test 1 nt segments?
- if len(existing_beads) == 0 or (existing_beads[0][0].get_nt_position(self,0) >= 0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__):
- # if len(existing_beads) > 0:
- # assert(existing_beads[0].get_nt_position(self) >= 0.5)
- c = self.nt_pos_to_contour(0)
- if self.num_nt == 1: c -= 0.4
- b = self._generate_one_bead(c, 0)
- existing_beads = [(b,0)] + existing_beads
-
- if (existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__) or len(existing_beads)==1:
- c = self.nt_pos_to_contour(self.num_nt-1)
- if self.num_nt == 1: c += 0.4
- b = self._generate_one_bead(c, 0)
- existing_beads.append( (b,1) )
- assert(len(existing_beads) > 1)
-
- ## Walk through existing_beads, add beads between
- tmp_children = [] # build list of children in nice order
- last = None
+ if last is not None:
+ last.make_intrahelical_neighbor(b)
+ last = b
+ tmp_children.append(b)
+ else:
+ ## First find points between-which beads must be generated
+ # 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 == "S001":
+ # pdb.set_trace()
- for I in range(len(existing_beads)-1):
- eb1,eb2 = [existing_beads[i][0] for i in (I,I+1)]
- ec1,ec2 = [existing_beads[i][1] for i in (I,I+1)]
- assert( (eb1,ec1) is not (eb2,ec2) )
+ # pdb.set_trace()
+ existing_beads0 = { (l.particle, l.particle.get_contour_position(self,l.address))
+ for l in self.locations if l.particle is not None }
+ existing_beads = sorted( list(existing_beads0), key=lambda bc: bc[1] )
- # if np.isclose(eb1.position[2], eb2.position[2]):
- # import pdb
+ # if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
# pdb.set_trace()
+ # return
+
+ for b,c in existing_beads:
+ assert(b.parent is not None)
+
+ ## Add ends if they don't exist yet
+ ## TODOTODO: test 1 nt segments?
+ if len(existing_beads) == 0 or (existing_beads[0][0].get_nt_position(self,0) >= 0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__):
+ # if len(existing_beads) > 0:
+ # assert(existing_beads[0].get_nt_position(self) >= 0.5)
+ c = self.nt_pos_to_contour(0)
+ if self.num_nt == 1: c -= 0.4
+ b = self._generate_one_bead(c, 0)
+ existing_beads = [(b,0)] + existing_beads
+
+ if (existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__) or len(existing_beads)==1:
+ c = self.nt_pos_to_contour(self.num_nt-1)
+ if self.num_nt == 1: c += 0.4
+ b = self._generate_one_bead(c, 0)
+ existing_beads.append( (b,1) )
+ assert(len(existing_beads) > 1)
+
+ ## Walk through existing_beads, add beads between
+ last = None
+
+ for I in range(len(existing_beads)-1):
+ eb1,eb2 = [existing_beads[i][0] for i in (I,I+1)]
+ ec1,ec2 = [existing_beads[i][1] for i in (I,I+1)]
+ assert( (eb1,ec1) is not (eb2,ec2) )
+
+ # if np.isclose(eb1.position[2], eb2.position[2]):
+ # import pdb
+ # pdb.set_trace()
- # print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
- e_ds = ec2-ec1
- num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead )
-
- ## Ensure there is a ssDNA bead between dsDNA beads
- if num_beads == 0 and isinstance(self,SingleStrandedSegment) and isinstance(eb1.parent,DoubleStrandedSegment) and isinstance(eb2.parent,DoubleStrandedSegment):
- num_beads = 1
- ## TODO similarly ensure there is a dsDNA bead between ssDNA beads
-
- ds = e_ds / (num_beads+1)
- nts = ds*self.num_nt
- eb1.num_nt += 0.5*nts
- eb2.num_nt += 0.5*nts
-
- ## Add beads
- if eb1.parent == self:
- tmp_children.append(eb1)
-
- s0 = ec1
- if last is not None:
- last.make_intrahelical_neighbor(eb1)
- last = eb1
- for j in range(num_beads):
- s = ds*(j+1) + s0
- # if self.name in ("51-2","51-3"):
- # if self.name in ("31-2",):
- # print(" adding bead at {}".format(s))
- b = self._generate_one_bead(s,nts)
-
- last.make_intrahelical_neighbor(b)
- last = b
- tmp_children.append(b)
+ # print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
+ e_ds = ec2-ec1
+ num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead )
+
+ ## Ensure there is a ssDNA bead between dsDNA beads
+ if num_beads == 0 and isinstance(self,SingleStrandedSegment) and isinstance(eb1.parent,DoubleStrandedSegment) and isinstance(eb2.parent,DoubleStrandedSegment):
+ num_beads = 1
+ ## TODO similarly ensure there is a dsDNA bead between ssDNA beads
+
+ ds = e_ds / (num_beads+1)
+ nts = ds*self.num_nt
+ eb1.num_nt += 0.5*nts
+ eb2.num_nt += 0.5*nts
+
+ ## Add beads
+ if eb1.parent == self:
+ tmp_children.append(eb1)
- last.make_intrahelical_neighbor(eb2)
+ s0 = ec1
+ if last is not None:
+ last.make_intrahelical_neighbor(eb1)
+ last = eb1
+ for j in range(num_beads):
+ s = ds*(j+1) + s0
+ # if self.name in ("51-2","51-3"):
+ # if self.name in ("31-2",):
+ # print(" adding bead at {}".format(s))
+ b = self._generate_one_bead(s,nts)
+
+ last.make_intrahelical_neighbor(b)
+ last = b
+ tmp_children.append(b)
+
+ last.make_intrahelical_neighbor(eb2)
+
+ if eb2.parent == self:
+ tmp_children.append(eb2)
- if eb2.parent == self:
- tmp_children.append(eb2)
# if self.name in ("31-2",):
# pdb.set_trace()
self._rebuild_children(tmp_children)
@@ -2249,10 +2262,17 @@ class SegmentModel(ArbdModel):
escapable_twist=escapable_twist)
def generate_bead_model(self,
- max_basepairs_per_bead = 7,
- max_nucleotides_per_bead = 4,
- local_twist=False,
- escapable_twist=True):
+ max_basepairs_per_bead = 7,
+ max_nucleotides_per_bead = 4,
+ local_twist=False,
+ escapable_twist=True,
+ sequence_dependent = False,
+ one_bead_per_monomer = False
+ ):
+ if sequence_dependent and not local_twist:
+ raise ValueError("Incompatible options: if sequence_dependent==True, then local_twist must be True")
+ if sequence_dependent and not one_bead_per_monomer:
+ raise ValueError("Incompatible options: if sequence_dependent==True, then one_bead_per_monomer must be True")
self.children = self.segments # is this okay?
self.clear_beads()
@@ -2291,97 +2311,120 @@ class SegmentModel(ArbdModel):
""" Generate beads at intrahelical junctions """
if self.DEBUG: print( "Adding intrahelical beads at junctions" )
+ if not one_bead_per_monomer:
+ ## 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)]
- ## 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)]
+ assert( A.particle is None )
+ assert( B.particle is None )
- assert( A.particle is None )
- assert( B.particle is None )
+ ## TODO: offload the work here to s1
+ # TODOTODO
+ a1,a2 = [l.address for l in (A,B)]
+ for a in (a1,a2):
+ assert( np.isclose(a,0) or np.isclose(a,1) )
- ## TODO: offload the work here to s1
- # TODOTODO
- a1,a2 = [l.address for l in (A,B)]
- for a in (a1,a2):
- assert( np.isclose(a,0) or np.isclose(a,1) )
-
- ## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
-
- def _get_nearest(seg,address):
- b = seg.get_nearest_bead(address)
- if b is not None:
- assert( b.parent is seg )
- """ if above assertion is true, no problem here """
- if np.abs(b.get_nt_position(seg) - seg.contour_to_nt_pos(address)) > 0.5:
- b = None
- return b
-
- """ Search to see whether bead at location is already found """
- b = None
- if isinstance(s1,DoubleStrandedSegment):
- b = _get_nearest(s1,a1)
- if b is None and isinstance(s2,DoubleStrandedSegment):
- b = _get_nearest(s2,a2)
-
- if b is not None and b.parent not in (s1,s2):
- b = None
+ ## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
- if b is None:
- ## need to generate a bead
- if isinstance(s2,DoubleStrandedSegment):
- b = s2._generate_one_bead(a2,0)
- else:
- b = s1._generate_one_bead(a1,0)
- A.particle = B.particle = b
- b.is_intrahelical = True
- b.locations.extend([A,B])
- if s1 is s2:
- b.is_circular_bead = True
+ def _get_nearest(seg,address):
+ b = seg.get_nearest_bead(address)
+ if b is not None:
+ assert( b.parent is seg )
+ """ if above assertion is true, no problem here """
+ if np.abs(b.get_nt_position(seg) - seg.contour_to_nt_pos(address)) > 0.5:
+ b = None
+ return b
- # pdb.set_trace()
+ """ Search to see whether bead at location is already found """
+ b = None
+ if isinstance(s1,DoubleStrandedSegment):
+ b = _get_nearest(s1,a1)
+ if b is None and isinstance(s2,DoubleStrandedSegment):
+ b = _get_nearest(s2,a2)
+
+ if b is not None and b.parent not in (s1,s2):
+ b = None
+
+ if b is None:
+ ## need to generate a bead
+ if isinstance(s2,DoubleStrandedSegment):
+ b = s2._generate_one_bead(a2,0)
+ else:
+ b = s1._generate_one_bead(a1,0)
+ A.particle = B.particle = b
+ b.is_intrahelical = True
+ b.locations.extend([A,B])
+ if s1 is s2:
+ b.is_circular_bead = True
- """ Generate beads at other junctions """
- for c,A,B in self.get_connections(exclude="intrahelical"):
- s1,s2 = [l.container for l in (A,B)]
- if A.particle is not None and B.particle is not None:
- continue
- # assert( A.particle is None )
- # assert( B.particle is None )
+ # pdb.set_trace()
- ## TODO: offload the work here to s1/s2 (?)
- a1,a2 = [l.address for l in (A,B)]
+ """ Generate beads at other junctions """
+ for c,A,B in self.get_connections(exclude="intrahelical"):
+ s1,s2 = [l.container for l in (A,B)]
+ if A.particle is not None and B.particle is not None:
+ continue
+ # assert( A.particle is None )
+ # assert( B.particle is None )
+
+ ## TODO: offload the work here to s1/s2 (?)
+ a1,a2 = [l.address for l in (A,B)]
+
+ def maybe_add_bead(location, seg, address, ):
+ if location.particle is None:
+ b = seg.get_nearest_bead(address)
+ try:
+ distance = seg.contour_to_nt_pos(np.abs(b.contour_position-address))
+ max_distance = min(max_basepairs_per_bead, max_nucleotides_per_bead)*0.5
+ if "is_intrahelical" in b.__dict__:
+ max_distance = 0.5
+ if distance >= max_distance:
+ raise Exception("except")
+
+ ## combine beads
+ b.update_position( 0.5*(b.contour_position + address) ) # avg position
+ except:
+ b = seg._generate_one_bead(address,0)
+ location.particle = b
+ b.locations.append(location)
+
+ maybe_add_bead(A,s1,a1)
+ maybe_add_bead(B,s2,a2)
+
+ """ Some tests """
+ for c,A,B in self.get_connections("intrahelical"):
+ for l in (A,B):
+ if l.particle is None: continue
+ assert( l.particle.parent is not None )
- def maybe_add_bead(location, seg, address, ):
- if location.particle is None:
- b = seg.get_nearest_bead(address)
- try:
- distance = seg.contour_to_nt_pos(np.abs(b.contour_position-address))
- max_distance = min(max_basepairs_per_bead, max_nucleotides_per_bead)*0.5
- if "is_intrahelical" in b.__dict__:
- max_distance = 0.5
- if distance >= max_distance:
- raise Exception("except")
-
- ## combine beads
- b.update_position( 0.5*(b.contour_position + address) ) # avg position
- except:
- b = seg._generate_one_bead(address,0)
- location.particle = b
- b.locations.append(location)
-
- maybe_add_bead(A,s1,a1)
- maybe_add_bead(B,s2,a2)
-
- """ Some tests """
- for c,A,B in self.get_connections("intrahelical"):
- for l in (A,B):
- if l.particle is None: continue
- assert( l.particle.parent is not None )
+ """ Generate beads in between """
+ if self.DEBUG: print("Generating beads")
+ for s in segments:
+ s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
- """ Generate beads in between """
- if self.DEBUG: print("Generating beads")
- for s in segments:
- s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
+ else:
+ """ Generate beads in each segment """
+ if self.DEBUG: print("Generating beads")
+ for s in segments:
+ s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead, one_bead_per_monomer )
+
+ """ Associate beads with junctions """
+ # for c,A,B in self.get_connections("intrahelical"):
+ for c,A,B in self.get_connections():
+ for l in (A,B):
+ seg = l.container
+ a = l.address
+ b = seg.get_nearest_bead(a)
+
+ l.particle = b
+ b.locations.append(l)
+ """
+ b.is_intrahelical = True
+ b.locations.extend([A,B])
+ if s1 is s2:
+ b.is_circular_bead = True
+ """
# """ Combine beads at junctions as needed """
# for c,A,B in self.get_connections():