Modified calculateElectrostaticPotentialForAtomGridPoint_kernel to apply...

Modified calculateElectrostaticPotentialForAtomGridPoint_kernel to apply periodic boundary conditions
Unit test for multipole grid potential changed to system comprised of 216 water molecules; 
Tinker potential values now based on PBC

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

@@ -144,6 +144,10 @@ __device__ void calculateElectrostaticPotentialForAtomGridPoint_kernel( volatile
     float yr                 = atomI.y - gridPoint.y;
     float zr                 = atomI.z - gridPoint.z;
 
+    xr                      -= floorf(xr*cSim.invPeriodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
+    yr                      -= floorf(yr*cSim.invPeriodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
+    zr                      -= floorf(zr*cSim.invPeriodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
+
     float r2                 = xr*xr + yr*yr + zr*zr;
     float r                  = sqrtf( r2 );
 

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

@@ -941,8 +941,7 @@ static void testQuadrupoleValidation( FILE* log ){
 static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::AmoebaNonbondedMethod nonbondedMethod,
                                                           AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
                                                           double cutoff, int inputPmeGridDimension, std::string testName,
-                                                          std::vector<Vec3>& forces, double& energy,
-                                                          std::vector<Vec3>& inputGrid, std::vector< double >& outputGridPotential, FILE* log ){
+                                                          std::vector<Vec3>& forces, double& energy, FILE* log ){
 
     // beginning of Multipole setup
 
@@ -1123,13 +1122,9 @@ static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::
 
     context.setPositions(positions);
 
-    if( testName == "testMultipoleGridPotential" ){
-        amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
-    } else {
     State state               = context.getState(State::Forces | State::Energy);
     forces                    = state.getForces();
     energy                    = state.getPotentialEnergy();
-    }
 
     return;
 
@@ -1148,12 +1143,8 @@ static void testMultipoleIonsAndWaterPMEDirectPolarization( FILE* log ) {
     std::vector<Vec3> forces;
     double energy;
 
-    std::vector<Vec3> inputGrid;
-    std::vector<double> outputGridPotential;
-
     setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Direct, 
-                                                  cutoff, inputPmeGridDimension, testName, forces, energy, 
-                                                  inputGrid, outputGridPotential, log );
+                                                  cutoff, inputPmeGridDimension, testName, forces, energy, log );
 
     std::vector<Vec3> expectedForces(numberOfParticles);
 
@@ -1187,11 +1178,9 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) {
     double energy;
 
     std::vector<Vec3> inputGrid;
-    std::vector<double> outputGridPotential;
 
     setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, 
-                                                  cutoff, inputPmeGridDimension, testName, forces, energy,
-                                                  inputGrid, outputGridPotential, log );
+                                                  cutoff, inputPmeGridDimension, testName, forces, energy, log );
 
     std::vector<Vec3> expectedForces(numberOfParticles);
 
@@ -1211,101 +1200,6 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) {
     double tolerance                = 5.0e-04;
     compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
 
-}
-// test computation of multipole potential on a grid 
-
-static void testMultipoleGridPotential( FILE* log ) {
-
-    std::string testName      = "testMultipoleGridPotential";
-    
-    int numberOfParticles     = 8;
-    int inputPmeGridDimension = 64;
-    double cutoff             = 0.70;
-
-    std::vector<Vec3> forces;
-    double energy;
-
-    // initialize grid
-
-    int gridSize  = 27;
-    std::vector<Vec3> inputGrid(gridSize);
-    inputGrid[0]  = Vec3( -3.2894000e+00,  -1.3098000e+00,  -1.7064092e+00);
-    inputGrid[1]  = Vec3( -3.2894000e+00,  -1.3098000e+00,   6.3928525e-01);
-    inputGrid[2]  = Vec3( -3.2894000e+00,  -1.3098000e+00,   2.9849797e+00);
-    inputGrid[3]  = Vec3( -3.2894000e+00,   5.1360000e-01,  -1.7064092e+00);
-    inputGrid[4]  = Vec3( -3.2894000e+00,   5.1360000e-01,   6.3928525e-01);
-    inputGrid[5]  = Vec3( -3.2894000e+00,   5.1360000e-01,   2.9849797e+00);
-    inputGrid[6]  = Vec3( -3.2894000e+00,   2.3370000e+00,  -1.7064092e+00);
-    inputGrid[7]  = Vec3( -3.2894000e+00,   2.3370000e+00,   6.3928525e-01);
-    inputGrid[8]  = Vec3( -3.2894000e+00,   2.3370000e+00,   2.9849797e+00);
-    inputGrid[9]  = Vec3( -2.3754000e+00,  -1.3098000e+00,  -1.7064092e+00);
-    inputGrid[10] = Vec3( -2.3754000e+00,  -1.3098000e+00,   6.3928525e-01);
-    inputGrid[11] = Vec3( -2.3754000e+00,  -1.3098000e+00,   2.9849797e+00);
-    inputGrid[12] = Vec3( -2.3754000e+00,   5.1360000e-01,  -1.7064092e+00);
-    inputGrid[13] = Vec3( -2.3754000e+00,   5.1360000e-01,   6.3928525e-01);
-    inputGrid[14] = Vec3( -2.3754000e+00,   5.1360000e-01,   2.9849797e+00);
-    inputGrid[15] = Vec3( -2.3754000e+00,   2.3370000e+00,  -1.7064092e+00);
-    inputGrid[16] = Vec3( -2.3754000e+00,   2.3370000e+00,   6.3928525e-01);
-    inputGrid[17] = Vec3( -2.3754000e+00,   2.3370000e+00,   2.9849797e+00);
-    inputGrid[18] = Vec3( -1.4614000e+00,  -1.3098000e+00,  -1.7064092e+00);
-    inputGrid[19] = Vec3( -1.4614000e+00,  -1.3098000e+00,   6.3928525e-01);
-    inputGrid[20] = Vec3( -1.4614000e+00,  -1.3098000e+00,   2.9849797e+00);
-    inputGrid[21] = Vec3( -1.4614000e+00,   5.1360000e-01,  -1.7064092e+00);
-    inputGrid[22] = Vec3( -1.4614000e+00,   5.1360000e-01,   6.3928525e-01);
-    inputGrid[23] = Vec3( -1.4614000e+00,   5.1360000e-01,   2.9849797e+00);
-    inputGrid[24] = Vec3( -1.4614000e+00,   2.3370000e+00,  -1.7064092e+00);
-    inputGrid[25] = Vec3( -1.4614000e+00,   2.3370000e+00,   6.3928525e-01);
-    inputGrid[26] = Vec3( -1.4614000e+00,   2.3370000e+00,   2.9849797e+00);
-
-    std::vector<double> outputGridPotential;
-
-    setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, 
-                                                  cutoff, inputPmeGridDimension, testName, forces, energy,
-                                                  inputGrid, outputGridPotential, log );
-
-    // TINKER computed grid values
-
-    std::vector<double> tinkerGridPotential(gridSize);
-    tinkerGridPotential[0]  =  -1.8587681e+01;
-    tinkerGridPotential[1]  =  -3.7731481e+01;
-    tinkerGridPotential[2]  =  -1.4093518e+01;
-    tinkerGridPotential[3]  =  -8.6733284e-01;
-    tinkerGridPotential[4]  =   1.6378362e+01;
-    tinkerGridPotential[5]  =   1.7545816e+01;
-    tinkerGridPotential[6]  =   1.2115085e+01;
-    tinkerGridPotential[7]  =   1.5693689e+02;
-    tinkerGridPotential[8]  =   5.6517394e+01;
-    tinkerGridPotential[9]  =  -2.8489173e+01;
-    tinkerGridPotential[10] =  -1.1037348e+02;
-    tinkerGridPotential[11] =  -7.1050586e+01;
-    tinkerGridPotential[12] =  -5.9901289e+00;
-    tinkerGridPotential[13] =  -4.0533846e+00;
-    tinkerGridPotential[14] =   1.0506199e+01;
-    tinkerGridPotential[15] =  -8.6421838e+00;
-    tinkerGridPotential[16] =   8.3729402e+01;
-    tinkerGridPotential[17] =   4.4286652e+01;
-    tinkerGridPotential[18] =  -3.4746066e+01;
-    tinkerGridPotential[19] =  -1.0763056e+03;
-    tinkerGridPotential[20] =  -2.0034673e+01;
-    tinkerGridPotential[21] =  -1.2131063e+01;
-    tinkerGridPotential[22] =  -2.4662346e+01;
-    tinkerGridPotential[23] =   1.4630799e+00;
-    tinkerGridPotential[24] =   5.1623298e+00;
-    tinkerGridPotential[25] =   3.3114482e+01;
-    tinkerGridPotential[26] =   2.5828622e+01;
-
-    double tolerance = 1.0e-04;
-    for( unsigned int ii = 0; ii < gridSize; ii++ ){
-        double difference = fabs( outputGridPotential[ii] - tinkerGridPotential[ii] );
-        //fprintf( stderr, "%2d %15.7e %15.7e %15.7e %15.7e\n", ii, difference, difference/fabs(tinkerGridPotential[ii]), outputGridPotential[ii], tinkerGridPotential[ii] );
-        if( difference > tolerance ){
-            std::stringstream details;
-            details << testName << "Multipole moment " << ii << " does not agree w/ TINKER computed moments: OpenMM=" << outputGridPotential[ii];
-            details << " TINKER=" <<  tinkerGridPotential[ii]  << " difference=" << difference;
-            throwException(__FILE__, __LINE__, details.str());
-        }
-    }
-
 }
 
 // setup for box of 216 water molecules -- used to test PME
@@ -1314,7 +1208,9 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::Am
                                                         AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
                                                         double cutoff, int inputPmeGridDimension, std::string& testName, 
                                                         std::vector<Vec3>& forces, double& energy,
-                                                        std::vector< double >& outputMultipoleMoments, FILE* log ){
+                                                        std::vector< double >& outputMultipoleMoments,
+                                                        std::vector< Vec3 >& inputGrid,
+                                                        std::vector< double >& outputGridPotential, FILE* log ){
 
     // beginning of Multipole setup
 
@@ -2099,6 +1995,8 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::Am
     if( testName == "testSystemMultipoleMoments" ){
         Vec3 origin( 0.0, 0.0, 0.0 );
         amoebaMultipoleForce->getSystemMultipoleMoments( origin, context, outputMultipoleMoments );
+    } else if( testName == "testMultipoleGridPotential" ){
+        amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
     } else {
         State state               = context.getState(State::Forces | State::Energy);
         forces                    = state.getForces();
@@ -2120,10 +2018,12 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) {
     std::vector<Vec3> forces;
     double energy;
     std::vector< double > outputMultipoleMoments;
+    std::vector< Vec3 > inputGrid;
+    std::vector< double > outputGridPotential;
 
     setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, 
                                                 cutoff, inputPmeGridDimension, testName,
-                                                forces, energy, outputMultipoleMoments, log );
+                                                forces, energy, outputMultipoleMoments, inputGrid, outputGridPotential, log );
     std::vector<Vec3> expectedForces(numberOfParticles);
 
     double expectedEnergy     =  -1.3268930e+04;
@@ -2795,10 +2695,12 @@ static void testSystemMultipoleMoments( FILE* log ) {
     std::vector<Vec3> forces;
     double energy;
     std::vector<double> outputMultipoleMoments;
+    std::vector< Vec3 > inputGrid;
+    std::vector< double > outputGridPotential;
 
     setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, 
                                                 cutoff, inputPmeGridDimension, testName,
-                                                forces, energy, outputMultipoleMoments, log );
+                                                forces, energy, outputMultipoleMoments, inputGrid, outputGridPotential, log );
 
     std::vector<double> tinkerMoments(13);
 
@@ -2842,6 +2744,105 @@ static void testSystemMultipoleMoments( FILE* log ) {
 
 }
 
+// test computation of multipole potential on a grid 
+
+static void testMultipoleGridPotential( FILE* log ) {
+
+    std::string testName      = "testMultipoleGridPotential";
+    
+    int numberOfParticles     = 648;
+    int inputPmeGridDimension = 24;
+    double cutoff             = 0.70;
+
+    std::vector<Vec3> forces;
+    double energy;
+
+    // initialize grid
+
+    int gridSize  = 27;
+    std::vector<Vec3> inputGrid(gridSize);
+
+    inputGrid[0]  = Vec3( -1.0318535e+00,  -1.0224502e+00,  -1.0202836e+00);
+    inputGrid[1]  = Vec3( -1.0318535e+00,  -1.0224502e+00,  -3.4032700e-02);
+    inputGrid[2]  = Vec3( -1.0318535e+00,  -1.0224502e+00,   9.5221820e-01);
+    inputGrid[3]  = Vec3( -1.0318535e+00,  -3.0961200e-02,  -1.0202836e+00);
+    inputGrid[4]  = Vec3( -1.0318535e+00,  -3.0961200e-02,  -3.4032700e-02);
+    inputGrid[5]  = Vec3( -1.0318535e+00,  -3.0961200e-02,   9.5221820e-01);
+    inputGrid[6]  = Vec3( -1.0318535e+00,   9.6052780e-01,  -1.0202836e+00);
+    inputGrid[7]  = Vec3( -1.0318535e+00,   9.6052780e-01,  -3.4032700e-02);
+    inputGrid[8]  = Vec3( -1.0318535e+00,   9.6052780e-01,   9.5221820e-01);
+    inputGrid[9]  = Vec3( -3.3969000e-02,  -1.0224502e+00,  -1.0202836e+00);
+    inputGrid[10] = Vec3( -3.3969000e-02,  -1.0224502e+00,  -3.4032700e-02);
+    inputGrid[11] = Vec3( -3.3969000e-02,  -1.0224502e+00,   9.5221820e-01);
+    inputGrid[12] = Vec3( -3.3969000e-02,  -3.0961200e-02,  -1.0202836e+00);
+    inputGrid[13] = Vec3( -3.3969000e-02,  -3.0961200e-02,  -3.4032700e-02);
+    inputGrid[14] = Vec3( -3.3969000e-02,  -3.0961200e-02,   9.5221820e-01);
+    inputGrid[15] = Vec3( -3.3969000e-02,   9.6052780e-01,  -1.0202836e+00);
+    inputGrid[16] = Vec3( -3.3969000e-02,   9.6052780e-01,  -3.4032700e-02);
+    inputGrid[17] = Vec3( -3.3969000e-02,   9.6052780e-01,   9.5221820e-01);
+    inputGrid[18] = Vec3(  9.6391550e-01,  -1.0224502e+00,  -1.0202836e+00);
+    inputGrid[19] = Vec3(  9.6391550e-01,  -1.0224502e+00,  -3.4032700e-02);
+    inputGrid[20] = Vec3(  9.6391550e-01,  -1.0224502e+00,   9.5221820e-01);
+    inputGrid[21] = Vec3(  9.6391550e-01,  -3.0961200e-02,  -1.0202836e+00);
+    inputGrid[22] = Vec3(  9.6391550e-01,  -3.0961200e-02,  -3.4032700e-02);
+    inputGrid[23] = Vec3(  9.6391550e-01,  -3.0961200e-02,   9.5221820e-01);
+    inputGrid[24] = Vec3(  9.6391550e-01,   9.6052780e-01,  -1.0202836e+00);
+    inputGrid[25] = Vec3(  9.6391550e-01,   9.6052780e-01,  -3.4032700e-02);
+    inputGrid[26] = Vec3(  9.6391550e-01,   9.6052780e-01,   9.5221820e-01);
+
+    std::vector<double> outputGridPotential;
+    std::vector< double > outputMultipoleMoments; 
+
+    setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, 
+                                                cutoff, inputPmeGridDimension, testName, forces, energy,
+                                                outputMultipoleMoments, inputGrid, outputGridPotential, log );
+
+    // TINKER computed grid values
+
+    std::vector<double> tinkerGridPotential(gridSize);
+
+    tinkerGridPotential[0]  =   7.5696473e+01;
+    tinkerGridPotential[1]  =  -1.5382673e+01;
+    tinkerGridPotential[2]  =   6.9487135e+01;
+    tinkerGridPotential[3]  =   1.5459661e+03;
+    tinkerGridPotential[4]  =  -4.5138366e+02;
+    tinkerGridPotential[5]  =   4.6553857e+01;
+    tinkerGridPotential[6]  =   1.8486978e+01;
+    tinkerGridPotential[7]  =  -5.9235079e+01;
+    tinkerGridPotential[8]  =   1.0448125e+01;
+    tinkerGridPotential[9]  =   1.9529845e+00;
+    tinkerGridPotential[10] =   5.3587575e+00;
+    tinkerGridPotential[11] =  -3.3794408e+01;
+    tinkerGridPotential[12] =   1.5197977e+01;
+    tinkerGridPotential[13] =   4.3613381e+02;
+    tinkerGridPotential[14] =   1.6620315e+01;
+    tinkerGridPotential[15] =   5.3724507e+01;
+    tinkerGridPotential[16] =   6.9026934e+02;
+    tinkerGridPotential[17] =  -1.2708548e+01;
+    tinkerGridPotential[18] =  -9.8826500e+02;
+    tinkerGridPotential[19] =   3.2952143e+01;
+    tinkerGridPotential[20] =   9.1165100e+02;
+    tinkerGridPotential[21] =   4.0844598e+01;
+    tinkerGridPotential[22] =  -1.3509483e+01;
+    tinkerGridPotential[23] =   1.4499749e+00;
+    tinkerGridPotential[24] =   1.5011664e+02;
+    tinkerGridPotential[25] =  -2.8581208e+02;
+    tinkerGridPotential[26] =   1.3960144e+02;
+
+    double tolerance = 3.0e-04;
+    for( unsigned int ii = 0; ii < gridSize; ii++ ){
+        double difference = fabs( (outputGridPotential[ii] - tinkerGridPotential[ii] )/tinkerGridPotential[ii] );
+        //(void) fprintf( stderr, "Grid: %2d %15.7e %15.7e %15.7e\n", ii, difference, outputGridPotential[ii], tinkerGridPotential[ii] );
+        if( difference > tolerance ){
+            std::stringstream details;
+            details << testName << " potential for grid point " << ii << " does not agree w/ TINKER computed value: OpenMM=" << outputGridPotential[ii];
+            details << " TINKER=" <<  tinkerGridPotential[ii]  << " difference=" << difference;
+            throwException(__FILE__, __LINE__, details.str());
+        }
+    }
+
+}
+
 int main( int numberOfArguments, char* argv[] ) {
 
     try {