Merge pull request #1533 from peastman/amoeba

Made AMOEBA API more consistent

plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h

@@ -147,11 +147,34 @@ public:
      */
     void setCutoffDistance(double distance);
 
+    /**
+     * Get the parameters to use for PME calculations.  If alpha is 0 (the default), these parameters are
+     * ignored and instead their values are chosen based on the Ewald error tolerance.
+     *
+     * @param[out] alpha   the separation parameter
+     * @param[out] nx      the number of grid points along the X axis
+     * @param[out] ny      the number of grid points along the Y axis
+     * @param[out] nz      the number of grid points along the Z axis
+     */
+    void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+
+    /**
+     * Set the parameters to use for PME calculations.  If alpha is 0 (the default), these parameters are
+     * ignored and instead their values are chosen based on the Ewald error tolerance.
+     *
+     * @param alpha   the separation parameter
+     * @param nx      the number of grid points along the X axis
+     * @param ny      the number of grid points along the Y axis
+     * @param nz      the number of grid points along the Z axis
+     */
+    void setPMEParameters(double alpha, int nx, int ny, int nz);
+
     /**
      * Get the Ewald alpha parameter.  If this is 0 (the default), a value is chosen automatically
      * based on the Ewald error tolerance.
      *
      * @return the Ewald alpha parameter
+     * @deprecated This method exists only for backward compatibility.  Use getPMEParameters() instead.
      */
     double getAEwald() const;
 
@@ -160,6 +183,7 @@ public:
      * based on the Ewald error tolerance.
      *
      * @param aewald alpha parameter
+     * @deprecated This method exists only for backward compatibility.  Use setPMEParameters() instead.
      */
     void setAEwald(double aewald);
 
@@ -175,6 +199,7 @@ public:
      * are chosen automatically based on the Ewald error tolerance.
      *
      * @return the PME grid dimensions
+     * @deprecated This method exists only for backward compatibility.  Use getPMEParameters() instead.
      */
     void getPmeGridDimensions(std::vector<int>& gridDimension) const;
 
@@ -183,6 +208,7 @@ public:
      * are chosen automatically based on the Ewald error tolerance.
      *
      * @param gridDimension   the PME grid dimensions
+     * @deprecated This method exists only for backward compatibility.  Use setPMEParameters() instead.
      */
     void setPmeGridDimensions(const std::vector<int>& gridDimension);
 
@@ -408,9 +434,8 @@ private:
     NonbondedMethod nonbondedMethod;
     PolarizationType polarizationType;
     double cutoffDistance;
-    double aewald;
-    int pmeBSplineOrder;
-    std::vector<int> pmeGridDimension;
+    double alpha;
+    int pmeBSplineOrder, nx, ny, nz;
     int mutualInducedMaxIterations;
     std::vector<double> extrapolationCoefficients;
 

plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -183,13 +183,34 @@ public:
      */
     void getParticleExclusions(int particleIndex, std::vector<int>& exclusions) const;
 
+    /**
+     * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     *
+     * @return the cutoff distance, measured in nm
+     */
+
+    double getCutoffDistance() const;
+    
+    /**
+     * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     *
+     * @param distance    the cutoff distance, measured in nm
+     */
+    void setCutoffDistance(double distance);
+
     /**
      * Set the cutoff distance.
+     * 
+     * @deprecated This method exists only for backward compatibility.  Use setCutoffDistance() instead.
      */
     void setCutoff(double cutoff);
 
     /**
      * Get the cutoff distance.
+     * 
+     * @deprecated This method exists only for backward compatibility.  Use getCutoffDistance() instead.
      */
     double getCutoff() const;
 

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -40,9 +40,7 @@ using std::string;
 using std::vector;
 
 AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), polarizationType(Mutual), pmeBSplineOrder(5), cutoffDistance(1.0), ewaldErrorTol(1e-4), mutualInducedMaxIterations(60),
-                                               mutualInducedTargetEpsilon(1.0e-02), scalingDistanceCutoff(100.0), electricConstant(138.9354558456), aewald(0.0) {
-    pmeGridDimension.resize(3);
-    pmeGridDimension[0] = pmeGridDimension[1] = pmeGridDimension[2];
+                                               mutualInducedTargetEpsilon(1.0e-02), scalingDistanceCutoff(100.0), electricConstant(138.9354558456), alpha(0.0), nx(0), ny(0), nz(0) {
     extrapolationCoefficients.push_back(-0.154);
     extrapolationCoefficients.push_back(0.017);
     extrapolationCoefficients.push_back(0.658);
@@ -81,12 +79,26 @@ void AmoebaMultipoleForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
+void AmoebaMultipoleForce::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
+    alpha = this->alpha;
+    nx = this->nx;
+    ny = this->ny;
+    nz = this->nz;
+}
+
+void AmoebaMultipoleForce::setPMEParameters(double alpha, int nx, int ny, int nz) {
+    this->alpha = alpha;
+    this->nx = nx;
+    this->ny = ny;
+    this->nz = nz;
+}
+
 double AmoebaMultipoleForce::getAEwald() const { 
-    return aewald; 
+    return alpha; 
 } 
  
 void AmoebaMultipoleForce::setAEwald(double inputAewald) { 
-    aewald = inputAewald; 
+    alpha = inputAewald; 
 } 
  
 int AmoebaMultipoleForce::getPmeBSplineOrder() const { 
@@ -94,26 +106,18 @@ int AmoebaMultipoleForce::getPmeBSplineOrder() const {
 } 
  
 void AmoebaMultipoleForce::getPmeGridDimensions(std::vector<int>& gridDimension) const { 
-    if (gridDimension.size() < 3) {
+    if (gridDimension.size() < 3)
         gridDimension.resize(3);
-    }
-    if (pmeGridDimension.size() > 2) {
-        gridDimension[0] = pmeGridDimension[0];
-        gridDimension[1] = pmeGridDimension[1];
-        gridDimension[2] = pmeGridDimension[2];
-    } else {
-        gridDimension[0] = gridDimension[1] = gridDimension[2] = 0;
-    }
-    return;
+    gridDimension[0] = nx;
+    gridDimension[1] = ny;
+    gridDimension[2] = nz;
 } 
  
 void AmoebaMultipoleForce::setPmeGridDimensions(const std::vector<int>& gridDimension) {
-    pmeGridDimension.resize(3);
-    pmeGridDimension[0] = gridDimension[0];
-    pmeGridDimension[1] = gridDimension[1];
-    pmeGridDimension[2] = gridDimension[2];
-    return;
-} 
+    nx = gridDimension[0];
+    ny = gridDimension[1];
+    nz = gridDimension[2];
+}
 
 void AmoebaMultipoleForce::getPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const {
     dynamic_cast<const AmoebaMultipoleForceImpl&>(getImplInContext(context)).getPMEParameters(alpha, nx, ny, nz);

plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp

@@ -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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors:                                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -102,12 +102,20 @@ void AmoebaVdwForce::getParticleExclusions(int particleIndex, std::vector< int >
 
 }
 
-void AmoebaVdwForce::setCutoff(double inputCutoff) {
+double AmoebaVdwForce::getCutoffDistance() const {
+    return cutoff;
+}
+
+void AmoebaVdwForce::setCutoffDistance(double inputCutoff) {
     cutoff = inputCutoff;
 }
 
+void AmoebaVdwForce::setCutoff(double inputCutoff) {
+    setCutoffDistance(inputCutoff);
+}
+
 double AmoebaVdwForce::getCutoff() const {
-    return cutoff;
+    return getCutoffDistance();
 }
 
 AmoebaVdwForce::NonbondedMethod AmoebaVdwForce::getNonbondedMethod() const {

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -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:                                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -62,7 +62,7 @@ void AmoebaVdwForceImpl::initialize(ContextImpl& context) {
     if (owner.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic) {
         Vec3 boxVectors[3];
         system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-        double cutoff = owner.getCutoff();
+        double cutoff = owner.getCutoffDistance();
         if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
             throw OpenMMException("AmoebaVdwForce: The cutoff distance cannot be greater than half the periodic box size.");
     }   
@@ -103,7 +103,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     }
 
     // Compute the VdW tapering coefficients.  Mostly copied from amoebaCudaGpu.cpp.
-    double cutoff = force.getCutoff();
+    double cutoff = force.getCutoffDistance();
     double vdwTaper = 0.90; // vdwTaper is a scaling factor, it is not a distance.
     double c0 = 0.0;
     double c1 = 0.0;

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

@@ -1146,11 +1146,10 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
         coefficients << cu.doubleToString(sum);
     }
     defines["EXTRAPOLATION_COEFFICIENTS_SUM"] = coefficients.str();
-    alpha = force.getAEwald();
     if (usePME) {
-        vector<int> pmeGridDimension;
-        force.getPmeGridDimensions(pmeGridDimension);
-        if (pmeGridDimension[0] == 0 || alpha == 0.0) {
+        int nx, ny, nz;
+        force.getPMEParameters(alpha, nx, ny, nz);
+        if (nx == 0 || alpha == 0.0) {
             NonbondedForce nb;
             nb.setEwaldErrorTolerance(force.getEwaldErrorTolerance());
             nb.setCutoffDistance(force.getCutoffDistance());
@@ -1159,9 +1158,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
             gridSizeY = CudaFFT3D::findLegalDimension(gridSizeY);
             gridSizeZ = CudaFFT3D::findLegalDimension(gridSizeZ);
         } else {
-            gridSizeX = CudaFFT3D::findLegalDimension(pmeGridDimension[0]);
-            gridSizeY = CudaFFT3D::findLegalDimension(pmeGridDimension[1]);
-            gridSizeZ = CudaFFT3D::findLegalDimension(pmeGridDimension[2]);
+            gridSizeX = CudaFFT3D::findLegalDimension(nx);
+            gridSizeY = CudaFFT3D::findLegalDimension(ny);
+            gridSizeZ = CudaFFT3D::findLegalDimension(nz);
         }
         defines["EWALD_ALPHA"] = cu.doubleToString(alpha);
         defines["SQRT_PI"] = cu.doubleToString(sqrt(M_PI));
@@ -2601,15 +2600,15 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
         replacements["EPSILON_COMBINING_RULE"] = "4";
     else
         throw OpenMMException("Illegal combining rule for sigma: "+sigmaCombiningRule);
-    double cutoff = force.getCutoff();
+    double cutoff = force.getCutoffDistance();
     double taperCutoff = cutoff*0.9;
-    replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoff());
+    replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoffDistance());
     replacements["TAPER_CUTOFF"] = cu.doubleToString(taperCutoff);
     replacements["TAPER_C3"] = cu.doubleToString(10/pow(taperCutoff-cutoff, 3.0));
     replacements["TAPER_C4"] = cu.doubleToString(15/pow(taperCutoff-cutoff, 4.0));
     replacements["TAPER_C5"] = cu.doubleToString(6/pow(taperCutoff-cutoff, 5.0));
     bool useCutoff = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff);
-    nonbonded->addInteraction(useCutoff, useCutoff, true, force.getCutoff(), exclusions,
+    nonbonded->addInteraction(useCutoff, useCutoff, true, force.getCutoffDistance(), exclusions,
         cu.replaceStrings(CudaAmoebaKernelSources::amoebaVdwForce2, replacements), 0);
     
     // Create the other kernels.

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors:                                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -594,9 +594,9 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
     nonbondedMethod  = force.getNonbondedMethod();
     if (nonbondedMethod == AmoebaMultipoleForce::PME) {
         usePme     = true;
-        alphaEwald = force.getAEwald();
+        pmeGridDimension.resize(3);
+        force.getPMEParameters(alphaEwald, pmeGridDimension[0], pmeGridDimension[1], pmeGridDimension[2]);
         cutoffDistance = force.getCutoffDistance();
-        force.getPmeGridDimensions(pmeGridDimension);
         if (pmeGridDimension[0] == 0 || alphaEwald == 0.0) {
             NonbondedForce nb;
             nb.setEwaldErrorTolerance(force.getEwaldErrorTolerance());
@@ -996,7 +996,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
     epsilonCombiningRule   = force.getEpsilonCombiningRule();
     useCutoff              = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff);
     usePBC                 = (force.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic);
-    cutoff                 = force.getCutoff();
+    cutoff                 = force.getCutoffDistance();
     neighborList           = useCutoff ? new NeighborList() : NULL;
     dispersionCoefficient  = force.getUseDispersionCorrection() ?  AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
 

plugins/amoeba/serialization/src/AmoebaMultipoleForceProxy.cpp

@@ -74,20 +74,18 @@ void AmoebaMultipoleForceProxy::serialize(const void* object, SerializationNode&
     node.setIntProperty("forceGroup", force.getForceGroup());
     node.setIntProperty("nonbondedMethod",                  force.getNonbondedMethod());
     node.setIntProperty("polarizationType",                 force.getPolarizationType());
-    //node.setIntProperty("pmeBSplineOrder",                  force.getPmeBSplineOrder());
-    //node.setIntProperty("mutualInducedIterationMethod",     force.getMutualInducedIterationMethod());
     node.setIntProperty("mutualInducedMaxIterations",       force.getMutualInducedMaxIterations());
 
     node.setDoubleProperty("cutoffDistance",                force.getCutoffDistance());
-    node.setDoubleProperty("aEwald",                        force.getAEwald());
+    double alpha;
+    int nx, ny, nz;
+    force.getPMEParameters(alpha, nx, ny, nz);
+    node.setDoubleProperty("aEwald",                        alpha);
     node.setDoubleProperty("mutualInducedTargetEpsilon",    force.getMutualInducedTargetEpsilon());
-    //node.setDoubleProperty("electricConstant",              force.getElectricConstant());
     node.setDoubleProperty("ewaldErrorTolerance",           force.getEwaldErrorTolerance());
 
-    std::vector<int> gridDimensions;
-    force.getPmeGridDimensions(gridDimensions);
     SerializationNode& gridDimensionsNode  = node.createChildNode("MultipoleParticleGridDimension");
-    gridDimensionsNode.setIntProperty("d0", gridDimensions[0]).setIntProperty("d1", gridDimensions[1]).setIntProperty("d2", gridDimensions[2]); 
+    gridDimensionsNode.setIntProperty("d0", nx).setIntProperty("d1", ny).setIntProperty("d2", nz); 
     
     SerializationNode& coefficients = node.createChildNode("ExtrapolationCoefficients");
     vector<double> coeff = force.getExtrapolationCoefficients();
@@ -144,22 +142,14 @@ void* AmoebaMultipoleForceProxy::deserialize(const SerializationNode& node) cons
         force->setNonbondedMethod(static_cast<AmoebaMultipoleForce::NonbondedMethod>(node.getIntProperty("nonbondedMethod")));
         if (version >= 2)
             force->setPolarizationType(static_cast<AmoebaMultipoleForce::PolarizationType>(node.getIntProperty("polarizationType")));
-        //force->setPmeBSplineOrder(node.getIntProperty("pmeBSplineOrder"));
-        //force->setMutualInducedIterationMethod(static_cast<AmoebaMultipoleForce::MutualInducedIterationMethod>(node.getIntProperty("mutualInducedIterationMethod")));
         force->setMutualInducedMaxIterations(node.getIntProperty("mutualInducedMaxIterations"));
 
         force->setCutoffDistance(node.getDoubleProperty("cutoffDistance"));
-        force->setAEwald(node.getDoubleProperty("aEwald"));
         force->setMutualInducedTargetEpsilon(node.getDoubleProperty("mutualInducedTargetEpsilon"));
-        //force->setElectricConstant(node.getDoubleProperty("electricConstant"));
         force->setEwaldErrorTolerance(node.getDoubleProperty("ewaldErrorTolerance"));
 
-        std::vector<int> gridDimensions;
         const SerializationNode& gridDimensionsNode  = node.getChildNode("MultipoleParticleGridDimension");
-        gridDimensions.push_back(gridDimensionsNode.getIntProperty("d0"));
-        gridDimensions.push_back(gridDimensionsNode.getIntProperty("d1"));
-        gridDimensions.push_back(gridDimensionsNode.getIntProperty("d2"));
-        force->setPmeGridDimensions(gridDimensions);
+        force->setPMEParameters(node.getDoubleProperty("aEwald"), gridDimensionsNode.getIntProperty("d0"), gridDimensionsNode.getIntProperty("d1"), gridDimensionsNode.getIntProperty("d2"));
     
         if (version >= 3) {
             const SerializationNode& coefficients = node.getChildNode("ExtrapolationCoefficients");

plugins/amoeba/serialization/src/AmoebaVdwForceProxy.cpp

@@ -48,7 +48,7 @@ void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node)
     node.setIntProperty("forceGroup", force.getForceGroup());
     node.setStringProperty("SigmaCombiningRule", force.getSigmaCombiningRule());
     node.setStringProperty("EpsilonCombiningRule", force.getEpsilonCombiningRule());
-    node.setDoubleProperty("VdwCutoff", force.getCutoff());
+    node.setDoubleProperty("VdwCutoff", force.getCutoffDistance());
 
     node.setIntProperty("method", (int) force.getNonbondedMethod());
 
@@ -82,7 +82,7 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
             force->setForceGroup(node.getIntProperty("forceGroup", 0));
         force->setSigmaCombiningRule(node.getStringProperty("SigmaCombiningRule"));
         force->setEpsilonCombiningRule(node.getStringProperty("EpsilonCombiningRule"));
-        force->setCutoff(node.getDoubleProperty("VdwCutoff"));
+        force->setCutoffDistance(node.getDoubleProperty("VdwCutoff"));
         force->setNonbondedMethod((AmoebaVdwForce::NonbondedMethod) node.getIntProperty("method"));
 
         const SerializationNode& particles = node.getChildNode("VdwParticles");

plugins/amoeba/serialization/tests/TestSerializeAmoebaMultipoleForce.cpp

@@ -116,11 +116,8 @@ void testSerialization() {
     ASSERT_EQUAL(force1.getCutoffDistance(),                force2.getCutoffDistance());
     ASSERT_EQUAL(force1.getNonbondedMethod(),               force2.getNonbondedMethod());
     ASSERT_EQUAL(force1.getAEwald(),                        force2.getAEwald());
-    //ASSERT_EQUAL(force1.getPmeBSplineOrder(),               force2.getPmeBSplineOrder());
-    //ASSERT_EQUAL(force1.getMutualInducedIterationMethod(),  force2.getMutualInducedIterationMethod());
     ASSERT_EQUAL(force1.getMutualInducedMaxIterations(),    force2.getMutualInducedMaxIterations());
     ASSERT_EQUAL(force1.getMutualInducedTargetEpsilon(),    force2.getMutualInducedTargetEpsilon());
-    //ASSERT_EQUAL(force1.getElectricConstant(),              force2.getElectricConstant());
     ASSERT_EQUAL(force1.getEwaldErrorTolerance(),           force2.getEwaldErrorTolerance());