Completed adding softcore Vdw support to the AMOEBA Cuda platform, but further testing is needed

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -91,7 +91,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     for (int i = 0; i < force.getNumParticles(); i++) {
         double sigma, epsilon, reduction;
         bool isAlchemical;
-        // The variables reduction, ivindex are not used.
+        // The variables reduction, ivindex and isAlchemical are not used.
         int ivindex;
         // Get the sigma and epsilon parameters, ignoring everything else.
         force.getParticleParameters(i, ivindex, sigma, epsilon, reduction, isAlchemical);
@@ -152,8 +152,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     int ndelta = int(double(nstep) * (range - cut));
     double rdelta = (range - cut) / double(ndelta);
     double offset = cut - 0.5 * rdelta;
-    double dhal = 0.07; // This magic number also appears in kCalculateAmoebaCudaVdw14_7.cu
-    double ghal = 0.12; // This magic number also appears in kCalculateAmoebaCudaVdw14_7.cu
+
+    // Buffered-14-7 buffering constants
+    double dhal = 0.07; 
+    double ghal = 0.12;
+
     double elrc = 0.0; // This number is incremented and passed out at the end
     double e = 0.0;
     double sigma, epsilon; // The pairwise sigma and epsilon parameters.

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

@@ -2365,35 +2365,45 @@ CudaCalcAmoebaVdwForceKernel::~CudaCalcAmoebaVdwForceKernel() {
 void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
     cu.setAsCurrent();
     sigmaEpsilon.initialize<float2>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon");
+    isAlchemical.initialize<int>(cu, cu.getPaddedNumAtoms(), "isAlchemical");
+    vdwLambda.initialize<float>(cu, 1, "vdwLambda");
     bondReductionAtoms.initialize<int>(cu, cu.getPaddedNumAtoms(), "bondReductionAtoms");
-    bondReductionFactors.initialize<float>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon");
+    bondReductionFactors.initialize<float>(cu, cu.getPaddedNumAtoms(), "bondReductionFactors");
     tempPosq.initialize(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "tempPosq");
     tempForces.initialize<long long>(cu, 3*cu.getPaddedNumAtoms(), "tempForces");
     
     // Record atom parameters.
     
     vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
+    vector<int> isAlchemicalVec(cu.getPaddedNumAtoms(), 0);
     vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
     vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
     vector<vector<int> > exclusions(cu.getNumAtoms());
     for (int i = 0; i < force.getNumParticles(); i++) {
         int ivIndex;
         double sigma, epsilon, reductionFactor;
-        bool isAlchemical;
-        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, isAlchemical);
+        bool alchemical;
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical);
         sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
+        isAlchemicalVec[i] = (alchemical) ? 1 : 0;
         bondReductionAtomsVec[i] = ivIndex;
         bondReductionFactorsVec[i] = (float) reductionFactor;
         force.getParticleExclusions(i, exclusions[i]);
         exclusions[i].push_back(i);
     }
     sigmaEpsilon.upload(sigmaEpsilonVec);
+    isAlchemical.upload(isAlchemicalVec);
     bondReductionAtoms.upload(bondReductionAtomsVec);
     bondReductionFactors.upload(bondReductionFactorsVec);
     if (force.getUseDispersionCorrection())
         dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force);
     else
         dispersionCoefficient = 0.0;               
+
+    // A single lambda parameter to define the alchemical vdW state.
+    vector<float> vdwLambdaVec(1, 0);
+    vdwLambdaVec[0] = 1.0;
+    vdwLambda.upload(vdwLambdaVec);
  
     // This force is applied based on modified atom positions, where hydrogens have been moved slightly
     // closer to their parent atoms.  We therefore create a separate CudaNonbondedUtilities just for
@@ -2401,6 +2411,8 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
     
     nonbonded = new CudaNonbondedUtilities(cu);
     nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float", 2, sizeof(float2), sigmaEpsilon.getDevicePointer()));
+    nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("isAlchemical", "int", 1, sizeof(int), isAlchemical.getDevicePointer()));
+    nonbonded->addArgument(CudaNonbondedUtilities::ParameterInfo("vdwLambda", "float", 1, sizeof(float), vdwLambda.getDevicePointer()));
     
     // Create the interaction kernel.
     
@@ -2425,6 +2437,11 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
         replacements["EPSILON_COMBINING_RULE"] = "4";
     else
         throw OpenMMException("Illegal combining rule for sigma: "+sigmaCombiningRule);
+
+    replacements["VDW_SOFTCORE_POWER"] = cu.intToString(force.getSoftcorePower());
+    replacements["VDW_SOFTCORE_ALPHA"] = cu.doubleToString(force.getSoftcoreAlpha()); 
+    replacements["VDW_ALCHEMICAL_METHOD"] = cu.intToString(force.getAlchemicalMethod()); 
+
     double cutoff = force.getCutoffDistance();
     double taperCutoff = cutoff*0.9;
     replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoffDistance());
@@ -2451,6 +2468,12 @@ double CudaCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeF
         hasInitializedNonbonded = true;
         nonbonded->initialize(system);
     }
+
+    // Not sure about passing a Context parameter to the Kernal?
+    vector<float> vdwLambdaVec(1, 0);
+    vdwLambdaVec[0] = context.getParameter(AmoebaVdwForce::Lambda());
+    vdwLambda.upload(vdwLambdaVec);
+
     cu.getPosq().copyTo(tempPosq);
     cu.getForce().copyTo(tempForces);
     void* prepareArgs[] = {&cu.getForce().getDevicePointer(), &cu.getPosq().getDevicePointer(), &tempPosq.getDevicePointer(),
@@ -2474,20 +2497,22 @@ void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context,
         throw OpenMMException("updateParametersInContext: The number of particles has changed");
     
     // Record the per-particle parameters.
-    
     vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
+    vector<int> isAlchemicalVec(cu.getPaddedNumAtoms(), 0);
     vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
     vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
     for (int i = 0; i < force.getNumParticles(); i++) {
         int ivIndex;
         double sigma, epsilon, reductionFactor;
-        bool isAlchemical;
-        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, isAlchemical);
+        bool alchemical;
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical);
         sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
+        isAlchemicalVec[i] = (alchemical) ? 1 : 0;
         bondReductionAtomsVec[i] = ivIndex;
         bondReductionFactorsVec[i] = (float) reductionFactor;
     }
     sigmaEpsilon.upload(sigmaEpsilonVec);
+    isAlchemical.upload(isAlchemicalVec);
     bondReductionAtoms.upload(bondReductionAtomsVec);
     bondReductionFactors.upload(bondReductionFactorsVec);
     if (force.getUseDispersionCorrection())

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

@@ -575,7 +575,8 @@ private:
     CudaArray sigmaEpsilon;
     CudaArray bondReductionAtoms;
     CudaArray bondReductionFactors;
-    CudaArray isAlchemical
+    CudaArray isAlchemical;
+    CudaArray vdwLambda;
     CudaArray tempPosq;
     CudaArray tempForces;
     CudaNonbondedUtilities* nonbonded;

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

@@ -27,22 +27,18 @@
     real epsilon = (epsilon_s == 0.0f ? (real) 0 : 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s));
 #endif
 
-    real softcore = 0.0;
-#if USE_SOFTCORE == 1
+    real softcore = 0.0f;
+#if VDW_ALCHEMICAL_METHOD == 1
     // Decouple
-    bool both = isAlchemicalI && isAlchemicalJ;
-    bool either = isAlchemicalI || isAlchemicalJ;
-    bool soft = !both && either;
-    if (soft) { 
+    if (isAlchemical1 != isAlchemical2) { 
        epsilon = epsilon * pow(VDW_LAMBDA, VDW_SOFTCORE_POWER);
-       softcore = VDW_SOFTCORE_ALPHA * (1.0 - VDW_LAMBDA) * (1.0 - VDW_LAMBDA);
+       softcore = VDW_SOFTCORE_ALPHA * (1.0 - vdwLambda) * (1.0 - vdwLambda);
     }
-#else if USE_SOFTCORE == 2 
+#else if VDW_ALCHEMICAL_METHOD == 2 
     // Annihilate
-    bool soft = isAlchemicalI || isAlchemicalJ;
-    if (soft) {
+    if (isAlchemical1 || isAlchemical2) {
        epsilon = epsilon * pow(VDW_LAMBDA, VDW_SOFTCORE_POWER);
-       softcore = VDW_SOFTCORE_ALPHA * (1.0 - VDW_LAMBDA) * (1.0 - VDW_LAMBDA);
+       softcore = VDW_SOFTCORE_ALPHA * (1.0 - vdwLambda) * (1.0 - vdwLambda);
     }
 @endif
 

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp

@@ -313,24 +313,18 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
                 double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]);
                 double softcore = 0.0;
 
-                // Apply softcore modifications using an alchemical method.
-                if (this->_alchemicalMethod != None) {
-                   if (isAlchemicalI || isAlchemical[jj]) {
-                      // Decouple maintains interactions between two softcore atoms.
-                      if (this->_alchemicalMethod == Decouple && isAlchemicalI && isAlchemical[jj]) {
-                         // No softcore 
-                      } else {
-                         combinedEpsilon *= pow(lambda, this->_n); 
-                         softcore = this->_alpha * pow(1.0 - lambda, 2);
-                      } 
-                   }
+                if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemical[jj])) {
+                   combinedEpsilon *= pow(lambda, this->_n);
+                   softcore = this->_alpha * pow(1.0 - lambda, 2);
+                } else if (this->_alchemicalMethod == Annihilate && (isAlchemicalI || isAlchemical[jj])) {
+                   combinedEpsilon *= pow(lambda, this->_n);
+                   softcore = this->_alpha * pow(1.0 - lambda, 2);
                 }
 
                 Vec3 force;
-                energy                     += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
-                                                               reducedPositions[ii], reducedPositions[jj],
-                                                               force);
-                
+                energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
+                                           reducedPositions[ii], reducedPositions[jj], force);
+
                 if (indexIVs[ii] == ii) {
                     forces[ii][0] -= force[0];
                     forces[ii][1] -= force[1];
@@ -390,22 +384,17 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
         int isAlchemicalI      = isAlchemical[siteI];
         int isAlchemicalJ      = isAlchemical[siteJ];
 
-        // Apply softcore modifications using an alchemical method.
-        if (this->_alchemicalMethod != None) {
-           if (isAlchemicalI || isAlchemicalJ) {
-              // Decouple maintains interactions between two softcore atoms.
-              if (this->_alchemicalMethod == Decouple && isAlchemicalI && isAlchemicalJ) {
-                 // No softcore
-              } else {
-                 combinedEpsilon *= pow(lambda, this->_n);
-                 softcore = this->_alpha * pow(1.0 - lambda, 2);
-              }
-           }
+        if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemicalJ)) {
+           combinedEpsilon *= pow(lambda, this->_n);
+           softcore = this->_alpha * pow(1.0 - lambda, 2);
+        } else if (this->_alchemicalMethod == Annihilate && (isAlchemicalI || isAlchemicalJ)) {
+           combinedEpsilon *= pow(lambda, this->_n);
+           softcore = this->_alpha * pow(1.0 - lambda, 2);
         }
 
         Vec3 force;
-        energy                     += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
-                                                       reducedPositions[siteI], reducedPositions[siteJ], force);
+        energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
+                                   reducedPositions[siteI], reducedPositions[siteJ], force);
                 
         if (indexIVs[siteI] == siteI) {
             forces[siteI][0] -= force[0];