DNA-RNA hybrids

.github/workflows/test-linux.yml

@@ -22,4 +22,5 @@ jobs:
     - name: Run Simple Test Script
       run: |
         pushd build
+        make test_run
         make test_quick

DRH_sequence_dependent_parameters.txt

@@ -0,0 +1,4 @@
+HYDR_A_U = 1.21
+HYDR_A_T = 1.37
+HYDR_rC_dG = 1.61
+HYDR_rG_dC = 1.77

README.md

@@ -30,10 +30,6 @@ The `analysis/paper_examples` folder contains examples for `oxDNA_analysis_tools
 
 **A:** No, oxDNA can run simulations on single cores or single GPUs only.
 
-**Q: Can I simulate systems containing both DNA and RNA?**
-
-**A:** Unfortunately not: at the moment there is no force field for that.
-
 ## Citing oxDNA
 
 [![DOI](https://joss.theoj.org/papers/10.21105/joss.04693/status.svg)](https://doi.org/10.21105/joss.04693)

docs/source/input.md

@@ -36,7 +36,7 @@ These are the options that control the overall behaviour of the simulation and o
 * `time_scale = linear/log_lin`: a linear time_scale will make oxDNA print linearly-spaced configurations. a log_lin will make it print linearly-spaced cycles of logarithmically-spaced configurations.
 * `print_conf_interval = <int>`: if the time scale is linear, this is the number of time steps between the outputing of configurations, otherwise this is just the first point of the logarithmic part of the log_lin time scale.
 * `print_energy_every = <int>`: number of time steps between the outputing of the energy (and of the other default observables such as acceptance ratios in Monte Carlo simulations).
-* `[interaction_type = DNA|DNA2|RNA|RNA2|LJ|...]`: Particle-particle interaction of choice. Check the documentation relative to the specific interaction for more details and the supported interaction-specific options. Defaults to `DNA`.
+* `[interaction_type = DNA|DNA2|RNA|RNA2|DRH|LJ|...]`: Particle-particle interaction of choice. Check the documentation relative to the specific interaction for more details and the supported interaction-specific options. Defaults to `DNA`.
 * `[max_io = <float>]`: the maximum rate at which the output is printed, in MB/s. This is a useful option to avoid filling up the disk too quickly. Increase the default value (1 MB/s) at your own risk! 
 * `[fix_diffusion = <bool>]`: if true, particles that leave the simulation box are brought back in via periodic boundary conditions. Defaults to `true`.
 * `[fix_diffusion_every = <int>]`: number of time steps every which the diffusion is fixed. Used only if `fix_diffusion = true`, defaults to 100000 ({math}`10^5`).
@@ -169,7 +169,15 @@ The following options control the behaviour of MC simulations.
 
 * `[salt = <float>]`: the salt concentration in molar (M). Defaults to 1.
 * `[mismatch_repulsion = <bool>]`: if `true`, add a repulsion between mismatches. Defaults to `false`.
-* `[mismatch_repulsion_strength = <float>]`: set the strength of the repulsion between mismatches. Used only if `mismatch_repulsion = true`, defaults to `1`. 
+* `[mismatch_repulsion_strength = <float>]`: set the strength of the repulsion between mismatches. Used only if `mismatch_repulsion = true`, defaults to `1`.
+
+## Options for `DRH` simulations
+
+* `[use_average_seq = <bool>]`: use the average-sequence parameters. Defaults to `true`.
+* `[seq_dep_file_DNA = <path>]`: path to the location of the file containing the DNA sequence-dependent parameters. Mandatory if `use_average_seq = false`.
+* `[seq_dep_file_RNA = <path>]`: path to the location of the file containing the DNA sequence-dependent parameters. Mandatory if `use_average_seq = false`.
+* `[seq_dep_file_DRH = <path>]`: path to the location of the file containing the DNA-RNA sequence-dependent parameters. Mandatory if `use_average_seq = false`.
+
 
 ## Options for `LJ` (Lennard-Jones) simulations
 

docs/source/install.md

@@ -118,7 +118,7 @@ cmake -DPython=1 -DPYTHON_EXECUTABLE=$HOME/miniconda3/bin/python -DPYTHON_INCLUD
 
 ## Testing
 
-* `make test_run` runs quick tests to check whether oxDNA has been correctly compiled or not.	
+* `make test_run` runs quick tests to check that oxDNA has been correctly compiled, and that the core interactions (`DNA2`, `RNA2` and `DRH`) compute the energy contributions in a nicked double strand correctly.
 * `make test_quick` runs longer tests to check that oxDNA works.
 * `make test_oxpy` checks that the Python bindings work.
 * `make test` runs all sets of tests above.

src/Backends/FFS_MD_CPUBackend.cpp

@@ -218,6 +218,7 @@ number FFS_MD_CPUBackend::pair_interaction_nonbonded_DNA_with_op(BaseParticle *p
 
 	_interaction->set_computed_r(computed_r);
 	number energy = _interaction->pair_interaction_term(DNAInteraction::HYDROGEN_BONDING, p, q, false, update_forces);
+
 	if(energy <= MAX_BOND_CUTOFF) {
 		_op.add_hb(q->index, p->index, energy);
 	}
@@ -235,6 +236,8 @@ number FFS_MD_CPUBackend::pair_interaction_nonbonded_DNA_with_op(BaseParticle *p
 		energy += _interaction->pair_interaction_term(RNA2Interaction::DEBYE_HUCKEL, p, q, false, update_forces);
 	}
 
+
+
 	return energy;
 }
 

src/Backends/VMMC_CPUBackend.cpp

@@ -7,6 +7,7 @@
 #include "../Utilities/Utils.h"
 #include "../Interactions/RNAInteraction2.h"
 #include "../Interactions/DNA2Interaction.h"
+#include "../Interactions/DRHInteraction.h"
 
 #include <set>
 #include <sstream>
@@ -192,6 +193,7 @@ void VMMC_CPUBackend::get_settings(input_file & inp) {
 	ok_interactions.push_back("DNA2_nomesh");
 	ok_interactions.push_back("RNA");
 	ok_interactions.push_back("RNA2");
+	ok_interactions.push_back("DRH");
 	if(getInputString(&inp, "interaction_type", inter, 0) == KEY_FOUND) {
 		// std::find points is equal to ok_interactions.end() if it can't find inter in ok_interactions.
 		if(std::find(ok_interactions.begin(), ok_interactions.end(), inter) == ok_interactions.end()) {
@@ -424,17 +426,26 @@ 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
-	if(dynamic_cast<DNA2Interaction *>(_interaction.get()) == NULL)
-	energy += _interaction->pair_interaction_term(DNAInteraction::COAXIAL_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) ) {
+		energy += _interaction->pair_interaction_term(DNAInteraction::COAXIAL_STACKING, p, q, false, false);
+	}
 
-	if(dynamic_cast<DNA2Interaction *>(_interaction.get()) != NULL) {
+	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);
 	}
 	else if(dynamic_cast<RNA2Interaction *>(_interaction.get()) != NULL) {
 		energy += _interaction->pair_interaction_term(RNA2Interaction::DEBYE_HUCKEL, p, q, false, false);
 	}
+
+
 
 	return energy;
 }

src/CMakeLists.txt

@@ -122,6 +122,8 @@ SET(interactions_SOURCES
 	Interactions/PatchyInteractionDan.cpp
 	Interactions/KFInteraction.cpp
 	Interactions/CustomInteraction.cpp
+	Interactions/DRHInteraction.cpp
+	Interactions/DRHInteraction_relax.cpp
 )
 
 SET(move_SOURCES

src/Interactions/DNA2Interaction.h

@@ -25,7 +25,7 @@
 
 #include "DNAInteraction.h"
 
-class DNA2Interaction: public DNAInteraction {
+class DNA2Interaction: virtual public DNAInteraction {
 
 protected:
 	float _salt_concentration;

src/Interactions/DNAInteraction.cpp

@@ -233,7 +233,15 @@ void DNAInteraction::get_settings(input_file &inp) {
 
 	if(getInputBool(&inp, "use_average_seq", &_average, 0) == KEY_FOUND) {
 		if(!_average) {
-			getInputString(&inp, "seq_dep_file", _seq_filename, 1);
+			std::string inter_type;
+			getInputString(&inp, "interaction_type", inter_type, 0);
+
+			if(inter_type.compare("DRH") || inter_type.compare("DRH_relax")) {
+				getInputString(&inp, "seq_dep_file_DNA", _seq_filename, 1);
+			} else {
+				getInputString(&inp, "seq_dep_file", _seq_filename, 1);
+			}
+
 			OX_LOG(Logger::LOG_INFO, "Using '%s' as the input for sequence-dependent values", _seq_filename.c_str());
 		}
 	}

src/Interactions/DNAInteraction.h

@@ -12,7 +12,7 @@
  * This is the default interaction.
  */
 
-class DNAInteraction : public BaseInteraction {
+class DNAInteraction : virtual public BaseInteraction {
 protected:
 	bool _average;
 	std::string _seq_filename;