Added easy commands for ignoring long bonds during coarse-simulation

bin/mrdna

@@ -53,7 +53,8 @@ parser.add_argument('--oxdna-steps', type=float, default=None,
                     help='Perform an oxDNA simulation instead of creating atomic model')
 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')
 
 parser.add_argument('--backbone-scale', type=float, default=1.0,
                     help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
@@ -155,6 +156,7 @@ def main():
             minimization_output_period = int(args.output_period),
             coarse_local_twist = args.coarse_local_twist,
             fix_linking_number = args.fix_linking_number,
+            bond_cutoff = args.coarse_bond_cutoff,
             coarse_output_period = int(args.output_period),
             fine_output_period = int(args.output_period),
             minimization_steps = 0, # int(args.minimization_steps),

mrdna/simulate.py

@@ -4,8 +4,22 @@ from .arbdmodel.coords import readArbdCoords, readAvgArbdCoords
 import shutil
 from . import get_resource_path
 
+import numpy as np
+
 ## TODO: implement replicas, initial conditions specified through some restart, and a custom simulation schedule
 
+def _remove_bonds_beyond(model, cutoff):
+    bonds = model.get_bonds()
+    new_bonds = []
+    for bond in bonds:
+        n1,n2,b,ex = bond
+        r2 = np.sum( (n1.position - n2.position)**2 )
+        if r2 > cutoff**2:
+            for l in [model] + model.children:
+                try:
+                    l.bonds.remove(bond)
+                except:
+                    ...
 
 def minimize_and_simulate_oxdna( model, 
                                  output_name,
@@ -61,7 +75,7 @@ def multiresolution_simulation( model, output_name,
                                 fine_output_period = 1e5,
                                 coarse_local_twist = False,
                                 fix_linking_number = False,
-                                remove_long_bonds=False,
+                                bond_cutoff = 0,
                                 backbone_scale=1.0,
                                 oxdna_steps = None,
                                 oxdna_output_period = None,
@@ -123,10 +137,11 @@ def multiresolution_simulation( model, output_name,
         # TODO: add the output directory in to the readArbdCoords functions or make it an attribute of the model object 
 
         """ Prepare coarse model for minimization and coarse simulation """
-        if remove_long_bonds:
-            raise NotImplementedError("TODO: remove long bonds")
         model.clear_beads()
         model.generate_bead_model( **coarse_model_args )
+        
+        if bond_cutoff > 0:
+            _remove_bonds_beyond(model, bond_cutoff)
 
         """ Minimization """
         if minimization_steps > 0:
@@ -137,11 +152,11 @@ def multiresolution_simulation( model, output_name,
 
         """ Coarse simulation """
         simargs = dict(timestep=150e-6 if coarse_local_twist else 200e-6, output_period=coarse_output_period, gpu=gpu)
-        if remove_long_bonds:
-            run_step(num_steps=0.05*coarse_steps, simargs=simargs)
+        if bond_cutoff > 0 and minimization_steps <= 0:
+            run_step(num_steps=0.25*coarse_steps, simargs=simargs)
             model.clear_beads()
             model.generate_bead_model( **coarse_model_args )
-            run_step(num_steps=0.95*coarse_steps, simargs=simargs)
+            run_step(num_steps=0.75*coarse_steps, simargs=simargs)
         else:
             run_step(num_steps=coarse_steps, simargs=simargs)