Add AttractivePlane Force (#114)

* Add a new Force: Attraction Plane

Implements an attraction plane that constrains particles to stay on one side of it while attracting them with a constant force (`type = attraction_plane`).  Documentation is also included.

analysis/src/oxDNA_analysis_tools/external_force_utils/forces.py

@@ -112,6 +112,25 @@ def repulsion_plane(particle:int, stiff:float, direction:List[float], position:L
     })
 
 
+def attraction_plane(particle:int, stiff:float, direction:List[float], position:List[float]) -> Dict:
+    """
+    A plane that pulls the affected particle towards it while staying on one side.
+
+    Parameters:
+        particle (int): the particle that the force acts upon.  -1 will act on whole system.
+        stiff (float): the stiffness of the trap (force = stiff * distance below plane)
+        direction ([float, float, float]): the normal vector to the plane
+        position ([float, float, float]): position of the plane (plane is d0*x + d1*y + d2*z + position = 0)
+    """
+    return({
+        "type" : "attraction_plane",
+        "particle" : particle,
+        "stiff" : stiff,
+        "dir" : direction,
+        "position" : position
+    })
+
+
 def repulsion_sphere(particle:int, center:List[float], stiff:float, r0:float, rate:float) -> Dict:
     """
     A sphere that encloses the particle

docs/source/forces.md

@@ -205,11 +205,45 @@ The following snippet defines a plane that acts on the whole system and will not
 
 ````
 
+## Attraction plane
+
+This kind of external force implements an attraction plane that attracts particles towards the plane with a constant force component while constraining them to stay on one side using a repulsive force component (compare repulsion plane).
+The repulsion is implemented as a harmonic interaction, but the stiffness can be made arbitrarily high to mimic a hard repulsion.
+The attraction is implemented as a constant force defined as: stiffness * a, where a = 1 units of length.
+Both attractive and repulsive component use the same stiffness value.
+
+A force of this kind is specified with `type = attraction_plane`. 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.
+* `stiff = <float>`: stiffness of the repulsion and strength of the attraction.
+* `dir = <float>,<float>,<float>`: the vector normal to the plane: it should point towards the half-plane where the attraction is not acting.
+* `position = <float>`: defines the position of the plane along the direction identified by the plane normal.
+
+If direction is `dir`{math}`=(u,v,w)`, then the plane contains all the points {math}`(x,y,z)` that satisfy the equation: {math}`u x + v y + w z + {\rm position} = 0`.
+Nucleotides with coordinates {math}`(x,y,z)` that satisfy {math}`u x + v y + w z + {\rm position} \ge 0` will feel a constant attraction force towards the plane equal to `stiff`.
+Nucleotides with coordinates {math}`(x,y,z)` that satisfy {math}`u x + v y + w z + {\rm position} \lt 0` will feel a repulsion force equal to `stiff` \* *D*, where {math}`D = | ux + vy + wz + \mbox{position}| / \sqrt{u^2 + v^2 + w^2 }` is the distance of the nucleotide from the plane.
+
+````{admonition} Example
+
+The following snippet defines an attraction plane that acts on on all nucleotides, pulling them towards the YZ plane (x-direction at origin).
+Particles with a positive x-coordinate are pulled towards the YZ plane by a constant force of 1 simulation unit \* x (48.6 pN \* x for DNA)
+A restoring force proportional to 1 simulation unit \* x (48.6 pN \* x for DNA) is exerted on all particles with a negative x-coordinate, constraining them at the plane.
+
+    {
+    type = attraction_plane
+    particle = -1
+    stiff = 1.00
+    dir = 1, 0, 0
+    position = 0
+    }
+
+````
+
 ## Repulsive sphere
 
 This force encloses particle(s) in a repulsive sphere. The repulsion force is harmonic.
 
-A force of this kind is specified with `type = sphere`. The relevant keys are: 
+A force of this kind is specified with `type = sphere`. 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.
 * `stiff = <float>`: stiffness of the repulsion.

docs/source/oat/forces.md

@@ -17,13 +17,15 @@ Python data structures for all forces defined by oxDNA
     oxDNA_analysis_tools.external_force_utils.forces.harmonic_trap
     oxDNA_analysis_tools.external_force_utils.forces.rotating_harmonic_trap
     oxDNA_analysis_tools.external_force_utils.forces.repulsion_plane
+    oxDNA_analysis_tools.external_force_utils.forces.attraction_plane
     oxDNA_analysis_tools.external_force_utils.forces.repulsion_sphere
 
 .. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.mutual_trap
 .. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.string
 .. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.harmonic_trap
 .. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.rotating_harmonic_trap
 .. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.repulsion_plane
+.. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.attraction_plane
 .. autofunction:: oxDNA_analysis_tools.external_force_utils.forces.repulsion_sphere
 ```
 

input_options.md

@@ -29,8 +29,8 @@ Core options:
         forces' configuration has to be provided (see external_forces_file)
     [external_forces_file = <path>]
         specifies the file containing all the external forces' configurations.
-        Currently there are six supported force types: string, twist, trap,
-        repulsion_plane, repulsion_plane_moving and mutual_trap (see
+        Currently there are several supported force types: string, twist, trap,
+        repulsion_plane, attraction_plane, repulsion_plane_moving and mutual_trap (see
         EXAMPLES/TRAPS for some examples)
     [back_in_box = <bool>]
         whether particles should be brought back into the box when a

src/Backends/SimBackend.h

@@ -43,7 +43,7 @@ class Timer;
  restart_step_counter = <boolean> (false means that the step counter will start from the value read in the configuration file, true means that the step counter will start from 0)
 
  [external_forces = <bool> (specifies whether there are external forces acting on the nucleotides or not. If it is set to 1, then a file which specifies the external forces' configuration has to be provided (see external_forces_file))]
- [external_forces_file = <path> (specifies the file containing all the external forces' configurations. Currently there are six supported force types: string, twist, trap, repulsion_plane, repulsion_plane_moving and mutual_trap (see EXAMPLES/TRAPS for some examples))]
+ [external_forces_file = <path> (specifies the file containing all the external forces' configurations. Currently there are several supported force types: string, twist, trap, repulsion_plane, attraction_plane, repulsion_plane_moving and mutual_trap (see EXAMPLES/TRAPS for some examples))]
 
  [back_in_box = <bool> (whether particles should be brought back into the box when a configuration is printed or not, defaults to false)]
 

src/CMakeLists.txt

@@ -167,6 +167,8 @@ SET(forces_SOURCES
 	Forces/Metadynamics/LT2DCOMTrap.cpp
 	Forces/Metadynamics/LTAtanCOMTrap.cpp
 	Forces/Metadynamics/meta_utils.cpp
+	Forces/AttractionPlane.cpp
+
 )
 
 SET(box_SOURCES

src/CUDA/Backends/CUDA_MD.cuh

@@ -170,6 +170,20 @@ __global__ void set_external_forces(c_number4 *poss, GPU_quat *orientations, CUD
 				}
 				break;
 			}
+			case CUDA_ATTRACTION_PLANE: {
+				number distance = extF.attractionplane.dir.x*ppos.x + extF.attractionplane.dir.y*ppos.y + extF.attractionplane.dir.z*ppos.z + extF.attractionplane.position;
+				if(distance < 0.0f) {
+					F.x += -distance * extF.attractionplane.stiff * extF.attractionplane.dir.x;
+					F.y += -distance * extF.attractionplane.stiff * extF.attractionplane.dir.y;
+					F.z += -distance * extF.attractionplane.stiff * extF.attractionplane.dir.z;
+				}
+				else{
+					F.x += -extF.attractionplane.stiff * extF.attractionplane.dir.x;
+					F.y += -extF.attractionplane.stiff * extF.attractionplane.dir.y;
+					F.z += -extF.attractionplane.stiff * extF.attractionplane.dir.z;
+				}
+				break;
+			}
 			case CUDA_REPULSION_PLANE_MOVING: {
 				for(int idx = extF.repulsionplanemoving.low_idx; idx <= extF.repulsionplanemoving.high_idx; idx++) {
 					c_number4 qpos = poss[idx];

src/CUDA/Backends/MD_CUDABackend.cu

@@ -159,6 +159,10 @@ void MD_CUDABackend::_apply_external_forces_changes() {
 					RepulsionPlane *p_force = (RepulsionPlane *) p->ext_forces[j];
 					init_RepulsionPlane_from_CPU(&cuda_force->repulsionplane, p_force);
 				}
+				else if(force_type == typeid(AttractionPlane) ) {
+					AttractionPlane *p_force = (AttractionPlane *) p->ext_forces[j];
+					init_AttractionPlane_from_CPU(&cuda_force->attractionplane, p_force);
+				}
 				else if(force_type == typeid(RepulsionPlaneMoving)) {
 					RepulsionPlaneMoving *p_force = (RepulsionPlaneMoving *) p->ext_forces[j];
 					init_RepulsionPlaneMoving_from_CPU(&cuda_force->repulsionplanemoving, p_force);

src/CUDA/CUDAForces.h

@@ -25,6 +25,7 @@
 #include "../Forces/RepulsiveEllipsoid.h"
 #include "../Forces/YukawaSphere.h"
 #include "../Forces/Metadynamics/LTCOMTrap.h"
+#include "../Forces/AttractionPlane.h"
 
 #include "CUDAUtils.h"
 
@@ -45,6 +46,8 @@
 #define CUDA_COM_FORCE 13
 #define CUDA_LR_COM_TRAP 14
 #define CUDA_YUKAWA_SPHERE 15
+#define CUDA_ATTRACTION_PLANE 16
+
 
 /**
  * @brief CUDA version of a ConstantRateForce.
@@ -151,6 +154,25 @@ void init_RepulsionPlane_from_CPU(repulsion_plane *cuda_force, RepulsionPlane *c
 	cuda_force->dir = make_float3(cpu_force->_direction.x, cpu_force->_direction.y, cpu_force->_direction.z);
 }
 
+/**
+ * @brief CUDA version of an AttractionPlane.
+ */
+struct attraction_plane {
+	int type;
+	c_number stiff;
+	c_number position;
+	float3 dir;
+};
+
+void init_AttractionPlane_from_CPU(attraction_plane *cuda_force, AttractionPlane *cpu_force) {
+	cuda_force->type = CUDA_ATTRACTION_PLANE;
+	cuda_force->stiff = cpu_force->_stiff;
+	cuda_force->position = cpu_force->_position;
+	cuda_force->dir = make_float3(cpu_force->_direction.x, cpu_force->_direction.y, cpu_force->_direction.z);
+}
+
+
+
 /**
  * @brief CUDA version of a RepulsionPlaneMoving.
  */
@@ -474,6 +496,7 @@ union CUDA_trap {
 	COM_force comforce;
 	lt_com_trap ltcomtrap;
 	Yukawa_sphere yukawasphere;
+	attraction_plane attractionplane;
 };
 
 #endif /* CUDAFORCES_H_ */

src/Forces/AttractionPlane.cpp

@@ -0,0 +1,63 @@
+/**
+ * @file    AttractionPlane.cpp
+ * @date    01/aug/2024
+ * @author  Matthies, Tilibit 
+ *
+ */
+
+#include "AttractionPlane.h"
+#include "../Utilities/oxDNAException.h"
+#include "../Particles/BaseParticle.h"
+
+AttractionPlane::AttractionPlane() :
+				BaseForce() {
+	_position = -1.;
+}
+
+
+std::tuple<std::vector<int>, std::string> AttractionPlane::init(input_file &inp) {
+	BaseForce::init(inp);
+
+	std::string particles_string;
+	getInputString(&inp, "particle", particles_string, 1);
+
+	getInputNumber(&inp, "stiff", &_stiff, 1);
+	getInputNumber(&inp, "position", &_position, 1);
+
+	int tmpi;
+	double tmpf[3];
+	std::string strdir;
+	getInputString(&inp, "dir", strdir, 1);
+	tmpi = sscanf(strdir.c_str(), "%lf,%lf,%lf", tmpf, tmpf + 1, tmpf + 2);
+	if(tmpi != 3) {
+		throw oxDNAException("Could not parse dir %s in external forces file. Aborting", strdir.c_str());
+	}
+	_direction = LR_vector((number) tmpf[0], (number) tmpf[1], (number) tmpf[2]);
+	_direction.normalize();
+
+	auto particle_ids = Utils::get_particles_from_string(CONFIG_INFO->particles(), particles_string, "AttractionPlane");
+	std::string description = Utils::sformat("AttractionPlane (stiff=%g, position=%g, dir=%g,%g,%g", _stiff, _position, _direction.x, _direction.y, _direction.z);
+
+	return std::make_tuple(particle_ids, description);
+}
+
+LR_vector AttractionPlane::value(llint step, LR_vector &pos) {
+	number distance_from_plane = this->_direction*pos + this->_position;
+
+	if(distance_from_plane >=  0.) {
+		number force = -this->_stiff * 1.0;  // constant times * unit of length
+		return force * this->_direction;
+	}
+	else
+		return -(distance_from_plane*this->_stiff)*this->_direction;
+}
+
+number AttractionPlane::potential(llint step, LR_vector &pos) {
+	number distance_from_plane = this->_direction*pos + this->_position;
+
+	if(distance_from_plane >= 0.) 
+		return  (number) (this->_stiff*distance_from_plane);
+	else 
+		return (number) (0.5*this->_stiff*SQR(distance_from_plane));
+}
+

src/Forces/AttractionPlane.h

@@ -0,0 +1,50 @@
+/**
+ * @file    AttractionPlane.h
+ * @date    01/aug/2024
+ * @author  Matthies, Tilibit 
+ *
+ */
+
+#ifndef ATTRACTIONPLANE_H_
+#define ATTRACTIONPLANE_H_
+
+#include "BaseForce.h"
+
+/**
+ * @brief External force field that attracts particles towards a plane.
+ *
+ * The definition of the plane is given by dir * (x,y,z) + position = 0.
+ * Example section in the external forces file:
+
+\n
+{\n
+  particle = -1    # acts on all particles\n
+  type = attraction_plane\n
+  position = 7.    # in simulation unit lengths\n
+  stiff = 50.      # quite stiff. Good for MC, not for MD\n
+  dir = 0.,0.,1.   # points towards positive z\n
+}\n\n
+
+ * @verbatim
+stiff = <float> (stiffness of the attraction.)
+dir = <float>,<float>,<float> (the vector normal to the plane: it should point towards the half-plane where the repulsion is not acting.)
+position = <float> (defines the position of the plane along the direction identified by the plane normal.)
+particle = <int> (index of the particle on which the force shall be applied. If -1, the force will be exerted on all the particles.)
+@endverbatim
+*/
+
+class AttractionPlane : public BaseForce {
+public:
+	number _position;
+
+	AttractionPlane ();
+	virtual ~AttractionPlane() {
+	}
+
+	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 // ATTRACTIONPLANE_H

src/Forces/ForceFactory.cpp

@@ -26,6 +26,8 @@
 #include "LJCone.h"
 #include "RepulsiveEllipsoid.h"
 #include "YukawaSphere.h"
+#include "AttractionPlane.h"
+
 
 // metadynamics-related forces
 #include "Metadynamics/LT2DCOMTrap.h"
@@ -70,6 +72,7 @@ void ForceFactory::add_force(input_file &inp, std::vector<BaseParticle *> &parti
 	else if(type_str.compare("twist") == 0) extF = std::make_shared<ConstantRateTorque>();
 	else if(type_str.compare("trap") == 0) extF = std::make_shared<MovingTrap>();
 	else if(type_str.compare("repulsion_plane") == 0) extF = std::make_shared<RepulsionPlane>();
+	else if(type_str.compare("attraction_plane") == 0) extF = std::make_shared<AttractionPlane>();
 	else if(type_str.compare("repulsion_plane_moving") == 0) extF = std::make_shared<RepulsionPlaneMoving>();
 	else if(type_str.compare("mutual_trap") == 0) extF = std::make_shared<MutualTrap>();
 	else if(type_str.compare("lowdim_trap") == 0) extF = std::make_shared<LowdimMovingTrap>();