Rough draft AmoebaVdwForce Cuda code; next steps are to add CUDA kernal...

Rough draft AmoebaVdwForce Cuda code; next steps are to add CUDA kernal compilation support for USE_SOFTCORE, VDW_LAMBDA, VDW_SOFTCORE_POWER, VDW_SOFTCORE_ALPHA, and per particle isAlchemicalI / isAlchemicalJ flags

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h

@@ -547,7 +547,7 @@ public:
      * Initialize the kernel.
      * 
      * @param system     the System this kernel will be applied to
-     * @param force      the AmoebaMultipoleForce this kernel will be used for
+     * @param force      the AmoebaVdwForce this kernel will be used for
      */
     void initialize(const System& system, const AmoebaVdwForce& force);
     /**
@@ -575,6 +575,7 @@ private:
     CudaArray sigmaEpsilon;
     CudaArray bondReductionAtoms;
     CudaArray bondReductionFactors;
+    CudaArray isAlchemical
     CudaArray tempPosq;
     CudaArray tempForces;
     CudaNonbondedUtilities* nonbonded;

plugins/amoeba/platforms/cuda/src/kernels/amoebaVdwForce2.cu

-{
+
 #ifdef USE_CUTOFF
     unsigned int includeInteraction = (!isExcluded && r2 < CUTOFF_SQUARED);
 #else
@@ -26,20 +26,50 @@
     real epsilon_s = SQRT(sigmaEpsilon1.y) + SQRT(sigmaEpsilon2.y);
     real epsilon = (epsilon_s == 0.0f ? (real) 0 : 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s));
 #endif
-    real r6 = r2*r2*r2;
-    real r7 = r6*r;
-    real sigma7 = sigma*sigma;
-    sigma7 = sigma7*sigma7*sigma7*sigma;
-    real rho = r7 + sigma7*0.12f;
-    real invRho = RECIP(rho);
-    real tau = 1.07f/(r + 0.07f*sigma);
-    real tau7 = tau*tau*tau;
-    tau7 = tau7*tau7*tau;
-    real dTau = tau/1.07f;
-    real tmp = sigma7*invRho;
-    real gTau = epsilon*tau7*r6*1.12f*tmp*tmp;
-    real termEnergy = epsilon*sigma7*tau7*((sigma7*1.12f*invRho)-2.0f);
-    real deltaE = -7.0f*(dTau*termEnergy+gTau);
+
+    real softcore = 0.0;
+#if USE_SOFTCORE == 1
+    // Decouple
+    bool both = isAlchemicalI && isAlchemicalJ;
+    bool either = isAlchemicalI || isAlchemicalJ;
+    bool soft = !both && either;
+    if (soft) { 
+       epsilon = epsilon * pow(VDW_LAMBDA, VDW_SOFTCORE_POWER);
+       softcore = VDW_SOFTCORE_ALPHA * (1.0 - VDW_LAMBDA) * (1.0 - VDW_LAMBDA);
+    }
+#else if USE_SOFTCORE == 2 
+    // Annihilate
+    bool soft = isAlchemicalI || isAlchemicalJ;
+    if (soft) {
+       epsilon = epsilon * pow(VDW_LAMBDA, VDW_SOFTCORE_POWER);
+       softcore = VDW_SOFTCORE_ALPHA * (1.0 - VDW_LAMBDA) * (1.0 - VDW_LAMBDA);
+    }
+@endif
+
+    // The Buffered-14-7 buffering constants.
+    real dhal = 0.07;
+    real ghal = 0.12;
+    // The buffering constant plus one.
+    real dhal1 = 1.07;
+    real ghal1 = 1.12;
+
+    real rho = r / sigma;
+    real rho2 = rho * rho;
+    real rho6 = rho2 * rho2 * rho2;
+    real rhoplus = rho + dhal;
+    real rhodec2 = rhoplus * rhoplus;
+    real rhodec = rhodec2 * rhodec2 * rhodec2;
+    real s1 = 1.0 / (softcore + rhodec * rhoplus);
+    real s2 = 1.0 / (softcore + rho6 * rho + 0.12);
+    real point72 = dhal1 * dhal1;
+    real t1 = dhal1 * point72 * point72 * point72 * s1;
+    real t2 = ghal1 * s2;
+    real t2min = t2 - 2;
+    real dt1 = -7.0 * rhodec * t1 * s1;
+    real dt2 = -7.0 * rho6 * t2 * s2;
+    real termEnergy = epsilon * t1 * t2min;
+    real deltaE = epsilon * (dt1 * t2min + t1 * dt2) / sigma;
+
 #ifdef USE_CUTOFF
     if (r > TAPER_CUTOFF) {
         real x = r-TAPER_CUTOFF;