From 78f1a225c3208c23a9cb07411e4e816cb9935e80 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 1 Sep 2020 16:44:09 -0500
Subject: [PATCH] Various bugfixes
---
src/ComputeForce.cuh | 10 ++++-----
src/TabulatedPotential.h | 47 +++++++++++++++++++++++++++++++++++++---
2 files changed, 49 insertions(+), 8 deletions(-)
diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh
index 275764f..c14c8b7 100644
--- a/src/ComputeForce.cuh
+++ b/src/ComputeForce.cuh
@@ -943,12 +943,12 @@ void computeProductPotentials(Vector3* force,
if (j == k) continue;
tmp_force *= energy_and_deriv[k].x;
}
- if (tmp_force == 0) continue;
-
- // TODO add energy
- // TODO make it work with replicas
SimplePotential& p = potentialList[ productPotential_list[i].x + j ];
- p.apply_force(pos,sys, force, &productPotentialParticles[part_idx], tmp_force);
+ if (tmp_force != 0) {
+ // TODO add energy
+ // TODO make it work with replicas
+ p.apply_force(pos,sys, force, &productPotentialParticles[part_idx], tmp_force);
+ }
part_idx += p.type==BOND? 2: p.type==ANGLE? 3: 4;
}
}
diff --git a/src/TabulatedPotential.h b/src/TabulatedPotential.h
index bebe922..c7933aa 100644
--- a/src/TabulatedPotential.h
+++ b/src/TabulatedPotential.h
@@ -258,7 +258,7 @@ public:
HOST DEVICE inline float2 compute_energy_and_deriv(float value) {
float2 ret;
if (type == DIHEDRAL) {
- ret = linearly_interpolate<true>(value, BD_PI);
+ ret = linearly_interpolate<true>(value, -BD_PI);
} else {
ret = linearly_interpolate<false>(value);
}
@@ -327,7 +327,14 @@ public:
int home = int( floorf(w) );
w = w - home;
// if (home < 0) return EnergyForce(v0[0], Vector3(0.0f));
- if (home < 0 || (home >= size && !is_periodic)) return make_float2(0.0f,0.0f);
+ if (home < 0) {
+ if (is_periodic) home += size;
+ else return make_float2(pot[0],0.0f);
+ }
+ else if (home >= size) {
+ if (is_periodic) home -= size;
+ else return make_float2(pot[size-1],0.0f);
+ }
float u0 = pot[home];
float du = home+1 < size ? pot[home+1]-u0 : is_periodic ? pot[0]-u0 : 0;
@@ -390,6 +397,9 @@ public:
float sin = sqrtf(1.0f - cos*cos);
energy_deriv /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
+ if (abs(sin) < 1e-3) {
+ printf("BAD ANGLE: sin, cos, energy_deriv, distab, distbc, distac2: (%f %f %f %f %f)\n",
+ sin,cos,energy_deriv,distab,distbc); }
// Calculate the forces
TwoVector3 force;
@@ -398,6 +408,37 @@ public:
return force;
}
+ // DEVICE inline TwoVector3 get_angle_force(const Vector3& ab,
+ // const Vector3& bc,
+ // float energy_deriv) const {
+ // // Find the distance between each pair of particles
+ // float distab = ab.length2();
+ // float distbc = bc.length2();
+
+ // float pre = distab*distbc - pow(ab.dot(bc),2);
+ // // if (pre < 1e-6) {
+ // // pre = 1e-3;
+ // // printf("BAD ANGLE: pre, energy_deriv, distab, distbc, ab.dot(bc): (%f %f %f %f %f)\n",
+ // // pre,energy_deriv,distab,distbc,ab.dot(bc));
+ // // } else pre = sqrt(pre);
+ // // if (distab == distbc) {
+ // // printf("GOOD ANGLE: pre, energy_deriv, distab, distbc, ab.dot(bc): (%f %f %f %f %f)\n",
+ // // pre,energy_deriv,distab,distbc,ab.dot(bc));
+ // // }
+
+ // pre = pre > 1e-6 ? sqrt(pre) : 1e-3;
+ // energy_deriv /= pre;
+
+ // TwoVector3 force;
+ // //force.v1 = energy_deriv * Vector3::element_mult( 1-Vector3::element_mult(ab,ab)/distab, bc);
+ // //force.v2 = energy_deriv * Vector3::element_mult( 1-Vector3::element_mult(bc,bc)/distbc, ab);
+
+ // Vector3 abbc = Vector3::element_mult(-ab,bc);
+ // force.v1 = -energy_deriv * (bc-Vector3::element_mult(abbc, -ab/distab));
+ // force.v2 = -energy_deriv * (-ab-Vector3::element_mult(abbc, bc/distbc));
+ // return force;
+ // }
+
DEVICE inline void apply_angle_force(const Vector3* __restrict__ pos,
const BaseGrid* __restrict__ sys,
Vector3* __restrict__ force,
@@ -434,7 +475,7 @@ public:
// return -atan2(sin_phi,cos_phi);
Vector3 f1, f2, f3; // forces
- float distbc = bc.length2();
+ float distbc = bc.length();
f1 = -distbc * crossABC.rLength2() * crossABC;
f3 = -distbc * crossBCD.rLength2() * crossBCD;
--
GitLab