diff --git a/docs/source/forces.md b/docs/source/forces.md index dd7e38515..cc4ee84ba 100644 --- a/docs/source/forces.md +++ b/docs/source/forces.md @@ -232,6 +232,35 @@ The following snippet defines a repulsive sphere that acts on the first 50 nucle ```` +## Yukawa sphere + +This force encloses particle(s) in a sphere that interacts through a purely repulsive [WCA](http://www.sklogwiki.org/SklogWiki/index.php/Weeks-Chandler-Andersen_reference_system_model) potential complemented by a Yukawa potential of the form {math}`A \exp(-r / \lambda)`, where {math}`A` is the interaction strength (which can be either positive or negative), {math}`r` is the distance between a particle and the sphere surface, {math}`\lambda` is the Debye length. The relevant keys are: + +* `particle = `: comma-separated list of indices of particles to apply the force to. A value of `-1` or `all` applies it to all particles. Entries separated by a dash "-" get expanded in a list of all the particles on a same strand comprised between the two indices. For instance, `particle = 1,2,5-7` applies the force to 1,2,5,6,7 if 5 and 7 are on the same strand. +* `radius = `: radius of the sphere. +* `debye_A = `: Strength the Yukawa interaction. +* `debye_length = `: Debye length of the Yukawa interaction. +* `WCA_epsilon = `: strength of the WCA repulsion, defaults to `1`. +* `WCA_sigma = `: diameter of the WCA repulsion, defaults to `1`. +* `WCA_n = `: exponent of the WCA repulsion, defaults to `6`. +* `[cutoff = ]`: cutoff of the interaction, defaults to four times the Debye length. +* `[center = ,,]`: centre of the sphere, defaults to `0,0,0`. + +````{admonition} Example + +The following snippet defines a Yukawa sphere that acts on all nucleotides, confining them within a sphere of radius 5 centred in {math}`(0, 0, 10)`. The WCA parameters are set to their defaults, while the Debye length and strength are `5` and `-10`, respectively. + + { + type = yukawa_sphere + radius = 5 + center = 0,0,10 + debye_length = 5 + debye_A = -10 + particle = all + } + +```` + ## `type = com` stiff = diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f7aaccfa8..c249c74f4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -161,6 +161,7 @@ SET(forces_SOURCES Forces/AlignmentField.cpp Forces/GenericCentralForce.cpp Forces/LJCone.cpp + Forces/YukawaSphere.cpp Forces/Metadynamics/LTCOMTrap.cpp Forces/Metadynamics/LTCOMAngleTrap.cpp Forces/Metadynamics/LT2DCOMTrap.cpp diff --git a/src/CUDA/Backends/CUDA_MD.cuh b/src/CUDA/Backends/CUDA_MD.cuh index 6b5031d70..00918139c 100644 --- a/src/CUDA/Backends/CUDA_MD.cuh +++ b/src/CUDA/Backends/CUDA_MD.cuh @@ -417,6 +417,36 @@ __global__ void set_external_forces(c_number4 *poss, GPU_quat *orientations, CUD } break; } + case CUDA_YUKAWA_SPHERE: { + c_number4 centre = make_c_number4(extF.yukawasphere.center.x, extF.yukawasphere.center.y, extF.yukawasphere.center.z, 0.); + + c_number radius = extF.yukawasphere.radius; + c_number epsilon = extF.yukawasphere.epsilon; + c_number sigma = extF.yukawasphere.sigma; + int WCA_n = extF.yukawasphere.WCA_n; + c_number WCA_cutoff = extF.yukawasphere.WCA_cutoff; + c_number debye_length = extF.yukawasphere.debye_length; + c_number debye_A = extF.yukawasphere.debye_A; + c_number cutoff = extF.yukawasphere.cutoff; + + c_number4 dist = box->minimum_image(centre, ppos); + c_number dist_surface = radius - _module(dist); + + if(dist_surface < cutoff) { + c_number dist_surface_sqr = SQR(dist_surface); + c_number4 direction = -dist / _module(dist); + + c_number4 force = direction * ((debye_A * exp(-dist_surface / debye_length)) * (1.f / (dist_surface * debye_length) + 1.f / (dist_surface_sqr))); + + if(dist_surface < WCA_cutoff) { + number WCA_part = pow(sigma / dist_surface, 6); + force += direction * (4 * epsilon * WCA_n * (2 * SQR(WCA_part) - WCA_part) / dist_surface); + } + F.x += force.x; + F.y += force.y; + F.z += force.z; + } + } default: { break; } diff --git a/src/CUDA/Backends/MD_CUDABackend.cu b/src/CUDA/Backends/MD_CUDABackend.cu index c270bc91e..1c2b30358 100644 --- a/src/CUDA/Backends/MD_CUDABackend.cu +++ b/src/CUDA/Backends/MD_CUDABackend.cu @@ -199,6 +199,10 @@ void MD_CUDABackend::_apply_external_forces_changes() { LTCOMTrap *p_force = (LTCOMTrap *) p->ext_forces[j]; init_LTCOMTrap_from_CPU(&cuda_force->ltcomtrap, p_force, first_time); } + else if(force_type == typeid(YukawaSphere)) { + YukawaSphere *p_force = (YukawaSphere *) p->ext_forces[j]; + init_YukawaSphere_from_CPU(&cuda_force->yukawasphere, p_force); + } else { throw oxDNAException("Only ConstantRate, MutualTrap, MovingTrap, LowdimMovingTrap, RepulsionPlane, " "RepulsionPlaneMoving, RepulsiveSphere, LJWall, ConstantRateTorque, GenericCentralForce, " diff --git a/src/CUDA/CUDAForces.h b/src/CUDA/CUDAForces.h index eae13c3f7..e5e682f80 100644 --- a/src/CUDA/CUDAForces.h +++ b/src/CUDA/CUDAForces.h @@ -23,6 +23,7 @@ #include "../Forces/GenericCentralForce.h" #include "../Forces/LJCone.h" #include "../Forces/RepulsiveEllipsoid.h" +#include "../Forces/YukawaSphere.h" #include "../Forces/Metadynamics/LTCOMTrap.h" #include "CUDAUtils.h" @@ -43,6 +44,7 @@ #define CUDA_REPULSIVE_ELLIPSOID 12 #define CUDA_COM_FORCE 13 #define CUDA_LR_COM_TRAP 14 +#define CUDA_YUKAWA_SPHERE 15 /** * @brief CUDA version of a ConstantRateForce. @@ -422,6 +424,35 @@ void init_LTCOMTrap_from_CPU(lt_com_trap *cuda_force, LTCOMTrap *cpu_force, bool CUDA_SAFE_CALL(cudaMemcpy(cuda_force->potential_grid, local_grid.data(), sizeof(c_number) * local_grid.size(), cudaMemcpyHostToDevice)); } +/** + * @brief CUDA version of a YukawaSphere. + */ +struct Yukawa_sphere { + int type; + float3 center; + c_number radius; + c_number epsilon; + c_number sigma; + int WCA_n; + c_number WCA_cutoff; + c_number debye_length; + c_number debye_A; + c_number cutoff; +}; + +void init_YukawaSphere_from_CPU(Yukawa_sphere *cuda_force, YukawaSphere *cpu_force) { + cuda_force->type = CUDA_YUKAWA_SPHERE; + cuda_force->center = make_float3(cpu_force->_center.x, cpu_force->_center.y, cpu_force->_center.z); + cuda_force->radius = cpu_force->_radius; + cuda_force->epsilon = cpu_force->_epsilon; + cuda_force->sigma = cpu_force->_sigma; + cuda_force->WCA_n = cpu_force->_WCA_n; + cuda_force->WCA_cutoff = cpu_force->_WCA_cutoff; + cuda_force->debye_length = cpu_force->_debye_length; + cuda_force->debye_A = cpu_force->_debye_A; + cuda_force->cutoff = cpu_force->_cutoff; +} + /** * @brief Used internally by CUDA classes to provide an inheritance-like mechanism for external forces. */ @@ -442,6 +473,7 @@ union CUDA_trap { repulsive_ellipsoid repulsiveellipsoid; COM_force comforce; lt_com_trap ltcomtrap; + Yukawa_sphere yukawasphere; }; #endif /* CUDAFORCES_H_ */ diff --git a/src/Forces/ForceFactory.cpp b/src/Forces/ForceFactory.cpp index 3e5606715..d31d43d43 100644 --- a/src/Forces/ForceFactory.cpp +++ b/src/Forces/ForceFactory.cpp @@ -25,6 +25,7 @@ #include "GenericCentralForce.h" #include "LJCone.h" #include "RepulsiveEllipsoid.h" +#include "YukawaSphere.h" // metadynamics-related forces #include "Metadynamics/LT2DCOMTrap.h" @@ -82,10 +83,11 @@ void ForceFactory::add_force(input_file &inp, std::vector &parti else if(type_str.compare("generic_central_force") == 0) extF = std::make_shared(); else if(type_str.compare("LJ_cone") == 0) extF = std::make_shared(); else if(type_str.compare("ellipsoid") == 0) extF = std::make_shared(); - else if (type_str.compare("meta_com_trap") == 0) extF = std::make_shared(); - else if (type_str.compare("meta_2D_com_trap") == 0) extF = std::make_shared(); - else if (type_str.compare("meta_atan_com_trap") == 0) extF = std::make_shared(); - else if (type_str.compare("meta_com_angle_trap") == 0) extF = std::make_shared(); + else if(type_str.compare("yukawa_sphere") == 0) extF = std::make_shared(); + else if(type_str.compare("meta_com_trap") == 0) extF = std::make_shared(); + else if(type_str.compare("meta_2D_com_trap") == 0) extF = std::make_shared(); + else if(type_str.compare("meta_atan_com_trap") == 0) extF = std::make_shared(); + else if(type_str.compare("meta_com_angle_trap") == 0) extF = std::make_shared(); else throw oxDNAException("Invalid force type `%s\'", type_str.c_str()); string group = string("default"); diff --git a/src/Forces/YukawaSphere.cpp b/src/Forces/YukawaSphere.cpp new file mode 100644 index 000000000..538bf42d1 --- /dev/null +++ b/src/Forces/YukawaSphere.cpp @@ -0,0 +1,93 @@ +/* + * YukawaSphere.cpp + * + * Created on: 31/may/2024 + * Author: Lorenzo + */ + +#include "YukawaSphere.h" +#include "../Utilities/oxDNAException.h" +#include "../Particles/BaseParticle.h" +#include "../Boxes/BaseBox.h" + +YukawaSphere::YukawaSphere() : + BaseForce() { + +} + +std::tuple, std::string> YukawaSphere::init(input_file &inp) { + BaseForce::init(inp); + + getInputNumber(&inp, "radius", &_radius, 1); + getInputNumber(&inp, "WCA_epsilon", &_epsilon, 0); + getInputNumber(&inp, "WCA_sigma", &_sigma, 0); + getInputNumber(&inp, "WCA_n", &_sigma, 0); + _WCA_cutoff = _sigma * pow(2., 1. / _WCA_n); + + getInputNumber(&inp, "debye_length", &_debye_length, 1); + getInputNumber(&inp, "debye_A", &_debye_A, 1); + if(getInputNumber(&inp, "cutoff", &_cutoff, 0) == KEY_NOT_FOUND) { + _cutoff = 4.0 * _debye_length; + } + + std::string particles_string; + getInputString(&inp, "particle", particles_string, 1); + + std::string strdir; + if(getInputString(&inp, "center", strdir, 0) == KEY_FOUND) { + auto center_as_vector = Utils::split_to_numbers(strdir, ","); + if(center_as_vector.size() == 3) { + _center = LR_vector(center_as_vector[0], center_as_vector[1], center_as_vector[2]); + } + else { + throw oxDNAException("Could not parse center %s in external forces file. Aborting", strdir.c_str()); + } + } + + auto particle_ids = Utils::get_particles_from_string(CONFIG_INFO->particles(), particles_string, "YukawaSphere"); + std::string description = Utils::sformat("YukawaSphere force (radius=%g, center=%g,%g,%g, eps=%g, sigma=%g, debye_length=%g, debye_A=%g, cutoff=%g)", _radius, _center.x, _center.y, _center.z, _epsilon, _sigma, _debye_length, _debye_A, _cutoff); + + return std::make_tuple(particle_ids, description); +} + +LR_vector YukawaSphere::value(llint step, LR_vector &pos) { + LR_vector dist = CONFIG_INFO->box->min_image(_center, pos); + number dist_surface = _radius - dist.module(); + if(dist_surface < 0) { + throw oxDNAException("Found a particle beyond the Yukawa sphere radius %lf. Check your initial configuration or your timestep.", _radius); + } + + LR_vector force; + if(dist_surface < _cutoff) { + number dist_surface_sqr = SQR(dist_surface); + LR_vector direction = -dist / dist.module(); + + force = direction * ((_debye_A * exp(-dist_surface / _debye_length)) * (1.0 / (dist_surface * _debye_length) + 1.0 / (dist_surface_sqr))); + + if(dist_surface < _WCA_cutoff) { + number WCA_part = pow(_sigma / dist_surface, 6); + force += direction * (4 * _epsilon * _WCA_n * (2 * SQR(WCA_part) - WCA_part) / dist_surface); + } + } + + return force; +} + +number YukawaSphere::potential(llint step, LR_vector &pos) { + number sqr_r = CONFIG_INFO->box->sqr_min_image_distance(_center, pos); + number dist_surface = _radius - sqrt(sqr_r); + if(dist_surface < 0) { + return 1e10; + } + + number energy = 0.; + if(dist_surface < _cutoff) { + energy = exp(-dist_surface / _debye_length) * _debye_A / dist_surface; + if(dist_surface < _WCA_cutoff) { + number WCA_part = pow(_sigma / dist_surface, 6); + energy += 4 * _epsilon * (SQR(WCA_part) - WCA_part) + _epsilon; + + } + } + return energy; +} diff --git a/src/Forces/YukawaSphere.h b/src/Forces/YukawaSphere.h new file mode 100644 index 000000000..49858f1d2 --- /dev/null +++ b/src/Forces/YukawaSphere.h @@ -0,0 +1,36 @@ +/** + * @file YukawaSphere.h + * @date 31/may/2024 + * @author Lorenzo + * + */ + +#ifndef YUKAWASPHERE_H_ +#define YUKAWASPHERE_H_ + +#include "BaseForce.h" + +class YukawaSphere: public BaseForce { +public: + /// center of the sphere + LR_vector _center = LR_vector(0., 0., 0.); + + number _radius; + number _epsilon = 1.0; + number _sigma = 1.0; + int _WCA_n = 6; + number _WCA_cutoff; + number _debye_length, _debye_A; + number _cutoff; + + YukawaSphere(); + virtual ~YukawaSphere() { + } + + std::tuple, std::string> init(input_file &inp) override; + + virtual LR_vector value(llint step, LR_vector &pos); + virtual number potential(llint step, LR_vector &pos); +}; + +#endif // YUKAWASPHERE_H_