Bond cutoff

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -19,59 +19,86 @@
 #include "bonded_interaction_data.hpp"
 #include "communication.hpp"
 
+#include <boost/range/numeric.hpp>
 #include <utils/constants.hpp>
 
 std::vector<Bonded_ia_parameters> bonded_ia_params;
 
+auto cutoff(int type, Bond_parameters const &bp) {
+  switch (type) {
+  case BONDED_IA_NONE:
+    return -1.;
+  case BONDED_IA_FENE:
+    return bp.fene.cutoff();
+  case BONDED_IA_HARMONIC:
+    return bp.harmonic.cutoff();
+  case BONDED_IA_HARMONIC_DUMBBELL:
+    return bp.harmonic_dumbbell.cutoff();
+  case BONDED_IA_QUARTIC:
+    return bp.quartic.cutoff();
+  case BONDED_IA_BONDED_COULOMB:
+    return bp.bonded_coulomb.cutoff();
+  case BONDED_IA_BONDED_COULOMB_SR:
+    return bp.bonded_coulomb_sr.cutoff();
+  case BONDED_IA_DIHEDRAL:
+    return bp.dihedral.cutoff();
+  case BONDED_IA_TABULATED_DISTANCE:
+  case BONDED_IA_TABULATED_ANGLE:
+  case BONDED_IA_TABULATED_DIHEDRAL:
+    return bp.tab.cutoff();
+  case BONDED_IA_SUBT_LJ:
+    return bp.subt_lj.cutoff();
+  case BONDED_IA_RIGID_BOND:
+    return bp.rigid_bond.cutoff();
+  case BONDED_IA_VIRTUAL_BOND:
+    return bp.virt.cutoff();
+  case BONDED_IA_ANGLE_HARMONIC:
+    return bp.angle_harmonic.cutoff();
+  case BONDED_IA_ANGLE_COSINE:
+    return bp.angle_cosine.cutoff();
+  case BONDED_IA_ANGLE_COSSQUARE:
+    return bp.angle_cossquare.cutoff();
+  case BONDED_IA_OIF_LOCAL_FORCES:
+    return bp.oif_local_forces.cutoff();
+  case BONDED_IA_OIF_GLOBAL_FORCES:
+    return bp.oif_global_forces.cutoff();
+  case BONDED_IA_IBM_TRIEL:
+    return bp.ibm_triel.cutoff();
+  case BONDED_IA_IBM_VOLUME_CONSERVATION:
+    return bp.ibmVolConsParameters.cutoff();
+  case BONDED_IA_IBM_TRIBEND:
+    return bp.ibm_tribend.cutoff();
+  case BONDED_IA_UMBRELLA:
+    return bp.umbrella.cutoff();
+  case BONDED_IA_THERMALIZED_DIST:
+    return bp.thermalized_bond.cutoff();
+  default:
+    throw std::runtime_error("Unknown bond type.");
+  }
+}
+
 double recalc_maximal_cutoff_bonded() {
-  auto max_cut_bonded = -1.;
+  auto const max_cut_bonded =
+      boost::accumulate(bonded_ia_params, -1.,
+                        [](auto max_cut, Bonded_ia_parameters const &bond) {
+                          return std::max(max_cut, cutoff(bond.type, bond.p));
+                        });
 
-  for (auto const &bonded_ia_param : bonded_ia_params) {
-    switch (bonded_ia_param.type) {
-    case BONDED_IA_FENE:
-      max_cut_bonded =
-          std::max(max_cut_bonded,
-                   bonded_ia_param.p.fene.r0 + bonded_ia_param.p.fene.drmax);
-      break;
-    case BONDED_IA_HARMONIC:
-      max_cut_bonded =
-          std::max(max_cut_bonded, bonded_ia_param.p.harmonic.r_cut);
-      break;
-    case BONDED_IA_THERMALIZED_DIST:
-      max_cut_bonded =
-          std::max(max_cut_bonded, bonded_ia_param.p.thermalized_bond.r_cut);
-      break;
-    case BONDED_IA_RIGID_BOND:
-      max_cut_bonded =
-          std::max(max_cut_bonded, std::sqrt(bonded_ia_param.p.rigid_bond.d2));
-      break;
-    case BONDED_IA_TABULATED_DISTANCE:
-      max_cut_bonded =
-          std::max(max_cut_bonded, bonded_ia_param.p.tab.pot->cutoff());
-      break;
-    case BONDED_IA_IBM_TRIEL:
-      max_cut_bonded =
-          std::max(max_cut_bonded, bonded_ia_param.p.ibm_triel.maxDist);
-      break;
-    default:
-      break;
-    }
-  }
+  /* Check if there are dihedrals */
+  auto const any_dihedrals = std::any_of(
+      bonded_ia_params.begin(), bonded_ia_params.end(), [](auto const &bond) {
+        switch (bond.type) {
+        case BONDED_IA_DIHEDRAL:
+        case BONDED_IA_TABULATED_DIHEDRAL:
+          return true;
+        default:
+          return false;
+        }
+      });
 
   /* dihedrals: the central particle is indirectly connected to the fourth
    * particle via the third particle, so we have to double the cutoff */
-  for (auto const &bonded_ia_param : bonded_ia_params) {
-    switch (bonded_ia_param.type) {
-    case BONDED_IA_DIHEDRAL:
-    case BONDED_IA_TABULATED_DIHEDRAL:
-      max_cut_bonded *= 2;
-      break;
-    default:
-      break;
-    }
-  }
-
-  return max_cut_bonded;
+  return (any_dihedrals) ? 2 * max_cut_bonded : max_cut_bonded;
 }
 
 void make_bond_type_exist(int type) {
@@ -95,6 +122,7 @@ int virtual_set_params(int bond_type) {
 
   bonded_ia_params[bond_type].type = BONDED_IA_VIRTUAL_BOND;
   bonded_ia_params[bond_type].num = 1;
+  bonded_ia_params[bond_type].p.virt = VirtualBond_Parameters{};
 
   /* broadcast interaction parameters */
   mpi_bcast_ia_params(bond_type, -1);

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -73,10 +73,6 @@ enum BondedInteraction {
   BONDED_IA_OIF_LOCAL_FORCES,
   /** Type of bonded interaction: OIF global forces. */
   BONDED_IA_OIF_GLOBAL_FORCES,
-  /** Type of bonded interaction: determining outward direction of OIF membrane
-   *  (not associated to a parameter struct).
-   */
-  BONDED_IA_OIF_OUT_DIRECTION,
   /** Type of bonded interaction is a wall repulsion (immersed boundary). */
   BONDED_IA_IBM_TRIEL,
   /** Type of bonded interaction is volume conservation force (immersed
@@ -110,6 +106,8 @@ struct Fene_bond_parameters {
   double drmax2;
   /** inverse square of @p drmax (internal parameter) */
   double drmax2i;
+
+  double cutoff() const { return r0 + drmax; }
 };
 
 /** Parameters for OIF global forces
@@ -126,6 +124,8 @@ struct Oif_global_forces_bond_parameters {
   double V0;
   /** Volume coefficient */
   double kv;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for OIF local forces
@@ -151,6 +151,8 @@ struct Oif_local_forces_bond_parameters {
   double kal;
   /** Viscous coefficient of the triangle vertices */
   double kvisc;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for harmonic bond Potential */
@@ -161,6 +163,8 @@ struct Harmonic_bond_parameters {
   double r;
   /** cutoff bond length */
   double r_cut;
+
+  double cutoff() const { return r_cut; }
 };
 
 /** Parameters for Thermalized bond **/
@@ -174,9 +178,10 @@ struct Thermalized_bond_parameters {
   double pref2_com;
   double pref1_dist;
   double pref2_dist;
+
+  double cutoff() const { return r_cut; }
 };
 
-#ifdef ROTATION
 /** Parameters for harmonic dumbbell bond Potential */
 struct Harmonic_dumbbell_bond_parameters {
   /** spring constant */
@@ -187,26 +192,33 @@ struct Harmonic_dumbbell_bond_parameters {
   double r;
   /** cutoff bond length */
   double r_cut;
+
+  double cutoff() const { return r_cut; }
 };
-#endif
 
 /** Parameters for quartic bond Potential */
 struct Quartic_bond_parameters {
   double k0, k1;
   double r;
   double r_cut;
+
+  double cutoff() const { return r_cut; }
 };
 
 /** Parameters for %Coulomb bond Potential */
 struct Bonded_coulomb_bond_parameters {
   /** %Coulomb prefactor */
   double prefactor;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for %Coulomb bond short-range Potential */
 struct Bonded_coulomb_sr_bond_parameters {
   /** charge factor */
   double q1q2;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for three-body angular potential (harmonic). */
@@ -215,6 +227,8 @@ struct Angle_harmonic_bond_parameters {
   double bend;
   /** equilibrium angle (default is 180 degrees) */
   double phi0;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for three-body angular potential (cosine). */
@@ -227,6 +241,8 @@ struct Angle_cosine_bond_parameters {
   double cos_phi0;
   /** sine of @p phi0 (internal parameter) */
   double sin_phi0;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for three-body angular potential (cossquare). */
@@ -237,31 +253,47 @@ struct Angle_cossquare_bond_parameters {
   double phi0;
   /** cosine of @p phi0 (internal parameter) */
   double cos_phi0;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for four-body angular potential (dihedral-angle potentials). */
 struct Dihedral_bond_parameters {
   double mult;
   double bend;
   double phase;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for n-body tabulated potential (n=2,3,4). */
 struct Tabulated_bond_parameters {
+  TabulatedBondedInteraction type;
   TabulatedPotential *pot;
+
+  double cutoff() const {
+    switch (type) {
+    case TAB_BOND_LENGTH:
+      return assert(pot), pot->cutoff();
+    default:
+      return -1.;
+    };
+  }
 };
 
-#ifdef UMBRELLA
 /** Parameters for umbrella potential */
 struct Umbrella_bond_parameters {
   double k;
   int dir;
   double r;
+
+  double cutoff() const { return std::numeric_limits<double>::infinity(); }
 };
-#endif
 
 /** Dummy parameters for subtracted-LJ Potential */
-struct Subt_lj_bond_parameters {};
+struct Subt_lj_bond_parameters {
+  double cutoff() const { return -1.; }
+};
 
 /** Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM */
 struct Rigid_bond_parameters {
@@ -273,6 +305,8 @@ struct Rigid_bond_parameters {
   /**Velocity Tolerance/Accuracy for termination of RATTLE/SHAKE iterations
    * during velocity corrections */
   double v_tol;
+
+  double cutoff() const { return std::sqrt(d2); }
 };
 
 enum class tElasticLaw { NeoHookean, Skalak };
@@ -299,6 +333,8 @@ struct IBM_Triel_Parameters {
   tElasticLaw elasticLaw;
   double k1;
   double k2;
+
+  double cutoff() const { return maxDist; }
 };
 
 /** Parameters for IBM volume conservation bond **/
@@ -309,10 +345,8 @@ struct IBM_VolCons_Parameters {
   double volRef;
   /** Spring constant for volume force */
   double kappaV;
-  // Whether to write out center-of-mass at each time step
-  // Actually this is more of an analysis function and does not strictly belong
-  // to volume conservation
-  //  bool writeCOM;
+
+  double cutoff() const { return -1.; }
 };
 
 /** Parameters for IBM tribend **/
@@ -322,6 +356,12 @@ struct IBM_Tribend_Parameters {
 
   /** Reference angle */
   double theta0;
+
+  double cutoff() const { return -1.; }
+};
+
+struct VirtualBond_Parameters {
+  double cutoff() const { return -1.; }
 };
 
 /** Union in which to store the parameters of an individual bonded interaction
@@ -331,9 +371,7 @@ union Bond_parameters {
   Oif_global_forces_bond_parameters oif_global_forces;
   Oif_local_forces_bond_parameters oif_local_forces;
   Harmonic_bond_parameters harmonic;
-#ifdef ROTATION
   Harmonic_dumbbell_bond_parameters harmonic_dumbbell;
-#endif
   Quartic_bond_parameters quartic;
   Bonded_coulomb_bond_parameters bonded_coulomb;
   Bonded_coulomb_sr_bond_parameters bonded_coulomb_sr;
@@ -342,15 +380,14 @@ union Bond_parameters {
   Angle_cossquare_bond_parameters angle_cossquare;
   Dihedral_bond_parameters dihedral;
   Tabulated_bond_parameters tab;
-#ifdef UMBRELLA
   Umbrella_bond_parameters umbrella;
-#endif
   Thermalized_bond_parameters thermalized_bond;
   Subt_lj_bond_parameters subt_lj;
   Rigid_bond_parameters rigid_bond;
   IBM_Triel_Parameters ibm_triel;
   IBM_VolCons_Parameters ibmVolConsParameters;
   IBM_Tribend_Parameters ibm_tribend;
+  VirtualBond_Parameters virt;
 };
 
 /** Defines parameters for a bonded interaction. */

src/core/bonded_interactions/bonded_tab.cpp

@@ -38,6 +38,7 @@ int tabulated_bonded_set_params(int bond_type,
 
   /* set types */
   auto tab_pot = bonded_ia_params[bond_type].p.tab.pot = new TabulatedPotential;
+  bonded_ia_params[bond_type].p.tab.type = tab_type;
 
   /* set number of interaction partners */
   switch (tab_type) {