From 0811c391df58da11fe344797f5fda6447b030252 Mon Sep 17 00:00:00 2001 From: Erik Poppleton Date: Wed, 20 Sep 2023 17:18:08 +0200 Subject: [PATCH] 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 | 4 ++-- .../external_force_utils/forces.py | 8 ++++++-- docs/source/forces.md | 11 +++++++++-- src/CUDA/Backends/CUDA_MD.cuh | 2 +- src/CUDA/CUDAForces.h | 2 ++ src/Forces/MutualTrap.cpp | 10 +++++++--- src/Forces/MutualTrap.h | 2 ++ 7 files changed, 29 insertions(+), 10 deletions(-) diff --git a/analysis/src/oxDNA_analysis_tools/db_to_force.py b/analysis/src/oxDNA_analysis_tools/db_to_force.py index 5f13d165c..b3921966b 100755 --- a/analysis/src/oxDNA_analysis_tools/db_to_force.py +++ b/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 diff --git a/analysis/src/oxDNA_analysis_tools/external_force_utils/forces.py b/analysis/src/oxDNA_analysis_tools/external_force_utils/forces.py index 123ef6b84..4771e5175 100644 --- a/analysis/src/oxDNA_analysis_tools/external_force_utils/forces.py +++ b/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) }) diff --git a/docs/source/forces.md b/docs/source/forces.md index 28f31523f..8c613dcc0 100644 --- a/docs/source/forces.md +++ b/docs/source/forces.md @@ -77,7 +77,7 @@ 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: @@ -85,7 +85,14 @@ A force of this kind is specified with `type = mutual_trap`. The relevant keys a * `particle = `: the particle on which to exert the force. * `ref_particle = `: particle to pull towards. * `stiff = `: stiffness of the trap. -* `r0 = `: equilibrium distance of the trap. +* `PBC = `: (default: 1) If 0, calculate the distance between particles without considering periodic boundary conditions. +* `rate = `: change `r0` by this much every time step. +* `stiff_rate = `: 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 diff --git a/src/CUDA/Backends/CUDA_MD.cuh b/src/CUDA/Backends/CUDA_MD.cuh index c1e9e8d67..b79ac7e2b 100644 --- a/src/CUDA/Backends/CUDA_MD.cuh +++ b/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; diff --git a/src/CUDA/CUDAForces.h b/src/CUDA/CUDAForces.h index 61c1825b9..eae13c3f7 100644 --- a/src/CUDA/CUDAForces.h +++ b/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; diff --git a/src/Forces/MutualTrap.cpp b/src/Forces/MutualTrap.cpp index ae1c4edcd..2f776fc1f 100644 --- a/src/Forces/MutualTrap.cpp +++ b/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::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::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{_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)); } diff --git a/src/Forces/MutualTrap.h b/src/Forces/MutualTrap.h index 9665a52b2..056df24a5 100644 --- a/src/Forces/MutualTrap.h +++ b/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();