Using oxDNA from Python

# clean up all files from the previous run
!rm trajectory.dat
!rm energy.dat
!rm *.pyidx
!rm relax_energy.dat
!rm bonds*
p = None

# all functions required to read a configuration using the new RyeReader
from oxDNA_analysis_tools.UTILS.RyeReader import describe, get_confs, inbox
# the function used to visualize a configuration in oxView
from oxDNA_analysis_tools.UTILS.oxview import oxdna_conf

top = "./hairpin.top"
traj = "./hairpin.conf"
input_file = "./oxDNA_input"
# RyeReader uses indexing allows for random access in the trajectory
top_info, traj_info = describe(top, traj)

# This is the new way to read configurations, as it returns a list we need [0]
# to access a single conf
ref_conf = get_confs(top_info, traj_info, 0, 1)[0]
# inbox the configuration
ref_conf = inbox(ref_conf)

# show an iframe displaying the configuration
# inbox_settings = ["None", "None"] makes sure oxview does no post processing of the
# loaded configuration
oxdna_conf(top_info, ref_conf, inbox_settings=["None","None"])

import oxpy
import multiprocessing
# define the number of simulation steps we run
STEPS = 2000000
INTERVAL = 10000
dt = 0.003

# helper function offloading the work to a new process
def spawn(f, args = ()):
    p = multiprocessing.Process(target=f, args= args )
    p.start()
    return p

# one simulation instance
# if you want to have multiple make sure the output files are different 
def relax_replica():
    with oxpy.Context():
        input = oxpy.InputFile()
        input.init_from_filename(input_file)
        # all input parameters provided to oxDNA need to be strings 
        input["steps"] = str(STEPS)
        input["print_conf_interval"] = str(INTERVAL)
        input["print_energy_every"] = str(INTERVAL)
        input["dt"] = str(dt)
        # make sure our output is not overcrowded
        input["no_stdout_energy"] = "true"
        # these settings turn a regular simulation into a relaxation one  
        input["T"] = "10C"
        input["max_backbone_force"] = "5"
        input["max_backbone_force_far"] = "10"
        input["energy_file"] = "relax_energy.dat"
        # discard the relax data and log info
        input["trajectory"] = "/dev/null"
        input["log_file"] = "/dev/null"

        # init the manager with the given input file
        manager = oxpy.OxpyManager(input)
        #run complete run's it till the number steps specified are reached 
        manager.run_complete()

# make sure we don't have anything running
try:
    p.terminate()
except:
    pass
# run a simulation and obtain a reference to the background process
p = spawn(relax_replica)

import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv("relax_energy.dat", delimiter="\s+",names=['time', 'U','P','K'])

# make sure our figure is bigger
plt.figure(figsize=(15,3)) 
# plot the energy
plt.plot(df.time/dt,df.U)
# and the line indicating the complete run
plt.ylim([-2,1])
plt.plot([STEPS,STEPS],[df.U.max(),df.U.min()-2], color="r")
plt.ylabel("Energy")
plt.xlabel("Steps")
# and the relaxed state line
plt.plot([0, STEPS], [-1,-1], color="g")
print("Relaxation is running:", p.is_alive())

#load and inbox
ti, di = describe(top, "./last_conf.dat")
conf = inbox(get_confs(ti,di, 0,1)[0])
# display
oxdna_conf(ti, conf, inbox_settings=["None","None"])

STEPS = 60000000
INTERVAL = 50000
# define our boiler plate code to run the simulation
def replica():
    with oxpy.Context():
        input = oxpy.InputFile()
        input.init_from_filename(input_file)
        # all input parameters provided to oxDNA need to be strings 
        input["conf_file"] = "last_conf.dat"
        input["steps"] = str(STEPS)
        input["print_conf_interval"] = str(INTERVAL)
        input["print_energy_every"] = str(INTERVAL)
        input["dt"] = str(dt)
        # make sure our output is not overcrowded
        input["no_stdout_energy"] = "true"
        input["log_file"] = "log_production.dat"
        
        # init the manager with the given input file
        manager = oxpy.OxpyManager(input)
        #run complete run's it till the number steps specified are reached 
        manager.run_complete()

#make sure we don't have anything running
try:
    p.terminate()
except:
    pass
#note that for more than one replica running you will have to change the guard code 
spawn(replica)

import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv("energy.dat", delimiter="\s+", names=['time', 'U','P','K'])

# make sure our figure is bigger
plt.figure(figsize=(15,3)) 
# plot the energy
plt.plot(df.time/dt, df.U)
# and the line indicating the complete run
plt.ylim([-1.5,0])
plt.plot([STEPS,STEPS], [df.U.max(), df.U.min()-2], color="r")
plt.ylabel("Energy")
plt.xlabel("Steps")
# and the relaxed state line
plt.plot([0, STEPS], [-1, -1], color="b")
plt.plot([0, STEPS], [-1.2, -1.2], color="g")
print("Simulation is running:", p.is_alive())
plt.show()


#load last conf
ti, di = describe(top, "./last_conf.dat")
conf = get_confs(ti,di, 0,1)[0]
print("Last configuration is:")

#produce the hydrogen bond overlay
!oat output_bonds -v bonds -p 2 oxDNA_input last_conf.dat

from json import loads
with open("bonds_total.json") as file:
    overlay = loads(file.read())

# display
oxdna_conf(ti, conf, overlay=overlay)