Skip to content

GitLab

Commit 0ca4dc50 authored by cmaffeo2's avatar cmaffeo2
Browse files

Initial implementation of export to oxDNA

parent 5fa68395
Hide whitespace changes
Inline Side-by-side
bin/mrdna
View file @ 0ca4dc50
......@@ -33,6 +33,12 @@ 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,
help='Simulation steps for fine model (50 fs/step)')
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('--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')
parser.add_argument('--debug', action='store_true',
......@@ -115,7 +121,9 @@ def main():
minimization_steps = 0, # int(args.minimization_steps),
coarse_steps = int(args.coarse_steps),
fine_steps = int(args.fine_steps),
backbone_scale = args.backbone_scale
backbone_scale = args.backbone_scale,
oxdna_steps = args.oxdna_steps,
oxdna_output_period = args.oxdna_output_period
)
simulate( **run_args )
......
mrdna/__init__.py
View file @ 0ca4dc50
......@@ -10,6 +10,7 @@ def get_resource_path(relative_path):
from .segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from .simulate import multiresolution_simulation
from .model.dna_sequence import read_sequence_file
from .reporting import report
......
mrdna/model/arbdmodel.py
View file @ 0ca4dc50
# -*- coding: utf-8 -*-
import abc
from pathlib import Path
import numpy as np
from copy import copy, deepcopy
from inspect import ismethod
......@@ -1175,3 +1175,253 @@ run {num_steps:d}
self.writePsf( output_name + ".psf" )
self.write_namd_configuration( output_name, output_directory = output_directory )
os.sync()
""" oxDNA """
## https://dna.physics.ox.ac.uk/index.php/Documentation#Input_file
def _write_oxdna_configuration(self, filename):
_angstroms_to_oxdna = 0.11739845 ## units "AA" "8.518e-10 m"
with open(filename,'w') as fh:
fh.write("""t = {temperature}
b = {dimX} {dimY} {dimZ}
E = 0 0 0
""".format(temperature=self.temperature,
dimX = self.dimensions[0]*_angstroms_to_oxdna,
dimY = self.dimensions[1]*_angstroms_to_oxdna,
dimZ = self.dimensions[2]*_angstroms_to_oxdna))
base_vec = np.array((1,0,0)) # TODO
norm_vec = np.array((0,0,1)) # TODO
## Traverse 3'-to-5'
for nt in [nt for s in self.strands for nt in s.oxdna_nt[::-1]]:
data = dict()
# o = nt.collapsedOrientation()
o = nt.orientation
for k,v in zip('x y z'.split(),nt.collapsedPosition()):
data[k] = v * _angstroms_to_oxdna
for k,v in zip('x y z'.split(),o.dot(base_vec)):
data['b'+k] = v
for k,v in zip('x y z'.split(),o.dot(norm_vec)):
data['n'+k] = v
fh.write("{x} {y} {z} {bx} {by} {bz} {nx} {ny} {nz} 0 0 0 0 0 0\n".format(**data)
)
def _write_oxdna_topology(self,filename):
strands = [s for s in self.strands if s.num_nt > 0]
with open(filename,'w') as fh:
fh.write("{num_nts} {num_strands}\n".format(
num_nts = sum([s.num_nt for s in strands]),
num_strands = len(self.strands)))
idx = 0
sidx = 1
for strand in strands:
prev = idx+len(strand.num_nt) if strand.is_circular else -1
last = idx if strand.is_circular else -1
## Traverse 3'-to-5'
sequence = [seq for s in strand.strand_segments
for seq in s.get_sequence()][::-1]
for seq in sequence[:-1]:
## strand seq 3' 5'
fh.write("{} {} {} {}\n".format(sidx, seq, prev, idx+1))
prev = idx
idx += 1
seq = sequence[-1]
fh.write("{} {} {} {}\n".format(sidx, seq, prev, last))
idx += 1
sidx += 1
def _write_oxdna_input(self, filename,
topology,
conf_file,
trajectory_file,
last_conf_file,
log_file,
num_steps = 1e6,
interaction_type = 'DNA2',
salt_concentration = 0.5,
print_conf_interval = None,
print_energy_every = None,
timestep = 0.003,
sim_type = "MD",
backend = None,
backend_precision = None,
seed = None,
newtonian_steps = 103,
diff_coeff = 2.50,
thermostat = "john",
list_type = "cells",
ensemble = "nvt",
delta_translation = 0.22,
delta_rotation = 0.22,
verlet_skin = 0.5,
max_backbone_force = 5,
):
if seed is None:
import random
seed = random.randint(1,99999)
temperature = self.temperature
num_steps = int(num_steps)
newtonian_steps = int(newtonian_steps)
if print_conf_interval is None:
# print_conf_interval = min(num_steps//100)
print_conf_interval = 10000
print_conf_interval = int(print_conf_interval)
if print_energy_every is None:
print_energy_every = print_conf_interval
print_energy_every = int(print_energy_every)
if max_backbone_force is not None:
max_backbone_force = 'max_backbone_force = {}'.format(max_backbone_force)
if interaction_type in ('DNA2',):
salt_concentration = 'salt_concentration = {}'.format(salt_concentration)
else:
salt_concentration = ''
if backend is None:
backend = 'CUDA' if sim_type == 'MD' else 'CPU'
if backend_precision is None:
backend_precision = 'mixed' if backend == 'CUDA' else 'double'
if sim_type == 'VMMC':
ensemble = 'ensemble = {}'.format(ensemble)
delta_translation = 'delta_translation = {}'.format(delta_translation)
delta_rotation = 'delta_rotation = {}'.format(delta_rotation)
else:
ensemble = ''
delta_translation = ''
delta_rotation = ''
with open(filename,'w') as fh:
fh.write("""##############################
#### PROGRAM PARAMETERS ####
##############################
interaction_type = {interaction_type}
{salt_concentration}
sim_type = {sim_type}
backend = {backend}
backend_precision = {backend_precision}
#debug = 1
seed = {seed}
##############################
#### SIM PARAMETERS ####
##############################
steps = {num_steps:d}
newtonian_steps = {newtonian_steps:d}
diff_coeff = {diff_coeff}
thermostat = {thermostat}
list_type = {list_type}
{ensemble}
{delta_translation}
{delta_rotation}
T = {temperature:f} K
dt = {timestep}
verlet_skin = {verlet_skin}
{max_backbone_force}
##############################
#### INPUT / OUTPUT ####
##############################
topology = {topology}
conf_file = {conf_file}
lastconf_file = {last_conf_file}
trajectory_file = {trajectory_file}
refresh_vel = 1
log_file = {log_file}
no_stdout_energy = 1
restart_step_counter = 1
energy_file = {log_file}.energy.dat
print_conf_interval = {print_conf_interval}
print_energy_every = {print_energy_every}
time_scale = linear
external_forces = 0
""".format( **locals() ))
def simulate_oxdna(self, output_name, directory='.', output_directory='output', topology=None, configuration=None, oxDNA=None, **oxdna_args):
if output_directory == '': output_directory='.'
d_orig = os.getcwd()
if not os.path.exists(directory):
os.makedirs(directory)
os.chdir(directory)
try:
if oxDNA is None:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
fname = os.path.join(path, "oxDNA")
if os.path.isfile(fname) and os.access(fname, os.X_OK):
oxDNA = fname
break
if oxDNA is None: raise Exception("oxDNA was not found")
if not os.path.exists(oxDNA):
raise Exception("oxDNA was not found")
if not os.path.isfile(oxDNA):
raise Exception("oxDNA was not found")
if not os.access(oxDNA, os.X_OK):
raise Exception("oxDNA is not executable")
if not os.path.exists(output_directory):
os.makedirs(output_directory)
elif not os.path.isdir(output_directory):
raise Exception("output_directory '%s' is not a directory!" % output_directory)
if configuration is None:
configuration = "{}.conf".format(output_name)
self._write_oxdna_configuration(configuration)
# elif not Path(configuration).exists():
# raise Exception("Unable to find oxDNA configuration file '{}'.".format(configuration))
if topology is None:
topology = "{}-topology.dat".format(output_name)
self._write_oxdna_topology(topology)
elif not Path(topology).exists():
raise Exception("Unable to find oxDNA topology file '{}'.".format(topology))
last_conf_file = "{}/{}.last.conf".format(output_directory,output_name)
input_file = "{}-input".format(output_name)
self._write_oxdna_input(input_file,
topology = topology,
conf_file = configuration,
trajectory_file = "{}/{}.dat".format(output_directory,output_name),
last_conf_file = last_conf_file,
log_file="{}/{}.log".format(output_directory,output_name),
**oxdna_args)
os.sync()
## TODO: call oxdna
cmd = [oxDNA, input_file]
cmd = tuple(str(x) for x in cmd)
print("Running oxDNA with: %s" % " ".join(cmd))
process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
for line in process.stdout:
sys.stdout.write(line)
sys.stdout.flush()
return topology,last_conf_file
finally:
os.chdir(d_orig)
mrdna/segmentmodel.py
View file @ 0ca4dc50
......@@ -38,6 +38,8 @@ TODO:
- remove recursive calls
- document
- develop unit test suite
- refactor parts of Segment into an abstract_polymer class
- make each call generate_bead_model, generate_atomic_model, generate_oxdna_model return an object with only have a reference to original object
"""
class CircularDnaError(Exception):
pass
......@@ -704,6 +706,28 @@ class Segment(ConnectableElement, Group):
return atoms
def _generate_oxdna_nucleotide(self, contour_position, is_fwd, seq):
bp_center = self.contour_to_position(contour_position)
orientation = self.contour_to_orientation(contour_position)
DefaultOrientation = rotationAboutAxis([0,0,1], 90)
if is_fwd:
DefaultOrientation = rotationAboutAxis([1,0,0], 180).dot(DefaultOrientation)
o = orientation.dot(DefaultOrientation)
if isinstance(self, SingleStrandedSegment):
pos = bp_center
else:
pos = bp_center - 5*o.dot(np.array((1,0,0)))
nt = PointParticle("oxdna_nt", position= pos,
orientation = o)
nt.contour_position = contour_position
return nt
def add_location(self, nt, type_, on_fwd_strand=True):
## Create location if needed, add to segment
c = self.nt_pos_to_contour(nt)
......@@ -1340,6 +1364,7 @@ class Strand(Group):
Group.__init__(self)
self.num_nt = 0
self.children = self.strand_segments = []
self.oxdna_nt = []
self.segname = segname
self.is_circular = is_circular
self.debug = False
......@@ -1452,26 +1477,54 @@ class Strand(Group):
return first_atomic_index
def update_atomic_orientations(self,default_orientation):
def generate_atomic_model(self, scale, first_atomic_index):
last = None
resid = 1
## TODO relabel "strand_segment"
strand_segment_count = 0
for s in self.strand_segments:
strand_segment_count += 1
seg = s.segment
contour = s.get_contour_points()
for c,seq,nt in zip(contour,s.get_sequence(),s.children):
# if s.end == s.start:
# pdb.set_trace()
# assert(s.end != s.start)
assert( s.num_nt == 1 or (np.linalg.norm( seg.contour_to_position(contour[-1]) - seg.contour_to_position(contour[0]) ) > 0.1) )
nucleotide_count = 0
for c,seq in zip(contour,s.get_sequence()):
nucleotide_count += 1
if last is None and not self.is_circular:
seq = "5"+seq
if strand_segment_count == len(s.strand_segments) and nucleotide_count == s.num_nt and not self.is_circular:
seq = seq+"3"
nt = seg._generate_atomic_nucleotide( c, s.is_fwd, seq, scale, s )
## Join last basepairs
if last is not None:
self.link_nucleotides(last,nt)
nt.__dict__['resid'] = resid
resid += 1
last = nt
nt._first_atomic_index = first_atomic_index
first_atomic_index += len(nt.children)
if self.is_circular:
self.link_nucleotides(last,self.strand_segments[0].children[0])
return first_atomic_index
def generate_oxdna_model(self):
for s in self.strand_segments:
seg = s.segment
contour = s.get_contour_points()
assert( s.num_nt == 1 or (np.linalg.norm( seg.contour_to_position(contour[-1]) - seg.contour_to_position(contour[0]) ) > 0.1) )
for c,seq in zip(contour,s.get_sequence()):
nt = seg._generate_oxdna_nucleotide( c, s.is_fwd, seq )
self.oxdna_nt.append(nt)
orientation = seg.contour_to_orientation(c)
## TODO: move this code (?)
if orientation is None:
axis = seg.contour_to_tangent(c)
angleVec = np.array([1,0,0])
if axis.dot(angleVec) > 0.9: angleVec = np.array([0,1,0])
angleVec = angleVec - angleVec.dot(axis)*axis
angleVec = angleVec/np.linalg.norm(angleVec)
y = np.cross(axis,angleVec)
orientation = np.array([angleVec,y,axis]).T
nt.orientation = orientation.dot(default_orientation) # this one should be correct
class SegmentModel(ArbdModel):
def __init__(self, segments=[], local_twist=True, escapable_twist=True,
......@@ -1619,6 +1672,8 @@ class SegmentModel(ArbdModel):
for strand in self.strands:
for s in strand.children:
s.clear_all()
s.oxdna_nt = []
for seg in self.segments:
for d in ('fwd','rev'):
seg.strand_pieces[d] = []
......@@ -2994,6 +3049,12 @@ proc calcforces {} {
first_atomic_index = s.generate_atomic_model(scale,first_atomic_index)
self._assign_basepairs()
def generate_oxdna_model(self, scale=1):
self.clear_beads()
self.children = self.strands
for s in self.strands:
s.generate_oxdna_model()
def vmd_tube_tcl(self, file_name="drawTubes.tcl"):
with open(file_name, 'w') as tclFile:
tclFile.write("## beginning TCL script \n")
......
mrdna/simulate.py
View file @ 0ca4dc50
......@@ -6,7 +6,46 @@ from . import get_resource_path
## TODO: implement replicas, initial conditions specified through some restart, and a custom simulation schedule
def minimize_and_simulate_oxdna( model,
output_name,
directory = None,
num_min_steps = 1e2,
num_steps = 1e7,
output_period = None,
**oxdna_args
):
d_orig = os.getcwd()
try:
if directory is not None:
if not os.path.exists(directory):
os.makedirs(directory)
os.chdir(directory)
model.generate_oxdna_model()
top = "{}-oxdna.top".format(output_name)
model._write_oxdna_topology(top)
top,conf = model.simulate_oxdna(output_name = "{}-oxdna-min".format(output_name),
topology = top,
sim_type = "min",
num_steps = num_min_steps,
print_conf_interval = 100,
print_energy_every = 10,
)
model.simulate_oxdna(output_name = "{}-oxdna".format(output_name),
topology = top,
configuration = conf,
num_steps = num_steps,
print_conf_interval = output_period,
print_energy_every = output_period,
**oxdna_args)
finally:
os.chdir(d_orig)
def multiresolution_simulation( model, output_name,
job_id=None,
......@@ -21,6 +60,8 @@ def multiresolution_simulation( model, output_name,
coarse_local_twist = False,
remove_long_bonds=False,
backbone_scale=1.0,
oxdna_steps = None,
oxdna_output_period = None,
):
## Round steps up to nearest multiple of output_period, plus 1
......@@ -32,8 +73,6 @@ def multiresolution_simulation( model, output_name,
if fine_steps == 1: fine_steps += 1
fine_steps = fine_steps*fine_output_period+1
ret = None
d_orig = os.getcwd()
try:
......@@ -113,20 +152,28 @@ def multiresolution_simulation( model, output_name,
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 """
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)
if oxdna_steps is not None:
minimize_and_simulate_oxdna( model,
output_name = output_name,
num_min_steps = 1e3,
num_steps = oxdna_steps,
output_period = oxdna_output_period )
else:
""" Atomic simulation """
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 )
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:
...
try:
shutil.copytree( get_resource_path("charmm36.nbfix"), "charmm36.nbfix" )
except FileExistsError:
...
ret = directory
......
Markdown is supported
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