Commit b349e427 authored by cmaffeo2's avatar cmaffeo2
Browse files

Prepared code structure in potential section so that dsDNA and ssDNA are...

Prepared code structure in potential section so that dsDNA and ssDNA are treated using separate defs
parent 3006dcba
......@@ -2210,6 +2210,29 @@ class SegmentModel(ArbdModel):
self._countParticleTypes() # probably not needed here
self._updateParticleOrder()
def k_dsdna_bond(d):
conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
elastic_modulus_times_area = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
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
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)
""" Add intrahelical bond potentials """
if self.DEBUG: print("Adding intrahelical bond potentials")
dists = dict() # intrahelical distances built for later use
......@@ -2226,22 +2249,16 @@ class SegmentModel(ArbdModel):
## TODO: could be sligtly smarter about sep
sep = 0.5*(b1.num_nt+b2.num_nt)
conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
elastic_modulus_times_area = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
d = 3.4*sep
k = conversion*elastic_modulus_times_area/d
k = k_dsdna_bond(d)
else:
## TODO: get better numbers our ssDNA model
elastic_modulus_times_area = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
d = 5*sep
if b1.type_.name[0] != b2.type_.name[0]:
""" Add a small extra distance to junction """
d += 3
k = k_ssdna_bond(d)
k = conversion*elastic_modulus_times_area/d
# print(sep,d,k)
if b1 not in dists:
dists[b1] = dict()
if b2 not in dists:
......@@ -2286,8 +2303,6 @@ class SegmentModel(ArbdModel):
# factor = bead_per_bp * (0.954-0.8944
# return Lp0 * bead_per_bp
empirical_compensation_factor = max_basepairs_per_bead
if self.DEBUG: print("Adding intrahelical angle potentials")
for b1,b2,b3 in self._get_intrahelical_angle_beads():
## TODO: could be slightly smarter about sep
......@@ -2296,8 +2311,8 @@ class SegmentModel(ArbdModel):
kT = 0.58622522 # kcal/mol
if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
Lp = get_effective_dsDNA_Lp(sep)
k = angle_spring_from_lp(sep,Lp)
k = k_dsdna_angle(sep)
if local_twist:
k_dihed = 0.25*k
k *= 0.75 # reduce because orientation beads impose similar springs
......@@ -2307,7 +2322,7 @@ class SegmentModel(ArbdModel):
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 = angle_spring_from_lp(sep,2)
k = k_ssdna_angle(sep)
else:
## Considered as a sscrossover below
continue
......@@ -2371,8 +2386,7 @@ class SegmentModel(ArbdModel):
parent = self._getParent( b1, b2 )
""" Add heuristic 90 degree potential to keep orientation bead orthogonal """
Lp = get_effective_dsDNA_Lp(sep)
k = 0.5*angle_spring_from_lp(sep,Lp)
k = 0.5*k_dsdna_angle(sep)
pot = self.get_angle_potential(k,90)
parent.add_angle(o1,b1,b2, pot)
parent.add_angle(b1,b2,o2, pot)
......@@ -2391,11 +2405,6 @@ class SegmentModel(ArbdModel):
pot = self.get_dihedral_potential(k,angle)
parent.add_dihedral(o1,b1,b2,o2, pot)
def k_angle(sep):
return 0.5*angle_spring_from_lp(sep,147)
def k_xover_angle(sep):
return 0.5 * k_angle(sep)
def add_local_crossover_strand_orientation_potential(b1,b2, b1_on_fwd_strand):
......
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