diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py index 40ab067a42792d1fa8fb0dcb6182b2284346c55d..be162ec19fc2ec59f2604731d15d316da9942b4b 100644 --- a/mrdna/segmentmodel.py +++ b/mrdna/segmentmodel.py @@ -2072,6 +2072,22 @@ class SegmentModel(ArbdModel): # print("Helix {}: bps {}, beads {}, separation {}".format(s.name, s.num_nt, beadsum, sepsum)) """ Add intrahelical angle potentials """ + def get_effective_dsDNA_Lp(sep): + + """ The persistence length of our model was found to be a + little off (probably due to NB interactions). This + attempts to compensate """ + + ## For 1 bp, Lp=559, for 25 Lp = 524 + beads_per_bp = sep/2 + Lp0 = 147 + # return 0.93457944*Lp0 ;# factor1 + return 0.97*Lp0 ;# factor2 + # 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 @@ -2080,7 +2096,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": - k = angle_spring_from_lp(sep,147) + Lp = get_effective_dsDNA_Lp(sep) + k = angle_spring_from_lp(sep,Lp) if local_twist: k_dihed = 0.25*k k *= 0.75 # reduce because orientation beads impose similar springs @@ -2134,7 +2151,8 @@ class SegmentModel(ArbdModel): parent = self._getParent( b1, b2 ) """ Add heuristic 90 degree potential to keep orientation bead orthogonal """ - k = 0.5*angle_spring_from_lp(sep,147) + Lp = get_effective_dsDNA_Lp(sep) + k = 0.5*angle_spring_from_lp(sep,Lp) pot = self.get_angle_potential(k,90) parent.add_angle(o1,b1,b2, pot) parent.add_angle(b1,b2,o2, pot) diff --git a/mrdna/version.py b/mrdna/version.py index 80ef677eddc2299361ad673ce725d573702ee2ad..b12c13951fab4d793851ebc8994e02d6b4937923 100644 --- a/mrdna/version.py +++ b/mrdna/version.py @@ -133,7 +133,7 @@ def get_version(abbrev=7): if __name__ == "__main__": - print( get_git_version() ) + print( get_version() ) def read_version_file(filename): with open(filename) as fh: