Skip to content

GitLab

Commit 713cc0f9 authored by Peter Eastman's avatar Peter Eastman
Browse files

CUDA implementation of periodic boundary conditions for AMOEBA bonded forces

parent 94fbbe9c
Hide whitespace changes
Inline Side-by-side
Showing with 532 additions and 256 deletions
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
......@@ -117,7 +117,7 @@ void CudaCalcAmoebaBondForceKernel::initialize(const System& system, const Amoeb
}
params->upload(paramVector);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = "0";
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = CudaAmoebaKernelSources::amoebaBondForce;
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params->getDevicePointer(), "float2");
replacements["CUBIC_K"] = cu.doubleToString(force.getAmoebaGlobalBondCubic());
......@@ -215,7 +215,7 @@ void CudaCalcAmoebaAngleForceKernel::initialize(const System& system, const Amoe
}
params->upload(paramVector);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = "0";
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = CudaAmoebaKernelSources::amoebaAngleForce;
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params->getDevicePointer(), "float2");
replacements["CUBIC_K"] = cu.doubleToString(force.getAmoebaGlobalAngleCubic());
......@@ -317,6 +317,7 @@ void CudaCalcAmoebaInPlaneAngleForceKernel::initialize(const System& system, con
}
params->upload(paramVector);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params->getDevicePointer(), "float2");
replacements["CUBIC_K"] = cu.doubleToString(force.getAmoebaGlobalInPlaneAngleCubic());
replacements["QUARTIC_K"] = cu.doubleToString(force.getAmoebaGlobalInPlaneAngleQuartic());
......@@ -419,6 +420,7 @@ void CudaCalcAmoebaPiTorsionForceKernel::initialize(const System& system, const
}
params->upload(paramVector);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params->getDevicePointer(), "float");
cu.getBondedUtilities().addInteraction(atoms, cu.replaceStrings(CudaAmoebaKernelSources::amoebaPiTorsionForce, replacements), force.getForceGroup());
cu.addForce(new ForceInfo(force));
......@@ -519,6 +521,7 @@ void CudaCalcAmoebaStretchBendForceKernel::initialize(const System& system, cons
params1->upload(paramVector);
params2->upload(paramVectorK);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params1->getDevicePointer(), "float3");
replacements["FORCE_CONSTANTS"] = cu.getBondedUtilities().addArgument(params2->getDevicePointer(), "float2");
replacements["RAD_TO_DEG"] = cu.doubleToString(180/M_PI);
......@@ -619,6 +622,7 @@ void CudaCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& system, c
}
params->upload(paramVector);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params->getDevicePointer(), "float");
replacements["CUBIC_K"] = cu.doubleToString(force.getAmoebaGlobalOutOfPlaneBendCubic());
replacements["QUARTIC_K"] = cu.doubleToString(force.getAmoebaGlobalOutOfPlaneBendQuartic());
......@@ -750,6 +754,7 @@ void CudaCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, c
gridValues->upload(gridValuesVec);
gridParams->upload(gridParamsVec);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["GRID_VALUES"] = cu.getBondedUtilities().addArgument(gridValues->getDevicePointer(), "float4");
replacements["GRID_PARAMS"] = cu.getBondedUtilities().addArgument(gridParams->getDevicePointer(), "float4");
replacements["TORSION_PARAMS"] = cu.getBondedUtilities().addArgument(torsionParams->getDevicePointer(), "int2");
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaInPlaneForce.cu
View file @ 713cc0f9
float2 angleParams = PARAMS[index];
real xad = pos1.x - pos4.x;
real yad = pos1.y - pos4.y;
real zad = pos1.z - pos4.z;
real3 ad = make_real3(pos1.x-pos4.x, pos1.y-pos4.y, pos1.z-pos4.z);
real3 bd = make_real3(pos2.x-pos4.x, pos2.y-pos4.y, pos2.z-pos4.z);
real3 cd = make_real3(pos3.x-pos4.x, pos3.y-pos4.y, pos3.z-pos4.z);
real xbd = pos2.x - pos4.x;
real ybd = pos2.y - pos4.y;
real zbd = pos2.z - pos4.z;
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ad)
APPLY_PERIODIC_TO_DELTA(bd)
APPLY_PERIODIC_TO_DELTA(cd)
#endif
real xcd = pos3.x - pos4.x;
real ycd = pos3.y - pos4.y;
real zcd = pos3.z - pos4.z;
real xt = yad*zcd - zad*ycd;
real yt = zad*xcd - xad*zcd;
real zt = xad*ycd - yad*xcd;
real xt = ad.y*cd.z - ad.z*cd.y;
real yt = ad.z*cd.x - ad.x*cd.z;
real zt = ad.x*cd.y - ad.y*cd.x;
real rt2 = xt*xt + yt*yt + zt*zt;
real delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2;
real delta = -(xt*bd.x + yt*bd.y + zt*bd.z) / rt2;
real xip = pos2.x + xt*delta;
real yip = pos2.y + yt*delta;
real zip = pos2.z + zt*delta;
real xap = pos1.x - xip;
real yap = pos1.y - yip;
real zap = pos1.z - zip;
real3 ap = make_real3(pos1.x-xip, pos1.y-yip, pos1.z-zip);
real3 cp = make_real3(pos3.x-xip, pos3.y-yip, pos3.z-zip);
real xcp = pos3.x - xip;
real ycp = pos3.y - yip;
real zcp = pos3.z - zip;
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ap)
APPLY_PERIODIC_TO_DELTA(cp)
#endif
real rap2 = xap*xap + yap*yap + zap*zap;
real rcp2 = xcp*xcp + ycp*ycp + zcp*zcp;
real rap2 = ap.x*ap.x + ap.y*ap.y + ap.z*ap.z;
real rcp2 = cp.x*cp.x + cp.y*cp.y + cp.z*cp.z;
real xm = ycp*zap - zcp*yap;
real ym = zcp*xap - xcp*zap;
real zm = xcp*yap - ycp*xap;
real xm = cp.y*ap.z - cp.z*ap.y;
real ym = cp.z*ap.x - cp.x*ap.z;
real zm = cp.x*ap.y - cp.y*ap.x;
real rm = max(SQRT(xm*xm + ym*ym + zm*zm), (real) 1e-6f);
real dotp = xap*xcp + yap*ycp + zap*zcp;
real dotp = ap.x*cp.x + ap.y*cp.y + ap.z*cp.z;
real product = SQRT(rap2*rcp2);
real cosine = (product > 0 ? (dotp/product) : 0);
cosine = max(min(cosine, (real) 1), (real) -1);
real angle;
if (cosine > 0.99f || cosine < -0.99f) {
real3 cross_prod = cross(make_real3(xap, yap, zap), make_real3(xcp, ycp, zcp));
real3 cross_prod = cross(ap, cp);
angle = ASIN(SQRT(dot(cross_prod, cross_prod)/(rap2*rcp2)))*RAD_TO_DEG;
if (cosine < 0.0f)
angle = 180-angle;
......@@ -67,13 +65,13 @@ dEdAngle *= RAD_TO_DEG;
real terma = -dEdAngle/(rap2*rm);
real termc = dEdAngle/(rcp2*rm);
real dedxia = terma * (yap*zm-zap*ym);
real dedyia = terma * (zap*xm-xap*zm);
real dedzia = terma * (xap*ym-yap*xm);
real dedxia = terma * (ap.y*zm-ap.z*ym);
real dedyia = terma * (ap.z*xm-ap.x*zm);
real dedzia = terma * (ap.x*ym-ap.y*xm);
real dedxic = termc * (ycp*zm-zcp*ym);
real dedyic = termc * (zcp*xm-xcp*zm);
real dedzic = termc * (xcp*ym-ycp*xm);
real dedxic = termc * (cp.y*zm-cp.z*ym);
real dedyic = termc * (cp.z*xm-cp.x*zm);
real dedzic = termc * (cp.x*ym-cp.y*xm);
real dedxip = -dedxia - dedxic;
real dedyip = -dedyia - dedyic;
......@@ -82,23 +80,23 @@ real dedzip = -dedzia - dedzic;
real delta2 = 2.0f*delta;
real ptrt2 = (dedxip*xt + dedyip*yt + dedzip*zt) / rt2;
real term = (zcd*ybd-ycd*zbd) + delta2*(yt*zcd-zt*ycd);
real dpdxia = delta*(ycd*dedzip-zcd*dedyip) + term*ptrt2;
real term = (cd.z*bd.y-cd.y*bd.z) + delta2*(yt*cd.z-zt*cd.y);
real dpdxia = delta*(cd.y*dedzip-cd.z*dedyip) + term*ptrt2;
term = (xcd*zbd-zcd*xbd) + delta2*(zt*xcd-xt*zcd);
real dpdyia = delta*(zcd*dedxip-xcd*dedzip) + term*ptrt2;
term = (cd.x*bd.z-cd.z*bd.x) + delta2*(zt*cd.x-xt*cd.z);
real dpdyia = delta*(cd.z*dedxip-cd.x*dedzip) + term*ptrt2;
term = (ycd*xbd-xcd*ybd) + delta2*(xt*ycd-yt*xcd);
real dpdzia = delta*(xcd*dedyip-ycd*dedxip) + term*ptrt2;
term = (cd.y*bd.x-cd.x*bd.y) + delta2*(xt*cd.y-yt*cd.x);
real dpdzia = delta*(cd.x*dedyip-cd.y*dedxip) + term*ptrt2;
term = (yad*zbd-zad*ybd) + delta2*(zt*yad-yt*zad);
real dpdxic = delta*(zad*dedyip-yad*dedzip) + term*ptrt2;
term = (ad.y*bd.z-ad.z*bd.y) + delta2*(zt*ad.y-yt*ad.z);
real dpdxic = delta*(ad.z*dedyip-ad.y*dedzip) + term*ptrt2;
term = (zad*xbd-xad*zbd) + delta2*(xt*zad-zt*xad);
real dpdyic = delta*(xad*dedzip-zad*dedxip) + term*ptrt2;
term = (ad.z*bd.x-ad.x*bd.z) + delta2*(xt*ad.z-zt*ad.x);
real dpdyic = delta*(ad.x*dedzip-ad.z*dedxip) + term*ptrt2;
term = (xad*ybd-yad*xbd) + delta2*(yt*xad-xt*yad);
real dpdzic = delta*(yad*dedxip-xad*dedyip) + term*ptrt2;
term = (ad.x*bd.y-ad.y*bd.x) + delta2*(yt*ad.x-xt*ad.y);
real dpdzic = delta*(ad.y*dedxip-ad.x*dedyip) + term*ptrt2;
dedxia = dedxia + dpdxia;
dedyia = dedyia + dpdyia;
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaOutOfPlaneBendForce.cu
View file @ 713cc0f9
// compute the value of the bond angle
real xab = pos1.x - pos2.x;
real yab = pos1.y - pos2.y;
real zab = pos1.z - pos2.z;
real xcb = pos3.x - pos2.x;
real ycb = pos3.y - pos2.y;
real zcb = pos3.z - pos2.z;
// compute the out-of-plane bending angle
real xdb = pos4.x - pos2.x;
real ydb = pos4.y - pos2.y;
real zdb = pos4.z - pos2.z;
real xad = pos1.x - pos4.x;
real yad = pos1.y - pos4.y;
real zad = pos1.z - pos4.z;
real xcd = pos3.x - pos4.x;
real ycd = pos3.y - pos4.y;
real zcd = pos3.z - pos4.z;
real rdb2 = xdb*xdb + ydb*ydb + zdb*zdb;
real rad2 = xad*xad + yad*yad + zad*zad;
real rcd2 = xcd*xcd + ycd*ycd + zcd*zcd;
real ee = xab*(ycb*zdb-zcb*ydb) + yab*(zcb*xdb-xcb*zdb) + zab*(xcb*ydb-ycb*xdb);
real dot = xad*xcd + yad*ycd + zad*zcd;
real3 ab = make_real3(pos1.x-pos2.x, pos1.y-pos2.y, pos1.z-pos2.z);
real3 cb = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
real3 db = make_real3(pos4.x-pos2.x, pos4.y-pos2.y, pos4.z-pos2.z);
real3 ad = make_real3(pos1.x-pos4.x, pos1.y-pos4.y, pos1.z-pos4.z);
real3 cd = make_real3(pos3.x-pos4.x, pos3.y-pos4.y, pos3.z-pos4.z);
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ab)
APPLY_PERIODIC_TO_DELTA(cb)
APPLY_PERIODIC_TO_DELTA(db)
APPLY_PERIODIC_TO_DELTA(ad)
APPLY_PERIODIC_TO_DELTA(cd)
#endif
real rdb2 = db.x*db.x + db.y*db.y + db.z*db.z;
real rad2 = ad.x*ad.x + ad.y*ad.y + ad.z*ad.z;
real rcd2 = cd.x*cd.x + cd.y*cd.y + cd.z*cd.z;
real ee = ab.x*(cb.y*db.z-cb.z*db.y) + ab.y*(cb.z*db.x-cb.x*db.z) + ab.z*(cb.x*db.y-cb.y*db.x);
real dot = ad.x*cd.x + ad.y*cd.y + ad.z*cd.z;
real cc = rad2*rcd2 - dot*dot;
real bkk2 = (cc != 0 ? (ee*ee)/(cc) : (real) 0);
bkk2 = rdb2 - bkk2;
real adXcd_0 = yad*zcd - zad*ycd;
real adXcd_1 = zad*xcd - xad*zcd;
real adXcd_2 = xad*ycd - yad*xcd;
real adXcd_0 = ad.y*cd.z - ad.z*cd.y;
real adXcd_1 = ad.z*cd.x - ad.x*cd.z;
real adXcd_2 = ad.x*cd.y - ad.y*cd.x;
real adXcd_nrm2 = adXcd_0*adXcd_0 + adXcd_1*adXcd_1 + adXcd_2*adXcd_2;
real adXcd_dot_db = xdb*adXcd_0 + ydb*adXcd_1 + zdb*adXcd_2;
real adXcd_dot_db = db.x*adXcd_0 + db.y*adXcd_1 + db.z*adXcd_2;
adXcd_dot_db /= SQRT(rdb2*adXcd_nrm2);
real angle = abs(ASIN(adXcd_dot_db));
......@@ -62,13 +54,13 @@ real dedcos = -deddt*eeSign/SQRT(cc*bkk2);
real term = ee / cc;
real dccdxia = (xad*rcd2-xcd*dot) * term;
real dccdyia = (yad*rcd2-ycd*dot) * term;
real dccdzia = (zad*rcd2-zcd*dot) * term;
real dccdxia = (ad.x*rcd2-cd.x*dot) * term;
real dccdyia = (ad.y*rcd2-cd.y*dot) * term;
real dccdzia = (ad.z*rcd2-cd.z*dot) * term;
real dccdxic = (xcd*rad2-xad*dot) * term;
real dccdyic = (ycd*rad2-yad*dot) * term;
real dccdzic = (zcd*rad2-zad*dot) * term;
real dccdxic = (cd.x*rad2-ad.x*dot) * term;
real dccdyic = (cd.y*rad2-ad.y*dot) * term;
real dccdzic = (cd.z*rad2-ad.z*dot) * term;
real dccdxid = -dccdxia - dccdxic;
real dccdyid = -dccdyia - dccdyic;
......@@ -76,17 +68,17 @@ real dccdzid = -dccdzia - dccdzic;
term = ee / rdb2;
real deedxia = ydb*zcb - zdb*ycb;
real deedyia = zdb*xcb - xdb*zcb;
real deedzia = xdb*ycb - ydb*xcb;
real deedxia = db.y*cb.z - db.z*cb.y;
real deedyia = db.z*cb.x - db.x*cb.z;
real deedzia = db.x*cb.y - db.y*cb.x;
real deedxic = yab*zdb - zab*ydb;
real deedyic = zab*xdb - xab*zdb;
real deedzic = xab*ydb - yab*xdb;
real deedxic = ab.y*db.z - ab.z*db.y;
real deedyic = ab.z*db.x - ab.x*db.z;
real deedzic = ab.x*db.y - ab.y*db.x;
real deedxid = ycb*zab - zcb*yab + xdb*term;
real deedyid = zcb*xab - xcb*zab + ydb*term;
real deedzid = xcb*yab - ycb*xab + zdb*term;
real deedxid = cb.y*ab.z - cb.z*ab.y + db.x*term;
real deedyid = cb.z*ab.x - cb.x*ab.z + db.y*term;
real deedzid = cb.x*ab.y - cb.y*ab.x + db.z*term;
// compute first derivative components for this angle
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaPiTorsionForce.cu
View file @ 713cc0f9
// compute the value of the pi-orbital torsion angle
real xad = pos1.x - pos4.x;
real yad = pos1.y - pos4.y;
real zad = pos1.z - pos4.z;
real xbd = pos2.x - pos4.x;
real ybd = pos2.y - pos4.y;
real zbd = pos2.z - pos4.z;
real xec = pos5.x - pos3.x;
real yec = pos5.y - pos3.y;
real zec = pos5.z - pos3.z;
real xgc = pos6.x - pos3.x;
real ygc = pos6.y - pos3.y;
real zgc = pos6.z - pos3.z;
real xip = yad*zbd - ybd*zad + pos3.x;
real yip = zad*xbd - zbd*xad + pos3.y;
real zip = xad*ybd - xbd*yad + pos3.z;
real xiq = yec*zgc - ygc*zec + pos4.x;
real yiq = zec*xgc - zgc*xec + pos4.y;
real ziq = xec*ygc - xgc*yec + pos4.z;
real3 ad = make_real3(pos1.x-pos4.x, pos1.y-pos4.y, pos1.z-pos4.z);
real3 bd = make_real3(pos2.x-pos4.x, pos2.y-pos4.y, pos2.z-pos4.z);
real3 ec = make_real3(pos5.x-pos3.x, pos5.y-pos3.y, pos5.z-pos3.z);
real3 gc = make_real3(pos6.x-pos3.x, pos6.y-pos3.y, pos6.z-pos3.z);
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ad)
APPLY_PERIODIC_TO_DELTA(bd)
APPLY_PERIODIC_TO_DELTA(ec)
APPLY_PERIODIC_TO_DELTA(gc)
#endif
real xip = ad.y*bd.z - bd.y*ad.z + pos3.x;
real yip = ad.z*bd.x - bd.z*ad.x + pos3.y;
real zip = ad.x*bd.y - bd.x*ad.y + pos3.z;
real xiq = ec.y*gc.z - gc.y*ec.z + pos4.x;
real yiq = ec.z*gc.x - gc.z*ec.x + pos4.y;
real ziq = ec.x*gc.y - gc.x*ec.y + pos4.z;
real xcp = pos3.x - xip;
real ycp = pos3.y - yip;
......@@ -112,21 +108,21 @@ real dedziq = ydc*dedxu - xdc*dedyu;
// compute first derivative components for individual atoms
real dedxia = ybd*dedzip - zbd*dedyip;
real dedyia = zbd*dedxip - xbd*dedzip;
real dedzia = xbd*dedyip - ybd*dedxip;
real dedxia = bd.y*dedzip - bd.z*dedyip;
real dedyia = bd.z*dedxip - bd.x*dedzip;
real dedzia = bd.x*dedyip - bd.y*dedxip;
real dedxib = zad*dedyip - yad*dedzip;
real dedyib = xad*dedzip - zad*dedxip;
real dedzib = yad*dedxip - xad*dedyip;
real dedxib = ad.z*dedyip - ad.y*dedzip;
real dedyib = ad.x*dedzip - ad.z*dedxip;
real dedzib = ad.y*dedxip - ad.x*dedyip;
real dedxie = ygc*dedziq - zgc*dedyiq;
real dedyie = zgc*dedxiq - xgc*dedziq;
real dedzie = xgc*dedyiq - ygc*dedxiq;
real dedxie = gc.y*dedziq - gc.z*dedyiq;
real dedyie = gc.z*dedxiq - gc.x*dedziq;
real dedzie = gc.x*dedyiq - gc.y*dedxiq;
real dedxig = zec*dedyiq - yec*dedziq;
real dedyig = xec*dedziq - zec*dedxiq;
real dedzig = yec*dedxiq - xec*dedyiq;
real dedxig = ec.z*dedyiq - ec.y*dedziq;
real dedyig = ec.x*dedziq - ec.z*dedxiq;
real dedzig = ec.y*dedxiq - ec.x*dedyiq;
dedxic = dedxic + dedxip - dedxie - dedxig;
dedyic = dedyic + dedyip - dedyie - dedyig;
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaStretchBendForce.cu
View file @ 713cc0f9
// compute the value of the bond angle
real xab = pos1.x - pos2.x;
real yab = pos1.y - pos2.y;
real zab = pos1.z - pos2.z;
real3 ab = make_real3(pos1.x-pos2.x, pos1.y-pos2.y, pos1.z-pos2.z);
real3 cb = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
real xcb = pos3.x - pos2.x;
real ycb = pos3.y - pos2.y;
real zcb = pos3.z - pos2.z;
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ab)
APPLY_PERIODIC_TO_DELTA(cb)
#endif
real rab = SQRT(xab*xab + yab*yab + zab*zab);
real rcb = SQRT(xcb*xcb + ycb*ycb + zcb*zcb);
real rab = SQRT(ab.x*ab.x + ab.y*ab.y + ab.z*ab.z);
real rcb = SQRT(cb.x*cb.x + cb.y*cb.y + cb.z*cb.z);
real xp = ycb*zab - zcb*yab;
real yp = zcb*xab - xcb*zab;
real zp = xcb*yab - ycb*xab;
real xp = cb.y*ab.z - cb.z*ab.y;
real yp = cb.z*ab.x - cb.x*ab.z;
real zp = cb.x*ab.y - cb.y*ab.x;
real rp = SQRT(xp*xp + yp*yp + zp*zp);
real dotp = xab*xcb + yab*ycb + zab*zcb;
real dotp = ab.x*cb.x + ab.y*cb.y + ab.z*cb.z;
real cosine = rab*rcb > 0 ? (dotp / (rab*rcb)) : (real) 1;
cosine = (cosine > 1 ? (real) 1 : cosine);
cosine = (cosine < -1 ? -(real) 1 : cosine);
real angle;
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(ab.x, ab.y, ab.z), make_real3(cb.x, cb.y, cb.z));
angle = ASIN(SQRT(dot(cross_prod, cross_prod))/(rab*rcb))*RAD_TO_DEG;
if (cosine < 0.0f)
angle = 180-angle;
......@@ -41,13 +41,13 @@ real dt = angle - RAD_TO_DEG*parameters.z;
real terma = rab*rp != 0 ? (-RAD_TO_DEG/(rab*rab*rp)) : (real) 0;
real termc = rcb*rp != 0 ? (RAD_TO_DEG/(rcb*rcb*rp)) : (real) 0;
real ddtdxia = terma * (yab*zp-zab*yp);
real ddtdyia = terma * (zab*xp-xab*zp);
real ddtdzia = terma * (xab*yp-yab*xp);
real ddtdxia = terma * (ab.y*zp-ab.z*yp);
real ddtdyia = terma * (ab.z*xp-ab.x*zp);
real ddtdzia = terma * (ab.x*yp-ab.y*xp);
real ddtdxic = termc * (ycb*zp-zcb*yp);
real ddtdyic = termc * (zcb*xp-xcb*zp);
real ddtdzic = termc * (xcb*yp-ycb*xp);
real ddtdxic = termc * (cb.y*zp-cb.z*yp);
real ddtdyic = termc * (cb.z*xp-cb.x*zp);
real ddtdzic = termc * (cb.x*yp-cb.y*xp);
// find chain rule terms for the bond length deviations
......@@ -61,13 +61,13 @@ real frc2 = ((rp != 0) ? force_constants.y : (real) 0);
real drkk = dr1*frc1 + dr2*frc2;
real ddrdxia = terma * xab;
real ddrdyia = terma * yab;
real ddrdzia = terma * zab;
real ddrdxia = terma * ab.x;
real ddrdyia = terma * ab.y;
real ddrdzia = terma * ab.z;
real ddrdxic = termc * xcb;
real ddrdyic = termc * ycb;
real ddrdzic = termc * zcb;
real ddrdxic = termc * cb.x;
real ddrdyic = termc * cb.y;
real ddrdzic = termc * cb.z;
// get the energy and master chain rule terms for derivatives
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaTorsionTorsionForce.cu
View file @ 713cc0f9
int2 torsionParams = TORSION_PARAMS[index];
real xba = pos2.x - pos1.x;
real yba = pos2.y - pos1.y;
real zba = pos2.z - pos1.z;
real xcb = pos3.x - pos2.x;
real ycb = pos3.y - pos2.y;
real zcb = pos3.z - pos2.z;
real xdc = pos4.x - pos3.x;
real ydc = pos4.y - pos3.y;
real zdc = pos4.z - pos3.z;
real xed = pos5.x - pos4.x;
real yed = pos5.y - pos4.y;
real zed = pos5.z - pos4.z;
real xt = yba*zcb - ycb*zba;
real yt = zba*xcb - zcb*xba;
real zt = xba*ycb - xcb*yba;
real xu = ycb*zdc - ydc*zcb;
real yu = zcb*xdc - zdc*xcb;
real zu = xcb*ydc - xdc*ycb;
real3 ba = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
real3 cb = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
real3 dc = make_real3(pos4.x-pos3.x, pos4.y-pos3.y, pos4.z-pos3.z);
real3 ed = make_real3(pos5.x-pos4.x, pos5.y-pos4.y, pos5.z-pos4.z);
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ba)
APPLY_PERIODIC_TO_DELTA(cb)
APPLY_PERIODIC_TO_DELTA(dc)
APPLY_PERIODIC_TO_DELTA(ed)
#endif
real xt = ba.y*cb.z - cb.y*ba.z;
real yt = ba.z*cb.x - cb.z*ba.x;
real zt = ba.x*cb.y - cb.x*ba.y;
real xu = cb.y*dc.z - dc.y*cb.z;
real yu = cb.z*dc.x - dc.z*cb.x;
real zu = cb.x*dc.y - dc.x*cb.y;
real rt2 = xt*xt + yt*yt + zt*zt;
real ru2 = xu*xu + yu*yu + zu*zu;
real rtru = SQRT(rt2 * ru2);
real xv = ydc*zed - yed*zdc;
real yv = zdc*xed - zed*xdc;
real zv = xdc*yed - xed*ydc;
real xv = dc.y*ed.z - ed.y*dc.z;
real yv = dc.z*ed.x - ed.z*dc.x;
real zv = dc.x*ed.y - ed.x*dc.y;
real rv2 = xv*xv + yv*yv + zv*zv;
real rurv = SQRT(ru2 * rv2);
real rcb = SQRT(xcb*xcb + ycb*ycb + zcb*zcb);
real rcb = SQRT(cb.x*cb.x + cb.y*cb.y + cb.z*cb.z);
real cosine1 = (rtru != 0 ? (xt*xu+yt*yu+zt*zu)/rtru : (real) 0);
cosine1 = (cosine1 > 1 ? (real) 1 : cosine1);
cosine1 = (cosine1 < -1 ? (real) -1 : cosine1);
......@@ -51,11 +47,11 @@ if (cosine1 > 0.99f || cosine1 < -0.99f) {
}
else
angle1 = RAD_TO_DEG*ACOS(cosine1);
real sign = xba*xu + yba*yu + zba*zu;
real sign = ba.x*xu + ba.y*yu + ba.z*zu;
angle1 = (sign < 0 ? -angle1 : angle1);
real value1 = angle1;
real rdc = SQRT(xdc*xdc + ydc*ydc + zdc*zdc);
real rdc = SQRT(dc.x*dc.x + dc.y*dc.y + dc.z*dc.z);
real cosine2 = (xu*xv + yu*yv + zu*zv) / rurv;
cosine2 = (cosine2 > 1 ? (real) 1 : cosine2);
cosine2 = (cosine2 < -1 ? (real) -1 : cosine2);
......@@ -70,7 +66,7 @@ if (cosine2 > 0.99f || cosine2 < -0.99f) {
}
else
angle2 = RAD_TO_DEG*ACOS(cosine2);
sign = xcb*xv + ycb*yv + zcb*zv;
sign = cb.x*xv + cb.y*yv + cb.z*zv;
angle2 = (sign < 0 ? -angle2 : angle2);
real value2 = angle2;
......@@ -83,24 +79,20 @@ real value2 = angle2;
int chiralAtomIndex = (torsionParams.x > -1 ? torsionParams.x : atom5);
real4 pos6 = posq[chiralAtomIndex];
real xac = pos6.x - pos3.x;
real yac = pos6.y - pos3.y;
real zac = pos6.z - pos3.z;
real xbc = pos2.x - pos3.x;
real ybc = pos2.y - pos3.y;
real zbc = pos2.z - pos3.z;
real3 ac = make_real3(pos6.x-pos3.x, pos6.y-pos3.y, pos6.z-pos3.z);
real3 bc = make_real3(pos2.x-pos3.x, pos2.y-pos3.y, pos2.z-pos3.z);
real3 dc1 = make_real3(pos4.x-pos3.x, pos4.y-pos3.y, pos4.z-pos3.z);
// xdc, ydc, zdc appear above
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ac)
APPLY_PERIODIC_TO_DELTA(bc)
APPLY_PERIODIC_TO_DELTA(dc1)
#endif
real xdc1 = pos4.x - pos3.x;
real ydc1 = pos4.y - pos3.y;
real zdc1 = pos4.z - pos3.z;
real c1 = ybc*zdc1 - zbc*ydc1;
real c2 = ydc1*zac - zdc1*yac;
real c3 = yac*zbc - zac*ybc;
real vol = xac*c1 + xbc*c2 + xdc1*c3;
real c1 = bc.y*dc1.z - bc.z*dc1.y;
real c2 = dc1.y*ac.z - dc1.z*ac.y;
real c3 = ac.y*bc.z - ac.z*bc.y;
real vol = ac.x*c1 + bc.x*c2 + dc1.x*c3;
sign = (vol > 0 ? (real) 1 : (real) -1);
sign = (torsionParams.x < 0 ? (real) 1 : sign);
value1 *= sign;
......@@ -170,66 +162,68 @@ dedang2 *= sign * RAD_TO_DEG;
// chain rule terms for first angle derivative components
real xca = pos3.x - pos1.x;
real yca = pos3.y - pos1.y;
real zca = pos3.z - pos1.z;
real3 ca = make_real3(pos3.x-pos1.x, pos3.y-pos1.y, pos3.z-pos1.z);
real3 db = make_real3(pos4.x-pos2.x, pos4.y-pos2.y, pos4.z-pos2.z);
real xdb = pos4.x - pos2.x;
real ydb = pos4.y - pos2.y;
real zdb = pos4.z - pos2.z;
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ca)
APPLY_PERIODIC_TO_DELTA(db)
#endif
real dedxt = dedang1 * (yt*zcb - ycb*zt) / (rt2*rcb);
real dedyt = dedang1 * (zt*xcb - zcb*xt) / (rt2*rcb);
real dedzt = dedang1 * (xt*ycb - xcb*yt) / (rt2*rcb);
real dedxu = -dedang1 * (yu*zcb - ycb*zu) / (ru2*rcb);
real dedyu = -dedang1 * (zu*xcb - zcb*xu) / (ru2*rcb);
real dedzu = -dedang1 * (xu*ycb - xcb*yu) / (ru2*rcb);
real dedxt = dedang1 * (yt*cb.z - cb.y*zt) / (rt2*rcb);
real dedyt = dedang1 * (zt*cb.x - cb.z*xt) / (rt2*rcb);
real dedzt = dedang1 * (xt*cb.y - cb.x*yt) / (rt2*rcb);
real dedxu = -dedang1 * (yu*cb.z - cb.y*zu) / (ru2*rcb);
real dedyu = -dedang1 * (zu*cb.x - cb.z*xu) / (ru2*rcb);
real dedzu = -dedang1 * (xu*cb.y - cb.x*yu) / (ru2*rcb);
// compute first derivative components for first angle
real dedxia = zcb*dedyt - ycb*dedzt;
real dedyia = xcb*dedzt - zcb*dedxt;
real dedzia = ycb*dedxt - xcb*dedyt;
real dedxia = cb.z*dedyt - cb.y*dedzt;
real dedyia = cb.x*dedzt - cb.z*dedxt;
real dedzia = cb.y*dedxt - cb.x*dedyt;
real dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu;
real dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu;
real dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu;
real dedxib = ca.y*dedzt - ca.z*dedyt + dc.z*dedyu - dc.y*dedzu;
real dedyib = ca.z*dedxt - ca.x*dedzt + dc.x*dedzu - dc.z*dedxu;
real dedzib = ca.x*dedyt - ca.y*dedxt + dc.y*dedxu - dc.x*dedyu;
real dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu;
real dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu;
real dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu;
real dedxic = ba.z*dedyt - ba.y*dedzt + db.y*dedzu - db.z*dedyu;
real dedyic = ba.x*dedzt - ba.z*dedxt + db.z*dedxu - db.x*dedzu;
real dedzic = ba.y*dedxt - ba.x*dedyt + db.x*dedyu - db.y*dedxu;
real dedxid = zcb*dedyu - ycb*dedzu;
real dedyid = xcb*dedzu - zcb*dedxu;
real dedzid = ycb*dedxu - xcb*dedyu;
real dedxid = cb.z*dedyu - cb.y*dedzu;
real dedyid = cb.x*dedzu - cb.z*dedxu;
real dedzid = cb.y*dedxu - cb.x*dedyu;
// chain rule terms for second angle derivative components
real xec = pos5.x - pos3.x;
real yec = pos5.y - pos3.y;
real zec = pos5.z - pos3.z;
real3 ec = make_real3(pos5.x-pos3.x, pos5.y-pos3.y, pos5.z-pos3.z);
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(ec)
#endif
real dedxu2 = dedang2 * (yu*zdc - ydc*zu) / (ru2*rdc);
real dedyu2 = dedang2 * (zu*xdc - zdc*xu) / (ru2*rdc);
real dedzu2 = dedang2 * (xu*ydc - xdc*yu) / (ru2*rdc);
real dedxv2 = -dedang2 * (yv*zdc - ydc*zv) / (rv2*rdc);
real dedyv2 = -dedang2 * (zv*xdc - zdc*xv) / (rv2*rdc);
real dedzv2 = -dedang2 * (xv*ydc - xdc*yv) / (rv2*rdc);
real dedxu2 = dedang2 * (yu*dc.z - dc.y*zu) / (ru2*rdc);
real dedyu2 = dedang2 * (zu*dc.x - dc.z*xu) / (ru2*rdc);
real dedzu2 = dedang2 * (xu*dc.y - dc.x*yu) / (ru2*rdc);
real dedxv2 = -dedang2 * (yv*dc.z - dc.y*zv) / (rv2*rdc);
real dedyv2 = -dedang2 * (zv*dc.x - dc.z*xv) / (rv2*rdc);
real dedzv2 = -dedang2 * (xv*dc.y - dc.x*yv) / (rv2*rdc);
// compute first derivative components for second angle
real dedxib2 = zdc*dedyu2 - ydc*dedzu2;
real dedyib2 = xdc*dedzu2 - zdc*dedxu2;
real dedzib2 = ydc*dedxu2 - xdc*dedyu2;
real dedxic2 = ydb*dedzu2 - zdb*dedyu2 + zed*dedyv2 - yed*dedzv2;
real dedyic2 = zdb*dedxu2 - xdb*dedzu2 + xed*dedzv2 - zed*dedxv2;
real dedzic2 = xdb*dedyu2 - ydb*dedxu2 + yed*dedxv2 - xed*dedyv2;
real dedxid2 = zcb*dedyu2 - ycb*dedzu2 + yec*dedzv2 - zec*dedyv2;
real dedyid2 = xcb*dedzu2 - zcb*dedxu2 + zec*dedxv2 - xec*dedzv2;
real dedzid2 = ycb*dedxu2 - xcb*dedyu2 + xec*dedyv2 - yec*dedxv2;
real dedxie2 = zdc*dedyv2 - ydc*dedzv2;
real dedyie2 = xdc*dedzv2 - zdc*dedxv2;
real dedzie2 = ydc*dedxv2 - xdc*dedyv2;
real dedxib2 = dc.z*dedyu2 - dc.y*dedzu2;
real dedyib2 = dc.x*dedzu2 - dc.z*dedxu2;
real dedzib2 = dc.y*dedxu2 - dc.x*dedyu2;
real dedxic2 = db.y*dedzu2 - db.z*dedyu2 + ed.z*dedyv2 - ed.y*dedzv2;
real dedyic2 = db.z*dedxu2 - db.x*dedzu2 + ed.x*dedzv2 - ed.z*dedxv2;
real dedzic2 = db.x*dedyu2 - db.y*dedxu2 + ed.y*dedxv2 - ed.x*dedyv2;
real dedxid2 = cb.z*dedyu2 - cb.y*dedzu2 + ec.y*dedzv2 - ec.z*dedyv2;
real dedyid2 = cb.x*dedzu2 - cb.z*dedxu2 + ec.z*dedxv2 - ec.x*dedzv2;
real dedzid2 = cb.y*dedxu2 - cb.x*dedyu2 + ec.x*dedyv2 - ec.y*dedxv2;
real dedxie2 = dc.z*dedyv2 - dc.y*dedzv2;
real dedyie2 = dc.x*dedzv2 - dc.z*dedxv2;
real dedzie2 = dc.y*dedxv2 - dc.x*dedyv2;
real3 force1 = make_real3(-dedxia, -dedyia, -dedzia);
real3 force2 = make_real3(-dedxib-dedxib2, -dedyib-dedyib2, -dedzib-dedzib2);
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaAngleForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -35,6 +35,7 @@
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomAngleForce.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
......@@ -273,6 +274,56 @@ void testOneAngle() {
compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions, then compare to an identical custom force.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numParticles = 3;
for (int ii = 0; ii < numParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaAngleForce* amoebaAngleForce = new AmoebaAngleForce();
double angle = 100.0;
double quadraticK = 1.0;
double cubicK = 1.0e-01;
double quarticK = 1.0e-02;
double penticK = 1.0e-03;
double sexticK = 1.0e-04;
amoebaAngleForce->addAngle(0, 1, 2, angle, quadraticK);
amoebaAngleForce->setAmoebaGlobalAngleCubic(cubicK);
amoebaAngleForce->setAmoebaGlobalAngleQuartic(quarticK);
amoebaAngleForce->setAmoebaGlobalAnglePentic(penticK);
amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK);
amoebaAngleForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaAngleForce);
CustomAngleForce* customForce = new CustomAngleForce("k2*delta^2 + k3*delta^3 + k4*delta^4 + k5*delta^5 + k6*delta^6; delta=theta-theta0");
customForce->addGlobalParameter("theta0", angle*M_PI/180);
customForce->addGlobalParameter("k2", quadraticK*pow(180/M_PI, 2.0));
customForce->addGlobalParameter("k3", cubicK*pow(180/M_PI, 3.0));
customForce->addGlobalParameter("k4", quarticK*pow(180/M_PI, 4.0));
customForce->addGlobalParameter("k5", penticK*pow(180/M_PI, 5.0));
customForce->addGlobalParameter("k6", sexticK*pow(180/M_PI, 6.0));
customForce->addAngle(0, 1, 2);
customForce->setUsesPeriodicBoundaryConditions(true);
customForce->setForceGroup(1);
system.addForce(customForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 2);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy, true, 1);
State s2 = context.getState(State::Forces | State::Energy, true, 2);
ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaAngleForce running test..." << std::endl;
......@@ -280,6 +331,7 @@ int main(int argc, char* argv[]) {
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testOneAngle();
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaBondForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -36,6 +36,7 @@
#include "openmm/internal/AssertionUtilities.h"
#include "CudaPlatform.h"
#include "openmm/Context.h"
#include "openmm/CustomBondForce.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
......@@ -199,6 +200,49 @@ void testTwoBond() {
compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions, then compare to an identical custom force.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numParticles = 2;
for (int ii = 0; ii < numParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaBondForce* amoebaBondForce = new AmoebaBondForce();
double bondLength = 1.5;
double quadraticK = 1.0;
double cubicK = 2.0;
double quarticK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic(quarticK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
amoebaBondForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaBondForce);
CustomBondForce* customForce = new CustomBondForce("k2*delta^2 + k3*delta^3 + k4*delta^4; delta=r-r0");
customForce->addGlobalParameter("r0", bondLength);
customForce->addGlobalParameter("k2", quadraticK);
customForce->addGlobalParameter("k3", cubicK);
customForce->addGlobalParameter("k4", quarticK);
customForce->addBond(0, 1);
customForce->setUsesPeriodicBoundaryConditions(true);
customForce->setForceGroup(1);
system.addForce(customForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 0);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy, true, 1);
State s2 = context.getState(State::Forces | State::Energy, true, 2);
ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaBondForce running test..." << std::endl;
......@@ -206,6 +250,7 @@ int main(int argc, char* argv[]) {
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testTwoBond();
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaInPlaneAngleForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -346,6 +346,50 @@ void testOneAngle() {
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 4;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaInPlaneAngleForce* amoebaInPlaneAngleForce = new AmoebaInPlaneAngleForce();
double angle = 65.0;
double quadraticK = 1.0;
double cubicK = 0.0e-01;
double quarticK = 0.0e-02;
double penticK = 0.0e-03;
double sexticK = 0.0e-04;
amoebaInPlaneAngleForce->addAngle(0, 1, 2, 3, angle, quadraticK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleCubic(cubicK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleQuartic(quarticK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAnglePentic(penticK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK);
amoebaInPlaneAngleForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaInPlaneAngleForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 1);
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[2] = Vec3(0, 0, -2);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaInPlaneAngleForce running test..." << std::endl;
......@@ -353,6 +397,7 @@ int main(int argc, char* argv[]) {
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testOneAngle();
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -448,6 +448,43 @@ void testOneOutOfPlaneBend2(int setId) {
}
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 4;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic(0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic(-0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
double kOutOfPlaneBend = 0.328682196E-01;
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
amoebaOutOfPlaneBendForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaOutOfPlaneBendForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
positions[3] = Vec3(0, 0, 1);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[3] = Vec3(0, 0, -2);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaOutOfPlaneBendForce running test..." << std::endl;
......@@ -456,6 +493,7 @@ int main(int argc, char* argv[]) {
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testOneOutOfPlaneBend();
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPiTorsionForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -288,6 +288,41 @@ void testOnePiTorsion() {
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 6;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
double kTorsion = 6.85;
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion);
amoebaPiTorsionForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaPiTorsionForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 0.5);
positions[3] = Vec3(0.4, 0.4, 0.4);
positions[4] = Vec3(1, 0, 1);
positions[5] = Vec3(1, 1, 0);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[0] = Vec3(0, -2, 0);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaPiTorsionForce running test..." << std::endl;
......@@ -295,6 +330,7 @@ int main(int argc, char* argv[]) {
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testOnePiTorsion();
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -266,6 +266,41 @@ void testOneStretchBend() {
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 3;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaStretchBendForce* amoebaStretchBendForce = new AmoebaStretchBendForce();
double abLength = 0.144800000E+01;
double cbLength = 0.101500000E+01;
double angleStretchBend = 0.108500000E+03*DegreesToRadians;
double kStretchBend = 1.0;
amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
amoebaStretchBendForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaStretchBendForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 1);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[2] = Vec3(0, 0, -2);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaStretchBendForce running test..." << std::endl;
......@@ -273,6 +308,7 @@ int main(int argc, char* argv[]) {
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testOneStretchBend();
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -2667,6 +2667,44 @@ void testTorsionTorsion(int systemId) {
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 6;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaTorsionTorsionForce* amoebaTorsionTorsionForce = new AmoebaTorsionTorsionForce();
int chiralCheckAtomIndex;
int gridIndex;
chiralCheckAtomIndex = 5;
gridIndex = 2;
amoebaTorsionTorsionForce->addTorsionTorsion(0, 1, 2, 3, 4, chiralCheckAtomIndex, 0);
amoebaTorsionTorsionForce->setTorsionTorsionGrid(0, getTorsionGrid(gridIndex));
amoebaTorsionTorsionForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaTorsionTorsionForce);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 0.5);
positions[3] = Vec3(0.4, 0.4, 0.4);
positions[4] = Vec3(1, 0, 1);
positions[5] = Vec3(1, 1, 0);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[0] = Vec3(0, -2, 0);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
......@@ -2675,6 +2713,7 @@ int main(int argc, char* argv[]) {
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
testTorsionTorsion(1);
testPeriodic();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaAngleForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaBondForce.cpp
View file @ 713cc0f9
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......
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