Commit 103618fb authored by Jason Swails's avatar Jason Swails
Browse files

Ugh... stop shadowing dot so I can use it as a function. I hate obscure compiler

errors.
parent 9d6aa7ae
...@@ -39,9 +39,9 @@ real ym = zcp*xap - xcp*zap; ...@@ -39,9 +39,9 @@ real ym = zcp*xap - xcp*zap;
real zm = xcp*yap - ycp*xap; real zm = xcp*yap - ycp*xap;
real rm = max(SQRT(xm*xm + ym*ym + zm*zm), (real) 1e-6f); real rm = max(SQRT(xm*xm + ym*ym + zm*zm), (real) 1e-6f);
real dot = xap*xcp + yap*ycp + zap*zcp; real dotp = xap*xcp + yap*ycp + zap*zcp;
real product = SQRT(rap2*rcp2); real product = SQRT(rap2*rcp2);
real cosine = (product > 0 ? (dot/product) : 0); real cosine = (product > 0 ? (dotp/product) : 0);
cosine = max(min(cosine, (real) 1), (real) -1); cosine = max(min(cosine, (real) 1), (real) -1);
real angle; real angle;
if (cosine > 0.99f || cosine < -0.99f) { if (cosine > 0.99f || cosine < -0.99f) {
...@@ -51,7 +51,7 @@ if (cosine > 0.99f || cosine < -0.99f) { ...@@ -51,7 +51,7 @@ if (cosine > 0.99f || cosine < -0.99f) {
angle = 180-angle; angle = 180-angle;
} }
else else
real angle = ACOS(cosine)*RAD_TO_DEG; angle = ACOS(cosine)*RAD_TO_DEG;
// if product == 0, set force/energy to 0 // if product == 0, set force/energy to 0
......
...@@ -17,19 +17,20 @@ real zp = xcb*yab - ycb*xab; ...@@ -17,19 +17,20 @@ real zp = xcb*yab - ycb*xab;
real rp = SQRT(xp*xp + yp*yp + zp*zp); real rp = SQRT(xp*xp + yp*yp + zp*zp);
real dot = xab*xcb + yab*ycb + zab*zcb; real dotp = xab*xcb + yab*ycb + zab*zcb;
real cosine = rab*rcb > 0 ? (dot / (rab*rcb)) : (real) 1; real cosine = rab*rcb > 0 ? (dotp / (rab*rcb)) : (real) 1;
cosine = (cosine > 1 ? (real) 1 : cosine); cosine = (cosine > 1 ? (real) 1 : cosine);
cosine = (cosine < -1 ? -(real) 1 : cosine); cosine = (cosine < -1 ? -(real) 1 : cosine);
real angle; real angle;
if (cosine > 0.99f || cosine < -0.99f) { if (cosine > 0.99f || cosine < -0.99f) {
// Highly unlikely a stretch-bend angle will be near 0 or 180, but just in case...
real3 cross_prod = cross(make_real3(xab, yab, zab), make_real3(xcb, ycb, zcb)); real3 cross_prod = cross(make_real3(xab, yab, zab), make_real3(xcb, ycb, zcb));
angle = ASIN(SQRT(dot(cross_prod, cross_prod)/(rap2*rcp2)))*RAD_TO_DEG; angle = ASIN(SQRT(dot(cross_prod, cross_prod))/(rab*rcb))*RAD_TO_DEG;
if (cosine < 0.0f) if (cosine < 0.0f)
angle = 180-angle; angle = 180-angle;
} }
else else
real angle = ACOS(cosine)*RAD_TO_DEG; angle = ACOS(cosine)*RAD_TO_DEG;
// find chain rule terms for the bond angle deviation // find chain rule terms for the bond angle deviation
......
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