diff --git a/analysis/tests/cleanup.sh b/analysis/tests/cleanup.sh
index 194ef0f80..8722bfb56 100755
--- a/analysis/tests/cleanup.sh
+++ b/analysis/tests/cleanup.sh
@@ -1 +1 @@
-rm align* angle* *.json *.pyidx *.png *.npy *.mp4 centroid* cluster* forces.txt gen_pairs.txt mean* min.dat sub*
+rm -f align* angle* *.json *.pyidx *.png *.npy *.mp4 centroid* cluster* forces.txt gen_pairs.txt mean* min.dat sub*
diff --git a/analysis/tests/test.sh b/analysis/tests/test.sh
index 5d4044356..eb6904a61 100755
--- a/analysis/tests/test.sh
+++ b/analysis/tests/test.sh
@@ -1,175 +1,78 @@
 #!/bin/bash
 
-echo "Testing align..."
-if
-    python ../src/oxDNA_analysis_tools/align.py minitraj.dat aligned.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
-
-echo ""
+echo "Testing oat"
+echo "-----------"
 
-echo "Testing backbone_flexibility..."
-if
-    python ../src/oxDNA_analysis_tools/backbone_flexibility.py rna_tile.top minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
+if [ $# -gt 0 ];
 then
-    echo "AN ERROR OCCURED"
+    PYTHON=$1
 else
-    echo "OK"
+    PYTHON="python"
 fi
 
-echo ""
+test() {
+    tmpfile=$(mktemp)
+    $PYTHON $1 &> $tmpfile
+    local code_res=$?
+    grep_out=$(grep -i "ERROR" $tmpfile)
+    local grep_res=$?
 
-echo "Testing bond_analysis..."
-if
-    python ../src/oxDNA_analysis_tools/bond_analysis.py -p2 input_rna minitraj.dat pairs.txt 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+    if [ $code_res -eq 0 ] && [ $grep_res -eq 1 ]
+    then
+        echo " OK"
+    else
+        echo ""
+        echo $grep_out
+        echo "--> AN ERROR OCCURRED <--"
+    fi
+}
 
-echo ""
+echo -n "Testing align..."
+test "../src/oxDNA_analysis_tools/align.py minitraj.dat aligned.dat"
 
-echo "Testing mean and deviations..."
-if
-    python  ../src/oxDNA_analysis_tools/mean.py -p 2 -d devs.json minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing backbone_flexibility..."
+test "../src/oxDNA_analysis_tools/backbone_flexibility.py rna_tile.top minitraj.dat"
 
-echo ""
+echo -n "Testing bond_analysis..."
+test "../src/oxDNA_analysis_tools/bond_analysis.py -p2 input_rna minitraj.dat pairs.txt"
 
-echo "Testing centroid with indexing..."
-if
-    python ../src/oxDNA_analysis_tools/centroid.py -i index.txt mean.dat minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing mean and deviations..."
+test "../src/oxDNA_analysis_tools/mean.py -p 2 -d devs.json minitraj.dat"
 
-echo ""
+echo -n "Testing centroid with indexing..."
+test "../src/oxDNA_analysis_tools/centroid.py -i index.txt mean.dat minitraj.dat"
 
-echo "Testing contact_map..."
-if
-    python ../src/oxDNA_analysis_tools/contact_map.py minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing contact_map..."
+test "../src/oxDNA_analysis_tools/contact_map.py minitraj.dat"
 
-echo ""
+echo -n "Testing distance and clustering (this one takes a while because of the plot)..."
+test "../src/oxDNA_analysis_tools/distance.py -c -i minitraj.dat 1 3 5 67 34 56"
 
-echo "Testing distance and clustering (this one takes a while because of the plot)..."
-if
-    python ../src/oxDNA_analysis_tools/distance.py -c -i minitraj.dat 1 3 5 67 34 56 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing duplex_finder..."
+test "../src/oxDNA_analysis_tools/duplex_finder.py input_rna minitraj.dat"
 
-echo ""
+echo -n "Testing duplex_angle_plotter..."
+test "../src/oxDNA_analysis_tools/duplex_angle_plotter.py -o angle.png -i angles.txt 7 37"
 
-echo "Testing duplex_finder..."
-if
-    python ../src/oxDNA_analysis_tools/duplex_finder.py input_rna minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing generate_force..."
+test "../src/oxDNA_analysis_tools/generate_force.py -f gen_pairs.txt input_rna minitraj.dat"
 
-echo ""
+echo -n "Testing minify..."
+test "../src/oxDNA_analysis_tools/minify.py -a -d2 minitraj.dat min.dat"
 
-echo "Testing duplex_angle_plotter..."
-if
-    python ../src/oxDNA_analysis_tools/duplex_angle_plotter.py -o angle.png -i angles.txt 7 37 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing multidimensional_scaling_mean..."
+test "../src/oxDNA_analysis_tools/multidimensional_scaling_mean.py minitraj.dat -o meanM -d devsM"
 
-echo ""
+echo -n "Testing output_bonds..."
+test "../src/oxDNA_analysis_tools/output_bonds.py -v energy.json input_rna minitraj.dat"
 
-echo "Testing generate_force..."
-if
-    python ../src/oxDNA_analysis_tools/generate_force.py -f gen_pairs.txt input_rna minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing pca"
+test "../src/oxDNA_analysis_tools/pca.py minitraj.dat mean.dat pca.json"
 
-echo ""
+echo -n "Testing subset_trajectory..."
+test "../src/oxDNA_analysis_tools/subset_trajectory.py minitraj.dat rna_tile.top -i index.txt sub.dat"
 
-echo "Testing minify..."
-if
-    python ../src/oxDNA_analysis_tools/minify.py -a -d2 minitraj.dat min.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
-
-echo ""
-
-echo "Testing multidimensional_scaling_mean..."
-if
-    python ../src/oxDNA_analysis_tools/multidimensional_scaling_mean.py minitraj.dat -o meanM -d devsM 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
-
-echo ""
-
-echo "Testing output_bonds..."
-if
-    python ../src/oxDNA_analysis_tools/output_bonds.py -v energy.json input_rna minitraj.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
-
-echo ""
-
-echo "Testing pca"
-if
-    python ../src/oxDNA_analysis_tools/pca.py minitraj.dat mean.dat pca.json 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
-
-echo "Testing subset_trajectory..."
-if
-    python ../src/oxDNA_analysis_tools/subset_trajectory.py minitraj.dat rna_tile.top -i index.txt sub.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
-
-echo ""
-
-echo "Testing superimpose..."
-if
-    python ../src/oxDNA_analysis_tools/superimpose.py mean.dat centroid.dat 2>&1 >/dev/null | grep -y "ERROR"
-then
-    echo "AN ERROR OCCURED"
-else
-    echo "OK"
-fi
+echo -n "Testing superimpose..."
+test "../src/oxDNA_analysis_tools/superimpose.py mean.dat centroid.dat"
 
-echo ""
+echo ""
\ No newline at end of file
diff --git a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
index 1dbc28956..2be3ab33f 100644
--- a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
+++ b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.cpp
@@ -31,7 +31,7 @@ void CGNucleicAcidsInteraction::get_settings(input_file &inp) {
 	getInputNumber(&inp, "DPS_3b_sigma", &_3b_sigma, 0);
 	getInputNumber(&inp, "DPS_3b_range", &_3b_range, 0);
 	getInputNumber(&inp, "DPS_3b_lambda", &_3b_lambda, 0);
-	getInputNumber(&inp, "DPS_mu", &_mu, 1.0);
+	getInputNumber(&inp, "DPS_tC", &_tC, 37.0);
 	getInputNumber(&inp, "DPS_dS_mod", &dS_mod, 1.0);
 	getInputNumber(&inp, "DPS_alpha_mod", &alpha_mod, 1.0);
 	getInputNumber(&inp, "DPS_bdG_threshold", &bdG_threshold, 1.0);
@@ -80,7 +80,8 @@ void CGNucleicAcidsInteraction::get_settings(input_file &inp) {
 void CGNucleicAcidsInteraction::init() {
 	_sqr_rfene = SQR(_rfene);
 	_PS_sqr_rep_rcut = pow(2. * _WCA_sigma, 2. / _PS_n);
-	_WCA_sigma_unbonded = _WCA_sigma * (6.0 / _bead_size - _3b_sigma) / 2.0;
+	// _WCA_sigma_unbonded = _WCA_sigma * (6.0 / _bead_size - _3b_sigma) / 2.0; // disabled for now
+	_WCA_sigma_unbonded = _WCA_sigma;
 
 	OX_LOG(Logger::LOG_INFO, "CGNA: WCA sigma = %lf, WCA sigma unbonded = %lf", _WCA_sigma, _WCA_sigma_unbonded);
 
@@ -631,38 +632,38 @@ void CGNucleicAcidsInteraction::allocate_particles(std::vector<BaseParticle*> &p
 }
 
 void CGNucleicAcidsInteraction::_parse_interaction_matrix() {
-	// parse the interaction matrix file
-	input_file inter_matrix_file;
-	inter_matrix_file.init_from_filename(_interaction_matrix_file);
-	const number _t37_ = 310.15;
-	const number _kB_ = 1.9872036;
-	//const number dS_mod = 1.87;
-	if(inter_matrix_file.state == ERROR) {
-		throw oxDNAException("Caught an error while opening the interaction matrix file '%s'", _interaction_matrix_file.c_str());
-	}
+        // parse the interaction matrix file
+        input_file inter_matrix_file;
+        inter_matrix_file.init_from_filename(_interaction_matrix_file);
+        const number _t37_ = 310.15;
+        const number _kB_ = 1.9872036;
+        number _mu = _t37_ / (_tC+273.15);
+        if(inter_matrix_file.state == ERROR) {
+                throw oxDNAException("Caught an error while opening the interaction matrix file '%s'", _interaction_matrix_file.c_str());
+        }
 
 	_interaction_matrix_size = _N_attractive_types + 1;
 	_3b_epsilon.resize(_interaction_matrix_size * _interaction_matrix_size, 0.);
 
-	ofstream myfile;
-	myfile.open("beta_eps_matrix.dat");
-	for(int i = 1; i <= _N_attractive_types; i++) {
-		for(int j = 1; j <= _N_attractive_types; j++) {
-			number valueH;
-			number valueS;
-			std::string keyH = Utils::sformat("dH[%d][%d]", i, j);
-			std::string keyS = Utils::sformat("dS[%d][%d]", i, j);
-			if(getInputNumber(&inter_matrix_file, keyH.c_str(), &valueH, 0) == KEY_FOUND && getInputNumber(&inter_matrix_file, keyS.c_str(), &valueS, 0) == KEY_FOUND) {
-				number beta_dG = (_mu * valueH * 1000 / _t37_ - valueS) / _kB_;
-				number beta_eps = -(beta_dG + dS_mod) / alpha_mod;
-				if(beta_dG < bdG_threshold && abs(i - j) > 2) {
-					_3b_epsilon[i + _interaction_matrix_size * j] = _3b_epsilon[j + _interaction_matrix_size * i] = beta_eps;
-					myfile << "beta_eps[" << i << "][" << j << "]=" << beta_eps << "\n";
-				}
-			}
-		}
-	}
-	myfile.close();
+        ofstream myfile;
+        myfile.open("beta_eps_matrix.dat");
+        for(int i = 1; i <= _N_attractive_types; i++) {
+                for(int j = 1; j <= _N_attractive_types; j++) {
+                        number valueH;
+                        number valueS;
+                        std::string keyH = Utils::sformat("dH[%d][%d]", i, j);
+                        std::string keyS = Utils::sformat("dS[%d][%d]", i, j);
+                        if(getInputNumber(&inter_matrix_file, keyH.c_str(), &valueH, 0) == KEY_FOUND && getInputNumber(&inter_matrix_file, keyS.c_str(), &valueS, 0) == KEY_FOUND) {
+                                number beta_dG = (_mu * valueH * 1000 / _t37_ - valueS)/_kB_;
+                                number beta_eps = -(beta_dG + dS_mod) / alpha_mod;
+                                if(beta_dG<bdG_threshold && abs(i - j) > 1) {
+                                        _3b_epsilon[i + _interaction_matrix_size * j] = _3b_epsilon[j + _interaction_matrix_size * i] = beta_eps;
+                                        myfile << "beta_eps[" << i << "][" << j << "]=" << beta_eps << "\n";
+                                }
+                        }
+                }
+        }
+        myfile.close();
 }
 
 void CGNucleicAcidsInteraction::read_topology(int *N_strands, std::vector<BaseParticle*> &particles) {
diff --git a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
index 3ec69c388..091d0e3b4 100644
--- a/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
+++ b/contrib/tostiguerra/src/CGNucleicAcidsInteraction.h
@@ -48,7 +48,7 @@ class CGNucleicAcidsInteraction: public BaseInteraction {
 	number _sqr_rfene;
 	number _WCA_sigma = 1.0, _WCA_sigma_unbonded;
 	number _PS_sqr_rep_rcut;
-	number _mu = 1.0;
+	number _tC = 37.0;
 	number dS_mod = 1.0;
 	number alpha_mod = 1.0;
 	number bdG_threshold = 1.0;
diff --git a/docs/docs_requirements.txt b/docs/docs_requirements.txt
index dbea5b018..4219b553a 100644
--- a/docs/docs_requirements.txt
+++ b/docs/docs_requirements.txt
@@ -1,7 +1,15 @@
-Sphinx==3.4.1
-sphinx_rtd_theme==0.5.0
-myst_parser==0.17.2
-docutils==0.16
-Jinja2<3.1
-sphinx-argparse>=0.3.1
-ipython>=7.26.0
\ No newline at end of file
+Sphinx==7.0.1
+sphinx-argparse==0.3.2
+sphinx-rtd-theme==2.0.0
+sphinxcontrib-applehelp==1.0.2
+sphinxcontrib-bibtex==2.5.0
+sphinxcontrib-devhelp==1.0.2
+sphinxcontrib-htmlhelp==2.0.0
+sphinxcontrib-jquery==4.1
+sphinxcontrib-jsmath==1.0.1
+sphinxcontrib-qthelp==1.0.3
+sphinxcontrib-serializinghtml==1.1.5
+myst_parser==2.0.0
+docutils==0.20.1
+Jinja2==3.0.3
+ipython>=7.26.0
diff --git a/docs/source/input.md b/docs/source/input.md
index 7ecc97849..c18a77735 100644
--- a/docs/source/input.md
+++ b/docs/source/input.md
@@ -71,7 +71,7 @@ These options control the behaviour of MD simulations.
 
 * `sim_type = MD|FFS_MD`: run either an MD or an FFS simulation.
 * `backend = CPU|CUDA`: MD simulations can be run either on single CPU cores or on single CUDA-enabled GPUs.
-* `backend_precision = <any>`: by default CPU simulations are run with `double` precision, CUDA with `mixed` precision (see [here](https://doi.org/10.1002/jcc.23763) for details). The CUDA backend also supports single precision (`backend_precision = float`), but we do not recommend to use it. Optionally, [by using CMake switches](install.md#CMake-options) it is possible to run CPU simulations in single precision or CUDA simulations in double precision.
+* `backend_precision = <any>`: by default CPU simulations are run with `double` precision, CUDA with `mixed` precision (see [here](https://doi.org/10.1002/jcc.23763) for details). The CUDA backend also supports single precision (`backend_precision = float`), but we do not recommend to use it. Optionally, [by using CMake switches](install.md#cmake-options) it is possible to run CPU simulations in single precision or CUDA simulations in double precision.
 * `dt`: the simulation time step. The higher this value, the longer time a simulation of a given number of time steps will correspond to. However, a value that is too large will result in numerical instabilities. Typical values range between 0.001 and 0.005.
 * `refresh_vel = <bool>`: if `true` the velocities of the particles in the initial configuration will be randomly sampled from a Boltzmann distribution corresponding to `T`. If `false`, the velocities in the `conf_file` will be used (or an error will be thrown if the `conf_file` doesn't include initialized velocities).
 * `[reset_initial_com_momentum = <bool>]`: if `true` the momentum of the centre of mass of the initial configuration will be set to 0. Defaults to `false` to enforce the reproducibility of the trajectory.
@@ -152,6 +152,7 @@ The following options control the behaviour of MC simulations.
 ## Common options for `DNA2` and `RNA2` simulations
 
 * `[dh_lambda = <float>]`: the value that lambda, which is a function of temperature (T) and salt concentration (I), should take when T = 300K and I = 1M. Defaults to the value from Debye-Huckel theory, 0.3616455.
+* `[debye_huckel_rhigh]`: the distance at which the smoothing of the Debye-Hucker repulsion begins. Defaults to three times the Debye screening length.
 * `[dh_strength = <float>]`: the value that scales the overall strength of the Debye-Huckel interaction. Defaults to 0.0543.
 * `[dh_half_charged_ends = <bool>]`: if `false`, nucleotides at the end of a strand carry a full charge, if `true` their charge is halved. Defaults to `true`.
 
diff --git a/docs/source/install.md b/docs/source/install.md
index a8182cc05..f848d193d 100644
--- a/docs/source/install.md
+++ b/docs/source/install.md
@@ -70,10 +70,6 @@ If you are on your own machine or you installed Python via Anaconda, the `-DOxpy
 * `make rovigatti` Compiles the observables and interactions in contrib/rovigatti
 * `make romano` Compiles the observables and interactions in contrib/romano
 
-### Known issues
-
-The list of known issues can be browsed online [here](https://github.com/lorenzo-rovigatti/oxDNA/issues).
-
 #### CMake compiler choice
 
 CMake searches your $PATH for compatible C and C++ compilers and uses the first ones it finds. If you want to use a different set than the default, you can override the compiler choice as follows:
@@ -105,6 +101,12 @@ cmake -DPython=1 -DPYTHON_EXECUTABLE=$HOME/miniconda3/bin/python -DPYTHON_INCLUD
   `python -m pip install ./analysis --no-build-isolation`
 * Sometimes installation will fail with `TypeError: expected string or bytes-like object`. This error is usually caused by older versions of either `oxpy` or `oxDNA-analysis-tools` somewhere in your `$PATH`. Remove them with `pip uninstall oxpy` and `pip uninstall oxDNA-analysis-tools` and try installing again.
 
+## Known issues
+
+* An `illegal instruction` is sometimes issued when the code is compiled on a CPU architecture and run on another, or when specific combinations of CPU architecture and compiler are used. Invoke CMake with `-DNATIVE_COMPILATION=Off` and re-compile the code to fix the issue.
+* A list of other known issues can be browsed online [here](https://github.com/lorenzo-rovigatti/oxDNA/issues).
+
+
 ## Using `oxpy` with old Python versions
 
 `oxpy` interfaces with oxDNA using [Pybind11](https://github.com/pybind/pybind11). In September 2023 we updated the version of pybind11 included with oxDNA from 2.2 to 2.11 due to changes to the Python C API which made older versions of Pybind11 incompatible with the current Python 3.11. This new version of pybind11 is only compatible with Python > 3.8. If, for some reason, you need to use an older version of Python 3 and cannot install a newer version in a virtual environment via, for example, [Conda](https://docs.conda.io/projects/miniconda/en/latest/miniconda-install.html), this can be done by using an older version of pybind11:
diff --git a/docs/source/observables.md b/docs/source/observables.md
index a73b5df8c..1b09d6389 100644
--- a/docs/source/observables.md
+++ b/docs/source/observables.md
@@ -106,7 +106,7 @@ Compute and print the hydrogen-bonding (HB) energy of all or selected nucleotide
 
 * `type = hb_energy`: the observable type.
 * `[pairs_file = <path>]`: an order parameter file containing the list of pairs whose total HB energy is to be computed.
-* `[base_file = <path>]`: file containing a list of nucleotides whose total HB energy is to be computed, one nucleotide per line. If both this option and `pairs_file` are set, the former is silently ignored.
+* `[bases_file = <path>]`: file containing a list of nucleotides whose total HB energy is to be computed, one nucleotide per line. If both this option and `pairs_file` are set, the former is silently ignored.
 
 ## Hydrogen bonds
 
diff --git a/examples/PARALLEL_TEMPERING/README b/examples/PARALLEL_TEMPERING/README
index c547d9476..16215caf3 100644
--- a/examples/PARALLEL_TEMPERING/README
+++ b/examples/PARALLEL_TEMPERING/README
@@ -9,7 +9,7 @@ parallel, sinchronize them, attempt exchanges, etc.
 
 
 COMPILING
-this backend requires MPI, so it only gets compiled when using "make mpi=1". In
+this backend requires MPI, so it only gets compiled when using "make MPI=1". In
 beetween compilations with and without MPI, a "make clean" is required.
 
 USAGE:
diff --git a/oxpy/CMakeLists.txt b/oxpy/CMakeLists.txt
index d91895d6c..e4a20798c 100644
--- a/oxpy/CMakeLists.txt
+++ b/oxpy/CMakeLists.txt
@@ -52,7 +52,13 @@ IF(Python)
 	    WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test
 	    COMMENT "Running oxpy tests" VERBATIM
 	)
-	add_dependencies(test test_oxpy)
+	# and one to test oat
+	add_custom_target(test_oat
+	    bash test.sh ${PYTHON_EXECUTABLE} && bash cleanup.sh
+	    WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/analysis/tests
+	    COMMENT "Running oxpy tests" VERBATIM
+	)
+	add_dependencies(test test_oxpy test_oat)
 	
 	CONFIGURE_FILE(make_install_setup.py
 		${OXPY_OUTPUT_DIR}/setup.py
diff --git a/oxpy/OxpyContext.cpp b/oxpy/OxpyContext.cpp
index d97987969..4922c1eab 100644
--- a/oxpy/OxpyContext.cpp
+++ b/oxpy/OxpyContext.cpp
@@ -38,16 +38,8 @@ void OxpyContext::exit(py::object exc_type, py::object exc_value, py::object tra
 
 	if(_print_coda) {
 		OX_LOG(Logger::LOG_NOTHING, "");
-		OX_LOG(Logger::LOG_NOTHING, R"c(Please cite these publications for any work that uses the oxDNA simulation package
-		- for the code:
-			* P. Šulc et al., J. Chem. Phys. 137, 135101 (2012)
-			* L. Rovigatti et al., J. Comput. Chem. 36, 1 (2015)
-		- for the oxDNA model:
-			* T. E. Ouldridge et al., J. Chem. Phys, 134, 085101 (2011)
-		- for the oxDNA2 model:
-			* B. E. K. Snodin et al., J. Chem. Phys. 142, 234901 (2015)
-		- for the oxRNA model:
-			* P. Šulc et al., J. Chem. Phys. 140, 235102 (2014))c");
+		OX_LOG(Logger::LOG_NOTHING, R"c(Please cite at least the following publication for any work that uses the oxDNA simulation package
+  * E. Poppleton et al., J. Open Source Softw. 8, 4693 (2023), DOI: 10.21105/joss.04693)c");
 	}
 }
 
@@ -63,7 +55,14 @@ void export_OxpyContext(py::module &m) {
                 # your code using oxpy goes here
 	)pbdoc");
 
-	context.def(py::init<bool>(), py::arg("print_coda") = true);
+	context.def(py::init<bool>(), py::arg("print_coda") = true, R"pbdoc(
+The default constructor takes an optional parameter that controls whether the final message gets printed or not.
+
+Parameters
+----------
+print_coda : bool
+    If True (default value) a message asking to cite the oxDNA JOSS paper will be printed at the end of the simulation.
+	)pbdoc");
 	context.def("__enter__", &OxpyContext::enter);
 	context.def("__exit__", &OxpyContext::exit);
 }
diff --git a/oxpy/bindings_includes/input_file.h b/oxpy/bindings_includes/input_file.h
index 58f5f7814..d75d8e188 100644
--- a/oxpy/bindings_includes/input_file.h
+++ b/oxpy/bindings_includes/input_file.h
@@ -21,6 +21,10 @@ Options can be read, set or overwritten by using square brackets::
 	my_input["sim_type"] = "VMMC"
 	print(my_input["sim_type"])
 
+Options can also be deleted::
+
+	del my_input["log_file"]
+
 To check if a specific option has been set you can use `in`::
 
 	if "debug" in my_input:
@@ -28,14 +32,19 @@ To check if a specific option has been set you can use `in`::
 	)pbdoc");
 
 	input.def(py::init<>());
+
 	input.def("__getitem__", &input_file::get_value,
 		py::arg("key") = '\0',
 		py::arg("mandatory") = 0,
 		py::arg("found") = true
 	);
+
 	input.def("__setitem__", &input_file::set_value);
+
 	input.def("__delitem__", &input_file::unset_value);
+
 	input.def("__str__", &input_file::to_string);
+
 	// this enables the use of "if 'something' in input_file"
 	input.def("__contains__", [](input_file &inp, std::string v) {
 		try {
@@ -47,6 +56,7 @@ To check if a specific option has been set you can use `in`::
 			return false;
 		}
 	});
+
 	input.def("init_from_filename", &input_file::init_from_filename, py::arg("filename"), R"pbdoc(
         Initialise the object from the input file passed as parameter.
 
@@ -55,6 +65,7 @@ To check if a specific option has been set you can use `in`::
         filename: str
             The name of the input file whence the options will be loaded.
 	)pbdoc");
+
 	input.def("get_bool", [](input_file &inp, std::string &key) {
 		bool found;
 		if (inp.true_values.find(inp.get_value(key, 0, found)) != inp.true_values.end()) {
@@ -76,6 +87,25 @@ To check if a specific option has been set you can use `in`::
 		bool
 			The boolean value of the option.
 	)pbdoc");
+
+	input.def("keys", [](input_file &inp) {
+		std::vector<std::string> keys;
+		for (auto &it : inp.keys) {
+			keys.push_back(it.first);
+		}
+		return keys;
+	}, R"pbdoc(
+		Return the list of keys stored in the input file. For instance, the following code prints all the options stored in an input file::
+
+			keys = my_input.keys()
+			for key in keys:
+				print(f"Option '{key}' is set to '{my_input[key]}'")
+
+		Returns
+		-------
+		List(str)
+			The list of keys stored in the input file.
+	)pbdoc");
 }
 
 #endif /* OXPY_BINDINGS_INCLUDES_INPUT_FILE_H_ */
diff --git a/src/Backends/MCBackend.cpp b/src/Backends/MCBackend.cpp
index 41872f670..94fa58504 100644
--- a/src/Backends/MCBackend.cpp
+++ b/src/Backends/MCBackend.cpp
@@ -77,7 +77,7 @@ void MCBackend::get_settings(input_file &inp) {
 			getInputNumber(&inp, "delta_volume", &_delta[MC_MOVE_VOLUME], 1);
 		}
 	}
-	else if(sim_type == "VMMC") {
+	else if(sim_type == "VMMC" || sim_type == "PT_VMMC") {
 		getInputNumber(&inp, "delta_translation", &_delta[MC_MOVE_TRANSLATION], 1);
 		getInputNumber(&inp, "delta_rotation", &_delta[MC_MOVE_ROTATION], 1);
 	}
diff --git a/src/Backends/PT_VMMC_CPUBackend.cpp b/src/Backends/PT_VMMC_CPUBackend.cpp
index 542a29a3c..63069f1b4 100644
--- a/src/Backends/PT_VMMC_CPUBackend.cpp
+++ b/src/Backends/PT_VMMC_CPUBackend.cpp
@@ -7,7 +7,7 @@
 
 #include "PT_VMMC_CPUBackend.h"
 
-#include "mpi.h"
+#include <mpi.h>
 
 PT_VMMC_CPUBackend::PT_VMMC_CPUBackend() :
 				VMMC_CPUBackend() {
@@ -28,8 +28,6 @@ PT_VMMC_CPUBackend::~PT_VMMC_CPUBackend() {
 }
 
 void PT_VMMC_CPUBackend::init() {
-	//VMMC_CPUBackend::init(conf_input);
-
 	if(_oxRNA_stacking) {
 		RNAInteraction *it = dynamic_cast<RNAInteraction*>(_interaction.get());
 		model = it->get_model();
@@ -43,11 +41,6 @@ void PT_VMMC_CPUBackend::init() {
 
 	_conf_filename = string(my_conf_filename);
 
-	fprintf(stderr, "REPLICA %d: reading configuration from %s\n", _my_mpi_id, my_conf_filename);
-	VMMC_CPUBackend::init();
-
-	_exchange_conf = new PT_serialized_particle_info[N()];
-
 	// check that temperatures are in order...
 	bool check2 = true;
 	if(_my_mpi_id == 0) {
@@ -63,23 +56,19 @@ void PT_VMMC_CPUBackend::init() {
 	}
 
 	// let's get our own temperature...
-	if(_npttemps != _mpi_nprocs)
+	if(_npttemps != _mpi_nprocs) {
 		throw oxDNAException("Number of PT temperatures does not match number of processes (%d != %d)", _npttemps, _mpi_nprocs);
+	}
 	_T = _pttemps[_my_mpi_id];
 	//OX_LOG(Logger::LOG_INFO, "Replica %d: Running at T=%g", _my_mpi_id, _T);
 	fprintf(stderr, "Replica %d: Running at T=%g\n", _my_mpi_id, _T);
 
-	OX_LOG(Logger::LOG_INFO, "Deleting previous interaction and creating one...");
-	//_interaction.init(_T);
-	//throw oxDNAException ("File %s, line %d: not implemented", __FILE__, __LINE__);
-	// here we should initialize the interaction to use the appropriate T
+	CONFIG_INFO->update_temperature(_T);
 
-	OX_LOG(Logger::LOG_INFO, "Deleting previous lists and creating one...");
-	//throw oxDNAException ("File %s, line %d: not implemented", __FILE__, __LINE__);
-	// here we should create our own lists...
-
-	//VMMC_CPUBackend::_compute_energy();
-	MC_CPUBackend::_compute_energy();
+	fprintf(stderr, "REPLICA %d: reading configuration from %s\n", _my_mpi_id, my_conf_filename);
+	VMMC_CPUBackend::init();
+	// N() returns no non-sense only after having called init()
+	_exchange_conf = new PT_serialized_particle_info[N()];
 
 	//fprintf (stderr, "REPLICA %d: Running at T=%g\n", _my_mpi_id, _T);
 
@@ -89,14 +78,15 @@ void PT_VMMC_CPUBackend::init() {
 	// changing filename
 	char extra[16];
 	sprintf(extra, "%d", _my_mpi_id);
-	strcat(_last_hist_file, extra);
-	strcat(_traj_hist_file, extra);
 
-	if(_reload_hist)
+	if(_reload_hist) {
 		strcat(_init_hist_file, extra);
+	}
 
 	// common weights file? if so, we have a single file
 	if(_have_us) {
+		strcat(_last_hist_file, extra);
+		strcat(_traj_hist_file, extra);
 
 		sprintf(_irresp_weights_file, "%s", _weights_file);
 
@@ -172,7 +162,6 @@ void PT_VMMC_CPUBackend::get_settings(input_file &inp) {
 			memcpy(_pttemps, tmpt, _npttemps * sizeof(double));
 		}
 		delete[] tmpt;
-		//abort ();
 	}
 
 	int tmpi = -1;
@@ -189,10 +178,12 @@ void PT_VMMC_CPUBackend::get_settings(input_file &inp) {
 
 	std::string inter_type("");
 	getInputString(&inp, "interaction_type", inter_type, 0);
-	if(inter_type.compare("DNA2") == 0)
+	if(inter_type.compare("DNA2") == 0) {
 		_oxDNA2_stacking = true;
-	if(inter_type.compare("RNA2") == 0 || inter_type.compare("RNA2") == 0)
+	}
+	if(inter_type.compare("RNA2") == 0 || inter_type.compare("RNA2") == 0) {
 		_oxRNA_stacking = true;
+	}
 }
 
 void PT_serialized_particle_info::read_from(BaseParticle *par) {
@@ -203,15 +194,6 @@ void PT_serialized_particle_info::read_from(BaseParticle *par) {
 void PT_serialized_particle_info::write_to(BaseParticle *par) {
 	par->pos = pos;
 	par->orientation = orientation;
-
-	/*
-	 par->en3 = 0; // = p->en3;
-	 par->en5 = 0; // p->en5;
-	 par->esn3 = 0; //p->esn3;
-	 par->esn5 = 0; // p->esn5;
-	 par->inclust = false;
-	 */
-
 }
 
 void PT_VMMC_CPUBackend::sim_step() {
@@ -351,8 +333,7 @@ void PT_VMMC_CPUBackend::sim_step() {
 					_exchange_energy.U_ext = buffer_energy.U_ext;
 					_exchange_energy.replica_id = buffer_energy.replica_id;
 
-					// rebuild my conf (which will now become the other guy's
-					// old one
+					// rebuild my conf (which will now become the other guy's old one
 					_rebuild_exchange_conf();
 					_rebuild_exchange_energy();
 
@@ -397,10 +378,10 @@ void PT_VMMC_CPUBackend::sim_step() {
 		}
 		fflush(stdout);
 		// debug
-		if(fabs(_T - _pttemps[_my_mpi_id]) > 1.e-7) {
+		if(fabs(_T - _pttemps[_my_mpi_id]) > 1e-7) {
 			fprintf(stderr, "DISASTRO\n\n\n");
 		}
-		// we should se the forces again if we have swapped conf
+		// we should set the forces again if we have swapped conf
 		if(_external_forces) {
 			BaseParticle *p;
 			for(int i = 0; i < N(); i++) {
@@ -410,8 +391,6 @@ void PT_VMMC_CPUBackend::sim_step() {
 		}
 
 	}
-
-	return;
 }
 
 void PT_VMMC_CPUBackend::_send_exchange_energy(int other_id) {
@@ -503,7 +482,7 @@ void PT_VMMC_CPUBackend::_rebuild_exchange_conf() {
 		}
 	}
 
-//here we reset order parameters
+	// here we reset order parameters
 	_op.reset();
 	int i, j;
 	number hpq;
@@ -522,7 +501,7 @@ void PT_VMMC_CPUBackend::_rebuild_exchange_conf() {
 
 	_op.fill_distance_parameters(_particles, _box.get());
 
-	//VMMC_CPUBackend::_update_metainfo ();
+	//VMMC_CPUBackend::_update_metainfo();
 	return;
 }
 
@@ -549,18 +528,3 @@ int PT_VMMC_CPUBackend::_MPI_receive_block_data(void *data, size_t size, int nod
 		return 1;
 	}
 }
-
-number PT_VMMC_CPUBackend::get_pt_acc() {
-	if(_my_mpi_id == (_mpi_nprocs - 1)) {
-		// I am never responsible, so by definition...
-		return 0.;
-	}
-	else {
-		return _pt_exchange_accepted / (number) _pt_exchange_tries;
-	}
-}
-
-char* PT_VMMC_CPUBackend::get_replica_info_str() {
-	sprintf(_replica_info, "%d %d", _my_mpi_id, _which_replica);
-	return _replica_info;
-}
diff --git a/src/Backends/PT_VMMC_CPUBackend.h b/src/Backends/PT_VMMC_CPUBackend.h
index 90b8ef8cc..3f076f993 100644
--- a/src/Backends/PT_VMMC_CPUBackend.h
+++ b/src/Backends/PT_VMMC_CPUBackend.h
@@ -50,8 +50,6 @@ class PT_VMMC_CPUBackend: public VMMC_CPUBackend {
 	/// keep track of who's got which replica
 	int _which_replica;
 
-	number _U_ext;
-
 	llint _pt_exchange_tries, _pt_exchange_accepted;
 
 	Weights _irresp_w;
@@ -76,19 +74,14 @@ class PT_VMMC_CPUBackend: public VMMC_CPUBackend {
 
 	bool _oxDNA2_stacking;
 	bool _oxRNA_stacking;
-	Model *model; //is necessary for temperature dependent potentials and exchange probability
+	Model *model; // necessary for temperature dependent potentials and exchange probability
 
 public:
 	PT_VMMC_CPUBackend();
 	virtual ~PT_VMMC_CPUBackend();
-	//void init(ifstream &conf_input);
-	void init();
-	void get_settings (input_file &inp);
 
-	int get_mpi_id () { return _my_mpi_id; }
-	int get_which_replica () { return _which_replica; }
-	char * get_replica_info_str ();
-	number get_pt_acc ();
+	void get_settings (input_file &inp);
+	void init();
 
 	void sim_step();
 };
diff --git a/src/Backends/SimBackend.cpp b/src/Backends/SimBackend.cpp
index 6d049461e..da4d3aece 100644
--- a/src/Backends/SimBackend.cpp
+++ b/src/Backends/SimBackend.cpp
@@ -334,6 +334,20 @@ void SimBackend::init() {
 	// initialise external forces
 	ForceFactory::instance()->make_forces(_particles, _box.get());
 
+	// now that molecules and external forces have been initialised we can check what strands can be shifted
+	// and undo the shifting done to particles that belong to those strands that shouldn't be shifted
+	std::vector<bool> printed(_molecules.size(), false); // used to print the message once per strand
+	for(auto p : _particles) {
+		if(!_molecules[p->strand_id]->shiftable()) {
+			p->pos = _box->get_abs_pos(p);
+			p->set_pos_shift(0, 0, 0);
+			if(!printed[p->strand_id] && _enable_fix_diffusion) {
+				OX_LOG(Logger::LOG_INFO, "Strand %d is not shiftable and therefore 'fix_diffusion' will not act on it", p->strand_id);
+				printed[p->strand_id] = true;
+			}
+		}
+	}
+
 	_interaction->set_box(_box.get());
 
 	_lists->init(_rcut);
@@ -362,6 +376,7 @@ LR_vector SimBackend::_read_next_binary_vector() {
 	return res;
 }
 
+// here we cannot use _molecules because it has not been initialised yet
 bool SimBackend::read_next_configuration(bool binary) {
 	double Lx, Ly, Lz;
 	// parse headers. Binary and ascii configurations have different headers, and hence
@@ -438,14 +453,8 @@ bool SimBackend::read_next_configuration(bool binary) {
 	// the following part is always carried out in double precision since we want to be able to restart from confs that have
 	// large numbers in the conf file and use float precision later
 	int k, i;
-	std::vector<int> nins(_N_strands);
-	std::vector<LR_vector> scdm(_N_strands);
-
-	// here we cannot use _molecules because it has not been initialised yet
-	for(k = 0; k < _N_strands; k++) {
-		nins[k] = 0;
-		scdm[k] = LR_vector((double) 0., (double) 0., (double) 0.);
-	}
+	std::vector<int> nins(_N_strands, 0);
+	std::vector<LR_vector> scdm(_N_strands, LR_vector((double) 0., (double) 0., (double) 0.));
 
 	i = 0;
 	std::string line;
@@ -529,18 +538,10 @@ bool SimBackend::read_next_configuration(bool binary) {
 
 	for(i = 0; i < N(); i++) {
 		BaseParticle *p = _particles[i];
-		k = p->strand_id;
-
-		LR_vector p_pos = p->pos;
 		if(_enable_fix_diffusion && !binary) {
 			// we need to manually set the particle shift so that the particle absolute position is the right one
-			LR_vector scdm_number(scdm[k].x, scdm[k].y, scdm[k].z);
-			_box->shift_particle(p, scdm_number);
-			p_pos.x -= _box->box_sides().x * (floor(scdm[k].x / _box->box_sides().x));
-			p_pos.y -= _box->box_sides().y * (floor(scdm[k].y / _box->box_sides().y));
-			p_pos.z -= _box->box_sides().z * (floor(scdm[k].z / _box->box_sides().z));
+			_box->shift_particle(p, scdm[p->strand_id]);
 		}
-		p->pos = LR_vector(p_pos.x, p_pos.y, p_pos.z);
 	}
 
 	_interaction->check_input_sanity(_particles);
@@ -664,10 +665,12 @@ void SimBackend::fix_diffusion() {
 	// change particle position and fix orientation matrix;
 	for(int i = 0; i < N(); i++) {
 		BaseParticle *p = _particles[i];
-		_box->shift_particle(p, _molecules[p->strand_id]->com);
-		p->orientation.orthonormalize();
-		p->orientationT = p->orientation.get_transpose();
-		p->set_positions();
+		if(_molecules[p->strand_id]->shiftable()) {
+			_box->shift_particle(p, _molecules[p->strand_id]->com);
+			p->orientation.orthonormalize();
+			p->orientationT = p->orientation.get_transpose();
+			p->set_positions();
+		}
 	}
 
 	number E_after = _interaction->get_system_energy(_particles, _lists.get());
diff --git a/src/Backends/VMMC_CPUBackend.cpp b/src/Backends/VMMC_CPUBackend.cpp
index 374deeea4..07766d100 100644
--- a/src/Backends/VMMC_CPUBackend.cpp
+++ b/src/Backends/VMMC_CPUBackend.cpp
@@ -97,15 +97,18 @@ void VMMC_CPUBackend::init() {
 	if(_have_us) {
 		_op.init_from_file(_op_file, _particles, N());
 		_w.init((const char *) _weights_file, &_op, _safe_weights, _default_weight);
-		if(_reload_hist)
-		_h.init(_init_hist_file, &_op, _etemps, _netemps);
-		else
-		_h.init(&_op, _etemps, _netemps);
+		if(_reload_hist) {
+			_h.init(_init_hist_file, &_op, _etemps, _netemps);
+		}
+		else {
+			_h.init(&_op, _etemps, _netemps);
+		}
 		_h.set_simtemp(_T);
 	}
 
-	if(_delta[MC_MOVE_TRANSLATION] * sqrt(3) > _verlet_skin)
-	throw oxDNAException("verlet_skin must be > delta_translation times sqrt(3) (the maximum displacement)");
+	if(_delta[MC_MOVE_TRANSLATION] * sqrt(3) > _verlet_skin) {
+		throw oxDNAException("verlet_skin must be > delta_translation times sqrt(3) (the maximum displacement)");
+	}
 
 	_vmmc_box_side = _box->box_sides()[0];
 
@@ -173,9 +176,6 @@ void VMMC_CPUBackend::init() {
 		_update_ops();
 		check_ops();
 	}
-
-	if(_overlap == true)
-	throw oxDNAException("There is an overlap in the initial configuration. Dying badly");
 }
 
 void VMMC_CPUBackend::get_settings(input_file & inp) {
@@ -219,14 +219,7 @@ void VMMC_CPUBackend::get_settings(input_file & inp) {
 		}
 	}
 
-	if(getInputBoolAsInt(&inp, "preserve_topology", &tmpi, 0) == KEY_FOUND) {
-		if(tmpi > 0) {
-			_preserve_topology = true;
-		}
-		else {
-			_preserve_topology = false;
-		}
-	}
+	getInputBool(&inp, "preserve_topology", &_preserve_topology, 0);
 
 	if(getInputBoolAsInt(&inp, "umbrella_sampling", &is_us, 0) != KEY_NOT_FOUND) {
 		if(is_us > 0) {
@@ -266,8 +259,9 @@ void VMMC_CPUBackend::get_settings(input_file & inp) {
 			// whether to print out zero entries in traj_hist and last_hist
 			if(getInputBoolAsInt(&inp, "skip_hist_zeros", &tmpi, 0) == KEY_FOUND) {
 				_skip_hist_zeros = tmpi > 0;
-				if(_skip_hist_zeros)
-				OX_LOG(Logger::LOG_INFO, "(VMMC_CPUBackend.cpp) Skipping zero entries in traj_hist and last_hist files");
+				if(_skip_hist_zeros) {
+					OX_LOG(Logger::LOG_INFO, "(VMMC_CPUBackend.cpp) Skipping zero entries in traj_hist and last_hist files");
+				}
 			}
 
 			// should we extrapolate the histogram at different
@@ -552,6 +546,7 @@ inline number VMMC_CPUBackend::build_cluster_small(movestr *moveptr, int maxsize
 				E_old = pp->en3;
 				E_pp_moved = _particle_particle_bonded_interaction_n3_VMMC(pp, qq, &stack_temp);
 				test1 = VMMC_link(E_pp_moved, E_old);
+				
 				if(test1 > _next_rand()) {
 					// prelink successful
 					store_particle(qq);
@@ -1056,9 +1051,6 @@ inline number VMMC_CPUBackend::build_cluster_cells(movestr *moveptr, int maxsize
 	}
 	*size = nclust;
 
-	if(_overlap)
-	printf("cera anche qui\n");
-
 	if(nclust > maxsize) {
 		pprime = 0.;
 		_dU = 0.;
@@ -1237,8 +1229,8 @@ inline void VMMC_CPUBackend::_move_particle(movestr *moveptr, BaseParticle *q, B
 		q->pos += moveptr->t;
 	}
 	else if(moveptr->type == MC_MOVE_ROTATION) {
-		//in this case, the translation vector is the point around which we rotate
-		BaseParticle * p_old = _particles_old[p->index];
+		// in this case, the translation vector is the point around which we rotate
+		BaseParticle *p_old = _particles_old[p->index];
 		LR_vector dr, drp;
 		if(p->strand_id == q->strand_id) {
 			dr = q->pos - p_old->pos;
@@ -1247,21 +1239,12 @@ inline void VMMC_CPUBackend::_move_particle(movestr *moveptr, BaseParticle *q, B
 			dr = _box->min_image(p_old->pos, q->pos);
 		}
 
-		// we move around the backbone site
-		//dr += moveptr->t;
-
 		drp = moveptr->R * dr;
-		//q->pos += (drp - dr); // accounting for PBC
 		q->pos = p->pos + drp;
 		q->orientation = moveptr->R * q->orientation;
 		q->orientationT = q->orientation.get_transpose();
 		q->set_positions();
 	}
-	else {
-		;
-	}
-
-	return;
 }
 
 inline void VMMC_CPUBackend::_fix_list(int p_index, int oldcell, int newcell) {
@@ -1329,14 +1312,15 @@ void VMMC_CPUBackend::sim_step() {
 
 	LR_vector tmp;
 
-	int * clust, nclust;
+	int *clust, nclust;
 	clust = new int[N()];
 
 	double oldweight, weight;
 	int windex, oldwindex;
 	oldweight = weight = 1.;
-	if(_have_us)
-	oldweight = _w.get_weight(_op.get_all_states(), &oldwindex);
+	if(_have_us) {
+		oldweight = _w.get_weight(_op.get_all_states(), &oldwindex);
+	}
 
 	// set the potential due to external forces
 	_U_ext = (number) 0.f;
@@ -1347,23 +1331,21 @@ void VMMC_CPUBackend::sim_step() {
 	}
 
 	for(int i = 0; i < N(); i++) {
-		if(_have_us)
-		_op.store();
+		if(_have_us) {
+			_op.store();
+		}
 		_dU_stack = 0.;
 
 		// seed particle;
 		int pi = (int) (drand48() * N());
 		BaseParticle *p = _particles[pi];
 
-		//select the move
-		//printf("generating move...\n");
+		// select the move
 		movestr move;
 		move.seed = pi;
 		move.type = (drand48() < 0.5) ? MC_MOVE_TRANSLATION : MC_MOVE_ROTATION;
 
-		//generate translation / rotataion
-		//LR_vector translation;
-		//LR_matrix rotation;
+		// generate translation / rotation
 		if(move.type == MC_MOVE_TRANSLATION) {
 			move.t = LR_vector(Utils::gaussian(), Utils::gaussian(), Utils::gaussian()) * _delta[MC_MOVE_TRANSLATION];
 			move.R = LR_matrix((number) 1., (number) 0., (number) 0., (number) 0., (number) 1., (number) 0., (number) 0., (number) 0., (number) 1.);
@@ -1376,8 +1358,6 @@ void VMMC_CPUBackend::sim_step() {
 			move.Rt = (move.R).get_transpose();
 			//move.t = _particles[move.seed]->int_centers[DNANucleotide::BACK] + _particles[move.seed]->pos;
 			move.t = _particles[move.seed]->int_centers[DNANucleotide::BACK];
-			if(fabs((move.t * move.t)) > 0.5)
-			printf("caca");
 		}
 		_last_move = move.type;
 
@@ -1385,10 +1365,12 @@ void VMMC_CPUBackend::sim_step() {
 		//number pprime = build_cluster(pi, _maxclust, clust, &nclust, tainted, &ntainted);
 		//printf("building cluster starting from %i...\n", move.seed);
 		number pprime;
-		if(_small_system)
-		pprime = build_cluster_small(&move, _maxclust, clust, &nclust);
-		else
-		pprime = build_cluster_cells(&move, _maxclust, clust, &nclust);
+		if(_small_system) {
+			pprime = build_cluster_small(&move, _maxclust, clust, &nclust);
+		}
+		else {
+			pprime = build_cluster_cells(&move, _maxclust, clust, &nclust);
+		}
 
 		assert(nclust >= 1);
 
@@ -1430,7 +1412,7 @@ void VMMC_CPUBackend::sim_step() {
 		//printf ("### %lf\n", _dU_stack);
 		_tries[_last_move]++;
 
-		//printf("## U: %lf dU: %lf, p': %lf, nclust: %d \n", _U, _dU, pprime, nclust);
+		// printf("## U: %lf dU: %lf, p': %lf, nclust: %d \n", _U, _dU, pprime, nclust);
 		if(_overlap == false && pprime > drand48()) {
 			if(nclust <= _maxclust)
 			_accepted[_last_move]++;
@@ -1496,18 +1478,11 @@ void VMMC_CPUBackend::sim_step() {
 
 			_overlap = false;
 
-			if(_have_us)
-			_op.restore();
-			//_op.print();
-			//printf ("rejected... checking metainfo...\n");
-			//if (_have_us) check_ops();
-			//_check_metainfo();
-			//printf ("rejected... checking metainfo: PASSED\n");
+			if(_have_us) {
+				_op.restore();
+			}
 		}
 
-		//check_overlaps();
-		//assert (_overlap == false);
-
 		/*
 		 // check ext potential
 		 number c_ext = 0.;
@@ -1525,8 +1500,9 @@ void VMMC_CPUBackend::sim_step() {
 		 // check ext potential done*/
 
 		// add to the histogram
-		if(_have_us && current_step() > _equilibration_steps)
-		_h.add(oldwindex, oldweight, _U, _U_stack, _U_ext);
+		if(_have_us && current_step() > _equilibration_steps) {
+			_h.add(oldwindex, oldweight, _U, _U_stack, _U_ext);
+		}
 
 		// reset the inclust property to the particles
 		for(int k = 0; k < nclust; k++) {
@@ -1537,29 +1513,26 @@ void VMMC_CPUBackend::sim_step() {
 	//check_ops();
 
 	delete[] clust;
-	//delete[] tainted;
 
 	// check energy for percolation
 	if(current_step() % (llint) _check_energy_every == 1) {
 		//printf ("checking energy for percolation..\n");
 		number U_from_tally = _U;
 		_compute_energy();
-		if((_U - U_from_tally) > 1.e-4)
-		throw oxDNAException("(VMMC_CPUBackend) Accumulated Energy (%g) and Energy computed from scratch (%g) don't match. Possibly percolating clusters. Your box is too small", U_from_tally, _U);
+		if((_U - U_from_tally) > 1.e-4) {
+			throw oxDNAException("(VMMC_CPUBackend) Accumulated Energy (%g) and Energy computed from scratch (%g) don't match. Possibly percolating clusters. Your box is too small", U_from_tally, _U);
+		}
 		//printf ("all ok (%g %g)... \n", U_from_tally, _U);
 	}
 
-	//get_time(&_timer, 1);
-	//process_times(&_timer);
 	_timer_move->pause();
 	_mytimer->pause();
 }
 
 void VMMC_CPUBackend::check_ops() {
-	if(!_have_us)
-	return;
-	//printf ("checking OP...\n");
-	assert(_have_us);
+	if(!_have_us) {
+		return;
+	}
 
 	int * state;
 
@@ -1610,7 +1583,6 @@ void VMMC_CPUBackend::check_ops() {
 	}
 
 	free(state);
-	return;
 }
 
 void VMMC_CPUBackend::_update_ops() {
@@ -1657,11 +1629,10 @@ void VMMC_CPUBackend::_print_pos(int id) {
 }
 
 inline void VMMC_CPUBackend::check_overlaps() {
-	int i, noverlaps;
+	int i, N_overlaps = 0;
 	//number epq;
 	BaseParticle *p, *q;
 
-	noverlaps = 0;
 	for(i = 0; i < N(); i++) {
 		p = _particles[i];
 		if(p->n3 != P_VIRTUAL) {
@@ -1669,7 +1640,7 @@ inline void VMMC_CPUBackend::check_overlaps() {
 			q = p->n3;
 			_particle_particle_bonded_interaction_n3_VMMC(p, q);
 			if(_overlap) {
-				noverlaps++;
+				N_overlaps++;
 				LR_vector rbb, r;
 				r = p->pos - q->pos;
 				rbb = r + p->int_centers[DNANucleotide::BACK] - q->int_centers[DNANucleotide::BACK];
@@ -1678,9 +1649,9 @@ inline void VMMC_CPUBackend::check_overlaps() {
 			}
 		}
 	}
-	assert(noverlaps == 0);
-	if(noverlaps > 0)
-	abort();
+	if(N_overlaps > 0) {
+		throw oxDNAException("VMMC: There is an overlap in the initial configuration.");
+	}
 }
 
 char * VMMC_CPUBackend::get_op_state_str() {
@@ -1704,8 +1675,9 @@ char * VMMC_CPUBackend::get_op_state_str() {
 void VMMC_CPUBackend::print_conf(bool reduced, bool only_last) {
 	SimBackend::print_conf(reduced, only_last);
 	if(_have_us) {
-		if(!only_last)
-		_h.print_to_file(_traj_hist_file, current_step(), false, _skip_hist_zeros);
+		if(!only_last) {
+			_h.print_to_file(_traj_hist_file, current_step(), false, _skip_hist_zeros);
+		}
 		_h.print_to_file(_last_hist_file, current_step(), true, _skip_hist_zeros);
 	}
 }
@@ -1713,8 +1685,9 @@ void VMMC_CPUBackend::print_conf(bool reduced, bool only_last) {
 void VMMC_CPUBackend::print_conf(bool only_last) {
 	SimBackend::print_conf(only_last);
 	if(_have_us) {
-		if(!only_last)
-		_h.print_to_file(_traj_hist_file, current_step(), false, _skip_hist_zeros);
+		if(!only_last) {
+			_h.print_to_file(_traj_hist_file, current_step(), false, _skip_hist_zeros);
+		}
 		_h.print_to_file(_last_hist_file, current_step(), true, _skip_hist_zeros);
 	}
 }
@@ -1764,7 +1737,6 @@ void VMMC_CPUBackend::_compute_energy() {
 	if(_overlap) {
 		throw oxDNAException("overlap found. Aborting..\n");
 	}
-
 }
 
 void VMMC_CPUBackend::_init_cells() {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0e1d8d325..f7aaccfa8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -53,6 +53,7 @@ SET(observables_SOURCES
 	Observables/AverageEnergy.cpp
 	Observables/ContactMap.cpp
 	Observables/AllVectors.cpp
+	Observables/ExternalForce.cpp
 	Observables/Configurations/Configuration.cpp
 	Observables/Configurations/BinaryConfiguration.cpp
 	Observables/Configurations/TclOutput.cpp
diff --git a/src/Interactions/DNA2Interaction.cpp b/src/Interactions/DNA2Interaction.cpp
index a811dc824..8c237cb9c 100644
--- a/src/Interactions/DNA2Interaction.cpp
+++ b/src/Interactions/DNA2Interaction.cpp
@@ -69,27 +69,25 @@ void DNA2Interaction::get_settings(input_file &inp) {
 	//OX_LOG(Logger::LOG_INFO,"dh_half_charged_ends = %s", _debye_huckel_half_charged_ends ? "true" : "false");
 
 	// lambda-factor (the dh length at T = 300K, I = 1.0)
-	float lambdafactor;
-	if(getInputFloat(&inp, "dh_lambda", &lambdafactor, 0) == KEY_FOUND) {
-		_debye_huckel_lambdafactor = (float) lambdafactor;
-	}
-	else {
+	if(getInputNumber(&inp, "dh_lambda", &_debye_huckel_lambdafactor, 0) != KEY_FOUND) {
 		_debye_huckel_lambdafactor = 0.3616455;
 	}
 
 	// the prefactor to the Debye-Huckel term
-	float prefactor;
-	if(getInputFloat(&inp, "dh_strength", &prefactor, 0) == KEY_FOUND) {
-		_debye_huckel_prefactor = (float) prefactor;
-	}
-	else {
+	if(getInputNumber(&inp, "dh_strength", &_debye_huckel_prefactor, 0) != KEY_FOUND) {
 		_debye_huckel_prefactor = 0.0543;
 	}
 
+	number lambda = _debye_huckel_lambdafactor * sqrt(_T / 0.1f) / sqrt(_salt_concentration);
+	// RHIGH gives the distance at which the smoothing begins
+	if(getInputNumber(&inp, "debye_huckel_rhigh", &_debye_huckel_RHIGH, 0) != KEY_FOUND) {
+		_debye_huckel_RHIGH = 3.0 * lambda;
+	}
+
 	// notify the user that major-minor grooving is switched on
 	// check whether it's set in the input file to avoid duplicate messages
-	int tmp;
-	if(_grooving && (getInputBoolAsInt(&inp, "major_minor_grooving", &tmp, 0) != KEY_FOUND)) {
+	bool tmp;
+	if(_grooving && (getInputBool(&inp, "major_minor_grooving", &tmp, 0) != KEY_FOUND)) {
 		OX_LOG(Logger::LOG_INFO, "Using different widths for major and minor grooves");
 	}
 }
@@ -118,8 +116,6 @@ void DNA2Interaction::init() {
 
 	// We wish to normalise with respect to T=300K, I=1M. 300K=0.1 s.u. so divide _T by 0.1
 	number lambda = _debye_huckel_lambdafactor * sqrt(_T / 0.1f) / sqrt(_salt_concentration);
-	// RHIGH gives the distance at which the smoothing begins
-	_debye_huckel_RHIGH = 3.0 * lambda;
 	_minus_kappa = -1.0 / lambda;
 
 	// these are just for convenience for the smoothing parameter computation
diff --git a/src/Interactions/DNAInteraction.cpp b/src/Interactions/DNAInteraction.cpp
index c4eecabaa..e2ca377ab 100644
--- a/src/Interactions/DNAInteraction.cpp
+++ b/src/Interactions/DNAInteraction.cpp
@@ -236,7 +236,7 @@ void DNAInteraction::get_settings(input_file &inp) {
 			std::string inter_type;
 			getInputString(&inp, "interaction_type", inter_type, 0);
 
-			if(inter_type.compare("NA") || inter_type.compare("NA_relax")) {
+			if(!inter_type.compare("NA") || !inter_type.compare("NA_relax")) {
 				getInputString(&inp, "seq_dep_file_DNA", _seq_filename, 1);
 			} else {
 				getInputString(&inp, "seq_dep_file", _seq_filename, 1);
diff --git a/src/Interactions/RNAInteraction.cpp b/src/Interactions/RNAInteraction.cpp
index 0e87dfff8..8fd05207a 100644
--- a/src/Interactions/RNAInteraction.cpp
+++ b/src/Interactions/RNAInteraction.cpp
@@ -54,7 +54,7 @@ void RNAInteraction::get_settings(input_file &inp) {
 		if(!_average) {
 			std::string inter_type;
 			getInputString(&inp, "interaction_type", inter_type, 0);
-			if(inter_type.compare("NA") || inter_type.compare("NA_relax")) {
+			if(!inter_type.compare("NA") || !inter_type.compare("NA_relax")) {
 				getInputString(&inp, "seq_dep_file_RNA", _seq_filename, 1);
 			} else {
 				getInputString(&inp, "seq_dep_file", _seq_filename, 1);
diff --git a/src/Interactions/RNAInteraction2.cpp b/src/Interactions/RNAInteraction2.cpp
index f3eda9d1e..f7a94b914 100644
--- a/src/Interactions/RNAInteraction2.cpp
+++ b/src/Interactions/RNAInteraction2.cpp
@@ -32,75 +32,42 @@ number RNA2Interaction::pair_interaction_nonbonded(BaseParticle *p, BaseParticle
 void RNA2Interaction::get_settings(input_file &inp) {
 	RNAInteraction::get_settings(inp);
 
-	float salt;
-	int mismatches;
-	float prefactor; // this is the strength of the interaction
-	float lambdafactor; //Lambda is _debye_huckel_LAMBDAFACTOR / salt^0.5
-	float rh;
-	int half_charged_ends = 0;
-	//getInputString(&inp, "topology", _topology_filename, 1);
-
 	// read salt concentration from file, or set it to the default value
-	if(getInputFloat(&inp, "salt_concentration", &salt, 0) == KEY_FOUND) {
-		_salt_concentration = (float) salt;
-	}
-	else {
+	if(getInputNumber(&inp, "salt_concentration", &_salt_concentration, 0) != KEY_FOUND) {
 		_salt_concentration = 1.0f;
 	}
 
 	// read lambda-factor (the dh length at I = 1.0), or set it to the default value
-	if(getInputFloat(&inp, "dh_lambda", &lambdafactor, 0) == KEY_FOUND) {
-		_debye_huckel_lambdafactor = (float) lambdafactor;
-	}
-	else {
+	if(getInputNumber(&inp, "dh_lambda", &_debye_huckel_lambdafactor, 0) != KEY_FOUND) {
 		_debye_huckel_lambdafactor = 0.3667258;
 	}
 
 	// read the prefactor to the potential, or set it to the default value
-	if(getInputFloat(&inp, "dh_strength", &prefactor, 0) == KEY_FOUND) {
-		_debye_huckel_prefactor = (float) prefactor;
-	}
-	else {
+	if(getInputNumber(&inp, "dh_strength", &_debye_huckel_prefactor, 0) != KEY_FOUND) {
 		_debye_huckel_prefactor = 0.0858; // 0.510473f;  Q = 0.0858
 	}
-	// read the cutoff distance (todo - set the default value to something 
-	// which makes sense, maybe that depends on lambdafactor and salt concentration
+
 	number lambda = _debye_huckel_lambdafactor * sqrt(_T / 0.1f) / sqrt(_salt_concentration);
-	if(getInputFloat(&inp, "debye_huckel_rhigh", &rh, 0) == KEY_FOUND) {
-		_debye_huckel_RHIGH = (float) rh;
-	}
-	else {
+	if(getInputNumber(&inp, "debye_huckel_rhigh", &_debye_huckel_RHIGH, 0) != KEY_FOUND) {
 		_debye_huckel_RHIGH = 3.0 * lambda;
 	}
 	// read the mismatch_repulsion flag (not implemented yet)
-	if(getInputBoolAsInt(&inp, "mismatch_repulsion", &mismatches, 0) == KEY_FOUND) {
-		_mismatch_repulsion = (bool) mismatches;
-	}
-	else {
+	if(getInputBool(&inp, "mismatch_repulsion", &_mismatch_repulsion, 0) != KEY_FOUND) {
 		_mismatch_repulsion = false;
 	}
 	// whether to halve the charge on the terminus of the strand or not (default true)
-	if(getInputBoolAsInt(&inp, "dh_half_charged_ends", &half_charged_ends, 0) == KEY_FOUND) {
-		_debye_huckel_half_charged_ends = (bool) half_charged_ends;
-	}
-	else {
+	if(getInputBool(&inp, "dh_half_charged_ends", &_debye_huckel_half_charged_ends, 0) != KEY_FOUND) {
 		_debye_huckel_half_charged_ends = true;
 	}
 
 	//log it 
-	OX_LOG(Logger::LOG_INFO,"Running Debye-Huckel at salt_concentration =  %g",_salt_concentration);
+	OX_LOG(Logger::LOG_INFO,"Running Debye-Huckel at salt_concentration =  %g", _salt_concentration);
 
 	if(_mismatch_repulsion) {
-		float temp;
-		if(getInputFloat(&inp, "mismatch_repulsion_strength", &temp, 0) == KEY_FOUND) {
-			_RNA_HYDR_MIS = temp;
-		}
-		else {
+		if(getInputNumber(&inp, "mismatch_repulsion_strength", &_RNA_HYDR_MIS, 0) != KEY_FOUND) {
 			_RNA_HYDR_MIS = 1;
 		}
-
 	}
-
 }
 
 //initialise the interaction
diff --git a/src/Observables/ExternalForce.cpp b/src/Observables/ExternalForce.cpp
new file mode 100644
index 000000000..63ac29c54
--- /dev/null
+++ b/src/Observables/ExternalForce.cpp
@@ -0,0 +1,47 @@
+/*
+ * ExternalForce.cpp
+ *
+ *  Created on: Feb 2, 2024
+ *      Author: rovigatti
+ */
+
+#include "ExternalForce.h"
+
+ExternalForce::ExternalForce() {
+
+}
+
+ExternalForce::~ExternalForce() {
+
+}
+
+void ExternalForce::get_settings(input_file &my_inp, input_file &sim_inp) {
+	BaseObservable::get_settings(my_inp, sim_inp);
+
+	std::string ids;
+	if(getInputString(&my_inp, "ids", ids, 0) == KEY_FOUND) {
+		_ids = Utils::split(ids, ',');
+	}
+}
+
+void ExternalForce::update_data(llint curr_step) {
+	static bool initialised = false;
+	if(!initialised) {
+		if(_ids.size() == 0) {
+			_force_averages.resize(_config_info->forces.size());
+			_ext_forces = _config_info->forces;
+		}
+		else {
+			_force_averages.resize(_ids.size());
+			for(auto id : _ids) {
+				_ext_forces.push_back(_config_info->get_force_by_id(id));
+			}
+			_ext_forces.resize(_ids.size());
+		}
+		initialised = true;
+	}
+}
+
+std::string ExternalForce::get_output_string(llint curr_step) {
+	return Utils::sformat("TO BE IMPLEMENTED");
+}
diff --git a/src/Observables/ExternalForce.h b/src/Observables/ExternalForce.h
new file mode 100644
index 000000000..e780b66f5
--- /dev/null
+++ b/src/Observables/ExternalForce.h
@@ -0,0 +1,32 @@
+/*
+ * ExternalForce.h
+ *
+ *  Created on: Feb 02, 2024
+ *      Author: rovigatti
+ */
+
+#ifndef EXTERNALFORCE_H_
+#define EXTERNALFORCE_H_
+
+#include "BaseObservable.h"
+
+/**
+ * @brief Outputs the force applied by one or more external forces.
+ */
+
+class ExternalForce: public BaseObservable {
+protected:
+	std::vector<std::string> _ids;
+	std::vector<std::shared_ptr<BaseForce>> _ext_forces;
+	std::vector<LR_vector> _force_averages;
+public:
+	ExternalForce();
+	virtual ~ExternalForce();
+
+	virtual void get_settings(input_file &my_inp, input_file &sim_inp);
+	void update_data(llint curr_step) override;
+
+	std::string get_output_string(llint curr_step);
+};
+
+#endif /* EXTERNALFORCE_H_ */
diff --git a/src/Observables/ObservableFactory.cpp b/src/Observables/ObservableFactory.cpp
index 2bdba9c1f..da6119025 100644
--- a/src/Observables/ObservableFactory.cpp
+++ b/src/Observables/ObservableFactory.cpp
@@ -49,6 +49,7 @@
 #include "ContactMap.h"
 #include "AllVectors.h"
 #include "StressAutocorrelation.h"
+#include "ExternalForce.h"
 
 #include "Configurations/PdbOutput.h"
 #include "Configurations/ChimeraOutput.h"
@@ -115,6 +116,7 @@ ObservablePtr ObservableFactory::make_observable(input_file &obs_inp) {
 	else if(!strncasecmp(obs_type, "contact_map", 512)) res = std::make_shared<ContactMap>();
 	else if(!strncasecmp(obs_type, "all_vectors", 512)) res = std::make_shared<AllVectors>();
 	else if(!strncasecmp(obs_type, "stress_autocorrelation", 512)) res = std::make_shared<StressAutocorrelation>();
+	else if(!strncasecmp(obs_type, "external_force", 512)) res = std::make_shared<ExternalForce>();
 	else {
 		res = PluginManager::instance()->get_observable(obs_type);
 		if(res == NULL) throw oxDNAException("Observable '%s' not found. Aborting", obs_type);
diff --git a/src/Particles/Molecule.cpp b/src/Particles/Molecule.cpp
index 85e86086e..bed144918 100644
--- a/src/Particles/Molecule.cpp
+++ b/src/Particles/Molecule.cpp
@@ -7,35 +7,55 @@
 
 #include "Molecule.h"
 
-#include "../Utilities/ConfigInfo.h"
 #include "../Boxes/BaseBox.h"
+#include "../Forces/ConstantRateTorque.h"
+#include "../Forces/ConstantRateForce.h"
+#include "../Forces/MovingTrap.h"
+#include "../Utilities/ConfigInfo.h"
 
 int Molecule::_current_id = 0;
 
 Molecule::Molecule() : _id(_next_id()) {
-
 }
 
 Molecule::~Molecule() {
-
 }
 
 void Molecule::add_particle(BaseParticle *p) {
-	particles.push_back(p);
+    particles.push_back(p);
+    _shiftable_dirty = true;
 }
 
 void Molecule::normalise() {
-	throw oxDNAException("The Molecule::normalise() method has not been implemented yet!");
+    throw oxDNAException("The Molecule::normalise() method has not been implemented yet!");
 }
 
 unsigned int Molecule::N() {
-	return particles.size();
+    return particles.size();
 }
 
 void Molecule::update_com() {
-	com = LR_vector(0., 0., 0.);
-	for(auto p : particles) {
-		com += p->pos;
-	}
-	com /= N();
+    com = LR_vector(0., 0., 0.);
+    for (auto p : particles) {
+        com += p->pos;
+    }
+    com /= N();
+}
+
+bool Molecule::shiftable() {
+    if (_shiftable_dirty) {
+        _is_shiftable = true;
+        for (auto p : particles) {
+            for (auto force : p->ext_forces) {
+                if (typeid(*force) == typeid(MovingTrap) ||
+                    typeid(*force) == typeid(ConstantRateTorque) || 
+					typeid(*force) == typeid(ConstantRateForce)) {
+                    _is_shiftable = false;
+                }
+            }
+        }
+        _shiftable_dirty = false;
+    }
+
+    return _is_shiftable;
 }
diff --git a/src/Particles/Molecule.h b/src/Particles/Molecule.h
index 9d55e2fa3..7de6adf43 100644
--- a/src/Particles/Molecule.h
+++ b/src/Particles/Molecule.h
@@ -53,12 +53,18 @@ struct Molecule {
 		return _id + 1;
 	}
 
+	bool shiftable();
+
 private:
 	static int _next_id() {
 	   return _current_id++;
 	}
 	const int _id;
 	static int _current_id;
+
+	/// @brief true if the shiftable conditions should be re-evaluated
+	bool _shiftable_dirty = false;
+	bool _is_shiftable = true;
 };
 
 #endif /* SRC_PARTICLES_MOLECULE_H_ */