Commit d150e38e authored by cmaffeo2's avatar cmaffeo2
Browse files

Added improved WLC model parameters for ssDNA, still imperfect

parent b349e427
......@@ -57,12 +57,11 @@ class TabulatedPotential(NonbondedScheme):
copyfile(self.tableFile, filename)
## Bonded potentials
class HarmonicPotential():
def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
self.k = k
class BasePotential():
def __init__(self, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
self.r0 = r0
self.rRange = rRange
self.resolution = 0.1
self.resolution = resolution
self.maxForce = maxForce
self.prefix = prefix
self.periodic = False
......@@ -71,9 +70,10 @@ class HarmonicPotential():
self.kscale_ = None # only used for
def filename(self):
# raise NotImplementedError("Not implemented")
return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
self.k*self.kscale_, self.r0)
raise NotImplementedError("Not implemented")
def potential(self,dr):
raise NotImplementedError("Not implemented")
def __str__(self):
return self.filename()
......@@ -89,7 +89,7 @@ class HarmonicPotential():
assert(rSpan > 0)
dr = np.mod( dr+0.5*rSpan, rSpan) - 0.5*rSpan
u = 0.5*self.k*dr**2
u = self.potential(dr)
maxForce = self.maxForce
if maxForce is not None:
......@@ -115,6 +115,18 @@ class HarmonicPotential():
np.savetxt( self.filename(), np.array([r, u]).T, fmt="%f" )
class HarmonicPotential(BasePotential):
def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix):
self.k = k
self.kscale_ = None
BasePotential.__init__(self, r0, rRange, resolution, maxForce, max_potential, prefix)
def filename(self):
return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
self.k*self.kscale_, self.r0)
def potential(self,dr):
return 0.5*self.k*dr**2
def __hash__(self):
assert(self.type_ != "None")
return hash((self.type_, self.k, self.r0, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
......@@ -149,3 +161,69 @@ class HarmonicDihedral(HarmonicPotential):
self.type_ = "dihedral"
self.kscale_ = (180.0/np.pi)**2
class WLCSKBond(BasePotential):
""" ## https://aip.scitation.org/doi/full/10.1063/1.4968020 """
def __init__(self, d, lp, kT, rRange=(0,50), resolution=0.02, maxForce=100, max_potential=None, prefix="potentials/"):
BasePotential.__init__(self, 0, rRange, resolution, maxForce, max_potential, prefix)
self.type_ = "wlcbond"
self.d = d # separation
self.lp = lp # persistence length
self.kT = kT
def filename(self):
return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
self.d, self.lp)
def potential(self, dr):
nk = self.d / (2*self.lp)
q2 = (dr / self.d)**2
a1,a2 = 1, -7.0/(2*nk)
a3 = 3.0/32 - 3.0/(8*nk) - 6.0/(4*nk**2)
p0,p1,p2,p3,p4 = 13.0/32, 3.4719,2.5064,-1.2906,0.6482
a4 = (p0 + p1/(2*nk) + p2*(2*nk)**-2) / (1+p3/(2*nk)+p4*(2*nk)**-2)
with np.errstate(divide='ignore',invalid='ignore'):
u = self.kT * nk * ( a1/(1-q2) - a2*np.log(1-q2) + a3*q2 - 0.5*a4*q2*(q2-2) )
max_force = np.diff(u[q2<1][-2:]) / np.diff(dr).mean()
max_u = u[q2<1][-1]
max_dr = dr[q2<1][-2]
assert( max_force >= 0 )
u[q2>=1] = (dr[q2>=1]-max_dr)*max_force + max_u
return u
def __hash__(self):
assert(self.type_ != "None")
return hash((self.type_, self.d, self.lp, self.kT, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
def __eq__(self, other):
for a in ("type_", "d", "lp", "kT" "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
if self.__dict__[a] != other.__dict__[a]:
return False
return True
class WLCSKAngle(BasePotential):
## https://aip.scitation.org/doi/full/10.1063/1.4968020
def __init__(self, d, lp, kT, rRange=(0,181), resolution=0.5, maxForce=None, max_potential=None, prefix="potentials/"):
BasePotential.__init__(self, 180, rRange, resolution, maxForce, max_potential, prefix)
self.type_ = "wlcangle"
self.d = d # separation
self.lp = lp # persistence length
self.kT = kT
def filename(self):
return "%s%s-%.3f-%.3f.dat" % (self.prefix, self.type_,
self.d, self.lp)
def potential(self,dr):
nk = self.d / (2*self.lp)
p1,p2,p3,p4 = -1.237, 0.8105, -1.0243, 0.4595
C = (1 + p1*(2*nk) + p2*(2*nk)**2) / (2*nk+p3*(2*nk)**2+p4*(2*nk)**3)
u = self.kT * C * (1-np.cos(dr * np.pi / 180))
return u
def __hash__(self):
assert(self.type_ != "None")
return hash((self.type_, self.d, self.lp, self.kT, self.rRange, self.resolution, self.maxForce, self.max_potential, self.prefix, self.periodic))
def __eq__(self, other):
for a in ("type_", "d", "lp", "kT" "rRange", "resolution", "maxForce", "max_potential", "prefix", "periodic"):
if self.__dict__[a] != other.__dict__[a]:
return False
return True
......@@ -1637,6 +1637,34 @@ class SegmentModel(ArbdModel):
while d < -180: d+=360
return self._get_potential("dihedral", kSpring, d, max_potential)
def _get_wlc_sk_bond_potential(self, d):
## https://aip.scitation.org/doi/full/10.1063/1.4968020
type_='wclsk-bond'
# lp = 2*5 # TODO place these parameters in a logical spot
lp = 15 # https://journals.aps.org/pre/abstract/10.1103/PhysRevE.86.021901
# https://www.pnas.org/content/pnas/109/3/799.full.pdf
lp = 12 # selected semi-empirically to make ssDNA force extension match
key = (type_, d, lp)
assert( d > 0.2 )
if key not in self._bonded_potential:
kT = self.temperature * 0.0019872065 # kcal/mol
self._bonded_potential[key] = WLCSKBond( d, lp, kT )
return self._bonded_potential[key]
def _get_wlc_sk_angle_potential(self, d):
## https://aip.scitation.org/doi/full/10.1063/1.4968020
type_='wclsk-angle'
## TODO move
lp = 15 # https://journals.aps.org/pre/abstract/10.1103/PhysRevE.86.021901
# https://www.pnas.org/content/pnas/109/3/799.full.pdf
lp = 12 # selected semi-empirically to make ssDNA force extension match
key = (type_, d, lp)
assert( d > 0.2 )
if key not in self._bonded_potential:
kT = self.temperature * 0.0019872065 # kcal/mol
self._bonded_potential[key] = WLCSKAngle( d, lp, kT )
return self._bonded_potential[key]
def _getParent(self, *beads ):
if np.all( [b1.parent == b2.parent
......@@ -1647,7 +1675,7 @@ class SegmentModel(ArbdModel):
def _get_twist_spring_constant(self, sep):
""" sep in nt """
kT = 0.58622522 # kcal/mol
kT = self.temperature * 0.0019872065 # kcal/mol
twist_persistence_length = 90 # set semi-arbitrarily as there is a large spread in literature
## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
## Assume A is small
......@@ -2217,18 +2245,16 @@ class SegmentModel(ArbdModel):
return conversion*elastic_modulus_times_area/d
def k_ssdna_bond(d):
conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
## TODO: get better numbers our ssDNA model
elastic_modulus_times_area = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
# conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
# ## TODO: get better numbers our ssDNA model
# elastic_modulus_times_area = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
return conversion*elastic_modulus_times_area/d
def k_dsdna_angle(sep):
Lp = get_effective_dsDNA_Lp(sep)
return angle_spring_from_lp(sep,Lp)
def k_ssdna_angle(sep):
return angle_spring_from_lp(sep,2)
def k_xover_angle(sep):
return 0.25 * angle_spring_from_lp(sep,147)
......@@ -2249,15 +2275,15 @@ class SegmentModel(ArbdModel):
## TODO: could be sligtly smarter about sep
sep = 0.5*(b1.num_nt+b2.num_nt)
if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
is_dsdna = b1.type_.name[0] == "D" and b2.type_.name[0] == "D"
if is_dsdna:
d = 3.4*sep
k = k_dsdna_bond(d)
else:
d = 5*sep
# d = 6.5*sep # https://journals.aps.org/pre/abstract/10.1103/PhysRevE.86.021901
d = 6.4*sep # https://journals.aps.org/pre/abstract/10.1103/PhysRevE.86.021901
if b1.type_.name[0] != b2.type_.name[0]:
""" Add a small extra distance to junction """
d += 3
k = k_ssdna_bond(d)
if b1 not in dists:
dists[b1] = dict()
......@@ -2271,7 +2297,11 @@ class SegmentModel(ArbdModel):
if b1 is b2: continue
# dists[[b1,b2]] = dists[[b2,b1]] = sep
bond = self.get_bond_potential(k,d)
if is_dsdna:
k = k_dsdna_bond(d)
bond = self.get_bond_potential(k,d)
else:
bond = self._get_wlc_sk_bond_potential(d)
parent.add_bond( b1, b2, bond, exclude=True )
# for s in self.segments:
......@@ -2309,7 +2339,6 @@ class SegmentModel(ArbdModel):
sep = 0.5*(0.5*b1.num_nt+b2.num_nt+0.5*b3.num_nt)
parent = self._getParent(b1,b2,b3)
kT = 0.58622522 # kcal/mol
if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
k = k_dsdna_angle(sep)
......@@ -2318,16 +2347,16 @@ class SegmentModel(ArbdModel):
k *= 0.75 # reduce because orientation beads impose similar springs
dihed = self.get_dihedral_potential(k_dihed,180)
parent.add_dihedral(b1,b2,b2.orientation_bead,b3, dihed)
angle = self.get_angle_potential(k,180)
elif b1.type_.name[0] == "S" and b2.type_.name[0] == "S" and b3.type_.name[0] == "S":
## TODO: get correct number from ssDNA model
k = k_ssdna_angle(sep)
# angle = self._get_wlc_sk_angle_potential(sep*6.5) # TODO move 6.5 somewhere sensible
angle = self._get_wlc_sk_angle_potential(sep*6.4) # TODO move 6.4 somewhere sensible
else:
## Considered as a sscrossover below
continue
angle = self.get_angle_potential(k,180)
parent.add_angle( b1, b2, b3, angle )
""" Add intrahelical exclusions """
......@@ -2643,6 +2672,7 @@ class SegmentModel(ArbdModel):
Lp = s1.twist_persistence_length/0.34 # set semi-arbitrarily as there is a large spread in literature
def get_spring(sep):
kT = self.temperature * 0.0019872065 # kcal/mol
fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - np.exp(-sep/Lp)
k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
return k[0][0] * 2*kT*0.00030461742
......
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