Skip to content
GitLab
Commit 4b50fbbc authored by cmaffeo2's avatar cmaffeo2
Browse files

Added direction to crossovers

parent d8689f42
Hide whitespace changes
Inline Side-by-side
segmentmodel.py
View file @ 4b50fbbc
......@@ -370,43 +370,22 @@ class DoubleStrandedSegment(Segment):
end3 = end3.end3
self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
def crossover(self, nt, other, other_nt, strands=None, other_strands=None):
def crossover(self, nt, other, other_nt, strands_fwd=[True,False]):
""" Add a crossover between two helices """
## Validate other, nt, other_nt
## TODO
## Fix input
if strands is None and other_strands is None:
## Regular double crossover
strands = ['fwd','rev']
other_strands = list(strands[::-1])
elif strands is None:
strands = []
for d in other_strands:
if d == 'fwd': strands.append('rev')
elif d == 'rev': strands.append('fwd')
elif other_strands == None:
other_strands = []
for d in strands:
if d == 'fwd': other_strands.append('rev')
elif d == 'rev': other_strands.append('fwd')
if isinstance(strands,str):
strands = [strands]
if isinstance(other_strands,str):
other_strands = [other_strands]
## Create locations, connections and add to segments
c = nt/self.num_nts
print(self.name,nt,c)
assert(c >= 0 and c <= 1)
loc = Location( self, address=c, type_ = strands )
loc = Location( self, address=c, type_ = strands_fwd[0] )
c = other_nt/other.num_nts
print(other.name,other_nt,c)
assert(c >= 0 and c <= 1)
other_loc = Location( other, address=c, type_ = other_strands )
self._connect(other, Connection( loc, other_loc, type_="crossover" )) # TODO: use multiple crossover types
other_loc = Location( other, address=c, type_ = strands_fwd[1] )
self._connect(other, Connection( loc, other_loc, type_="crossover" ))
## Real work
def _connect_ends(self, end1, end2, type_, force_connection):
......@@ -746,13 +725,14 @@ class SegmentModel(ArbdModel):
A.particle = B.particle = b
else:
## TODO fix this for ssDNA vs dsDNA
## TODO fix this for ssDNA vs dsDNA (maybe a1/a2 should be calculated differently)
a1,a2 = [l.address for l in (A,B)]
a1,a2 = [a - (0.5/s.num_nts) if a == 0 else a + (0.5/s.num_nts) for a,s in zip((a1,a2),(s1,s2))]
b = s1._generate_one_bead(a1,0)
if isinstance(s2,DoubleStrandedSegment):
b = s2._generate_one_bead(a2,0)
else:
b = s1._generate_one_bead(a1,0)
A.particle = B.particle = b
...
""" Generate beads at other junctions """
for c,A,B in self.get_connections(exclude="intrahelical"):
......@@ -858,11 +838,13 @@ class SegmentModel(ArbdModel):
# print(sep,d,k)
if b1 not in dists:
dists[b1] = []
dists[b1] = dict()
if b2 not in dists:
dists[b2] = []
dists[b1].append([b2,sep])
dists[b2].append([b1,sep])
dists[b2] = dict()
# dists[b1].append([b2,sep])
# dists[b2].append([b1,sep])
dists[b1][b2] = sep
dists[b2][b1] = sep
# dists[[b1,b2]] = dists[[b2,b1]] = sep
bond = self.get_bond_potential(k,d)
......@@ -895,7 +877,7 @@ class SegmentModel(ArbdModel):
beads = dists.keys()
def _recursively_get_beads_within(b1,d,done=[]):
ret = []
for b2,sep in dists[b1]:
for b2,sep in dists[b1].items():
if b2 in done: continue
if sep < d:
ret.append( b2 )
......@@ -926,7 +908,7 @@ class SegmentModel(ArbdModel):
for b1 in beads:
if "orientation_bead" not in b1.__dict__: continue
for b2,sep in dists[b1]:
for b2,sep in dists[b1].items():
if "orientation_bead" not in b2.__dict__: continue
p1,p2 = [b.parent for b in (b1,b2)]
......@@ -953,11 +935,13 @@ class SegmentModel(ArbdModel):
""" Add connection potentials """
for c,A,B in self.get_connections("terminal_crossover"):
## TODO: use a better description here
b1,b2 = [loc.particle for loc in (c.A,c.B)]
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
for c,A,B in self.get_connections("crossover"):
## TODO: avoid double-counting for double-crossovers
b1,b2 = [loc.particle for loc in (c.A,c.B)]
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
......@@ -966,57 +950,68 @@ class SegmentModel(ArbdModel):
## TODO
...
else:
## TODO: fwd/rev
k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
if 'orientation_bead' in b1.__dict__:
t0 = 90 + 60
# if f1: t0 -= 120
if A.type_: t0 -= 120
o = b1.orientation_bead
pot = self.get_angle_potential(k,t0)
self.add_angle( o,b1,b2, pot )
else:
t0 = 90 + 60
# if f2: t0 -= 120
if B.type_: t0 -= 120
o = b2.orientation_bead
pot = self.get_angle_potential(k,t0)
self.add_angle( b1,b2,o, pot )
## TODO: fix the following to ensure that neighbors are indeed above and below
### preferably with a function call to b
u1,u2 = [b.intrahelical_neighbors[-1] for b in (b1,b2)] # TODO: fix this
d1,d2 = [b.intrahelical_neighbors[0] for b in (b1,b2)] # TODO: fix this
## TODO make k depend on sep
k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
t0 = 90
pot = self.get_dihedral_potential(k,t0)
if 'orientation_bead' in b1.__dict__:
o1 = b1.orientation_bead
if u2 is not None:
k = k_fn( dists[b2][u2] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o1,b1,b2,u2, pot )
elif d2 is not None:
self.add_dihedral( o2,b2,b1,d2, pot )
k = k_fn( dists[b2][d2] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o1,b1,b2,d2, pot )
if 'orientation_bead' in b2.__dict__:
o2 = b2.orientation_bead
if u1 is not None:
k = k_fn( dists[b1][u1] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o2,b2,b1,u1, pot )
elif d1 is not None:
k = k_fn( dists[b1][d1] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o2,b2,b1,d1, pot )
if 'orientation_bead' in b1.__dict__ and 'orientation_bead' in b2.__dict__:
if u1 is not None and u2 is not None:
t0 = 0
k = k_fn( 0.5*(dists[b1][u1]+dists[b2][u2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( u1,b1,b2,u2, pot )
elif d1 is not None and d2 is not None:
t0 = 0
k = k_fn( 0.5*(dists[b1][d1]+dists[b2][d2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( d1,b1,b2,d2, pot )
elif d1 is not None and u2 is not None:
t0 = 180
k = k_fn( 0.5*(dists[b1][d1]+dists[b2][u2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( d1,b1,b2,u2, pot )
elif u1 is not None and d2 is not None:
t0 = 180
k = k_fn( 0.5*(dists[b1][u1]+dists[b2][d2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( u1,b1,b2,d2, pot )
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