diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst
index 967f11ec242..34876f44cb2 100644
--- a/doc/sphinx/installation.rst
+++ b/doc/sphinx/installation.rst
@@ -444,8 +444,6 @@ Fluid dynamics and fluid structure interaction
- ``LB_BOUNDARIES_GPU``
-- ``AFFINITY``
-
- ``LB_ELECTROHYDRODYNAMICS`` Enables the implicit calculation of electro-hydrodynamics for charged
particles and salt ions in an electric field.
diff --git a/maintainer/configs/maxset.hpp b/maintainer/configs/maxset.hpp
index c1cda0b13b4..28ac699bb2b 100644
--- a/maintainer/configs/maxset.hpp
+++ b/maintainer/configs/maxset.hpp
@@ -77,6 +77,5 @@
#define OIF_GLOBAL_FORCES
#define OIF_LOCAL_FORCES
#define MEMBRANE_COLLISION
-#define AFFINITY
#define ADDITIONAL_CHECKS
diff --git a/src/config/features.def b/src/config/features.def
index 787cf93b90f..e554af14f6d 100644
--- a/src/config/features.def
+++ b/src/config/features.def
@@ -92,7 +92,6 @@ THOLE requires ELECTROSTATICS
/* Fluid-Structure Interactions (object in fluid) */
OIF_LOCAL_FORCES
OIF_GLOBAL_FORCES
-AFFINITY requires EXPERIMENTAL_FEATURES
MEMBRANE_COLLISION requires EXPERIMENTAL_FEATURES
/* Immersed-Boundary Bayreuth version */
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index de03ad86f82..854ff2dae48 100644
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -115,7 +115,6 @@ set(EspressoCore_SRC
nonbonded_interactions/smooth_step.cpp
nonbonded_interactions/thole.cpp
nonbonded_interactions/wca.cpp
- object-in-fluid/affinity.cpp
object-in-fluid/membrane_collision.cpp
object-in-fluid/oif_global_forces.cpp
object-in-fluid/oif_local_forces.cpp
diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp
index 04f968c234c..c6fc53365d4 100644
--- a/src/core/forces_inline.hpp
+++ b/src/core/forces_inline.hpp
@@ -56,7 +56,6 @@
#include "nonbonded_interactions/thole.hpp"
#include "nonbonded_interactions/wca.hpp"
#include "npt.hpp"
-#include "object-in-fluid/affinity.hpp"
#include "object-in-fluid/membrane_collision.hpp"
#include "object-in-fluid/oif_global_forces.hpp"
#include "object-in-fluid/oif_local_forces.hpp"
@@ -248,18 +247,6 @@ inline void add_non_bonded_pair_force(Particle &p1, Particle &p2,
torque2 = &_torque2;
#endif
- /***********************************************/
- /* bond creation and breaking */
- /***********************************************/
-
-#ifdef AFFINITY
- /* affinity potential */
- // Prevent jump to non-inlined function
- if (dist < ia_params.affinity.cut) {
- force += affinity_pair_force(p1, p2, ia_params, d, dist);
- }
-#endif
-
/***********************************************/
/* non-bonded pair potentials */
/***********************************************/
diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp
index b8faa457897..72191efb755 100644
--- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp
+++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp
@@ -148,10 +148,6 @@ static double recalc_maximal_cutoff(const IA_parameters &data) {
(data.soft_sphere.cut + data.soft_sphere.offset));
#endif
-#ifdef AFFINITY
- max_cut_current = std::max(max_cut_current, data.affinity.cut);
-#endif
-
#ifdef MEMBRANE_COLLISION
max_cut_current = std::max(max_cut_current, data.membrane.cut);
#endif
diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
index 3a3b7d5cfd8..ae87235f4c1 100644
--- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
+++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
@@ -135,17 +135,6 @@ struct SoftSphere_Parameters {
double offset = 0.0;
};
-/** affinity potential */
-struct Affinity_Parameters {
- int type = -1;
- double kappa = INACTIVE_CUTOFF;
- double r0 = INACTIVE_CUTOFF;
- double Kon = INACTIVE_CUTOFF;
- double Koff = INACTIVE_CUTOFF;
- double maxBond = INACTIVE_CUTOFF;
- double cut = INACTIVE_CUTOFF;
-};
-
/** membrane collision potential */
struct Membrane_Parameters {
double a = 0.0;
@@ -252,10 +241,6 @@ struct IA_parameters {
SoftSphere_Parameters soft_sphere;
#endif
-#ifdef AFFINITY
- Affinity_Parameters affinity;
-#endif
-
#ifdef MEMBRANE_COLLISION
Membrane_Parameters membrane;
#endif
diff --git a/src/core/object-in-fluid/CMakeLists.txt b/src/core/object-in-fluid/CMakeLists.txt
deleted file mode 100644
index 8b137891791..00000000000
--- a/src/core/object-in-fluid/CMakeLists.txt
+++ /dev/null
@@ -1 +0,0 @@
-
diff --git a/src/core/object-in-fluid/affinity.cpp b/src/core/object-in-fluid/affinity.cpp
deleted file mode 100644
index ec949aead02..00000000000
--- a/src/core/object-in-fluid/affinity.cpp
+++ /dev/null
@@ -1,54 +0,0 @@
-/*
- * Copyright (C) 2010-2019 The ESPResSo project
- * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
- * Max-Planck-Institute for Polymer Research, Theory Group
- *
- * This file is part of ESPResSo.
- *
- * ESPResSo is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * ESPResSo is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see .
- */
-/** \file
- *
- * Implementation of \ref affinity.hpp
- */
-#include "affinity.hpp"
-
-#ifdef AFFINITY
-#include "communication.hpp"
-
-#include
-
-int affinity_set_params(int part_type_a, int part_type_b, int afftype,
- double kappa, double r0, double Kon, double Koff,
- double maxBond, double cut) {
- IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b);
-
- if (!data)
- return ES_ERROR;
-
- data->affinity.type = afftype;
- data->affinity.kappa = kappa;
- data->affinity.r0 = r0;
- data->affinity.Kon = Kon;
- data->affinity.Koff = Koff;
- data->affinity.maxBond = maxBond;
- data->affinity.cut = cut;
-
- /* broadcast interaction parameters */
- mpi_bcast_ia_params(part_type_a, part_type_b);
-
- return ES_OK;
-}
-
-#endif
diff --git a/src/core/object-in-fluid/affinity.hpp b/src/core/object-in-fluid/affinity.hpp
deleted file mode 100644
index a2b048f838a..00000000000
--- a/src/core/object-in-fluid/affinity.hpp
+++ /dev/null
@@ -1,675 +0,0 @@
-/*
- * Copyright (C) 2010-2019 The ESPResSo project
- * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
- * Max-Planck-Institute for Polymer Research, Theory Group
- *
- * This file is part of ESPResSo.
- *
- * ESPResSo is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * ESPResSo is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see .
- */
-#ifndef AFFINITY_H
-#define AFFINITY_H
-
-/** \file
- * Routines to calculate the affinity force for a particle pair.
- * \ref forces.cpp
- */
-
-#include "grid.hpp"
-#include "integrate.hpp"
-#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "particle_data.hpp"
-#include "random.hpp"
-
-#ifdef AFFINITY
-
-int affinity_set_params(int part_type_a, int part_type_b, int afftype,
- double kappa, double r0, double Kon, double Koff,
- double maxBond, double cut);
-
-/** Calculate soft-sphere potential force between particle p1 and p2 */
-inline Utils::Vector3d affinity_pair_force(Particle &p1, Particle &p2,
- IA_parameters const &ia_params,
- Utils::Vector3d const &d,
- double dist) {
-
- // The affinity potential has the first argument affinity.type. This is to
- // differentiate between different implementations. For example one
- // implementation can take into account the detachment force, another not.
- int aff_type_extracted = 0;
- int period_for_output = -1;
- if (ia_params.affinity.type > 10) {
- aff_type_extracted = ia_params.affinity.type % 10;
- period_for_output = ia_params.affinity.type - aff_type_extracted;
- } else {
- aff_type_extracted = ia_params.affinity.type;
- }
-
- auto const unfolded_pos = unfolded_position(p1.r.p, p1.l.i, box_geo.length());
- auto const vec = p1.p.bond_site - unfolded_pos;
- auto const len = vec.norm();
-
- Utils::Vector3d force{};
- if (aff_type_extracted == 1) {
- /************************
- *
- * Here I can implement the affinity force.
- * I have the position of the particle - p1, and under p1.p.bond_site I
- * have the coordinate of the bond_site. Also, under d[3] I have the vector
- * towards the constraint meaning that force on p1 should be in the
- * direction of d[3].
- *
- * Algorithm:
- * 1. First check is, whether I am in the cut-off radius: ?dist <
- * affinity.cut?.
- * 2. Then I check whether there exists a bond from the current particle:
- * ?bond_site != -1?
- * 3. If yes, then I maintain the bond. I put the forces and afterwards I
- * decide whether the bond will break or not.
- * 4. If no, I maintain the creation of a bond. First I check whether I am
- * in the area of possible bond creation: ?dist < affinity.r0?
- * 5. If yes, I run the decision algorithm for bond creation and I either
- * create or does not create the bond.
- * 6. If I am not in the area of possible bond creation I do nothing
- *
- * comments:
- * - strength of the force is proportional to the difference of
- * actual bond length and relaxed bond length
- * - bond is always created, no probability is involved
- * - if bond length reaches maxBond, the bond immediately ruptures,
- * no probability is involved
- *********************/
- double fac = 0.0;
- if (dist < ia_params.affinity.cut) { // Checking whether I am inside the
- // interaction cut-off radius.
- if (dist > 0.0) {
- // printf("bond_site: %f %f
- // %f\n",p1.p.bond_site[0],p1.p.bond_site[1],p1.p.bond_site[2]);
- if ((p1.p.bond_site[0] >= 0) && (p1.p.bond_site[1] >= 0) &&
- (p1.p.bond_site[2] >= 0)) // Checking whether any bond exists
- { // Bond exists
- if (len > ia_params.affinity.r0) {
- fac = ia_params.affinity.kappa * (len - ia_params.affinity.r0);
- // printf("len %f r0 %f\n",len, ia_params.affinity.r0);
- } else {
- fac = 0.0;
- }
- // double ftemp = 0;
- force += (fac / len) * vec;
- // printf("%f ",ftemp);
- // Decision whether I should break the bond: if the bond length is
- // greater than maxBond, it breaks.
- if (len > ia_params.affinity.maxBond) {
- p1.p.bond_site = {-1, -1, -1};
- }
- } else if (dist < ia_params.affinity
- .r0) { // Bond does not exist, we are inside
- // of possible bond creation area,
- // let's talk about creating a bond
- // This implementation creates bond always
- p1.p.bond_site = unfolded_pos - d;
- }
- }
- }
- }
- if (aff_type_extracted == 2) { // second implementation of affinity
- /************************
- *
- * Here I can implement the affinity force.
- * I have the position of the particle - p1, and under p1.p.bond_site I
- * have the coordinate of the bond_site. Also, under d[3] I have the vector
- * towards the constraint meaning that force on p1 should be in the
- * direction of d[3].
- *
- * Algorithm:
- * 1. First check is whether I am in the cut-off radius: ?dist <
- * affinity.cut?.
- * 2. Then I check whether there exists a bond from the current particle:
- * ?bond_site != -1?
- * 3. If yes, then I maintain the bond. I put the forces and afterwards I
- * decide whether the bond will break or not.
- * 4. If no, I maintain the creation of a bond. First I check whether I am
- * in the area of possible bond creation: ?dist < affinity.r0?
- * 5. If yes, I run the decision algorithm for bond creation and I either
- * create or does not create the bond.
- * 6. If I am not in the area of possible bond creation I do nothing
- *
- *
- * comments:
- * - strength of the force is proportional to the difference of actual
- * bond length and relaxed bond length
- * - bond is created with probability 1-exp(-Kon*timestep)
- * - maxBond is not used, we use probability 1-exp(-Koff*timestep) to
- * break the bond
- * - Koff depends on the bond length via Koff = K0*exp(F/Fd) =
- * K0*exp(kappa(r-r0)/Fd)
- * - here, ia_params.Koff gives us K_0, off rate when bond is relaxed
- * - here, maxBond is used as detachment force F_d
- * - the original check for ensuring, that particle flows out of the
- * cut-off radius and the bond remains active is replaced with fixed
- * check, that bond length must not be greater that 0.8 cut_off
- *********************/
- double fac = 0.0;
- if (dist < ia_params.affinity.cut) { // Checking whether I am inside the
- // interaction cut-off radius.
- if (dist > 0.0) {
- // printf("bond_site: %f %f
- // %f\n",p1.p.bond_site[0],p1.p.bond_site[1],p1.p.bond_site[2]);
- if ((p1.p.bond_site[0] >= 0) && (p1.p.bond_site[1] >= 0) &&
- (p1.p.bond_site[2] >= 0)) // Checking whether any bond exists
- { // Bond exists
- if (len > ia_params.affinity.r0) {
- fac = ia_params.affinity.kappa * (len - ia_params.affinity.r0);
- // printf("len %f r0 %f\n",len, ia_params.affinity.r0);
- } else {
- fac = 0.0;
- }
- // double ftemp = 0;
- force += (fac / len) * vec;
- // ftemp += abs((fac / len) * vec);
- // if (ftemp > 0.000000000000001) printf("%f ",fac);
- // Decision whether I should break the bond:
- // First, force exerted on bond is stored in fac
- double tmpF = fac;
- // Then, zero force off rate K_0 is stored at ia_params.affinity.Koff
- double tmpK0 = ia_params.affinity.Koff;
- // Then, detachment force is stored in ia_params.affinity.maxBond
- double tmpFd = ia_params.affinity.maxBond;
- // Then, compute Koff
- double tmpKoff = tmpK0 * exp(tmpF / tmpFd);
- // Finally, compute Poff
- double Poff = 1.0 - exp(-tmpKoff * time_step);
- // printf("%f ", Poff);
- if (len <
- 0.8 *
- ia_params.affinity.cut) { // in other implementation, maxBond
- // is used here. However, in this
- // implementation, we need maxBond
- // for setting detachment force F_d
- double decide = d_random();
- if (decide < Poff) {
- p1.p.bond_site = {-1, -1, -1};
- }
- } else {
- p1.p.bond_site = {-1, -1, -1};
- // printf("breaking: out of cut");
- }
- // Checkpoint output:
- if (period_for_output > 0)
- if (((int)floor(sim_time / time_step) % period_for_output == 0) &&
- (len > ia_params.affinity.r0)) {
- FILE *fp;
- double tmpPon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- fp = fopen("affinity_check.dat", "a");
- fprintf(fp, "sim_time %f, period_for_output %d aff type: %d ",
- sim_time, period_for_output, aff_type_extracted);
- fprintf(fp,
- "Pon %f, Kon %f, particle %d, Poff = %f, F = %f, Koff = "
- "%f, K0 = %f, len = %f \n",
- tmpPon, ia_params.affinity.Kon, p1.p.identity, Poff, tmpF,
- tmpKoff, tmpK0, len);
- fclose(fp);
- }
- } else if (dist < ia_params.affinity
- .r0) { // Bond does not exist, we are inside
- // of possible bond creation area,
- // let's talk about creating a bond
- double Pon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- // The probability is given by function Pon(x)= 1 - e^(-x) where x is
- // Kon*dt.
- double decide = d_random();
- if (decide <
- Pon) { // the bond will be created only with probability Pon.
- p1.p.bond_site = unfolded_pos - d;
- } else {
- // printf("In range, not creating: Pon = %f, decide = %f", Pon,
- // decide);
- }
- }
- }
- }
- }
- if (aff_type_extracted == 3) {
- /************************
- *
- * Here I can implement the affinity force.
- * I have the position of the particle - p1, and under p1.p.bond_site I
- * have the coordinate of the bond_site. Also, under d[3] I have the vector
- * towards the constraint meaning that force on p1 should be in the
- * direction of d[3].
- *
- * Algorithm:
- * 1. First check is whether I am in the cut-off radius: ?dist <
- * affinity.cut?.
- * 2. Then I check whether there exists a bond from the current particle:
- * ?bond_site != -1?
- * 3. If yes, then I maintain the bond. I put the forces and afterwards I
- * decide whether the bond will break or not.
- * 4. If no, I maintain the creation of a bond. First I check whether I am
- * in the area of possible bond creation: ?dist < affinity.r0?
- * 5. If yes, I run the decision algorithm for bond creation and I either
- * create or does not create the bond.
- * 6. If I am not in the area of possible bond creation I do nothing
- *
- *
- * comments:
- * - strength of the force is proportional to the difference of
- * actual bond length and relaxed bond length
- * - bond is created with probability 1-exp(-Kon*timestep)
- * - bond is ruptured with probability 1-exp(-Koff*timestep)
- * - to break the bond Koff is given as parameter, is not
- * dependent on the force nor the bond length
- * - here, maxBond stands for ensuring, that particle flows out of the
- * cut-off radius and the bond remains active
- * - maxBond should be always less than cut_off radius
- *********************/
- double fac = 0.0;
- if ((dist < ia_params.affinity.cut)) { // Checking whether I am inside
- // the interaction cut-off radius.
- if (dist > 0.0) {
- // printf("bond_site: %f %f
- // %f\n",p1.p.bond_site[0],p1.p.bond_site[1],p1.p.bond_site[2]);
- if ((p1.p.bond_site[0] >= 0) && (p1.p.bond_site[1] >= 0) &&
- (p1.p.bond_site[2] >= 0)) // Checking whether any bond exists
- { // Bond exists
- if (len > ia_params.affinity.r0) {
- fac =
- ia_params.affinity.kappa * (len - ia_params.affinity.r0) / len;
- // printf("len %f r0 %f\n",len, ia_params.affinity.r0);
- } else
- fac = 0.0;
- // double ftemp = 0;
- force += fac * vec;
- // printf("%f ",ftemp);
- // Decision whether I should break the bond:
- // The random decision algorithm is much more complicated with Fd
- // detachment force etc. Here, I use much simpler rule, the same as
- // with Kon, except that the probability of bond breakage increases
- // with prolongation of the bond. If the bond reaches
-
- double Poff = 1.0 - exp(-ia_params.affinity.Koff * time_step);
- if (len < ia_params.affinity.maxBond) {
- double decide = d_random();
- if (decide < Poff) {
- p1.p.bond_site = {-1, -1, -1};
- }
-
- } else {
- p1.p.bond_site = {-1, -1, -1};
- // printf("breaking: out of cut");
- }
- } else if (dist < ia_params.affinity
- .r0) { // Bond does not exist, we are inside
- // of possible bond creation area,
- // let's talk about creating a bond
- double Pon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- // The probability is given by function Pon(x)= 1 - e^(-x) where x is
- // Kon*dt.
- double decide = d_random();
- if (decide <
- Pon) { // the bond will be created only with probability Pon.
- p1.p.bond_site = unfolded_pos - d;
- }
- }
- }
- }
- }
- if (aff_type_extracted == 4) {
- /************************
- *
- * Here I can implement the affinity force.
- * I have the position of the particle - p1, and under p1.p.bond_site I
- * have the coordinate of the bond_site. Also, under d[3] I have the vector
- * towards the constraint meaning that force on p1 should be in the
- * direction of d[3].
- *
- * Algorithm:
- * 1. First check is whether I am in the cut-off radius: ?dist <
- * affinity.cut?.
- * 2. Then I check whether there exists a bond from the current particle:
- * ?bond_site != -1?
- * 3. If yes, then I maintain the bond. I put the forces and afterwards I
- * decide whether the bond will break or not.
- * 4. If no, I maintain the creation of a bond. First I check whether I am
- * in the area of possible bond creation: ?dist < affinity.r0?
- * 5. If yes, I run the decision algorithm for bond creation and I either
- * create or does not create the bond.
- * 6. If I am not in the area of possible bond creation I do nothing
- *
- *
- * comments:
- * - strength of the force is proportional to the actual bond length
- * - bond is created with probability 1-exp(-Kon*timestep)
- * - maxBond is not used, we use probability 1-exp(-Koff*timestep) to
- * break the bond
- * - Koff depends on the bond length via Koff = K0*exp(F/Fd) =
- * K0*exp(kappa*r/Fd)
- * - here, ia_params.Koff gives us K_0, off rate when bond is relaxed
- * - here, maxBond is used as detachment force F_d
- * - the original check for ensuring, that particle flows out of the
- * cut-off radius and the bond remains active is replaced with fixed
- * check, that bond length must not be greater that 0.8 cut_off
- *********************/
- double fac = 0.0;
- if (dist < ia_params.affinity.cut) { // Checking whether I am inside the
- // interaction cut-off radius.
- if (dist > 0.0) {
- // printf("bond_site: %f %f
- // %f\n",p1.p.bond_site[0],p1.p.bond_site[1],p1.p.bond_site[2]);
- if ((p1.p.bond_site[0] >= 0) && (p1.p.bond_site[1] >= 0) &&
- (p1.p.bond_site[2] >= 0)) // Checking whether any bond exists
- { // Bond exists
- fac = ia_params.affinity.kappa * len;
- // double ftemp = 0;
- force += (fac / len) * vec;
- // Decision whether I should break the bond:
- // First, force exerted on bond is stored in fac
- double tmpF = fac;
- // Then, zero force off rate K_0 is stored at ia_params.affinity.Koff
- double tmpK0 = ia_params.affinity.Koff;
- // Then, detachment force is stored in ia_params.affinity.maxBond
- double tmpFd = ia_params.affinity.maxBond;
- // Then, compute Koff
- double tmpKoff = tmpK0 * exp(tmpF / tmpFd);
- // Finally, compute Poff
- double Poff = 1.0 - exp(-tmpKoff * time_step);
- // printf("%f ", Poff);
- if (len <
- 0.8 *
- ia_params.affinity.cut) { // in other implementation, maxBond
- // is used here. However, in this
- // implementation, we need maxBond
- // for setting detachment force F_d
- double decide = d_random();
- if (decide < Poff) {
- p1.p.bond_site = {-1, -1, -1};
- }
-
- } else {
- p1.p.bond_site = {-1, -1, -1};
- // printf("breaking: out of cut");
- }
- // Checkpoint output:
- if (period_for_output > 0)
- if ((int)floor(sim_time / time_step) % period_for_output == 0) {
- FILE *fp;
- double tmpPon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- fp = fopen("affinity_check.dat", "a");
- fprintf(fp, "sim_time %f, period_for_output %d aff type: %d ",
- sim_time, period_for_output, aff_type_extracted);
- fprintf(fp,
- "Pon %f, Kon %f, particle %d, Poff = %f, F = %f, Koff = "
- "%f, K0 = %f, len = %f \n",
- tmpPon, ia_params.affinity.Kon, p1.p.identity, Poff, tmpF,
- tmpKoff, tmpK0, len);
- fclose(fp);
- }
- } else if (dist < ia_params.affinity
- .r0) { // Bond does not exist, we are inside
- // of possible bond creation area,
- // let's talk about creating a bond
- double Pon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- // The probability is given by function Pon(x)= 1 - e^(-x) where x is
- // Kon*dt.
- double decide = d_random();
- if (decide <
- Pon) { // the bond will be created only with probability Pon.
- p1.p.bond_site = unfolded_pos - d;
- } else {
- // printf("In range, not creating: Pon = %f, decide = %f", Pon,
- // decide);
- }
- }
- }
- }
- }
- if (aff_type_extracted == 5) { // second implementation of affinity
- /************************
- *
- * Here I can implement the affinity force.
- * I have the position of the particle - p1, and under p1.p.bond_site I
- * have the coordinate of the bond_site. Also, under d[3] I have the vector
- * towards the constraint meaning that force on p1 should be in the
- * direction of d[3].
- *
- * Algorithm:
- * 1. First check is whether I am in the cut-off radius: ?dist <
- * affinity.cut?.
- * 2. Then I check whether there exists a bond from the current particle:
- * ?bond_site != -1?
- * 3. If yes, then I maintain the bond. I put the forces and afterwards I
- * decide whether the bond will break or not.
- * 4. If no, I maintain the creation of a bond. First I check whether I am
- * in the area of possible bond creation: ?dist < affinity.r0?
- * 5. If yes, I run the decision algorithm for bond creation and I either
- * create or does not create the bond.
- * 6. If I am not in the area of possible bond creation I do nothing
- *
- *
- * comments:
- * - strength of the force is proportional to the difference of actual
- * bond length and 75% of the relaxed bond length
- * - bond is created with probability 1-exp(-Kon*timestep)
- * - maxBond is not used, we use probability 1-exp(-Koff*timestep) to
- * break the bond
- * - Koff depends on the bond length via Koff = K0*exp(F/Fd) =
- * K0*exp(kappa(r-0.75*r0)/Fd)
- * - here, ia_params.Koff gives us K_0, off rate when bond is relaxed
- * - here, maxBond is used as detachment force F_d
- * - the original check for ensuring, that particle flows out of the
- * cut-off radius and the bond remains active is replaced with fixed
- * check, that bond length must not be greater that 0.8 cut_off
- *********************/
- double fac = 0.0;
- if (dist < ia_params.affinity.cut) { // Checking whether I am inside the
- // interaction cut-off radius.
- if (dist > 0.0) {
- // printf("bond_site: %f %f
- // %f\n",p1.p.bond_site[0],p1.p.bond_site[1],p1.p.bond_site[2]);
- if ((p1.p.bond_site[0] >= 0) && (p1.p.bond_site[1] >= 0) &&
- (p1.p.bond_site[2] >= 0)) // Checking whether any bond exists
- { // Bond exists
- if (len > 0.75 * (ia_params.affinity.r0)) {
- fac = ia_params.affinity.kappa *
- (len - 0.75 * (ia_params.affinity.r0));
- // printf("len %f r0 %f\n",len, ia_params.affinity.r0);
- } else
- fac = 0.0;
- // double ftemp = 0;
- force += (fac / len) * vec;
- // if (ftemp > 0.000000000000001) printf("%f ",fac);
- // Decision whether I should break the bond:
- // First, force exerted on bond is stored in fac
- double tmpF = fac;
- // Then, zero force off rate K_0 is stored at ia_params.affinity.Koff
- double tmpK0 = ia_params.affinity.Koff;
- // Then, detachment force is stored in ia_params.affinity.maxBond
- double tmpFd = ia_params.affinity.maxBond;
- // Then, compute Koff
- double tmpKoff = tmpK0 * exp(tmpF / tmpFd);
- // Finally, compute Poff
- double Poff = 1.0 - exp(-tmpKoff * time_step);
- // printf("%f ", Poff);
- if (len <
- 0.8 *
- ia_params.affinity.cut) { // in other implementation, maxBond
- // is used here. However, in this
- // implementation, we need maxBond
- // for setting detachment force F_d
- double decide = d_random();
- if (decide < Poff) {
- p1.p.bond_site = {-1, -1, -1};
- }
-
- } else {
- p1.p.bond_site = {-1, -1, -1};
- // printf("breaking: out of cut");
- }
- // Checkpoint output:
- if (period_for_output > 0)
- if (((int)floor(sim_time / time_step) % period_for_output == 0) &&
- (len > ia_params.affinity.r0)) {
- FILE *fp;
- double tmpPon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- fp = fopen("affinity_check.dat", "a");
- fprintf(fp, "sim_time %f, period_for_output %d aff type: %d ",
- sim_time, period_for_output, aff_type_extracted);
- fprintf(fp,
- "Pon %f, Kon %f, particle %d, Poff = %f, F = %f, Koff = "
- "%f, K0 = %f, len = %f \n",
- tmpPon, ia_params.affinity.Kon, p1.p.identity, Poff, tmpF,
- tmpKoff, tmpK0, len);
- fclose(fp);
- }
- } else if (dist < ia_params.affinity
- .r0) { // Bond does not exist, we are inside
- // of possible bond creation area,
- // let's talk about creating a bond
- double Pon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- // The probability is given by function Pon(x)= 1 - e^(-x) where x is
- // Kon*dt.
- double decide = d_random();
- if (decide <
- Pon) { // the bond will be created only with probability Pon.
- p1.p.bond_site = unfolded_pos - d;
- } else {
- // printf("In range, not creating: Pon = %f, decide = %f", Pon,
- // decide);
- }
- }
- }
- }
- }
- if (aff_type_extracted == 6) { // second implementation of affinity
- /************************
- *
- * Here I can implement the affinity force.
- * I have the position of the particle - p1, and under p1.p.bond_site I
- * have the coordinate of the bond_site. Also, under d[3] I have the vector
- * towards the constraint meaning that force on p1 should be in the
- * direction of d[3].
- *
- * Algorithm:
- * 1. First check is whether I am in the cut-off radius: ?dist <
- * affinity.cut?.
- * 2. Then I check whether there exists a bond from the current particle:
- * ?bond_site != -1?
- * 3. If yes, then I maintain the bond. I put the forces and afterwards I
- * decide whether the bond will break or not.
- * 4. If no, I maintain the creation of a bond. First I check whether I am
- * in the area of possible bond creation: ?dist < affinity.r0?
- * 5. If yes, I run the decision algorithm for bond creation and I either
- * create or does not create the bond.
- * 6. If I am not in the area of possible bond creation I do nothing
- *
- *
- * comments:
- * - strength of the force is proportional to the difference of
- * actual bond length and the relaxed bond length
- * - bond is created with probability 1-exp(-Kon*timestep)
- * - maxBond is not used, we use probability 1-exp(-Koff*timestep) to
- * break the bond
- * - Koff depends on the bond length via Koff = K0*exp(F/Fd) =
- * K0*exp(kappa(r-0.75*r0)/Fd)
- * - here, ia_params.Koff gives us K_0, off rate when bond is relaxed
- * - here, maxBond is used as detachment force F_d
- * - the original check for ensuring, that particle flows out of the
- * cut-off radius and the bond remains active is replaced with fixed
- * check, that bond length must not be greater that 0.8 cut_off
- *********************/
- double fac = 0.0;
- if (dist < ia_params.affinity.cut) { // Checking whether I am inside the
- // interaction cut-off radius.
- if (dist > 0.0) {
- // printf("bond_site: %f %f
- // %f\n",p1.p.bond_site[0],p1.p.bond_site[1],p1.p.bond_site[2]);
- if ((p1.p.bond_site[0] >= 0) && (p1.p.bond_site[1] >= 0) &&
- (p1.p.bond_site[2] >= 0)) // Checking whether any bond exists
- { // Bond exists
- if (len > 1.0 * (ia_params.affinity.r0)) {
- fac = ia_params.affinity.kappa *
- (len - 1.0 * (ia_params.affinity.r0));
- // printf("len %f r0 %f\n",len, ia_params.affinity.r0);
- } else
- fac = 0.0;
- // double ftemp = 0;
- force += (fac / len) * vec;
- // if (ftemp > 0.000000000000001) printf("%f ",fac);
- // Decision whether I should break the bond:
- // First, force exerted on bond is stored in fac
- double tmpF = fac;
- // Then, zero force off rate K_0 is stored at ia_params.affinity.Koff
- double tmpK0 = ia_params.affinity.Koff;
- // Then, detachment force is stored in ia_params.affinity.maxBond
- double tmpFd = ia_params.affinity.maxBond;
- // Then, compute Koff
- double tmpKoff = tmpK0 * exp(tmpF / tmpFd);
- // Finally, compute Poff
- double Poff = 1.0 - exp(-tmpKoff * time_step);
- // printf("%f ", Poff);
- if (len <
- 0.8 *
- ia_params.affinity.cut) { // in other implementation, maxBond
- // is used here. However, in this
- // implementation, we need maxBond
- // for setting detachment force F_d
- double decide = d_random();
- if (decide < Poff) {
- p1.p.bond_site = {-1, -1, -1};
- }
-
- } else {
- p1.p.bond_site = {-1, -1, -1};
- // printf("breaking: out of cut");
- }
- // Checkpoint output:
- if (period_for_output > 0)
- if (((int)floor(sim_time / time_step) % period_for_output == 0) &&
- (len > ia_params.affinity.r0)) {
- FILE *fp;
- double tmpPon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- fp = fopen("affinity_check.dat", "a");
- fprintf(fp, "sim_time %f, period_for_output %d aff type: %d ",
- sim_time, period_for_output, aff_type_extracted);
- fprintf(fp,
- "Pon %f, Kon %f, particle %d, Poff = %f, F = %f, Koff = "
- "%f, K0 = %f, len = %f \n",
- tmpPon, ia_params.affinity.Kon, p1.p.identity, Poff, tmpF,
- tmpKoff, tmpK0, len);
- fclose(fp);
- }
- } else if (dist < ia_params.affinity
- .r0) { // Bond does not exist, we are inside
- // of possible bond creation area,
- // let's talk about creating a bond
- double Pon = 1.0 - exp(-ia_params.affinity.Kon * time_step);
- // The probability is given by function Pon(x)= 1 - e^(-x) where x is
- // Kon*dt.
- double decide = d_random();
- if (decide <
- Pon) { // the bond will be created only with probability Pon.
- p1.p.bond_site = unfolded_pos - d;
- } else {
- // printf("In range, not creating: Pon = %f, decide = %f", Pon,
- // decide);
- }
- }
- }
- }
- }
- return force;
-}
-
-#endif /* ifdef AFFINITY */
-#endif
diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp
index 594c216019d..e3a2b482669 100644
--- a/src/core/particle_data.cpp
+++ b/src/core/particle_data.cpp
@@ -137,9 +137,6 @@ using UpdatePropertyMessage = boost::variant
#ifdef ROTATIONAL_INERTIA
, UpdateProperty
#endif
-#ifdef AFFINITY
- , UpdateProperty
-#endif
#ifdef MEMBRANE_COLLISION
, UpdateProperty
#endif
@@ -884,13 +881,6 @@ void rotate_particle(int part, const Utils::Vector3d &axis, double angle) {
}
#endif
-#ifdef AFFINITY
-void set_particle_affinity(int part, double *bond_site) {
- mpi_update_particle_property(
- part, Utils::Vector3d(bond_site, bond_site + 3));
-}
-#endif
-
#ifdef MEMBRANE_COLLISION
void set_particle_out_direction(int part, double *out_direction) {
mpi_update_particle_propertyp.bond_site.data();
-}
-#endif
-
#ifdef MEMBRANE_COLLISION
void pointer_to_out_direction(const Particle *p, const double *&res) {
res = p->p.out_direction.data();
diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp
index e377247670d..b5183d16171 100644
--- a/src/core/particle_data.hpp
+++ b/src/core/particle_data.hpp
@@ -110,11 +110,6 @@ struct ParticleProperties {
static constexpr Utils::Vector3d rinertia = {1., 1., 1.};
#endif
-#ifdef AFFINITY
- /** parameters for affinity mechanisms */
- Utils::Vector3d bond_site = {-1., -1., -1.};
-#endif
-
#ifdef MEMBRANE_COLLISION
/** parameters for membrane collision mechanisms */
Utils::Vector3d out_direction = {0., 0., 0.};
@@ -638,14 +633,6 @@ void set_particle_rotation(int part, int rot);
*/
void rotate_particle(int part, const Utils::Vector3d &axis, double angle);
-#ifdef AFFINITY
-/** Call only on the master node: set particle affinity.
- * @param part the particle.
- * @param bond_site its new site of the affinity bond.
- */
-void set_particle_affinity(int part, double *bond_site);
-#endif
-
#ifdef MEMBRANE_COLLISION
/** Call only on the master node: set particle out_direction.
* @param part the particle.
@@ -990,9 +977,6 @@ void pointer_to_swimming(Particle const *p,
#ifdef ROTATIONAL_INERTIA
void pointer_to_rotational_inertia(Particle const *p, double const *&res);
#endif
-#ifdef AFFINITY
-void pointer_to_bond_site(Particle const *p, double const *&res);
-#endif
#ifdef MEMBRANE_COLLISION
void pointer_to_out_direction(const Particle *p, const double *&res);
diff --git a/src/python/espressomd/interactions.pxd b/src/python/espressomd/interactions.pxd
index 538142b0cd0..b0258eaf361 100644
--- a/src/python/espressomd/interactions.pxd
+++ b/src/python/espressomd/interactions.pxd
@@ -123,15 +123,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
double cut
double offset
- cdef struct Affinity_Parameters:
- int type
- double kappa
- double r0
- double Kon
- double Koff
- double maxBond
- double cut
-
cdef struct Membrane_Parameters:
double a
double n
@@ -178,8 +169,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
LJGen_Parameters ljgen
- Affinity_Parameters affinity
-
Membrane_Parameters membrane
SoftSphere_Parameters soft_sphere
@@ -324,11 +313,6 @@ IF SOFT_SPHERE:
double a, double n,
double cut, double offset)
-IF AFFINITY:
- cdef extern from "object-in-fluid/affinity.hpp":
- cdef int affinity_set_params(int part_type_a, int part_type_b,
- int afftype, double kappa, double r0,
- double Kon, double Koff, double maxBond, double cut)
IF TABULATED:
cdef extern from "nonbonded_interactions/nonbonded_tab.hpp":
int tabulated_set_params(int part_type_a, int part_type_b,
diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx
index 998a6861263..fc449528d84 100644
--- a/src/python/espressomd/interactions.pyx
+++ b/src/python/espressomd/interactions.pyx
@@ -1382,69 +1382,6 @@ IF SOFT_SPHERE == 1:
"""
return {"a", "n", "cutoff"}
-IF AFFINITY == 1:
-
- cdef class AffinityInteraction(NonBondedInteraction):
-
- def validate_params(self):
- if self._params["affinity_cut"] < 0:
- raise ValueError("Affinity cutoff has to be >=0")
- if self._params["affinity_type"] < 0:
- raise ValueError("Affinity type has to be >=0")
- if self._params["affinity_kappa"] < 0:
- raise ValueError("Affinity kappa has to be >=0")
- if self._params["affinity_r0"] < 0:
- raise ValueError("Affinity r0 has to be >=0")
- if self._params["affinity_Kon"] < 0:
- raise ValueError("Affinity Kon has to be >=0")
- if self._params["affinity_Koff"] < 0:
- raise ValueError("Affinity Koff has to be >=0")
- if self._params["affinity_maxBond"] < 0:
- raise ValueError("Affinity maxBond has to be >=0")
- return True
-
- def _get_params_from_es_core(self):
- cdef IA_parameters * ia_params
- ia_params = get_ia_param_safe(self._part_types[0],
- self._part_types[1])
- return {
- "affinity_type": ia_params.affinity.type,
- "affinity_kappa": ia_params.affinity.kappa,
- "affinity_r0": ia_params.affinity.r0,
- "affinity_Kon": ia_params.affinity.Kon,
- "affinity_Koff": ia_params.affinity.Koff,
- "affinity_maxBond": ia_params.affinity.maxBond,
- "affinity_cut": ia_params.affinity.cut}
-
- def is_active(self):
- return (self._params["affinity_kappa"] > 0)
-
- def _set_params_in_es_core(self):
- if affinity_set_params(self._part_types[0],
- self._part_types[1],
- self._params["affinity_type"],
- self._params["affinity_kappa"],
- self._params["affinity_r0"],
- self._params["affinity_Kon"],
- self._params["affinity_Koff"],
- self._params["affinity_maxBond"],
- self._params["affinity_cut"]):
- raise Exception("Could not set Affinity parameters")
-
- def default_params(self):
- return {}
-
- def type_name(self):
- return "Affinity"
-
- def valid_keys(self):
- return {"affinity_type", "affinity_kappa", "affinity_r0",
- "affinity_Kon", "affinity_Koff", "affinity_maxBond", "affinity_cut"}
-
- def required_keys(self):
- return {"affinity_type", "affinity_kappa", "affinity_r0",
- "affinity_Kon", "affinity_Koff", "affinity_maxBond", "affinity_cut"}
-
IF MEMBRANE_COLLISION == 1:
@@ -1693,8 +1630,6 @@ class NonBondedInteractionHandle:
_type1, _type2)
IF SOFT_SPHERE:
self.soft_sphere = SoftSphereInteraction(_type1, _type2)
- IF AFFINITY:
- self.affinity = AffinityInteraction(_type1, _type2)
IF LENNARD_JONES_GENERIC:
self.generic_lennard_jones = GenericLennardJonesInteraction(
_type1, _type2)
diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd
index 1586c89fbad..0822c2e22af 100644
--- a/src/python/espressomd/particle_data.pxd
+++ b/src/python/espressomd/particle_data.pxd
@@ -124,10 +124,6 @@ cdef extern from "particle_data.hpp":
void set_particle_out_direction(int part, double out_direction[3])
void pointer_to_out_direction(particle * p, double * & res)
- IF AFFINITY:
- void set_particle_affinity(int part, double bond_site[3])
- void pointer_to_bond_site(particle * p, double * & res)
-
IF MASS:
void pointer_to_mass(particle * p, double * & res)
diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx
index 571f98ea8f3..6f6c3cfeb67 100644
--- a/src/python/espressomd/particle_data.pyx
+++ b/src/python/espressomd/particle_data.pyx
@@ -570,25 +570,6 @@ cdef class ParticleHandle:
return np.array(
[out_direction[0], out_direction[1], out_direction[2]])
- IF AFFINITY:
- property bond_site:
- """OIF bond_site"""
-
- def __set__(self, _bond_site):
- cdef double bond_site[3]
- check_type_or_throw_except(
- _bond_site, 3, float, "bond_site has to be 3 floats")
- for i in range(3):
- bond_site[i] = _bond_site[i]
- set_particle_affinity(self.id, bond_site)
-
- def __get__(self):
- self.update_particle_data()
- cdef const double * bond_site = NULL
- pointer_to_bond_site(self.particle_data, bond_site)
- return array_locked([bond_site[0], bond_site[1], bond_site[2]])
-
-
# Charge
property q:
"""