Commit 73c13369 authored by cmaffeo2's avatar cmaffeo2
Browse files

Slightly reorganized simulation protocol; added options to mrdna script and...

Slightly reorganized simulation protocol; added options to mrdna script and simulate.py for using twist in the coarsest step
parent 8dafa9f9
#! /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()
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment