core: Removed random includes and random forces in non-bonded interac…

src/core/bonded_interactions/fene.hpp

@@ -27,10 +27,7 @@
  */
 
 #include "bonded_interaction_data.hpp"
-#include "debug.hpp"
-#include "errorhandling.hpp"
 #include "particle_data.hpp"
-#include "random.hpp"
 #include "utils.hpp"
 
 /************************************************************/
@@ -51,8 +48,6 @@ int fene_set_params(int bond_type, double k, double drmax, double r0);
 inline int calc_fene_pair_force(Particle *p1, Particle *p2,
                                 Bonded_ia_parameters *iaparams, double dx[3],
                                 double force[3]) {
-  int i;
-
   const double len2 = sqrlen(dx);
   const double len = sqrt(len2);
   const double dr = len - iaparams->p.fene.r0;
@@ -62,41 +57,15 @@ inline int calc_fene_pair_force(Particle *p1, Particle *p2,
 
   double fac =
       -iaparams->p.fene.k * dr / ((1.0 - dr * dr * iaparams->p.fene.drmax2i));
-  if (fabs(dr) > ROUND_ERROR_PREC) {
-    if (len > ROUND_ERROR_PREC) { /* Regular case */
-      fac /= len;
-    } else { /* dx[] == 0: the force is undefined. Let's use a random direction
-              */
-      for (int i = 0; i < 3; i++)
-        dx[i] = d_random() - 0.5;
-      fac /= sqrt(sqrlen(dx));
-    }
+  if (len > ROUND_ERROR_PREC) {
+    fac /= len;
   } else {
     fac = 0.0;
   }
 
-  FENE_TRACE(if (fac > 50)
-                 fprintf(stderr,
-                         "WARNING: FENE force factor between Pair (%d,%d) "
-                         "large: %f at distance %f\n",
-                         p1->p.identity, p2->p.identity, fac, sqrt(len2)));
-
-  for (i = 0; i < 3; i++)
+  for (int i = 0; i < 3; i++)
     force[i] = fac * dx[i];
 
-  ONEPART_TRACE(if (p1->p.identity == check_id)
-                    fprintf(stderr,
-                            "%d: OPT: FENE f = (%.3e,%.3e,%.3e) with part "
-                            "id=%d at dist %f fac %.3e\n",
-                            this_node, p1->f.f[0], p1->f.f[1], p1->f.f[2],
-                            p2->p.identity, sqrt(len2), fac));
-  ONEPART_TRACE(if (p2->p.identity == check_id)
-                    fprintf(stderr,
-                            "%d: OPT: FENE f = (%.3e,%.3e,%.3e) with part "
-                            "id=%d at dist %f fac %.3e\n",
-                            this_node, p2->f.f[0], p2->f.f[1], p2->f.f[2],
-                            p1->p.identity, sqrt(len2), fac));
-
   return 0;
 }
 
@@ -108,8 +77,6 @@ inline int fene_pair_energy(Particle *p1, Particle *p2,
 
   /* check bond stretching */
   if (dr >= iaparams->p.fene.drmax) {
-    runtimeErrorMsg() << "FENE bond broken between particles " << p1->p.identity
-                      << " and " << p2->p.identity;
     return 1;
   }
 

src/core/bonded_interactions/harmonic.hpp

@@ -29,9 +29,7 @@
 /************************************************************/
 
 #include "bonded_interaction_data.hpp"
-#include "debug.hpp"
 #include "particle_data.hpp"
-#include "random.hpp"
 #include "utils.hpp"
 
 /// set the parameters for the harmonic potential
@@ -50,44 +48,22 @@ int harmonic_set_params(int bond_type, double k, double r, double r_cut);
 inline int calc_harmonic_pair_force(Particle *p1, Particle *p2,
                                     Bonded_ia_parameters *iaparams,
                                     double dx[3], double force[3]) {
-  int i;
-  double fac;
   double dist2 = sqrlen(dx);
   double dist = sqrt(dist2);
-  double dr;
 
   if ((iaparams->p.harmonic.r_cut > 0.0) && (dist > iaparams->p.harmonic.r_cut))
     return 1;
 
-  dr = dist - iaparams->p.harmonic.r;
-  fac = -iaparams->p.harmonic.k * dr;
-  if (fabs(dr) > ROUND_ERROR_PREC) {
-    if (dist > ROUND_ERROR_PREC) { /* Regular case */
-      fac /= dist;
-    } else { /* dx[] == 0: the force is undefined. Let's use a random direction
-              */
-      for (i = 0; i < 3; i++)
-        dx[i] = d_random() - 0.5;
-      fac /= sqrt(sqrlen(dx));
-    }
+  auto const dr = dist - iaparams->p.harmonic.r;
+  auto fac = -iaparams->p.harmonic.k * dr;
+  if (dist > ROUND_ERROR_PREC) { /* Regular case */
+    fac /= dist;
   } else {
     fac = 0;
   }
 
-  for (i = 0; i < 3; i++)
+  for (int i = 0; i < 3; i++)
     force[i] = fac * dx[i];
-  ONEPART_TRACE(if (p1->p.identity == check_id)
-                    fprintf(stderr,
-                            "%d: OPT: HARMONIC f = (%.3e,%.3e,%.3e) with part "
-                            "id=%d at dist %f fac %.3e\n",
-                            this_node, p1->f.f[0], p1->f.f[1], p1->f.f[2],
-                            p2->p.identity, dist2, fac));
-  ONEPART_TRACE(if (p2->p.identity == check_id)
-                    fprintf(stderr,
-                            "%d: OPT: HARMONIC f = (%.3e,%.3e,%.3e) with part "
-                            "id=%d at dist %f fac %.3e\n",
-                            this_node, p2->f.f[0], p2->f.f[1], p2->f.f[2],
-                            p1->p.identity, dist2, fac));
 
   return 0;
 }

src/core/bonded_interactions/harmonic_dumbbell.hpp

@@ -29,9 +29,7 @@
 /************************************************************/
 
 #include "bonded_interaction_data.hpp"
-#include "debug.hpp"
 #include "particle_data.hpp"
-#include "random.hpp"
 #include "utils.hpp"
 
 #ifdef ROTATION
@@ -53,60 +51,24 @@ int harmonic_dumbbell_set_params(int bond_type, double k1, double k2, double r,
 inline int calc_harmonic_dumbbell_pair_force(Particle *p1, Particle *p2,
                                              Bonded_ia_parameters *iaparams,
                                              double dx[3], double force[3]) {
-  double fac;
   double dist2 = sqrlen(dx);
   double dist = sqrt(dist2);
-  double dr;
 
   if ((iaparams->p.harmonic_dumbbell.r_cut > 0.0) &&
       (dist > iaparams->p.harmonic_dumbbell.r_cut))
     return 1;
 
-  dr = dist - iaparams->p.harmonic_dumbbell.r;
-  fac = -iaparams->p.harmonic_dumbbell.k1 * dr;
-  if (fabs(dr) > ROUND_ERROR_PREC) {
-    if (dist > ROUND_ERROR_PREC) /* Regular case */
-      fac /= dist;
-    else { /* dx[] == 0: the force is undefined. Let's use a random direction */
-      for (int i = 0; i < 3; i++)
-        dx[i] = d_random() - 0.5;
-      fac /= sqrt(sqrlen(dx));
-    }
-  } else {
-    fac = 0;
-  }
+  auto const dr = dist - iaparams->p.harmonic_dumbbell.r;
+  auto const normalizer = (dist > ROUND_ERROR_PREC) ? 1. / dist : 0.0;
+  auto const fac = -iaparams->p.harmonic_dumbbell.k1 * dr * normalizer;
 
   for (int i = 0; i < 3; i++)
     force[i] = fac * dx[i];
 
-  double dhat[3];
-  dhat[0] = dx[0] / dist;
-  dhat[1] = dx[1] / dist;
-  dhat[2] = dx[2] / dist;
-
-  double da[3];
-  const Vector3d director1 = p1->r.calc_director();
-  da[0] = dhat[1] * director1[2] - dhat[2] * director1[1];
-  da[1] = dhat[2] * director1[0] - dhat[0] * director1[2];
-  da[2] = dhat[0] * director1[1] - dhat[1] * director1[0];
-
-  p1->f.torque[0] += iaparams->p.harmonic_dumbbell.k2 * da[0];
-  p1->f.torque[1] += iaparams->p.harmonic_dumbbell.k2 * da[1];
-  p1->f.torque[2] += iaparams->p.harmonic_dumbbell.k2 * da[2];
-
-  ONEPART_TRACE(if (p1->p.identity == check_id)
-                    fprintf(stderr,
-                            "%d: OPT: HARMONIC f = (%.3e,%.3e,%.3e) with part "
-                            "id=%d at dist %f fac %.3e\n",
-                            this_node, p1->f.f[0], p1->f.f[1], p1->f.f[2],
-                            p2->p.identity, dist2, fac));
-  ONEPART_TRACE(if (p2->p.identity == check_id)
-                    fprintf(stderr,
-                            "%d: OPT: HARMONIC f = (%.3e,%.3e,%.3e) with part "
-                            "id=%d at dist %f fac %.3e\n",
-                            this_node, p2->f.f[0], p2->f.f[1], p2->f.f[2],
-                            p1->p.identity, dist2, fac));
+  auto const dhat = Vector3d{dx[0], dx[1], dx[2]} * normalizer;
+  auto const da = dhat.cross(p1->r.calc_director());
 
+  p1->f.torque += iaparams->p.harmonic_dumbbell.k2 * da;
   return 0;
 }
 

src/core/bonded_interactions/quartic.hpp

@@ -31,7 +31,6 @@
 #include "bonded_interaction_data.hpp"
 #include "debug.hpp"
 #include "particle_data.hpp"
-#include "random.hpp"
 #include "utils.hpp"
 
 /// set the parameters for the quartic potential

src/core/cuda_interface.cpp

@@ -27,7 +27,6 @@
 #include "energy.hpp"
 #include "grid.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "random.hpp"
 #include "utils/mpi/gather_buffer.hpp"
 #include "utils/mpi/scatter_buffer.hpp"
 #include "utils/serialization/CUDA_particle_data.hpp"

src/core/object-in-fluid/affinity.hpp

@@ -226,10 +226,7 @@ inline void add_affinity_pair_force(Particle *p1, Particle *p2,
             if (decide < Poff) {
               for (j = 0; j < 3; j++)
                 p1->p.bond_site[j] = -1;
-              // printf("bond broken. Poff = %f, F = %f, Koff = %f, K0 = %f, len
-              // = %f\n", Poff, tmpF, tmpKoff, tmpK0, len);
             }
-
           } else {
             for (j = 0; j < 3; j++)
               p1->p.bond_site[j] = -1;

src/core/object-in-fluid/membrane_collision.hpp

@@ -35,7 +35,6 @@
 #include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "particle_data.hpp"
-#include "random.hpp"
 
 int membrane_collision_set_params(int part_type_a, int part_type_b, double a,
                                   double n, double cut, double offset);