diff --git a/mrdna/model/arbdmodel.py b/mrdna/model/arbdmodel.py
index 197bddbce8755f8e09fe87ec8beecd9988a8d1c6..de24a3a5675f60073a0cb629410740f0211d6083 100644
--- a/mrdna/model/arbdmodel.py
+++ b/mrdna/model/arbdmodel.py
@@ -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)
-
-
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 0273c5dbe183390f5b80110569f6b8cd97c54834..6f5efa809c6b7b118f87eb0ffbbe8732823200d0 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -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")