diff --git a/bin/mrdna b/bin/mrdna
index 2430aadeb170b2e17ffeab2c646e60f2cc51cd99..888142c2622838ff5776567e8c8bd297286e43c7 100755
--- a/bin/mrdna
+++ b/bin/mrdna
@@ -32,6 +32,8 @@ parser.add_argument('--output-period', type=float, default=1e4,
 #                     help='Simulation steps between DCD frames')
 parser.add_argument('--coarse-local-twist', action='store_true',
                     help='Use local twist for coarsest representation?')
+parser.add_argument('--fix-linking-number', action='store_true',
+                    help='Fix the linking number so that it cannot change once local twist is introduced?')
 parser.add_argument('--coarse-steps', type=float, default=1e7,
                     help='Simulation steps for coarse model (200 fs/step)')
 parser.add_argument('--fine-steps', type=float, default=1e7,
@@ -81,7 +83,6 @@ def main():
         from mrdna.readers import read_atomic_pdb as read_model
     else:
         raise Exception("Unrecognized input file '{}'".format(infile))
-                    
     model = read_model( str(infile), debye_length=args.debye_length )
 
     if args.output_prefix is not None:
@@ -120,6 +121,7 @@ def main():
             gpu = args.gpu,
             minimization_output_period = int(args.output_period),
             coarse_local_twist = args.coarse_local_twist,
+            fix_linking_number = args.fix_linking_number,
             coarse_output_period = int(args.output_period),
             fine_output_period = int(args.output_period),
             minimization_steps = 0, # int(args.minimization_steps),
@@ -148,7 +150,7 @@ if __name__ == '__main__':
         while True:
             try:
                 # pdb._runscript(mainpyfile)
-                self._user_requested_quit = False
+                pdb._user_requested_quit = False
                 pdb.runcall(main)
                 if pdb._user_requested_quit:
                     break
diff --git a/mrdna/simulate.py b/mrdna/simulate.py
index b7c034d6418c2b93dbc01e6d318d25a359330197..223d9eea5716686300a2984a8edf9f1a023af956 100644
--- a/mrdna/simulate.py
+++ b/mrdna/simulate.py
@@ -58,6 +58,7 @@ def multiresolution_simulation( model, output_name,
                                 coarse_output_period = 1e5,
                                 fine_output_period = 1e5,
                                 coarse_local_twist = False,
+                                fix_linking_number = False,
                                 remove_long_bonds=False,
                                 backbone_scale=1.0,
                                 oxdna_steps = None,
@@ -114,7 +115,7 @@ def multiresolution_simulation( model, output_name,
         coarse_model_args = dict( max_basepairs_per_bead = 5,
                                   max_nucleotides_per_bead = 5,
                                   local_twist = coarse_local_twist,
-                                  escapable_twist = True )
+                                  escapable_twist = fix_linking_number )
         
         # TODO: add the output directory in to the readArbdCoords functions or make it an attribute of the model object 
 
@@ -132,7 +133,7 @@ def multiresolution_simulation( model, output_name,
             model.generate_bead_model( **coarse_model_args )
 
         """ Coarse simulation """
-        simargs = dict(timestep=200e-6, output_period=coarse_output_period, gpu=gpu)
+        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)
             model.clear_beads()
@@ -143,9 +144,10 @@ def multiresolution_simulation( model, output_name,
 
         """ Fine simulation """ 
         simargs = dict(timestep=40e-6, output_period=fine_output_period, gpu=gpu)
-        model.clear_beads()
-        model.generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
-        run_step(num_steps=fine_steps, simargs=simargs, do_coordinate_average=(fine_steps > 2*fine_output_period+1))
+        if not fix_linking_number:
+            model.clear_beads()
+            model.generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
+            run_step(num_steps=fine_steps, simargs=simargs, do_coordinate_average=(fine_steps > 2*fine_output_period+1))
 
         """ Freeze twist """
         model.clear_beads()