From 2ea602db09351e00e8131113233bff3c83d63b37 Mon Sep 17 00:00:00 2001 From: Lorenzo Rovigatti Date: Fri, 15 Nov 2024 17:37:37 +0100 Subject: [PATCH 1/3] Fix a bug in the OXPY_REMD/reshuffle.py script, and strengthen the remd.py script --- examples/OXPY_REMD/README.md | 2 +- examples/OXPY_REMD/remd.py | 69 +++++++++++++++++++++------------ examples/OXPY_REMD/reshuffle.py | 36 +++++++++-------- 3 files changed, 66 insertions(+), 41 deletions(-) diff --git a/examples/OXPY_REMD/README.md b/examples/OXPY_REMD/README.md index c47647cff..f75013e74 100644 --- a/examples/OXPY_REMD/README.md +++ b/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 ``` diff --git a/examples/OXPY_REMD/remd.py b/examples/OXPY_REMD/remd.py index 81ed41740..072de6c2b 100644 --- a/examples/OXPY_REMD/remd.py +++ b/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,15 +148,15 @@ 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] @@ -143,16 +164,16 @@ def remd_log(pid, *args): # 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 diff --git a/examples/OXPY_REMD/reshuffle.py b/examples/OXPY_REMD/reshuffle.py index 0916e9a45..575e8064a 100644 --- a/examples/OXPY_REMD/reshuffle.py +++ b/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]) From 6096fdbbcbb8e30d2b2dd457667de71db791ca54 Mon Sep 17 00:00:00 2001 From: lorenzo-rovigatti Date: Sat, 16 Nov 2024 11:10:11 +0100 Subject: [PATCH 2/3] Fix a bug that made VMMC simulate the wrong model when used with the DNA2 interaction This commit also adds a new test to try to better catch future errors --- src/Backends/VMMC_CPUBackend.cpp | 7 ++--- test/DNA/DSDNA8/VMMC/quick_compare | 2 +- test/DNA/DSDNA8/VMMC_DNA2/quick_compare | 1 + test/DNA/DSDNA8/VMMC_DNA2/quick_input | 38 +++++++++++++++++++++++++ test/DNA/SSDNA15/VMMC/quick_compare | 2 +- 5 files changed, 44 insertions(+), 6 deletions(-) create mode 100644 test/DNA/DSDNA8/VMMC_DNA2/quick_compare create mode 100644 test/DNA/DSDNA8/VMMC_DNA2/quick_input diff --git a/src/Backends/VMMC_CPUBackend.cpp b/src/Backends/VMMC_CPUBackend.cpp index 07766d100..67b6b858a 100644 --- a/src/Backends/VMMC_CPUBackend.cpp +++ b/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(_interaction.get()) == NULL) || (dynamic_cast(_interaction.get()) == NULL) ) { + // the DNA and RNA/RNA2 interactions use the original coaxial stacking term + if(dynamic_cast(_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(_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(_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); diff --git a/test/DNA/DSDNA8/VMMC/quick_compare b/test/DNA/DSDNA8/VMMC/quick_compare index 74a088ec4..f34c38c12 100644 --- a/test/DNA/DSDNA8/VMMC/quick_compare +++ b/test/DNA/DSDNA8/VMMC/quick_compare @@ -1 +1 @@ -ColumnAverage::energy.dat::2::-1.37970256144::0.15 +ColumnAverage::energy.dat::2::-1.37970256144::0.01 diff --git a/test/DNA/DSDNA8/VMMC_DNA2/quick_compare b/test/DNA/DSDNA8/VMMC_DNA2/quick_compare new file mode 100644 index 000000000..7f912b2d8 --- /dev/null +++ b/test/DNA/DSDNA8/VMMC_DNA2/quick_compare @@ -0,0 +1 @@ +ColumnAverage::energy.dat::2::-1.379235::0.02 diff --git a/test/DNA/DSDNA8/VMMC_DNA2/quick_input b/test/DNA/DSDNA8/VMMC_DNA2/quick_input new file mode 100644 index 000000000..5b8f47cbe --- /dev/null +++ b/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 diff --git a/test/DNA/SSDNA15/VMMC/quick_compare b/test/DNA/SSDNA15/VMMC/quick_compare index 57987ea7e..b9f9cdfe8 100644 --- a/test/DNA/SSDNA15/VMMC/quick_compare +++ b/test/DNA/SSDNA15/VMMC/quick_compare @@ -1 +1 @@ -ColumnAverage::energy.dat::2::-0.700783241758::0.26 +ColumnAverage::energy.dat::2::-0.700783241758::0.15 From b78d23c73c26571eb36def3a2f59f586585461ae Mon Sep 17 00:00:00 2001 From: lorenzo-rovigatti Date: Sat, 16 Nov 2024 11:14:04 +0100 Subject: [PATCH 3/3] Prepare the repo for the new release --- CHANGELOG | 47 ++++++++++++++++++++++++++++++++++++++++++++++- README.md | 2 +- 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 3e2ee39cf..90bff4414 100644 --- a/CHANGELOG +++ b/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] diff --git a/README.md b/README.md index 44736d804..752c92a48 100644 --- a/README.md +++ b/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?**