Skip to content
Snippets Groups Projects
Commit c117923f authored by cmaffeo2's avatar cmaffeo2
Browse files

Added exclusions andnd angles to segmentmodel

parent d3361ce8
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
nonbonded.py
+ 8
− 7
View file @ c117923f
...@@ -67,7 +67,7 @@ class HarmonicPotential(): ...@@ -67,7 +67,7 @@ class HarmonicPotential():
self.prefix = prefix self.prefix = prefix
self.periodic = False self.periodic = False
self.type_ = "None" self.type_ = "None"
self.kscale_ = None self.kscale_ = None # only used for
def filename(self): def filename(self):
# raise NotImplementedError("Not implemented") # raise NotImplementedError("Not implemented")
...@@ -90,8 +90,9 @@ class HarmonicPotential(): ...@@ -90,8 +90,9 @@ class HarmonicPotential():
u = 0.5*self.k*dr**2 u = 0.5*self.k*dr**2
if self.maxForce is not None: maxForce = self.maxForce
assert(self.maxForce > 0) if maxForce is not None:
assert(maxForce > 0)
f = np.diff(u)/np.diff(r) f = np.diff(u)/np.diff(r)
f[f>maxForce] = maxForce f[f>maxForce] = maxForce
f[f<-maxForce] = -maxForce f[f<-maxForce] = -maxForce
...@@ -121,15 +122,15 @@ class HarmonicBond(HarmonicPotential): ...@@ -121,15 +122,15 @@ class HarmonicBond(HarmonicPotential):
self.kscale_ = 1.0 self.kscale_ = 1.0
class HarmonicAngle(HarmonicPotential): class HarmonicAngle(HarmonicPotential):
def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, prefix="potentials/"): def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, prefix="potentials/"):
HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix) HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix)
self.type_ = "angle" self.type_ = "angle"
self.kscale_ = (180.0/pi)**2 self.kscale_ = (180.0/np.pi)**2
class HarmonicDihedral(HarmonicPotential): class HarmonicDihedral(HarmonicPotential):
def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, prefix="potentials/"): def __init__(self, k, r0, rRange=(-180,180), resolution=0.1, maxForce=None, prefix="potentials/"):
HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix) HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, prefix)
self.periodic = True self.periodic = True
self.type_ = "dihedral" self.type_ = "dihedral"
self.kscale_ = (180.0/pi)**2 self.kscale_ = (180.0/np.pi)**2
segmentmodel.py
+ 111
− 50
View file @ c117923f
...@@ -94,7 +94,7 @@ class Segment(ConnectableElement, Group): ...@@ -94,7 +94,7 @@ class Segment(ConnectableElement, Group):
ret.append( tmp ) ret.append( tmp )
return ret return ret
def get_beads_before_bead(self, bead, number, inclusive=True): def get_beads_before_bead(self, bead, number, inclusive=False):
## Assume that consecutive beads in self.children are bonded ## Assume that consecutive beads in self.children are bonded
i = self.children.index(bead) i = self.children.index(bead)
l = len(self.children) l = len(self.children)
...@@ -105,7 +105,7 @@ class Segment(ConnectableElement, Group): ...@@ -105,7 +105,7 @@ class Segment(ConnectableElement, Group):
if inclusive: start = 0 if inclusive: start = 0
return [self.children[i-j] for j in range(start,number)] return [self.children[i-j] for j in range(start,number)]
def get_beads_after_bead(self, bead, number, inclusive=True): def get_beads_after_bead(self, bead, number, inclusive=False):
## Assume that consecutive beads in self.children are bonded ## Assume that consecutive beads in self.children are bonded
i = self.children.index(bead) i = self.children.index(bead)
l = len(self.children) l = len(self.children)
...@@ -116,6 +116,10 @@ class Segment(ConnectableElement, Group): ...@@ -116,6 +116,10 @@ class Segment(ConnectableElement, Group):
if inclusive: start = 0 if inclusive: start = 0
return [self.children[i+i] for j in range(start,number)] return [self.children[i+i] for j in range(start,number)]
# def get_bead_pairs_within(self, cutoff):
# for b1,b2 in self.get_all_consecutive_beads(self, number)
def _generate_beads(self, bead_model, max_nts_per_bead=4): def _generate_beads(self, bead_model, max_nts_per_bead=4):
""" Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """ """ Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
...@@ -334,39 +338,68 @@ class SegmentModel(ArbdModel): ...@@ -334,39 +338,68 @@ class SegmentModel(ArbdModel):
if c.type_ == "intrahelical": if c.type_ == "intrahelical":
b1,b2 = [loc.particle for loc in (c.A,c.B)] b1,b2 = [loc.particle for loc in (c.A,c.B)]
for b in (b1,b2): assert( b is not None ) for b in (b1,b2): assert( b is not None )
## TODO: make this code more robust
## TODO: find bug that makes things be out of order
try: try:
b0 = s1.get_beads_before_bead(b1,1)[0] b0 = s1.get_beads_before_bead(b1,1)
assert(len(b0) == 1)
b0 = b0[0]
assert( b0 is not None ) assert( b0 is not None )
ret.append( [b0,b1,b2] ) ret.append( [b0,b1,b2] )
except: except:
... ...
try: try:
b0 = s1.get_beads_after_bead(b1,1)[0] b0 = s1.get_beads_after_bead(b1,1)
assert(len(b0) == 1)
b0 = b0[0]
assert( b0 is not None ) assert( b0 is not None )
ret.append( [b2,b1,b0] ) ret.append( [b2,b1,b0] )
except: except:
... ...
try: try:
b3 = s2.get_beads_before_bead(b2,1)[0] b3 = s2.get_beads_before_bead(b2,1)
assert(len(b3) == 1)
b3 = b3[0]
assert( b3 is not None ) assert( b3 is not None )
ret.append( [b3,b2,b1] ) ret.append( [b3,b2,b1] )
except: except:
... ...
try: try:
b3 = s2.get_beads_after_bead(b2,1)[0] b3 = s2.get_beads_after_bead(b2,1)
assert(len(b3) == 1)
b3 = b3[0]
assert( b3 is not None ) assert( b3 is not None )
ret.append( [b1,b2,b3] ) ret.append( [b1,b2,b3] )
except: except:
... ...
return ret return ret
# def _get_intrahelical_bead_pairs_within(self, cutoff):
# dist = dict()
# for b1,b2 in self._get_intrahelical_beads:
# dist(b1,b2)
# ret = []
# for s in self.segments:
# ret.extend( s.get_bead_pairs_within(cutoff) )
# for s1 in self.segments:
# for c in s1.connections:
# if c.A.container != s1: continue
# s2 = c.B.container
# if c.type_ == "intrahelical":
# ret
def _get_potential(self, type_, kSpring, d): def _get_potential(self, type_, kSpring, d):
key = (type_,kSpring,d) key = (type_,kSpring,d)
if key not in self._bonded_potential: if key not in self._bonded_potential:
if type_ == "bond": if type_ == "bond":
self._bonded_potential[key] = HarmonicBond(kSpring,d) self._bonded_potential[key] = HarmonicBond(kSpring,d)
elif type_ == "angle": elif type_ == "angle":
self._bonded_potential[key] = HarmonicAngle(kSpring,d) self._bonded_potential[key] = HarmonicAngle(kSpring,d) # , resolution = 1, maxForce=0.1)
elif type_ == "dihedral": elif type_ == "dihedral":
self._bonded_potential[key] = HarmonicDihedral(kSpring,d) self._bonded_potential[key] = HarmonicDihedral(kSpring,d)
else: else:
...@@ -415,8 +448,8 @@ class SegmentModel(ArbdModel): ...@@ -415,8 +448,8 @@ class SegmentModel(ArbdModel):
beadtype_s[key] = b.type_ = t beadtype_s[key] = b.type_ = t
""" Add intrahelical potentials """ """ Add intrahelical bond potentials """
## First replace intrahelical bonds dists = dict() # built for later use
for b1,b2 in self._get_intrahelical_beads(): for b1,b2 in self._get_intrahelical_beads():
if b1.parent == b2.parent: if b1.parent == b2.parent:
sep = 0.5*(b1.num_nts+b2.num_nts) sep = 0.5*(b1.num_nts+b2.num_nts)
...@@ -425,17 +458,81 @@ class SegmentModel(ArbdModel): ...@@ -425,17 +458,81 @@ class SegmentModel(ArbdModel):
sep = 1 sep = 1
parent = self parent = self
conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
if b1.type_.name[0] == "D" and b1.type_.name[0] == "D": if b1.type_.name[0] == "D" and b1.type_.name[0] == "D":
k = 10.0/np.sqrt(sep) # TODO: determine from simulations elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
d = 3.4*sep d = 3.4*sep
k = conversion*elastic_modulus/d
else: else:
## TODO: get correct numbers from ssDNA model ## TODO: get better numbers our ssDNA model
k = 1.0/np.sqrt(sep) elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
d = 5*sep d = 5*sep
k = conversion*elastic_modulus/d
if b1 not in dists:
dists[b1] = []
if b2 not in dists:
dists[b2] = []
dists[b1].append([b2,sep])
dists[b2].append([b1,sep])
# dists[[b1,b2]] = dists[[b2,b1]] = sep
bond = self.get_bond_potential(k,d) bond = self.get_bond_potential(k,d)
parent.add_bond( b1, b2, bond, exclude=True ) parent.add_bond( b1, b2, bond, exclude=True )
""" Add intrahelical angle potentials """
for b1,b2,b3 in self._get_intrahelical_angle_beads():
sep = 0
if b1.parent == b2.parent:
sep += 0.5*(b1.num_nts+b2.num_nts)
else:
sep += 1
if b2.parent == b3.parent:
sep += 0.5*(b2.num_nts+b3.num_nts)
else:
sep += 1
if b1.parent == b2.parent and b2.parent == b3.parent:
parent = b1.parent
else:
parent = self
kT = 0.58622522 # kcal/mol
if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
## <cos(q)> = exp(-s/Lp) = integrate( x^4 exp(-A x^2) / 2, {x, 0, pi} ) / integrate( x^2 exp(-A x^2), {x, 0, pi} )
## <cos(q)> ~ 1 - 3/4A
## where A = k_spring / (2 kT)
k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
else:
## TODO: get correct number from ssDNA model
k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2
# k *= 1e-6
angle = self.get_angle_potential(k,180)
parent.add_angle( b1, b2, b3, angle )
""" Add intrahelical exclusions """
beads = dists.keys()
def _recursively_get_beads_within(b1,d,done=[]):
ret = []
for b2,sep in dists[b1]:
if b2 in done: continue
if sep < d:
ret.append( b2 )
done.append( b2 )
tmp = _recursively_get_beads_within(b2, d-sep, done)
if len(tmp) > 0: ret.extend(tmp)
return ret
exclusions = set()
for b1 in beads:
exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
for b1,b2 in exclusions:
if b1.parent == b2.parent:
parent = b1.parent
else:
parent = self
parent.add_exclusion( b1, b2 )
""" Add connection potentials """ """ Add connection potentials """
## TODO ## TODO
...@@ -449,39 +546,3 @@ class SegmentModel(ArbdModel): ...@@ -449,39 +546,3 @@ class SegmentModel(ArbdModel):
# r # r
# ... # ...
if __name__ == "__main__":
seg1 = DoubleStrandedSegment("strand", num_nts = 46)
seg2 = SingleStrandedSegment("strand",
start_position = seg1.end_position + np.array((1,0,1)),
num_nts = 12)
seg3 = SingleStrandedSegment("strand",
start_position = seg1.start_position + np.array((-1,0,-1)),
end_position = seg1.end_position + np.array((-1,0,1)),
num_nts = 128)
seg1.start3
seg1.start5
seg1.end3
seg1.end5
seg1.connect_end3(seg2)
seg1.connect_end5(seg3)
seg1.connect_start3(seg3)
model = SegmentModel( [seg1, seg2, seg3],
dimensions=(5000,5000,5000),
)
model.useNonbondedScheme( nbDnaScheme )
model.simulate( outputPrefix = 'strand-test', outputPeriod=1e4, numSteps=1e6, gpu=1 )
# seg = SingleStrandedSegment("strand", num_nts = 21)
# generate_bead_model( [seg] )
# for b in seg:
# print(b.num_nts, b.position)
segmentmodel_example.py
+ 4
− 3
View file @ c117923f
...@@ -17,12 +17,12 @@ if __name__ == "__main__": ...@@ -17,12 +17,12 @@ if __name__ == "__main__":
seg1 = DoubleStrandedSegment("strand", num_nts = 46) seg1 = DoubleStrandedSegment("strand", num_nts = 46)
seg2 = SingleStrandedSegment("strand", seg2 = SingleStrandedSegment("strand",
start_position = seg1.end_position + np.array((1,0,1)), start_position = seg1.end_position + np.array((5,0,5)),
num_nts = 12) num_nts = 12)
seg3 = SingleStrandedSegment("strand", seg3 = SingleStrandedSegment("strand",
start_position = seg1.start_position + np.array((-1,0,-1)), start_position = seg1.start_position + np.array((-5,0,-5)),
end_position = seg1.end_position + np.array((-1,0,1)), end_position = seg1.end_position + np.array((-5,0,5)),
num_nts = 128) num_nts = 128)
seg1.start3 seg1.start3
...@@ -37,5 +37,6 @@ if __name__ == "__main__": ...@@ -37,5 +37,6 @@ if __name__ == "__main__":
model = SegmentModel( [seg1, seg2, seg3], model = SegmentModel( [seg1, seg2, seg3],
dimensions=(5000,5000,5000), dimensions=(5000,5000,5000),
) )
model.useNonbondedScheme( nbDnaScheme ) model.useNonbondedScheme( nbDnaScheme )
model.simulate( outputPrefix = 'strand-test', outputPeriod=1e4, numSteps=1e6, gpu=1 ) model.simulate( outputPrefix = 'strand-test', outputPeriod=1e4, numSteps=1e6, gpu=1 )
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