From ed88faeb90c635493a02e35f2fa98213ab5d082d Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 19 Dec 2017 13:03:49 -0600
Subject: [PATCH] HYC: Updated bound checking for angle lookups; Capped maximum
force for dihedral angle (perhaps this should be reverted)
---
src/TabulatedMethods.cuh | 31 ++++++++++++++++++++++---------
1 file changed, 22 insertions(+), 9 deletions(-)
diff --git a/src/TabulatedMethods.cuh b/src/TabulatedMethods.cuh
index cae844c..bbd9b82 100644
--- a/src/TabulatedMethods.cuh
+++ b/src/TabulatedMethods.cuh
@@ -44,18 +44,19 @@ __device__ inline void computeAngle(const TabulatedAnglePotential* __restrict__
// tableAngle[1] stores the potential at angle_step * 2, etc.
// 'home' is the index after which 'convertedAngle' would appear if it were stored in the table
int home = int(floor(angle));
- assert(home >= 0);
- assert(home+1 < a->size);
+ home = (home >= a->size) ? (a->size)-1 : home;
+ //assert(home >= 0);
+ //assert(home+1 < a->size);
// // Make angle the distance from [0,1) from the first index in the potential array index
// angle -= home;
// Linearly interpolate the potential
float U0 = a->pot[home];
- float dUdx = (a->pot[(home+1)] - U0) * a->angle_step_inv;
+ float dUdx = (a->pot[(((home+1)==(a->size)) ? (a->size)-1 : home+1)] - U0) * a->angle_step_inv;
// float energy = (dUdx * angle) + U0;
float sin = sqrtf(1.0f - cos*cos);
- dUdx /= abs(sin) > 1e-2 ? sin : 1e-2; // avoid singularity
+ dUdx /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
// Calculate the forces
Vector3 force1 = -(dUdx*distab) * (ab * (cos*distab) + bc * distbc); // force on particle 1
@@ -111,14 +112,17 @@ __device__ inline void computeDihedral(const TabulatedDihedralPotential* __restr
float t = (angle + BD_PI) * d->angle_step_inv;
int home = (int) floorf(t);
t = t - home;
+ //home = home % (d->size);
+ home = (home < d->size) ? home : d->size-1;
+ int home1 = (home + 1) >= d->size ? (home+1-d->size) : home+1;
- assert(home >= 0);
- assert(home < d->size);
+ //assert(home >= 0);
+ //assert(home < d->size);
// home = home % size;
- int home1 = (home + 1) >= d->size ? (home+1-d->size) : home+1;
+ //int home1 = (home + 1) >= d->size ? (home+1-d->size) : home+1;
- assert(home1 > 0);
- assert(home1 < d->size);
+ //assert(home1 >= 0);
+ //assert(home1 < d->size);
//================================================
// Linear interpolation
@@ -136,6 +140,15 @@ __device__ inline void computeDihedral(const TabulatedDihedralPotential* __restr
// printf("%f %d %d (%.4f %.4f) %.2f %f\n",force,home,home1, d->pot[home], d->pot[home1], dU, d->angle_step_inv);
//assert( force < 10000.0f )
+ //if( force != force )
+ //force = 0.f;
+ assert(force == force);
+ if ( force > 1000.0f )
+ force = 1000.0f;
+ if ( force < -1000.0f )
+ force = -1000.0f;
+ //assert( force < 10000.0f );
+
f1 *= force;
f2 *= force;
f3 *= force;
--
GitLab