From a72a792fdae5234c2e570ae2ce9e666efdf94629 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Mon, 9 Mar 2020 13:01:37 -0500
Subject: [PATCH] Added option for assigning hj equilibrium angle
---
bin/mrdna | 5 ++++-
mrdna/segmentmodel.py | 18 ++++++++++++++----
2 files changed, 18 insertions(+), 5 deletions(-)
diff --git a/bin/mrdna b/bin/mrdna
index e43df23..a604fd7 100755
--- a/bin/mrdna
+++ b/bin/mrdna
@@ -73,6 +73,9 @@ parser.add_argument('--draw-cylinders', action='store_true',
parser.add_argument('--draw-tubes', action='store_true',
help='Whether or not to draw the tubes')
+parser.add_argument('--hj-equilibrium-angle', type=float, default=0.0,
+ help='Specify an equilibrium dihedral angle for the Holliday junction angle potential; the default value works well for origami')
+
parser.add_argument('--version', action='version',
version='%(prog)s {}'.format(__version__),
help='Print the version of mrdna')
@@ -113,7 +116,7 @@ def main():
from mrdna.readers import read_list as read_model
else:
raise Exception("Unrecognized input file '{}'".format(infile))
- model = read_model( str(infile), debye_length=args.debye_length, temperature=args.temperature )
+ model = read_model( str(infile), debye_length=args.debye_length, temperature=args.temperature, hj_equilibrium_angle = args.hj_equilibrium_angle )
if args.dimensions is None:
dim = [max(x,1000) for x in model.dimensions_from_structure( padding_factor=2.5 )]
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 0e52a44..99315a5 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -1600,6 +1600,7 @@ class SegmentModel(ArbdModel):
nonbondedResolution=0, DEBUG=0,
integrator='Brown',
debye_length = None,
+ hj_equilibrium_angle = 0,
extra_bd_file_lines = "",
):
self.DEBUG = DEBUG
@@ -1617,6 +1618,10 @@ class SegmentModel(ArbdModel):
# self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
self.children = self.segments = segments
+ if (hj_equilibrium_angle > 180) or (hj_equilibrium_angle < -180):
+ raise ValueError("Holliday junction equilibrium dihedral angle must be in the range from -180 to 180 degrees!")
+ self.hj_equilibrium_angle = hj_equilibrium_angle
+
self._generate_bead_callbacks = []
self._bonded_potential = dict() # cache for bonded potentials
@@ -2712,18 +2717,23 @@ class SegmentModel(ArbdModel):
a = None
if u1 is not None and u2 is not None:
- t0 = 0
+ t0 = self.hj_equilibrium_angle
a,b,c,d = (u1,b1,b2,u2)
elif d1 is not None and d2 is not None:
- t0 = 0
+ t0 = self.hj_equilibrium_angle
a,b,c,d = (d1,b1,b2,d2 )
elif d1 is not None and u2 is not None:
- t0 = 180
+ t0 = self.hj_equilibrium_angle-180
a,b,c,d = (d1,b1,b2,u2)
elif u1 is not None and d2 is not None:
- t0 = 180
+ t0 = self.hj_equilibrium_angle-180
a,b,c,d = (u1,b1,b2,d2)
+ if t0 > 180:
+ t0 = t0-360
+ elif t0 < -180:
+ t0 = t0+360
+
## TODO?: Check length-dependence of this potential
if a is not None:
k_intrinsic = 0.00086
--
GitLab