From 73c133699eaf5e32fbb9c9e7930c8327d0f611bf Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Sun, 2 Dec 2018 10:15:30 -0600
Subject: [PATCH] Slightly reorganized simulation protocol; added options to
mrdna script and simulate.py for using twist in the coarsest step
---
bin/mrdna | 23 ++++++--
mrdna/simulate.py | 141 ++++++++++++++++++++++++----------------------
2 files changed, 94 insertions(+), 70 deletions(-)
diff --git a/bin/mrdna b/bin/mrdna
index 4ff5a2e..789661c 100755
--- a/bin/mrdna
+++ b/bin/mrdna
@@ -1,14 +1,21 @@
#! /usr/bin/env python
# -*- coding: utf-8 -*-
-
+import sys
+if sys.version_info[0] != 3:
+ raise Exception("This script must be run with python3")
+
import argparse
import re
import pathlib
from mrdna import __version__ as __version__
from mrdna.simulate import multiresolution_simulation as simulate
+"""Easy multiresolution simulations of DNA nanotechnology objects using ARBD"""
+
parser = argparse.ArgumentParser(prog="mrdna",
- description="""Easy multiresolution simulations of DNA nanotechnology objects using ARBD""")
+ description=__doc__,
+ formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+)
parser.add_argument('-o','--output-prefix', type=str, default=None,
help="Name for your job's output")
@@ -18,6 +25,10 @@ parser.add_argument('-g','--gpu', type=int, default=0,
help='GPU for simulation; check nvidia-smi for availability')
parser.add_argument('--output-period', type=float, default=1e4,
help='Simulation steps between DCD frames')
+# parser.add_argument('--minimization-steps', type=float, default=0,
+# 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('--coarse-steps', type=float, default=1e7,
help='Simulation steps for coarse model (200 fs/step)')
parser.add_argument('--fine-steps', type=float, default=1e7,
@@ -97,17 +108,20 @@ def main():
job_id = "job-" + prefix,
directory = args.directory,
gpu = args.gpu,
+ minimization_output_period = int(args.output_period),
+ coarse_local_twist = args.coarse_local_twist,
coarse_output_period = int(args.output_period),
fine_output_period = int(args.output_period),
+ minimization_steps = 0, # int(args.minimization_steps),
coarse_steps = int(args.coarse_steps),
fine_steps = int(args.fine_steps),
backbone_scale = args.backbone_scale
)
-
simulate( **run_args )
if __name__ == '__main__':
+
if args.debug:
class Restart(Exception):
@@ -132,7 +146,8 @@ if __name__ == '__main__':
print("\t" + " ".join(args))
except SystemExit:
# In most cases SystemExit does not warrant a post-mortem session.
- print("The program exited via sys.exit(). Exit status:", end=' ')
+ # print("The program exited via sys.exit(). Exit status:", end=' ')
+ print("The program exited via sys.exit(). Exit status:")
print(sys.exc_info()[1])
except SyntaxError:
traceback.print_exc()
diff --git a/mrdna/simulate.py b/mrdna/simulate.py
index 2e04647..c069b1d 100644
--- a/mrdna/simulate.py
+++ b/mrdna/simulate.py
@@ -9,24 +9,34 @@ namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multi
## TODO: implement replicas, initial conditions specified through some restart, and a custom simulation schedule
+
+
def multiresolution_simulation( model, output_name,
job_id=None,
gpu = 0,
directory=None,
+ minimization_steps = 0,
coarse_steps = 5e7,
fine_steps = 5e7,
+ minimization_output_period = 1e5,
coarse_output_period = 1e5,
fine_output_period = 1e5,
+ coarse_local_twist = False,
remove_long_bonds=False,
backbone_scale=1.0,
):
## Round steps up to nearest multiple of output_period, plus 1
+ if minimization_steps > 0:
+ raise NotImplementedError
+ minimization_steps = ((minimization_steps+minimization_output_period-1)//minimization_output_period)*minimization_output_period+1
coarse_steps = ((coarse_steps+coarse_output_period-1)//coarse_output_period)*coarse_output_period+1
fine_steps = ((fine_steps+fine_output_period-1)//fine_output_period)
if fine_steps == 1: fine_steps += 1
fine_steps = fine_steps*fine_output_period+1
+
+
ret = None
d_orig = os.getcwd()
try:
@@ -38,89 +48,88 @@ def multiresolution_simulation( model, output_name,
os.chdir(directory)
output_directory = "output"
+
+ simulation_step = 0
+
+ def run_step(num_steps, simargs, do_coordinate_average=False):
+ """ Function runs one step of the simulation protocol process """
+ nonlocal simulation_step
+
+ output_prefix = "{}-{}".format(output_name, simulation_step)
+ full_output_prefix = "{}/{}".format(output_directory,output_prefix)
+
+ model.simulate( output_name = output_prefix, num_steps = num_steps, **simargs )
+
+ if do_coordinate_average:
+ try:
+ coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
+ except:
+ coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
+
+ else:
+ try:
+ coordinates = readArbdCoords('%s.restart' % full_output_prefix)
+ except:
+ coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
+
+ model.update_splines(coordinates)
+ simulation_step += 1
+
+ coarse_model_args = dict( max_basepairs_per_bead = 5,
+ max_nucleotides_per_bead = 5,
+ local_twist = coarse_local_twist,
+ escapable_twist = True )
# TODO: add the output directory in to the readArbdCoords functions or make it an attribute of the model object
- """ Coarse simulation """
- simargs = dict(timestep=200e-6, output_period=coarse_output_period, gpu=gpu)
+ """ 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 )
- model._clear_beads()
- model._generate_bead_model( 5, 5, local_twist=False )
+ """ Minimization """
+ if minimization_steps > 0:
+ simargs = dict(timestep=200e-6, output_period=minimization_output_period, gpu=gpu)
+ run_step(num_steps=minimization_steps, simargs=simargs)
+ model.clear_beads()
+ model.generate_bead_model( **coarse_model_args )
+ """ Coarse simulation """
+ simargs = dict(timestep=200e-6, output_period=coarse_output_period, gpu=gpu)
if remove_long_bonds:
- output_prefix = "%s-0" % output_name
- full_output_prefix = "%s/%s" % (output_directory,output_prefix)
- ## TODO Remove long bonds
- model.simulate( output_name = output_prefix, numSteps=0.05*coarse_steps, **simargs )
- try:
- coordinates = readArbdCoords('%s.restart' % full_output_prefix)
- except:
- coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
- output_prefix = "%s-1" % output_name
- model.update_splines(coordinates)
- model.simulate( output_prefix = output_prefix, num_steps=0.95*coarse_steps, **simargs )
+ run_step(num_steps=0.05*coarse_steps, simargs=simargs)
+ model.clear_beads()
+ model.generate_bead_model( **coarse_model_args )
+ run_step(num_steps=0.95*coarse_steps, simargs=simargs)
else:
- output_prefix = "%s-1" % output_name
- full_output_prefix = "%s/%s" % (output_directory,output_prefix)
- model.simulate( output_name = output_prefix, num_steps=coarse_steps, **simargs )
- try:
- coordinates = readArbdCoords('%s.restart' % full_output_prefix)
- except:
- coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
-
+ run_step(num_steps=coarse_steps, simargs=simargs)
""" Fine simulation """
- output_prefix = "%s-2" % output_name
- full_output_prefix = "%s/%s" % (output_directory,output_prefix)
- simargs['timestep'] = 50e-6
- simargs['output_period'] = fine_output_period
- model.update_splines(coordinates)
- model._clear_beads()
- model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
- model.simulate( output_name = output_prefix, num_steps=fine_steps, **simargs )
- if fine_steps <= 2*fine_output_period+1:
- try:
- coordinates = readArbdCoords('%s.restart' % full_output_prefix)
- except:
- coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
- else:
- try:
- coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
- except:
- coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
-
+ simargs = dict(timestep=50e-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))
""" Freeze twist """
- output_prefix = "%s-3" % output_name
- full_output_prefix = "%s/%s" % (output_directory,output_prefix)
- model.update_splines(coordinates)
- model._clear_beads()
- model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=False )
- model.simulate( output_name = output_prefix, num_steps=fine_steps, **simargs )
- if fine_steps <= 2*fine_output_period+1:
- try:
- coordinates = readArbdCoords('%s.restart' % full_output_prefix)
- except:
- coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
- else:
- try:
- coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
- except:
- coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
-
+ model.clear_beads()
+ model.generate_bead_model( 1, 1, local_twist=True, escapable_twist=False )
+ run_step(num_steps=fine_steps, simargs=simargs, do_coordinate_average=(fine_steps > 2*fine_output_period+1))
""" Atomic simulation """
- output_prefix = "%s-4" % output_name
- full_output_prefix = "%s/%s" % (output_directory,output_prefix)
- model.update_splines(coordinates)
- model._clear_beads()
- model._generate_atomic_model(scale=backbone_scale)
+ model.clear_beads()
+ model.generate_atomic_model(scale=backbone_scale)
+
+ output_prefix = "{}-{}".format(output_name, simulation_step)
+ # full_output_prefix = "%s/%s" % (output_directory,output_prefix)
+
+ model.write_atomic_ENM( output_prefix )
+ model.atomic_simulate( output_name = output_prefix )
+
try:
shutil.copytree( get_resource_path("charmm36.nbfix"), "charmm36.nbfix" )
except FileExistsError:
...
- model.write_atomic_ENM( output_prefix )
- model.atomic_simulate( output_name = output_prefix )
ret = directory
--
GitLab