[Standard Relaxation Procedure] Clarification Needed

[Zehui127]

Dear Team,

I am new to using oxDNA and am seeking assistance with the standard relaxation procedure.

I am constructing a pipeline for running oxDNA simulations with a given DNA sequence. The workflow is as follows:

1.	[seq] -> nabc -> tacoxDNA -> oxDNA -> tacoxDNA -> [configurations in PDB format]

I have used the script below to run the oxDNA simulation.

However, I am encountering issues with the relaxation step. Specifically, I am looking for the equivalent of ticking the relaxation box available on the oxDNA website (https://oxdna.org/).

From my understanding, it might be possible to adjust the simulation parameters for relaxation in the configuration file as described in the oxDNA documentation on relaxation. My question is: can I achieve relaxation by running two consecutive light simulations with modified parameters?

Any guidance or example configurations/scripts for properly implementing the relaxation procedure would be greatly appreciated.

# visualisation library
from oxDNA_analysis_tools.UTILS.oxview import from_path
# boilerplate code to setup and monitor simulations
from oxDNA_analysis_tools.UTILS.boilerplate import setup_simulation,\
get_default_input,\
Simulation
import os
import time
from tqdm import tqdm

NUM_OF_STEPS = 1000000
OUTPUT_DIR = 'sampled_dna_strut.dat'
PWD = os.path.dirname(os.path.abspath(__file__))
SEQ_DEP_FILE = os.path.join(PWD,"oxDNA2_average_sequence_parameters.txt")

def __compress_dat(conf):
    """Generate a dat file string with newlines and additional zeroes."""
    dat_string = f"t = {conf.time}\nb = {conf.box[0]} {conf.box[1]} {conf.box[2]}\nE = {conf.energy[0]} {conf.energy[1]} {conf.energy[2]}\n"
    dat_string += "\n".join([f"{p[0]:.3f} {p[1]:.3f} {p[2]:.3f} {a1[0]:.8f} {a1[1]:.8f} {a1[2]:.8f} {a3[0]:.8f} {a3[1]:.8f} {a3[2]:.8f} 0.0 0.0 0.0 0.0 0.0 0.0"
                             for p, a1, a3 in zip(conf.positions, conf.a1s, conf.a3s)])
    return dat_string

def save_dat(conf, file_path):
    """Save the config to a .dat file."""
    dat_string = __compress_dat(conf)

    with open(file_path, "w") as file:
        file.write(dat_string)

    print(f"Configuration saved to {file_path}")

# visualisation library
from oxDNA_analysis_tools.UTILS.oxview import from_path
# boilerplate code to setup and monitor simulations
from oxDNA_analysis_tools.UTILS.boilerplate import setup_simulation,\
get_default_input,\
Simulation
# get default settings for a CUDA simulation
parameters = get_default_input()
print(" Simulate with the following configuration")
print(parameters)
# setup sequence dependence
parameters["use_average_seq"] = "no"
parameters["seq_dep_file"] = SEQ_DEP_FILE
parameters['steps'] = NUM_OF_STEPS
# setup simulation
# we will simulate the 75 degree layered tile from
# https://doi.org/10.1021/jacs.8b07180
input_file = setup_simulation("/home/v-zehuili/repositories/HPC_Sim_DNA/ssf_exp1.pdb.top", # topology
"/home/v-zehuili/repositories/HPC_Sim_DNA/ssf_exp1.pdb.oxdna", # starting configuraiton
"./simulation_long", # output path
parameters, # simulation parameters
kill_out_dir = True)
# create a simulation class
s = Simulation(input_file)
# run the simulation
s.run()

# Calculate the total number of configurations to be generated
conf_interval = int(float(parameters['print_conf_interval']))
TOTAL_CONFIG_NUM = NUM_OF_STEPS // conf_interval

# Sleep for 10 seconds
print("Waiting for Simulate to Start")
time.sleep(10)


# Loop to monitor and save the configurations
pbar = tqdm(total=TOTAL_CONFIG_NUM, desc="Simulation progress")
num_configs = 0

while num_configs<TOTAL_CONFIG_NUM:
    num_configs = s.get_conf_count()
    print(num_configs)
    pbar.n = num_configs
    pbar.refresh()
    time.sleep(1)  # Add a small delay to avoid excessive checking
pbar.close()

print(s.get_conf_count(),TOTAL_CONFIG_NUM)
config = s.get_conf(5)[1]  # Assuming the configuration number is 11, adjust as needed
save_dat(config, OUTPUT_DIR)
s.terminate()

[lorenzo-rovigatti]

Hi! I can't remember on top of my head what oxdna.org does, but maybe @zoombya or @ErikPoppleton do.
About your general question, relaxation depends on your initial structure, and my experience is that it is hard to find a reliable protocol that

1. is as efficient as possible
2. works every time

Most of the times you'll need some sort of supervision (or to check some heuristics that let your code evaluate whether more relaxation is required).

[ErikPoppleton]

The relaxation protocol on oxdna.org uses the files from oxDNA/analysis/example_input_files (actually I think it runs the simulations a bit longer than those files for consistency). We have found that this procedure works for almost all structures. Depending on how good your initial configuration is, you might need to adjust the length of MD relax step or modify the max_backbone_force parameter to get good relaxation. One thing we strongly recommend is to generate an external force file with relatively strong forces (we arbitrarily use 3.1) holding all pre-existing base pairs together as relaxing stretched bonds can tear your initial structure apart.

The procedure was first described in section 3.2 of this protocol. We discuss it further in section 4 of this protocol.

[zoombya]

Erik's answer is extremely comprehensive, additionally I wanted to add that it makes little sense to use
save_dat(config, OUTPUT_DIR) \ since you already save while oxDNA is running.
OAT has the minify script to handle trajectories, and if you want to convert them to PDB afterwards you would need the raw data anyway.
Additionally the boilerplate Simulation class provides methods to visualize even trajectories, doing the minify step you copied from the library.