Merge branch 'master' of https://github.com/lorenzo-rovigatti/oxDNA

CHANGELOG

@@ -1,4 +1,49 @@
-v. 3.6.1
+v. 3.7.0
+	- Add two new options to manage the format of the numerical output of some observables [6d03f40]
+	- Deprecate many old plugins in the contrib/rovigatti folder [95bd2e0]
+	- Add compiler flags to support Ampere GPUs [21d7e48]
+	- Update generate-sa.py to support generation of dsDNA with sticky ends and add an example to show its functionalities [08dcc2c]
+	- Add an attraction plane that constrains particles to stay on one side of it while attracting them with a constant force (`type = attraction_plane`).  Documentation is also included [f0e0c9d]
+	- Port the ANNaMo interaction to CUDA [3568256]
+	- Make it possible to print the kinetic and total energies with any precision [f42d61d]
+	- Make the NA interaction support custom bases (fix #129) [306a276]
+	- Make it possible to simulate arbitrary LJ binary mixtures [b191c3d]
+	- DPS: Add an optional Yukawa repulsion to the model [dbdc726]
+	- DPS: Add an input file key to disable three body interaction [2d62ba7]
+	- DPS: Fix a bug whereby species with 0 particles would not be handled correctly [2d33bcc]
+	- Fix a bug whereby the CUDA code would crash for particular initial configurations made by nucleotides aligned along z (#120) [2b69edf]
+	- Fix an example (#123) [f4853f9]
+	- Fix a bug whereby the code would hang if bases_file (used by HBEnergy) was not readable [a20df96]
+	- Fix a bug that made VMMC simulate the wrong model when used with the DNA2 interaction [6096fdb]
+	- Improve the docs (see *e.g.* #116)
+	oxpy:
+	  - Fix #97 (FFS Simulation Type Fails to Log Correctly When Running Consecutive Processes) [2b93f63]
+	  - Add a workaround documentation to a known oxpy bug [960a8c5]
+	  - Fix #138: make oxpy's get_bool() to handle non-lowercase booleans [208049d]
+	  - Improve the docs
+	oat:
+	  - Duplex angle plotter now works with ANY nucleotide id in the duplex [eb17a42]
+	  - Fixed bug where output_bonds would fail instead of exiting normally when no plot was requested [c287226]
+	  - Fixed bug where duplexes were detected as longer than they should be if they became unpaired [e1580e6]
+	  - Fixed bug where inboxing would always center the structure, center=False works now [b23bdb7]
+	  - Fixed parser and data type errors in oxDNA_PDB [d951742]
+	  - Fixed bug where strands would sometimes be backwards in converted PDBs and fixed bug where the box center was incorrect during inboxing [2f179ad]
+	  - In oxview.py, now from_path supports multiple systems loading [0050492]
+	  - Updated bond_analysis to be able to drop data files [27d2fbb]
+	  - Added __len__ methods to System and Strand [ac5e632]
+	  - Changed read_force_file to interpret particle IDs as ints rather than floats [e7f6818]
+	  - Make oxDNA_PDB callable from scripting interface [a20207b]
+	  - Improved consistency in oat cli calls [3d628eb]
+	  - Updated oat superimpose to allow separate index files for each configuration [8a6ed9d]
+	  - Removal of now removed 'center' keywords from align function call [1cb3b40]
+	  - Fix error in input file parsing (#119) [bca41ab]
+	  - Created forces/pairs to dot-bracket converters and a skeleton script for trajectory analysis [cce8009]
+	  - Fixed bugs in pdb output file naming and improved RNA detection [e753007]
+	  - Fix parsing error for strand circularity in oat's new topology reader [a168f8d]
+	  - Faster linear version of decimate script [adbba3c]
+	  - Added PDB -> oxDNA converter to oat [e0084bc]
+
+v. 3.6.1 [6854eec]
 	- Fix some bugs and docs errors due to the DRH -> NA change (see #81) [79728d3 and 7cf4b94]
 	- Add and document an option to control the DH cut-off (see #83) [95a9f0d]
 	- Make the MPI VMMC backend runnable once again, after it had been broken for years (see #85) [8becd5f]

README.md

@@ -28,7 +28,7 @@ The `analysis/paper_examples` folder contains examples for `oxDNA_analysis_tools
 
 **Q: Can oxDNA be run on multiple CPU cores or GPUs?**
 
-**A:** No, oxDNA can run simulations on single cores or single GPUs only.
+**A:** No, oxDNA can run simulations on single cores or single GPUs only. However, the Python bindings (`oxpy`) can be used to implement parallel algorithms on top of `oxDNA`. See for instance `examples/METADYNAMICS` and `examples/OXPY_REMD`.
 
 **Q: How is the oxDNA force field defined?**
 

examples/OXPY_REMD/README.md

@@ -39,7 +39,7 @@ To run a REMD simulation execute the following commands:
 ```shell
 ./clean.sh
 
-mpirun -n 4  python3  remd.py  input_md
+mpirun -n 4 python3 remd.py input_md
 
 python3 reshuffle.py input_md history.json reshuffled
 ```

examples/OXPY_REMD/remd.py

@@ -8,6 +8,8 @@
 from mpi4py import MPI
 from random import random
 
+import argparse
+
 import oxpy
 
 try:
@@ -33,42 +35,61 @@ def get_order(i, rank, locations, mpi_nprocs):
         return im_responsible, im_rank, talk_to 
     else:
         return None    
+
+
+def remd_log(pid, *args):
+    if pid == rank:
+        print("REMD_LOG:", pid, ":", *args)
 
 
 comm = MPI.COMM_WORLD
-rank = comm.Get_rank()  # process id  
-mpi_nprocs = comm.Get_size()  # number of processes running 
+rank = comm.Get_rank()  # process id
+mpi_nprocs = comm.Get_size()  # number of processes running
 locations = list(range(mpi_nprocs))  # curent location of the exchange temperature
 exchange = np.zeros(mpi_nprocs)  # keep track of exchanges
 exchange_tries = np.zeros(mpi_nprocs)  # and attempts
 rates = np.zeros(mpi_nprocs)
 
-if rank == 0:
-    history = []
-
+parser = argparse.ArgumentParser(description="Run an REMD simulation.")
+parser.add_argument('input_file', type=str, nargs=1, help="The input file of the simulation")
+args = None
+try:
+    if rank == 0:
+        args = parser.parse_args()
+finally:
+    args = comm.bcast(args, root=0)
 
-def remd_log(pid, *args):
-    if pid == rank:
-        print("x", pid, ":", *args)
+if args == None:
+    exit(1)
+else:
+    input_file = args.input_file[0]
 
+if rank == 0:
+    history = []
 
 with oxpy.Context():
     input = oxpy.InputFile()  # load up conf and input
-    input.init_from_filename(sys.argv[1])  # as 1st is the name of the script
+    input.init_from_filename(input_file)  # as 1st is the name of the script
+
     # have separate files for every running instance
-
     input["energy_file"] = f"energy_{rank}.dat"
     input["lastconf_file"] = f"last_conf_{rank}.dat"
     input["CUDA_device"] = str(rank)
     # figure out steps to run and our temperature
     pt_move_every = int(input["pt_move_every"])
-    steps_to_run = int(int(input["steps"]) / pt_move_every) 
+    steps_to_run = int(int(input["steps"]) / pt_move_every)
     # setup temperatures
     temperatures = input["pt_temp_list"].split(",")
 
     input["T"] = temperatures[rank]  # set process temperature
     input["trajectory_file"] = f"trajectory_{rank}.dat"
 
+    # check that the number of processes is equal to the number of temperatures
+    if rank == 0:
+        if len(temperatures) != mpi_nprocs:
+            print(f"The number of temperatures ({len(temperatures)}) should match the number of MPI processes ({mpi_nprocs})", file=sys.stderr)
+            comm.Abort(errorcode=1)
+
     # convert temperaures to ox units
     if "K" in input["T"]:
         temperatures = [oxpy.utils.Kelvin_to_oxDNA(t.replace("K", "")) for t in temperatures]   
@@ -86,10 +107,10 @@ def remd_log(pid, *args):
         locations = comm.bcast(locations, root=0)
         # let's wait for everyone to get started 
         swap = False
-        resut = get_order(i, rank , locations, mpi_nprocs)
-        if resut:
+        result = get_order(i, rank , locations, mpi_nprocs)
+        if result:
             # unpack result
-            im_responsible, im_rank, talk_to = resut
+            im_responsible, im_rank, talk_to = result
             if im_responsible:
                 irresp_energy = comm.recv(source=talk_to)
 
@@ -127,32 +148,32 @@ def remd_log(pid, *args):
 
         # handle updates to locations
         if rank == 0:
-            swaped = [] 
+            swapped = [] 
             # handle 0
-            if  resut and im_responsible:
+            if  result and im_responsible:
                 exchange_tries[0] += 1
                 exchange_tries[talk_to] += 1 
                 if swap:
                     locations[0] , locations[talk_to] = locations[talk_to] , locations[0]
-                    swaped.append(0)
-                    swaped.append(talk_to)    
+                    swapped.append(0)
+                    swapped.append(talk_to)    
                     exchange[0] += 1
                     exchange[talk_to] += 1
                     rates [0] = exchange[0] / exchange_tries[0]
                     rates [talk_to] = exchange[talk_to] / exchange_tries[talk_to]
 
             # receive from everyone but me the updates
             for j in range(1, mpi_nprocs):
-                resut = comm.recv(source=j)
-                if resut:
-                    swap, im_responsible, im_rank, talk_to = resut
+                result = comm.recv(source=j)
+                if result:
+                    swap, im_responsible, im_rank, talk_to = result
                     if im_responsible:
                         exchange_tries[j] += 1
                         exchange_tries[talk_to] += 1
                         if swap: 
                             locations[j] , locations[talk_to] = locations[talk_to] , locations[j]
-                            swaped.append(j)
-                            swaped.append(talk_to)
+                            swapped.append(j)
+                            swapped.append(talk_to)
                             exchange[j] += 1
                             exchange[talk_to] += 1
                             rates [j] = exchange[0] / exchange_tries[j]
@@ -166,7 +187,7 @@ def remd_log(pid, *args):
             remd_log(0, "rates:", rates)
         else:
             # notify 0 of update
-            if resut:
+            if result:
                 comm.send((swap, im_responsible, im_rank, talk_to), 0)
             else:
                 comm.send(None, 0)  # nothing happened for this round

examples/OXPY_REMD/reshuffle.py

@@ -1,22 +1,22 @@
 # reintegrate trajectory
 from json import loads
 from oxDNA_analysis_tools.UTILS.RyeReader import describe, get_confs, get_input_parameter, write_conf
+import shutil, sys
 
 import argparse
 
 if __name__ == "__main__":
-    #handle commandline arguments
-    parser = argparse.ArgumentParser(description="Retrive temperature trajectories from an REMD run.")
-    parser.add_argument('input_file', type=str, nargs=1, help="The inputfile used to run the simulation")
-    parser.add_argument('history_file', type=str, nargs=1, help="The json history file keeping the remd run information.")
-    parser.add_argument('output_file', type=str, nargs=1, help="name of the file to save the output trajectories to")
+    # handle commandline arguments
+    parser = argparse.ArgumentParser(description="Retrieve temperature trajectories and last configurations from an REMD run.")
+    parser.add_argument('input_file', type=str, nargs=1, help="The input file used to run the simulation")
+    parser.add_argument('history_file', type=str, nargs=1, help="The json history file keeping the REMD run information.")
+    parser.add_argument('output_file', type=str, nargs=1, help="Prefix used to build the name of the file to save the output trajectories to")
 
     args = parser.parse_args()
     input_file = args.input_file[0]
     history_file = args.history_file[0]
     out_file = args.output_file[0]
 
-
     # parse the input file to retrieve the temperature settings
     pt_temp_list = None
     with open(input_file) as file:
@@ -29,21 +29,25 @@
         history = loads(file.read())
 
     top_name = get_input_parameter(input_file, 'topology')
-    traj_name = "trajectory"#get_input_parameter(input_file, 'trajectory_file')
-    #traj_base = traj_name.split('.')[:-1]
-    traj_ext = "dat"#traj_name.split('.')[-1]
+    traj_name = "trajectory"
+    traj_ext = "dat"
     tis = []
     dis = []
     for i in range(len(pt_temp_list)):
-        ti, di = describe(top_name, f'{traj_name}_{i}.{traj_ext}')
+        ti, di = describe(top_name, f"trajectory_{i}.dat")
         tis.append(ti)
         dis.append(di)
-
-    out_files = [f"./{out_file}_{T}.dat" for T in pt_temp_list]
 
+    # generate the ordered trajectories        
+    out_files = [f"./{out_file}_{T}.dat" for T in pt_temp_list]
     read_poses = [0 for _ in range(len(pt_temp_list))]
-    for i,locations in enumerate(history):
-        for j,l in enumerate(locations):
-            configuration = get_confs(tis[l], dis[l], read_poses[l], 1)[0]
-            read_poses[l] += 1
+    for i, locations in enumerate(history):
+        for j, l in enumerate(locations):
+            configuration = get_confs(tis[j], dis[j], read_poses[j], 1)[0]
+            read_poses[j] += 1
             write_conf(out_files[l], configuration, append=True, include_vel=dis[0].incl_v)
+
+    # print the last_conf files
+    last_files = [f"./last_conf_{T}.dat" for T in pt_temp_list]
+    for i, l in enumerate(history[-1]):
+        shutil.copy(f"last_conf_{i}.dat", last_files[l])

src/Backends/VMMC_CPUBackend.cpp

@@ -420,17 +420,16 @@ inline number VMMC_CPUBackend::_particle_particle_nonbonded_interaction_VMMC(Bas
 	energy += _interaction->pair_interaction_term(DNAInteraction::NONBONDED_EXCLUDED_VOLUME, p, q, false, false);
 	energy += _interaction->pair_interaction_term(DNAInteraction::CROSS_STACKING, p, q, false, false);
 
-	// all interactions except DNA2Interaction use the DNAInteraction coaxial stacking*
-	// *the hybrid interaction is a second exception
-	if( (dynamic_cast<DNA2Interaction *>(_interaction.get()) == NULL) || (dynamic_cast<DRHInteraction *>(_interaction.get()) == NULL) ) {
+	// the DNA and RNA/RNA2 interactions use the original coaxial stacking term
+	if(dynamic_cast<DNA2Interaction *>(_interaction.get()) == NULL) {
 		energy += _interaction->pair_interaction_term(DNAInteraction::COAXIAL_STACKING, p, q, false, false);
 	}
 
+	// the DNA2 and DRH interactions use DNA2's coaxial stacking term
 	if(dynamic_cast<DRHInteraction *>(_interaction.get()) != NULL) {
 		energy += _interaction->pair_interaction_term(DRHInteraction::COAXIAL_STACKING, p, q, false, false);
 		energy += _interaction->pair_interaction_term(DRHInteraction::DEBYE_HUCKEL, p, q, false, false);
 	}
-
 	else if(dynamic_cast<DNA2Interaction *>(_interaction.get()) != NULL) {
 		energy += _interaction->pair_interaction_term(DNA2Interaction::COAXIAL_STACKING, p, q, false, false);
 		energy += _interaction->pair_interaction_term(DNA2Interaction::DEBYE_HUCKEL, p, q, false, false);

test/DNA/DSDNA8/VMMC/quick_compare

@@ -1 +1 @@
-ColumnAverage::energy.dat::2::-1.37970256144::0.15
+ColumnAverage::energy.dat::2::-1.37970256144::0.01

test/DNA/DSDNA8/VMMC_DNA2/quick_compare

@@ -0,0 +1 @@
+ColumnAverage::energy.dat::2::-1.379235::0.02

test/DNA/DSDNA8/VMMC_DNA2/quick_input

@@ -0,0 +1,38 @@
+####  PROGRAM PARAMETERS  ####
+backend = CPU
+#debug = 1
+#seed = 4982
+
+####    SIM PARAMETERS    ####
+sim_type = VMMC
+ensemble = NVT
+
+interaction_type = DNA2
+salt_concentration = 0.5
+
+delta_translation = 0.22
+delta_rotation = 0.22
+
+steps = 5e5
+newtonian_steps = 103
+diff_coeff = 2.50
+#pt = 0.1
+thermostat = john
+
+T = 20C 
+dt = 0.005
+verlet_skin = 0.5
+
+####    INPUT / OUTPUT    ####
+topology = ../dsdna8.top
+conf_file = ../init.dat
+trajectory_file = trajectory.dat
+refresh_vel = 1
+log_file = log.dat
+no_stdout_energy = 1
+restart_step_counter = 1
+energy_file = energy.dat
+print_conf_interval = 1e5
+print_energy_every = 1e3 
+time_scale = linear
+external_forces = 0

test/DNA/SSDNA15/VMMC/quick_compare

@@ -1 +1 @@
-ColumnAverage::energy.dat::2::-0.700783241758::0.26
+ColumnAverage::energy.dat::2::-0.700783241758::0.15