Commit 14125aa9 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed bugs in OpenCL CustomHbondForce

parent f2012ad8
......@@ -2596,10 +2596,14 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
vector<mm_int4> acceptorBufferVector(numAcceptors);
vector<int> donorBufferCounter(numParticles, 0);
for (int i = 0; i < numDonors; i++)
donorBufferVector[i] = mm_int4(donorBufferCounter[donorVector[i].x]++, donorBufferCounter[donorVector[i].y]++, donorBufferCounter[donorVector[i].z]++, 0);
donorBufferVector[i] = mm_int4(donorVector[i].x > -1 ? donorBufferCounter[donorVector[i].x]++ : 0,
donorVector[i].y > -1 ? donorBufferCounter[donorVector[i].y]++ : 0,
donorVector[i].z > -1 ? donorBufferCounter[donorVector[i].z]++ : 0, 0);
vector<int> acceptorBufferCounter(numParticles, 0);
for (int i = 0; i < numAcceptors; i++)
acceptorBufferVector[i] = mm_int4(acceptorBufferCounter[acceptorVector[i].x]++, acceptorBufferCounter[acceptorVector[i].y]++, acceptorBufferCounter[acceptorVector[i].z]++, 0);
acceptorBufferVector[i] = mm_int4(acceptorVector[i].x > -1 ? acceptorBufferCounter[acceptorVector[i].x]++ : 0,
acceptorVector[i].y > -1 ? acceptorBufferCounter[acceptorVector[i].y]++ : 0,
acceptorVector[i].z > -1 ? acceptorBufferCounter[acceptorVector[i].z]++ : 0, 0);
donorBufferIndices->upload(donorBufferVector);
acceptorBufferIndices->upload(acceptorBufferVector);
int maxBuffers = 1;
......
......@@ -26,15 +26,15 @@ float4 deltaPeriodic(float4 vec1, float4 vec2) {
* Compute the angle between two vectors. The w component of each vector should contain the squared magnitude.
*/
float computeAngle(float4 vec1, float4 vec2) {
float dot = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
float cosine = dot/sqrt(vec1.w*vec2.w);
float dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
float cosine = dotProduct/sqrt(vec1.w*vec2.w);
float angle;
if (cosine > 0.99f || cosine < -0.99f) {
// We're close to the singularity in acos(), so take the cross product and use asin() instead.
float4 cross_prod = cross(vec1, vec2);
float4 crossProduct = cross(vec1, vec2);
float scale = vec1.w*vec2.w;
angle = asin(sqrt(dot(cross_prod, cross_prod)/scale));
angle = asin(sqrt(dot(crossProduct, crossProduct)/scale));
if (cosine < 0.0f)
angle = M_PI-angle;
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment