Add a new force

The new Yukawa sphere can be used to enclose particle(s) in a sphere that interacts through a purely repulsive WCA potential complemented by a Yukawa potential. See docs for details.

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 = <int>`: 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 = <float>`: radius of the sphere.
+* `debye_A = <float>`: Strength the Yukawa interaction.
+* `debye_length = <float>`: Debye length of the Yukawa interaction.
+* `WCA_epsilon = <float>`: strength of the WCA repulsion, defaults to `1`.
+* `WCA_sigma = <float>`: diameter of the WCA repulsion, defaults to `1`.
+* `WCA_n = <int>`: exponent of the WCA repulsion, defaults to `6`.
+* `[cutoff = <float>]`: cutoff of the interaction, defaults to four times the Debye length.
+* `[center = <float>,<float>,<float>]`: 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 = <float>

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

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;
 			}

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, "

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_ */

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<BaseParticle *> &parti
 	else if(type_str.compare("generic_central_force") == 0) extF = std::make_shared<GenericCentralForce>();
 	else if(type_str.compare("LJ_cone") == 0) extF = std::make_shared<LJCone>();
 	else if(type_str.compare("ellipsoid") == 0) extF = std::make_shared<RepulsiveEllipsoid>();
-	else if (type_str.compare("meta_com_trap") == 0) extF = std::make_shared<LTCOMTrap>();
-	else if (type_str.compare("meta_2D_com_trap") == 0) extF = std::make_shared<LT2DCOMTrap>();
-	else if (type_str.compare("meta_atan_com_trap") == 0) extF = std::make_shared<LTAtanCOMTrap>();
-	else if (type_str.compare("meta_com_angle_trap") == 0) extF = std::make_shared<LTCOMAngleTrap>();
+	else if(type_str.compare("yukawa_sphere") == 0) extF = std::make_shared<YukawaSphere>();
+	else if(type_str.compare("meta_com_trap") == 0) extF = std::make_shared<LTCOMTrap>();
+	else if(type_str.compare("meta_2D_com_trap") == 0) extF = std::make_shared<LT2DCOMTrap>();
+	else if(type_str.compare("meta_atan_com_trap") == 0) extF = std::make_shared<LTAtanCOMTrap>();
+	else if(type_str.compare("meta_com_angle_trap") == 0) extF = std::make_shared<LTCOMAngleTrap>();
 	else throw oxDNAException("Invalid force type `%s\'", type_str.c_str());
 
 	string group = string("default");

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::vector<int>, 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;
+}

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::vector<int>, 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_