added stiff_rate to mutual_trap and updated documentation to include …

…both rate and stiff_rate

analysis/src/oxDNA_analysis_tools/db_to_force.py

@@ -50,7 +50,7 @@ def parse_dot_bracket(input:str) -> np.ndarray:
 
     return output
 
-def db_to_forcelist(db_str:str, stiff:float, reverse:bool) -> List[Dict]:
+def db_to_forcelist(db_str:str, stiff:float, reverse:bool, r0:float=1.2, PBC:bool=True, rate:float=0, stiff_rate:float=0) -> List[Dict]:
     """
         Convert a dot-bracket string to oxDNA mutual traps
 
@@ -82,7 +82,7 @@ def db_to_forcelist(db_str:str, stiff:float, reverse:bool) -> List[Dict]:
     #p is particle id, q is paired particle id
     for p, q in enumerate(db_idx):
         if q != -1:
-            force_list.append(mutual_trap(p, q, stiff, 1.2, True))
+            force_list.append(mutual_trap(p, q, stiff, r0, True, rate=rate, stiff_rate=stiff_rate))
 
     return force_list
 

analysis/src/oxDNA_analysis_tools/external_force_utils/forces.py

@@ -1,7 +1,7 @@
 #see: https://dna.physics.ox.ac.uk/index.php/Documentation#External_Forces
 from typing import Dict, List, Literal
 
-def mutual_trap(particle:int, ref_particle:int, stiff:float, r0:float, PBC:bool) -> Dict:
+def mutual_trap(particle:int, ref_particle:int, stiff:float, r0:float, PBC:bool, rate:float=0, stiff_rate:float=0) -> Dict:
     """
     A spring force that pulls a particle towards the position of another particle
 
@@ -11,13 +11,17 @@ def mutual_trap(particle:int, ref_particle:int, stiff:float, r0:float, PBC:bool)
         stiff (float): the force constant of the spring (in simulation units)
         r0 (float): the equlibrium distance of the spring
         PBC (0 or 1): does the force calculation take PBC into account (almost always 1)
+        rate (float): changes r0 by this much every time step
+        stiff_rate (float): changes stiff by this much every time step
     """
     return({
         "type" : "mutual_trap",
         "particle" : particle,
         "ref_particle" : ref_particle,
-        "stiff" : stiff, 
+        "stiff" : stiff,
+        "stiff_rate" : stiff_rate, 
         "r0" : r0,
+        "rate" : rate,
         "PBC" : int(PBC)
     })
 

docs/source/forces.md

@@ -77,15 +77,22 @@ The following bit of code will create an external force on the first nucleotide
 This force is useful to form initial configurations. It is a harmonic force that at every moment pulls a particle towards a reference particle. It is possible to specify the separation at which the force will be 0.
 
 ````{warning}
-Please note that the reference particle (`ref_particle` below) will not feel any force, thus making the name mutual trap somewhat misleading. If you want to have an actual *mutual* trap you will need to add a force on the reference particle.
+Please note that the reference particle (`ref_particle` below) will not feel any force, thus making the name mutual trap somewhat misleading. If you want to have an actual *mutual* trap you will need to add a corresponding force on the reference particle.
 ````
 
 A force of this kind is specified with `type = mutual_trap`. The relevant keys are: 
 
 * `particle = <int>`: the particle on which to exert the force.
 * `ref_particle = <int>`: particle to pull towards.
 * `stiff = <float>`: stiffness of the trap.
-* `r0 = <float`: equilibrium distance of the trap.
+* `r0 = <float>`: equilibrium distance of the trap.
+* `PBC = <bool>`: (default: 1) If 0, calculate the distance between particles without considering periodic boundary conditions.
+* `rate = <float>`: change `r0` by this much every time step.
+* `stiff_rate = <float>`: change `stiff` by this much every time step.
+
+````{warning}
+PBC should almost always be 1. The only common exception is if you are simulating a single long strand where you want to pull the ends together.
+````
 
 ````{admonition} Example
 

src/CUDA/Backends/CUDA_MD.cuh

@@ -134,7 +134,7 @@ __global__ void set_external_forces(c_number4 *poss, GPU_quat *orientations, CUD
 				c_number4 dr = (extF.mutual.PBC) ? box->minimum_image(ppos, qpos) : qpos - ppos;
 				c_number dr_abs = _module(dr);
 
-				c_number4 force = dr * ((dr_abs - (extF.mutual.r0 + extF.mutual.rate * step)) * extF.mutual.stiff / dr_abs);
+				c_number4 force = dr * ((dr_abs - (extF.mutual.r0 + extF.mutual.rate * step)) * (extF.mutual.stiff + (step * extF.mutual.stiff_rate)) / dr_abs);
 
 				F.x += force.x;
 				F.y += force.y;

src/CUDA/CUDAForces.h

@@ -75,12 +75,14 @@ struct mutual_trap {
 	int p_ind;
 	bool PBC;
 	c_number rate;
+	c_number stiff_rate;
 };
 
 void init_MutualTrap_from_CPU(mutual_trap *cuda_force, MutualTrap *cpu_force) {
 	cuda_force->type = CUDA_TRAP_MUTUAL;
 	cuda_force->rate = cpu_force->_rate;
 	cuda_force->stiff = cpu_force->_stiff;
+	cuda_force->stiff_rate = cpu_force->_stiff_rate;
 	cuda_force->r0 = cpu_force->_r0;
 	cuda_force->p_ind = cpu_force->_p_ptr->index;
 	cuda_force->PBC = cpu_force->PBC;

src/Forces/MutualTrap.cpp

@@ -15,6 +15,8 @@ MutualTrap::MutualTrap() :
 	_particle = -2;
 	_p_ptr = NULL;
 	_r0 = -1.;
+	_rate = -1;
+	_stiff_rate = -1;
 	PBC = false;
 }
 
@@ -28,6 +30,8 @@ std::tuple<std::vector<int>, std::string> MutualTrap::init(input_file &inp) {
 	getInputBool(&inp, "PBC", &PBC, 0);
 	_rate = 0.f; //default rate is 0
 	getInputNumber(&inp, "rate", &_rate, 0);
+	_stiff_rate = 0.f; //default stiff_rate is 0
+	getInputNumber(&inp, "stiff_rate", &_stiff_rate, 0);
 
 	int N = CONFIG_INFO->particles().size();
 	if(_ref_id < 0 || _ref_id >= N) {
@@ -42,7 +46,7 @@ std::tuple<std::vector<int>, std::string> MutualTrap::init(input_file &inp) {
 		throw oxDNAException("Cannot apply MutualTrap to all particles. Aborting");
 	}
 
-	std::string description = Utils::sformat("MutualTrap (stiff=%g, rate=%g, r0=%g, ref_particle=%d, PBC=%d)", _stiff, _rate, _r0, _ref_id, PBC);
+	std::string description = Utils::sformat("MutualTrap (stiff=%g, stiff_rate=%g, r0=%g, rate=%g, ref_particle=%d, PBC=%d)", _stiff, _stiff_rate, _r0, _rate, _ref_id, PBC);
 
 	return std::make_tuple(std::vector<int>{_particle}, description);
 }
@@ -58,10 +62,10 @@ LR_vector MutualTrap::_distance(LR_vector u, LR_vector v) {
 
 LR_vector MutualTrap::value(llint step, LR_vector &pos) {
 	LR_vector dr = _distance(pos, CONFIG_INFO->box->get_abs_pos(_p_ptr));
-	return (dr / dr.module()) * (dr.module() - (_r0 + (_rate * step))) * _stiff;
+	return (dr / dr.module()) * (dr.module() - (_r0 + (_rate * step))) * (_stiff + (_stiff_rate * step));
 }
 
 number MutualTrap::potential(llint step, LR_vector &pos) {
 	LR_vector dr = _distance(pos, CONFIG_INFO->box->get_abs_pos(_p_ptr));
-	return pow(dr.module() - (_r0 + (_rate * step)), 2) * ((number) 0.5) * _stiff;
+	return pow(dr.module() - (_r0 + (_rate * step)), 2) * ((number) 0.5) * (_stiff + (_stiff_rate * step));
 }

src/Forces/MutualTrap.h

@@ -18,6 +18,8 @@ class MutualTrap: public BaseForce {
 public:
 	BaseParticle * _p_ptr;
 	number _r0;
+	number _rate;
+	number _stiff_rate;
 	bool PBC;
 
 	MutualTrap();