Merge pull request #352 from swails/master

Fix long-range vdW correction for large(ish) systems

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -203,11 +203,11 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
     
     // Count the total number of particle pairs for each pair of classes.
     
-    map<pair<int, int>, int> interactionCount;
+    map<pair<int, int>, long long int> interactionCount;
     if (force.getNumInteractionGroups() == 0) {
         // Count the particles of each class.
         
-        vector<int> classCounts(numClasses, 0);
+        vector<long long int> classCounts(numClasses, 0);
         for (int i = 0; i < numParticles; i++)
             classCounts[atomClass[i]]++;
         for (int i = 0; i < numClasses; i++) {
@@ -246,9 +246,10 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
     for (int i = 0; i < numClasses; i++)
         for (int j = i; j < numClasses; j++)
             sum += interactionCount[make_pair(i, j)]*integrateInteraction(expression, classes[i], classes[j], force, context);
-    int numInteractions = (numParticles*(numParticles+1))/2;
+    double nPart = (double) numParticles;
+    double numInteractions = (nPart*(nPart+1))/2;
     sum /= numInteractions;
-    return 2*M_PI*numParticles*numParticles*sum;
+    return 2*M_PI*nPart*nPart*sum;
 }
 
 double CustomNonbondedForceImpl::integrateInteraction(Lepton::CompiledExpression& expression, const vector<double>& params1, const vector<double>& params2,

openmmapi/src/NonbondedForceImpl.cpp

@@ -239,7 +239,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
     for (map<pair<double, double>, int>::const_iterator entry = classCounts.begin(); entry != classCounts.end(); ++entry) {
         double sigma = entry->first.first;
         double epsilon = entry->first.second;
-        int count = (entry->second*(entry->second+1))/2;
+        double count = (double) entry->second;
+        count *= (count + 1) / 2;
         double sigma2 = sigma*sigma;
         double sigma6 = sigma2*sigma2*sigma2;
         sum1 += count*epsilon*sigma6*sigma6;
@@ -251,7 +252,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
         for (map<pair<double, double>, int>::const_iterator class2 = classCounts.begin(); class2 != class1; ++class2) {
             double sigma = 0.5*(class1->first.first+class2->first.first);
             double epsilon = sqrt(class1->first.second*class2->first.second);
-            int count = class1->second*class2->second;
+            double count = (double) class1->second;
+            count *= (double) class2->second;
             double sigma2 = sigma*sigma;
             double sigma6 = sigma2*sigma2*sigma2;
             sum1 += count*epsilon*sigma6*sigma6;
@@ -259,8 +261,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
             if (useSwitch)
                 sum3 += count*epsilon*(evalIntegral(cutoff, switchDist, cutoff, sigma)-evalIntegral(switchDist, switchDist, cutoff, sigma));
         }
-    int numParticles = system.getNumParticles();
-    int numInteractions = (numParticles*(numParticles+1))/2;
+    double numParticles = (double) system.getNumParticles();
+    double numInteractions = (numParticles*(numParticles+1))/2;
     sum1 /= numInteractions;
     sum2 /= numInteractions;
     sum3 /= numInteractions;