Commit 205bb869 authored by cmaffeo2's avatar cmaffeo2
Browse files

Moved oxDNA from arbdmodel.py to segmentmodel.py

parent a8ed3cac
......@@ -1191,259 +1191,3 @@ run {num_steps:d}
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 = None,
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',):
if salt_concentration is None:
try:
## units "80 epsilon0 295 k K / (2 (AA)**2 e**2/particle)" mM
salt_concentration = 9331.3126/self.debye_length**2
except:
salt_concentration = 0.5
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)
......@@ -3248,12 +3248,265 @@ proc calcforces {} {
first_atomic_index = s.generate_atomic_model(scale,first_atomic_index)
self._assign_basepairs()
""" OxDNA """
## https://dna.physics.ox.ac.uk/index.php/Documentation#Input_file
def generate_oxdna_model(self, scale=1):
self.clear_beads()
self.children = self.strands
for s in self.strands:
s.generate_oxdna_model()
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 = None,
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',):
if salt_concentration is None:
try:
## units "80 epsilon0 295 k K / (2 (AA)**2 e**2/particle)" mM
salt_concentration = 9331.3126/self.debye_length**2
except:
salt_concentration = 0.5
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)
""" Visualization """
def vmd_tube_tcl(self, file_name="drawTubes.tcl"):
with open(file_name, 'w') as tclFile:
tclFile.write("## beginning TCL script \n")
......
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