diff --git a/mrdna/model/nonbonded.py b/mrdna/model/nonbonded.py
index 3dc4a6cbbd62930619685c551bb0ccb154e84634..53ce19fe4c16558acd284d0ce5d70bc0dee70124 100644
--- a/mrdna/model/nonbonded.py
+++ b/mrdna/model/nonbonded.py
@@ -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
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 8c8029337c843c85ded2f49c3467ed36a0c66de1..be15a4602c22ac178d2f15de953aa4c237a94f6d 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -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