Added option for assigning hj equilibrium angle

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 )]

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