Skip to content

Commit

Permalink
core: Removed random includes and random forces in non-bonded interac…
Browse files Browse the repository at this point in the history 
…tions.
  • Loading branch information
@fweik
fweik committed Nov 12, 2018
1 parent f80bd28 commit 5dd5a81
Show file tree
Hide file tree
Showing 7 changed files with 15 additions and 116 deletions.
39 changes: 3 additions & 36 deletions src/core/bonded_interactions/fene.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"

/************************************************************/
Expand All @@ -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;
Expand All @@ -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;
}

Expand All @@ -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;
}

Expand Down
34 changes: 5 additions & 29 deletions src/core/bonded_interactions/harmonic.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
}
Expand Down
50 changes: 6 additions & 44 deletions src/core/bonded_interactions/harmonic_dumbbell.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,7 @@
/************************************************************/

#include "bonded_interaction_data.hpp"
#include "debug.hpp"
#include "particle_data.hpp"
#include "random.hpp"
#include "utils.hpp"

#ifdef ROTATION
Expand All @@ -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;
}

Expand Down
1 change: 0 additions & 1 deletion src/core/bonded_interactions/quartic.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 0 additions & 1 deletion src/core/cuda_interface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
5 changes: 1 addition & 4 deletions src/core/object-in-fluid/affinity.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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;
Expand Down
1 change: 0 additions & 1 deletion src/core/object-in-fluid/membrane_collision.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down

0 comments on commit 5dd5a81

Please sign in to comment.