From 5dd5a81a90515027ee25c5a4c2ed7f006753d2de Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 12 Nov 2018 16:01:03 +0100 Subject: [PATCH] core: Removed random includes and random forces in non-bonded interactions. --- src/core/bonded_interactions/fene.hpp | 39 ++------------- src/core/bonded_interactions/harmonic.hpp | 34 ++----------- .../bonded_interactions/harmonic_dumbbell.hpp | 50 +++---------------- src/core/bonded_interactions/quartic.hpp | 1 - src/core/cuda_interface.cpp | 1 - src/core/object-in-fluid/affinity.hpp | 5 +- .../object-in-fluid/membrane_collision.hpp | 1 - 7 files changed, 15 insertions(+), 116 deletions(-) diff --git a/src/core/bonded_interactions/fene.hpp b/src/core/bonded_interactions/fene.hpp index ea90ac5b2af..3f63409b678 100644 --- a/src/core/bonded_interactions/fene.hpp +++ b/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; } diff --git a/src/core/bonded_interactions/harmonic.hpp b/src/core/bonded_interactions/harmonic.hpp index 5628c7f4417..a3a0bbfc010 100644 --- a/src/core/bonded_interactions/harmonic.hpp +++ b/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; } diff --git a/src/core/bonded_interactions/harmonic_dumbbell.hpp b/src/core/bonded_interactions/harmonic_dumbbell.hpp index 3a082d67768..a6ba498ef13 100644 --- a/src/core/bonded_interactions/harmonic_dumbbell.hpp +++ b/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; } diff --git a/src/core/bonded_interactions/quartic.hpp b/src/core/bonded_interactions/quartic.hpp index 343a707a905..6b5c294ae92 100644 --- a/src/core/bonded_interactions/quartic.hpp +++ b/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 diff --git a/src/core/cuda_interface.cpp b/src/core/cuda_interface.cpp index 31089e22b73..5812a4decc7 100644 --- a/src/core/cuda_interface.cpp +++ b/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" diff --git a/src/core/object-in-fluid/affinity.hpp b/src/core/object-in-fluid/affinity.hpp index 11239d98688..845e1f2073d 100644 --- a/src/core/object-in-fluid/affinity.hpp +++ b/src/core/object-in-fluid/affinity.hpp @@ -30,8 +30,8 @@ #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "particle_data.hpp" -#include "random.hpp" #include "utils.hpp" +#include "random.hpp" #ifdef AFFINITY @@ -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; diff --git a/src/core/object-in-fluid/membrane_collision.hpp b/src/core/object-in-fluid/membrane_collision.hpp index c6f11ee8d10..a6ef29551a8 100644 --- a/src/core/object-in-fluid/membrane_collision.hpp +++ b/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);