Continuing implementation of PME reciprocal space calculation

plugins/amoeba/openmmapi/include/AmoebaMultipoleForce.h

@@ -113,20 +113,6 @@ public:
      */
     void setCutoffDistance(double distance);
 
-    /**
-     * Get the aEwald parameter
-     *
-     * @return the Ewald parameter
-     */
-    double getAEwald() const;
-
-    /**
-     * Set the aEwald parameter
-     *
-     * @param Ewald parameter
-     */
-    void setAEwald(double aewald);
-
     /**
      * Add multipole-related info for a particle 
      *
@@ -272,6 +258,20 @@ public:
      * @param the electric constant
      */
     void setElectricConstant( double inputElectricConstant ); 
+    /**
+     * Get the error tolerance for Ewald summation.  This corresponds to the fractional error in the forces
+     * which is acceptable.  This value is used to select the reciprocal space cutoff and separation
+     * parameter so that the average error level will be less than the tolerance.  There is not a
+     * rigorous guarantee that all forces on all atoms will be less than the tolerance, however.
+     */
+    double getEwaldErrorTolerance() const;
+    /**
+     * Get the error tolerance for Ewald summation.  This corresponds to the fractional error in the forces
+     * which is acceptable.  This value is used to select the reciprocal space cutoff and separation
+     * parameter so that the average error level will be less than the tolerance.  There is not a
+     * rigorous guarantee that all forces on all atoms will be less than the tolerance, however.
+     */
+    void setEwaldErrorTolerance(double tol);
     
 protected:
     ForceImpl* createImpl();
@@ -279,27 +279,15 @@ private:
 
     AmoebaNonbondedMethod nonbondedMethod;
     double cutoffDistance;
-    double aewald;
     MutualInducedIterationMethod mutualInducedIterationMethod;
     int mutualInducedMaxIterations;
     double mutualInducedTargetEpsilon;
     double scalingDistanceCutoff;
     double electricConstant;
+    double ewaldErrorTol;
 
     class MultipoleInfo;
-
-// Retarded visual studio compiler complains about being unable to 
-// export private stl class members.
-// This stanza explains that it should temporarily shut up.
-#if defined(_MSC_VER)
-#pragma warning(push)
-#pragma warning(disable:4251)
-#endif
     std::vector<MultipoleInfo> multipoles;
-#if defined(_MSC_VER)
-#pragma warning(pop)
-#endif
-
 };
 
 class AmoebaMultipoleForce::MultipoleInfo {

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -36,7 +36,7 @@
 
 using namespace OpenMM;
 
-AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), aewald(0.0), mutualInducedIterationMethod(SOR), mutualInducedMaxIterations(60),
+AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), ewaldErrorTol(5e-4), mutualInducedIterationMethod(SOR), mutualInducedMaxIterations(60),
                                                mutualInducedTargetEpsilon(1.0e-05), scalingDistanceCutoff(100.0), electricConstant(138.9354558456) {
 }
 
@@ -56,14 +56,6 @@ void AmoebaMultipoleForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
-double AmoebaMultipoleForce::getAEwald() const {
-    return aewald;
-}
-
-void AmoebaMultipoleForce::setAEwald(double inputAewald ) {
-    aewald = inputAewald;
-}
-
 AmoebaMultipoleForce::MutualInducedIterationMethod AmoebaMultipoleForce::getMutualInducedIterationMethod( void ) const {
     return mutualInducedIterationMethod;
 }
@@ -104,6 +96,14 @@ void AmoebaMultipoleForce::setElectricConstant( double inputElectricConstant ) {
     electricConstant = inputElectricConstant;
 }
 
+double AmoebaMultipoleForce::getEwaldErrorTolerance() const {
+    return ewaldErrorTol;
+}
+
+void AmoebaMultipoleForce::setEwaldErrorTolerance(double tol) {
+    ewaldErrorTol = tol;
+}
+
 int AmoebaMultipoleForce::addParticle( double charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole, int axisType, 
                                         int multipoleAtomId1, int multipoleAtomId2, double thole, double dampingFactor, double polarity) {
     multipoles.push_back(MultipoleInfo( charge, molecularDipole, molecularQuadrupole,  axisType, multipoleAtomId1,  multipoleAtomId2, thole, dampingFactor, polarity));

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

@@ -31,6 +31,7 @@
 #include "kernels/amoebaCudaKernels.h"
 #include "internal/AmoebaMultipoleForceImpl.h"
 #include "internal/AmoebaWcaDispersionForceImpl.h"
+#include "openmm/internal/NonbondedForceImpl.h"
 
 #include <cmath>
 #ifdef _MSC_VER
@@ -588,8 +589,16 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
                                     static_cast<float>( force.getMutualInducedTargetEpsilon()),
                                     nonbondedMethod,
                                     static_cast<float>( force.getCutoffDistance()),
-                                    static_cast<float>( force.getAEwald()),
                                     static_cast<float>( force.getElectricConstant()) );
+    if (nonbondedMethod == AmoebaMultipoleForce::PME) {
+        double alpha;
+        int xsize, ysize, zsize;
+        NonbondedForce nb;
+        nb.setEwaldErrorTolerance(force.getEwaldErrorTolerance());
+        nb.setCutoffDistance(force.getCutoffDistance());
+        NonbondedForceImpl::calcPMEParameters(system, nb, alpha, xsize, ysize, zsize);
+        gpuSetAmoebaPMEParameters(data.getAmoebaGpu(), (float) alpha, xsize, ysize, zsize);
+    }
 
 }
 

plugins/amoeba/platforms/cuda/src/kernels/AmoebaGpu.cpp

@@ -351,7 +351,6 @@ void gpuPrintCudaAmoebaGmxSimulation(amoebaGpuContext amoebaGpu, FILE* log )
     (void) fprintf( log, "     pM_ScaleIndices                    %p\n",      amoebaGpu->amoebaSim.pM_ScaleIndices );
     (void) fprintf( log, "     sqrtPi                             %15.7e\n",  amoebaGpu->amoebaSim.sqrtPi );
     (void) fprintf( log, "     cutoffDistance2                    %15.7e\n",  amoebaGpu->amoebaSim.cutoffDistance2 );
-    (void) fprintf( log, "     aewald                             %15.7e\n",  amoebaGpu->amoebaSim.aewald );
     (void) fprintf( log, "     electric                           %15.7e\n",  amoebaGpu->amoebaSim.electric );
     (void) fprintf( log, "     box                                %15.7e %15.7e %15.7e\n", gpu->sim.periodicBoxSizeX, gpu->sim.periodicBoxSizeY, gpu->sim.periodicBoxSizeZ);
     (void) fprintf( log, "     gkc                                %15.7e\n",  amoebaGpu->amoebaSim.gkc );
@@ -1464,7 +1463,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
                                      const std::vector< std::vector< std::vector<int> > >& multipoleParticleCovalentInfo, const std::vector<int>& covalentDegree,
                                      const std::vector<int>& minCovalentIndices,  const std::vector<int>& minCovalentPolarizationIndices, int maxCovalentRange, 
                                      int mutualInducedIterativeMethod, int mutualInducedMaxIterations, float mutualInducedTargetEpsilon,
-                                     int nonbondedMethod, float cutoffDistance, float aewald, float electricConstant ){
+                                     int nonbondedMethod, float cutoffDistance, float electricConstant ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -1565,10 +1564,8 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
     }
     amoebaGpu->amoebaSim.cutoffDistance2             = cutoffDistance*cutoffDistance;
     amoebaGpu->amoebaSim.sqrtPi                      = sqrt( 3.1415926535897932384626433832795 );
-    amoebaGpu->amoebaSim.aewald                      = aewald;
     amoebaGpu->amoebaSim.electric                    = electricConstant;
 
-    amoebaGpu->gpuContext->sim.alphaEwald            = aewald;
     amoebaGpu->gpuContext->sim.nonbondedCutoff       = cutoffDistance;
     tabulateErfc(amoebaGpu->gpuContext); 
 

plugins/amoeba/platforms/cuda/src/kernels/amoebaCudaTypes.h

@@ -128,7 +128,6 @@ struct cudaAmoebaGmxSimulation {
     unsigned int paddedNumberOfAtoms;               // padded number of atoms
     float cutoffDistance2;                          // cutoff distance squared for PME
     float sqrtPi;                                   // sqrt(PI)
-    float aewald;                                   // aewald parameter
     float scalingDistanceCutoff;                    // scaling cutoff
     float2*         pDampingFactorAndThole;         // Thole & damping factors
 

plugins/amoeba/platforms/cuda/src/kernels/amoebaGpuTypes.h

@@ -154,7 +154,6 @@ struct _amoebaGpuContext {
 
     int multipoleNonbondedMethod;
     double cutoffDistance;
-    double aewald;
 
     // mutual induced field
 
@@ -308,7 +307,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
                                      const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo, const std::vector<int>& covalentDegree,
                                      const std::vector<int>& minCovalentIndices,  const std::vector<int>& minCovalentPolarizationIndices, int maxCovalentRange,
                                      int mutualInducedIterationMethod, int mutualInducedMaxIterations, float mutualInducedTargetEpsilon,
-                                     int nonbondedMethod, float cutoffDistance, float aewald, float electricConstant );
+                                     int nonbondedMethod, float cutoffDistance, float electricConstant );
 
 
 extern "C"
@@ -326,6 +325,9 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
                                 const std::string& sigmaCombiningRule,
                                 const std::string& epsilonCombiningRule,
                                 const std::vector< std::vector<int> >& allExclusions );
+extern "C"
+void gpuSetAmoebaPMEParameters(amoebaGpuContext amoebaGpu, float alpha, int gridSizeX, int gridSizeY, int gridSizeZ);
+
 extern "C"
 void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu,  const std::vector< std::vector<int> >& exclusions );
 

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

@@ -45,7 +45,7 @@ void kReduceDirectSelfFields_kernel( unsigned int fieldComponents, unsigned int
 
     // Reduce field
 
-    const float term = (4.0f/3.0f)*(cAmoebaSim.aewald*cAmoebaSim.aewald*cAmoebaSim.aewald)/cAmoebaSim.sqrtPi;
+    const float term = (4.0f/3.0f)*(cSim.alphaEwald*cSim.alphaEwald*cSim.alphaEwald)/cAmoebaSim.sqrtPi;
     while (pos < fieldComponents)
     {   
 
@@ -123,12 +123,12 @@ __device__ void calculateFixedFieldRealSpacePairIxn_kernel( FixedFieldParticle&
 
     // calculate the error function damping terms
 
-    float ralpha      = cAmoebaSim.aewald*r;
+    float ralpha      = cSim.alphaEwald*r;
     float bn[4];
 
     bn[0]             = erfc(ralpha)/r;
-    float alsq2       = 2.0f*cAmoebaSim.aewald*cAmoebaSim.aewald;
-    float alsq2n      = 1.0f/(cAmoebaSim.sqrtPi*cAmoebaSim.aewald);
+    float alsq2       = 2.0f*cSim.alphaEwald*cSim.alphaEwald;
+    float alsq2n      = 1.0f/(cAmoebaSim.sqrtPi*cSim.alphaEwald);
     float exp2a       = exp(-(ralpha*ralpha));
     alsq2n           *= alsq2;
     bn[1]             = (bn[0]+alsq2n*exp2a)/r2;

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

@@ -212,14 +212,14 @@ __device__ void calculateRealSpaceEwaldPairIxn_kernel( RealSpaceEwaldParticle& a
       
         // calculate the real space error function terms;
 
-        float ralpha = cAmoebaSim.aewald*r;
+        float ralpha = cSim.alphaEwald*r;
 
                bn[0] = fastErfc(ralpha)/r;
 
-        float alsq2  = 2.0f*cAmoebaSim.aewald*cAmoebaSim.aewald;
+        float alsq2  = 2.0f*cSim.alphaEwald*cSim.alphaEwald;
         float alsq2n = 0.0f;
-        if( cAmoebaSim.aewald > 0.0f){
-            alsq2n = 1.0f/(cAmoebaSim.sqrtPi*cAmoebaSim.aewald);
+        if( cSim.alphaEwald > 0.0f){
+            alsq2n = 1.0f/(cAmoebaSim.sqrtPi*cSim.alphaEwald);
         }
         float exp2a  = exp(-(ralpha*ralpha));
 

plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp

@@ -2007,7 +2007,6 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
     } else {
         multipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff );
     }
-    multipoleForce->setAEwald( aewald );
     multipoleForce->setCutoffDistance( cutoffDistance );
     system.setDefaultPeriodicBoxVectors( Vec3(box[0], 0.0, 0.0), Vec3(0.0, box[1], 0.0), Vec3(0.0, 0.0, box[2]) );
     if( log ){
@@ -2112,7 +2111,6 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
         double polarityConversion      = AngstromToNm*AngstromToNm*AngstromToNm;
         double dampingFactorConversion = sqrt( AngstromToNm );
       
-        multipoleForce->setAEwald(                multipoleForce->getAEwald()/AngstromToNm);
         multipoleForce->setCutoffDistance(        multipoleForce->getCutoffDistance()*AngstromToNm );
         multipoleForce->setScalingDistanceCutoff( multipoleForce->getScalingDistanceCutoff()*AngstromToNm );
 
@@ -2167,8 +2165,8 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
                         (useOpenMMUnits ? "OpenMM" : "Amoeba") );
 
         std::string nonbondedMethod = multipoleForce->getNonbondedMethod( ) == AmoebaMultipoleForce::PME ? "PME" : "NoCutoff";
-        (void) fprintf( log, "NonbondedMethod=%s aEwald=%15.7e cutoff=%15.7e.\n", nonbondedMethod.c_str(),
-                         multipoleForce->getAEwald(), multipoleForce->getCutoffDistance() );
+        (void) fprintf( log, "NonbondedMethod=%s cutoff=%15.7e.\n", nonbondedMethod.c_str(),
+                         multipoleForce->getCutoffDistance() );
 
         Vec3 a,b,c;
         system.getDefaultPeriodicBoxVectors( a, b, c );

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPME.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Ewald summation method cuda implementation of NonbondedForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "AmoebaMultipoleForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testPMEWater() {
+    Platform& platform = Platform::getPlatformByName("Cuda");
+    System system;
+    system.addParticle(16);
+    system.addParticle(1);
+    system.addParticle(1);
+    VerletIntegrator integrator(0.01);
+    AmoebaMultipoleForce* mp = new AmoebaMultipoleForce();
+    mp->setNonbondedMethod(AmoebaMultipoleForce::PME);
+    vector<double> dipole(3, 0.0);
+    dipole[2] = 7.556121361e-2;
+    vector<double> quadrupole(9, 0.0);
+    quadrupole[0] = 3.540307211e-2;
+    quadrupole[4] = -3.902570771e-2;
+    quadrupole[8] = 3.622635596e-3;
+    mp->addParticle(-0.51966, dipole, quadrupole, 1, 1, 2, 0.39, 9.707801995e-1, 0.837);
+    dipole[0] = -2.042094848e-2;
+    dipole[2] = -3.078753000e-2;
+    quadrupole[0] = -3.428482490e-3;
+    quadrupole[2] = -1.894859639e-4;
+    quadrupole[4] = -1.002408752e-2;
+    quadrupole[6] = -1.894859639e-4;
+    quadrupole[8] = 1.345257001e-2;
+    mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 2, 0.39, 8.897068742e-1, 0.496);
+    mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 1, 0.39, 8.897068742e-1, 0.496);
+    mp->setCutoffDistance(1.0);
+//    nonbonded->setEwaldErrorTolerance(TOL);
+    system.setDefaultPeriodicBoxVectors(Vec3(2, 0, 0), Vec3(0, 2, 0), Vec3(0, 0, 2));
+    system.addForce(mp);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    const double angle = 109.47*M_PI/180;
+    const double dOH = 0.1;
+    positions[0] = Vec3();
+    positions[1] = Vec3(dOH, 0, 0);
+    positions[2] = Vec3(dOH*std::cos(angle), dOH*std::sin(angle), 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+}
+
+int main() {
+    try {
+        Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
+        testPMEWater();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+