Added option to specify custom persistence length to mrdna script for first...

Added option to specify custom persistence length to mrdna script for first 25% of coarse simulation

bin/mrdna

@@ -57,7 +57,9 @@ parser.add_argument('--oxdna-steps', type=float, default=None,
 parser.add_argument('--oxdna-output-period', type=float, default=1e4,
                     help='Simulation steps between oxDNA configuration and energy output')
 parser.add_argument('--coarse-bond-cutoff', type=float, default = 0,
-                    help='Ignore bonds beyond this cutoff during first step of simulation; a value of 0 implies bonds are not ignored')
+                    help='Ignore bonds beyond this cutoff during 25% of coarse simulation, adding a step to the simulation. Default value of zero implies bonds are not ignored')
+parser.add_argument('--coarse-persistence-length', type=float, default = None,
+                    help='Specify a custom persistence length (in nm) for the first 25% of coarse simulation (after coarse-bond-cutoff simulation if specified), adding a step to the simulation')
 
 parser.add_argument('--crossover-to-intrahelical-cutoff', type=float, default=-1,
                     help='Set the distance (in Angstroms) beyond which crossovers between helix ends are converted to intrahelical connections; a negative value means no crossovers will be converted')
@@ -181,6 +183,7 @@ def main():
             coarse_local_twist = args.coarse_local_twist,
             fix_linking_number = args.fix_linking_number,
             bond_cutoff = args.coarse_bond_cutoff,
+            coarse_persistence_length = args.coarse_persistence_length,
             coarse_output_period = int(args.output_period),
             fine_output_period = int(args.output_period),
             minimization_steps = 0, # int(args.minimization_steps),

mrdna/segmentmodel.py

@@ -1304,6 +1304,7 @@ class DoubleStrandedSegment(Segment):
                  num_turns = None,
                  start_orientation = None,
                  twist_persistence_length = 90,
+                 persistence_length = 50,
                  **kwargs):
         
         self.helical_rise = 10.44
@@ -1324,7 +1325,8 @@ class DoubleStrandedSegment(Segment):
             start_orientation = np.eye(3) # np.array(((1,0,0),(0,1,0),(0,0,1)))
         self.start_orientation = start_orientation
         self.twist_persistence_length = twist_persistence_length
-
+        self.persistence_length = persistence_length
+        
         self.nicks = []
 
         self.start = self.start5 = Location( self, address=0, type_= "end5" )
@@ -2749,12 +2751,14 @@ class SegmentModel(ArbdModel):
             
             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_dsdna_angle(sep, Lp=None):
+            Lp_eff = get_effective_dsDNA_Lp(sep, Lp)
+            Lp_eff = Lp_eff/0.34 # convert from nm to bp
+            return angle_spring_from_lp(sep,Lp_eff)
 
-        def k_xover_angle(sep):
-            return 0.25 * angle_spring_from_lp(sep,147)
+        def k_xover_angle(sep, Lp=50):
+            Lp = Lp/0.34 # convert from nm to bp
+            return 0.25 * angle_spring_from_lp(sep,Lp)
 
 
         """ Add intrahelical bond potentials """
@@ -2821,15 +2825,13 @@ 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):
-
+        def get_effective_dsDNA_Lp(sep, Lp0=50):
             """ 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
@@ -2843,9 +2845,8 @@ class SegmentModel(ArbdModel):
             sep = dists[b2][b1] + dists[b2][b3]
 
 
-            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
-                
-                k = k_dsdna_angle(sep)
+            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":                
+                k = k_dsdna_angle(sep, seg.persistence_length)
                 if local_twist:
                     k_dihed = 0.25*k
                     k *= 0.75    # reduce because orientation beads impose similar springs
@@ -2933,7 +2934,7 @@ class SegmentModel(ArbdModel):
                     parent = self._getParent( b1, b2 )
 
                     """ Add heuristic 90 degree potential to keep orientation bead orthogonal """
-                    k = 0.25*k_dsdna_angle(sep)
+                    k = 0.25*k_dsdna_angle(sep, parent.persistence_length)
                     pot = self.get_angle_potential(k,90)
                     parent.add_angle(o1,b1,b2, pot)
                     parent.add_angle(b1,b2,o2, pot)

mrdna/simulate.py

@@ -2,7 +2,7 @@ import os
 import tempfile
 from .arbdmodel.coords import readArbdCoords, readAvgArbdCoords
 import shutil
-from . import get_resource_path
+from . import get_resource_path, DoubleStrandedSegment
 
 import numpy as np
 
@@ -85,6 +85,7 @@ def multiresolution_simulation( model, output_name,
                                 coarse_local_twist = False,
                                 fix_linking_number = False,
                                 bond_cutoff = 0,
+                                coarse_persistence_length = None,
                                 backbone_scale=1.0,
                                 oxdna_steps = None,
                                 oxdna_output_period = None,
@@ -166,13 +167,32 @@ def multiresolution_simulation( model, output_name,
         """ Coarse simulation """
         simargs = dict(timestep=100e-6 if coarse_local_twist else 200e-6, output_period=coarse_output_period,
                        write_pqr=write_pqr, gpu=gpu, arbd=arbd)
+
+        remaining_coarse_steps = coarse_steps
         if bond_cutoff > 0 and minimization_steps <= 0:
             run_step(num_steps=0.25*coarse_steps, simargs=simargs)
+            remaining_coarse_steps -= 0.25*coarse_steps
             model.clear_beads()
             model.generate_bead_model( **coarse_model_args )
-            run_step(num_steps=0.75*coarse_steps, simargs=simargs)
-        else:
-            run_step(num_steps=coarse_steps, simargs=simargs)
+        if coarse_persistence_length is not None:
+            Lps = []
+            segs = [seg for seg in model.segments if isinstance(seg,DoubleStrandedSegment)]
+            for seg in segs:
+                Lps.append(seg.persistence_length)
+                seg.persistence_length = coarse_persistence_length
+            model.clear_beads()
+            model.generate_bead_model( **coarse_model_args )
+            run_step(num_steps=0.25*coarse_steps, simargs=simargs)            
+            remaining_coarse_steps -= 0.25*coarse_steps
+
+            ## Restore original persistence length
+            for Lp,seg in zip(Lps,segs):
+                seg.persistence_length = Lp
+
+            model.clear_beads()
+            model.generate_bead_model( **coarse_model_args )
+        
+        run_step(num_steps=remaining_coarse_steps, simargs=simargs)
 
         """ Fine simulation """ 
         simargs = dict(timestep=40e-6, output_period=fine_output_period, write_pqr=write_pqr, gpu=gpu, arbd=arbd)