Continue PME

plugins/amoeba/openmmapi/include/AmoebaMultipoleForce.h

@@ -127,6 +127,34 @@ public:
      */		 
     void setAEwald(double aewald);		 
 		 
+   /**		   
+     * Get the B-spline order parameter		 
+     *		 
+     * @return the B-spline order parameter
+     */		 
+    int getPmeBSplineOrder( ) const;
+        
+   /**		   
+     * Set the B-spline order parameter		 
+     *		 
+     * @param the B-spline order parameter
+     */		 
+    void setPmeBSplineOrder(int inputBSplineOrder);
+     
+   /**		   
+     * Get the PME grid dimensions		 
+     *		 
+     * @return the PME grid dimensions
+     */		 
+    void getPmeGridDimensions( std::vector<int>& gridDimension ) const;
+         
+   /**		   
+     * Set the PME grid dimensions		 
+     *		 
+     * @param the PME grid dimensions
+     */		 
+    void setPmeGridDimensions( const std::vector<int>& gridDimension );
+    
     /**
      * Add multipole-related info for a particle 
      *
@@ -294,6 +322,8 @@ private:
     AmoebaNonbondedMethod nonbondedMethod;
     double cutoffDistance;
     double aewald;
+    int pmeBSplineOrder;
+    std::vector<int> pmeGridDimension;
     MutualInducedIterationMethod mutualInducedIterationMethod;
     int mutualInducedMaxIterations;
     double mutualInducedTargetEpsilon;

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -36,7 +36,7 @@
 
 using namespace OpenMM;
 
-AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), ewaldErrorTol(5e-4), mutualInducedIterationMethod(SOR), mutualInducedMaxIterations(60),
+AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), pmeBSplineOrder(5), cutoffDistance(1.0), ewaldErrorTol(5e-4), mutualInducedIterationMethod(SOR), mutualInducedMaxIterations(60),
                                                mutualInducedTargetEpsilon(1.0e-05), scalingDistanceCutoff(100.0), electricConstant(138.9354558456) {
 }
 
@@ -64,6 +64,36 @@ void AmoebaMultipoleForce::setAEwald(double inputAewald ) {
     aewald = inputAewald;		 
 }		 
  
+int AmoebaMultipoleForce::getPmeBSplineOrder( void ) const {		 
+    return pmeBSplineOrder;		 
+}		 
+		 
+void AmoebaMultipoleForce::setPmeBSplineOrder(int inputBSplineOrder) {			 
+    pmeBSplineOrder = inputBSplineOrder;		 
+}		 
+ 
+void AmoebaMultipoleForce::getPmeGridDimensions( std::vector<int>& gridDimension ) const {		 
+    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;
+}		 
+		 
+void AmoebaMultipoleForce::setPmeGridDimensions( const std::vector<int>& gridDimension ) {
+    pmeGridDimension.resize(3);
+    pmeGridDimension[0] = gridDimension[0];
+    pmeGridDimension[1] = gridDimension[1];
+    pmeGridDimension[2] = gridDimension[2];
+    return;
+}		 
+ 
 AmoebaMultipoleForce::MutualInducedIterationMethod AmoebaMultipoleForce::getMutualInducedIterationMethod( void ) const {
     return mutualInducedIterationMethod;
 }

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

@@ -597,7 +597,16 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
         NonbondedForce nb;
         nb.setEwaldErrorTolerance(force.getEwaldErrorTolerance());
         nb.setCutoffDistance(force.getCutoffDistance());
-        NonbondedForceImpl::calcPMEParameters(system, nb, alpha, xsize, ysize, zsize);
+        std::vector<int> pmeGridDimension;
+        force.getPmeGridDimensions( pmeGridDimension );
+        if( 1 || pmeGridDimension[0] == 0 ){
+            NonbondedForceImpl::calcPMEParameters(system, nb, alpha, xsize, ysize, zsize);
+        } else {
+            alpha = force.getAEwald();
+            xsize = pmeGridDimension[0];
+            ysize = pmeGridDimension[1];
+            zsize = pmeGridDimension[2];
+        }
         gpuSetAmoebaPMEParameters(data.getAmoebaGpu(), (float) alpha, xsize, ysize, zsize);
     }
 

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

@@ -351,6 +351,7 @@ 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, "     alpha Ewald                        %15.7e\n",  gpu->sim.alphaEwald );
+    (void) fprintf( log, "     PME grid dimensions                %6d %6d %6d\n",  gpu->sim.pmeGridSize.x, gpu->sim.pmeGridSize.y, gpu->sim.pmeGridSize.z);
     (void) fprintf( log, "     cutoffDistance2                    %15.7e\n",  amoebaGpu->amoebaSim.cutoffDistance2 );
     (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);

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

@@ -6,7 +6,7 @@
 #include "amoebaCudaKernels.h"
 #include "kCalculateAmoebaCudaUtilities.h"
 
-//#define AMOEBA_DEBUG
+#define AMOEBA_DEBUG
 
 static __constant__ cudaGmxSimulation cSim;
 static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;
@@ -39,7 +39,7 @@ __launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
 #else
 __launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
 #endif
-void kReduceDirectSelfFields_kernel( unsigned int fieldComponents, unsigned int outputBuffers,  float* fieldIn, float* fieldOut )
+static void kReducePmeEFieldPolar_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* EFieldReciprocal,  float* fieldIn, float* fieldOut )
 {
     unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
 
@@ -52,7 +52,7 @@ void kReduceDirectSelfFields_kernel( unsigned int fieldComponents, unsigned int
 
         // self-term included here
 
-        float totalField = term*cAmoebaSim.pLabFrameDipole[pos];
+        float totalField = EFieldReciprocal[pos] + term*cAmoebaSim.pLabFrameDipole[pos];
 
         float* pFt       = fieldIn + pos;
         unsigned int i   = outputBuffers;
@@ -80,19 +80,74 @@ void kReduceDirectSelfFields_kernel( unsigned int fieldComponents, unsigned int
     }   
 }
 
+
+__global__ 
+#if (__CUDA_ARCH__ >= 200)
+__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
+#elif (__CUDA_ARCH__ >= 130)
+__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
+#else
+__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
+#endif
+static void kReducePmeEField_kernel( unsigned int fieldComponents, unsigned int outputBuffers,  float* fieldIn, float* fieldOut )
+{
+    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
+
+    // Reduce field
+
+    const float term = (4.0f/3.0f)*(cSim.alphaEwald*cSim.alphaEwald*cSim.alphaEwald)/cAmoebaSim.sqrtPi;
+    //const float term = 0.0f;
+    while (pos < fieldComponents)
+    {   
+
+        // self-term included here
+
+        float totalField = term*cAmoebaSim.pLabFrameDipole[pos];
+
+        float* pFt       = fieldIn + pos;
+        unsigned int i   = outputBuffers;
+        while (i >= 4)
+        {   
+            totalField += pFt[0] + pFt[fieldComponents] + pFt[2*fieldComponents] + pFt[3*fieldComponents];
+            pFt        += fieldComponents*4;
+            i          -= 4;
+        }   
+
+        if (i >= 2)
+        {   
+            totalField += pFt[0] + pFt[fieldComponents];
+            pFt        += fieldComponents*2;
+            i          -= 2;
+        }   
+
+        if (i > 0)
+        {   
+            totalField += pFt[0];
+        }   
+
+        fieldOut[pos]  += totalField;
+        pos            += gridDim.x * blockDim.x;
+    }   
+}
+
 // reduce psWorkArray_3_1 -> EField
 // reduce psWorkArray_3_2 -> EFieldPolar
 
 static void kReducePmeDirectE_Fields(amoebaGpuContext amoebaGpu )
 {
-    kReduceDirectSelfFields_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
-                               amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->outputBuffers,
-                               amoebaGpu->psWorkArray_3_1->_pDevStream[0], amoebaGpu->psE_Field->_pDevStream[0] );
+
+    // E_FieldPolar = E_Field (reciprocal) + E_FieldPolar (direct) + self
+
+    kReducePmeEFieldPolar_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
+                                   amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->outputBuffers,
+                                   amoebaGpu->psE_Field->_pDevStream[0], amoebaGpu->psWorkArray_3_2->_pDevStream[0], amoebaGpu->psE_FieldPolar->_pDevStream[0] );
     LAUNCHERROR("kReducePmeE_Fields1");
 
-    kReduceDirectSelfFields_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
-                               amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->outputBuffers,
-                               amoebaGpu->psWorkArray_3_2->_pDevStream[0], amoebaGpu->psE_FieldPolar->_pDevStream[0] );
+    // E_Field = E_Field (reciprocal) + E_Field (direct) + self
+
+    kReducePmeEField_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
+                              amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->outputBuffers,
+                              amoebaGpu->psWorkArray_3_1->_pDevStream[0], amoebaGpu->psE_Field->_pDevStream[0] );
     LAUNCHERROR("kReducePmeE_Fields2");
 }
 
@@ -318,6 +373,20 @@ __device__ void calculateFixedFieldRealSpacePairIxn_kernel( FixedFieldParticle&
 #define METHOD_NAME(a, b) a##N2ByWarp##b
 #include "kCalculateAmoebaCudaPmeFixedEField.h"
 
+/**---------------------------------------------------------------------------------------
+
+   Report whether a number is a nan or infinity
+
+   @param number               number to test
+   @return 1 if number is  nan or infinity; else return 0
+
+   --------------------------------------------------------------------------------------- */
+
+static int isNanOrInfinity( double number ){
+    return (number != number || number == std::numeric_limits<double>::infinity() || number == -std::numeric_limits<double>::infinity()) ? 1 : 0; 
+}
+
+
 /**---------------------------------------------------------------------------------------
 
    Compute fixed electric field using PME
@@ -335,6 +404,7 @@ static void cudaComputeAmoebaPmeDirectFixedEField( amoebaGpuContext amoebaGpu )
     gpuContext gpu    = amoebaGpu->gpuContext;
 
 #ifdef AMOEBA_DEBUG
+
     static const char* methodName = "computeCudaAmoebaPmeFixedEField";
     if( amoebaGpu->log ){
         (void) fprintf( amoebaGpu->log, "\n%s\n", methodName ); (void) fflush( amoebaGpu->log );
@@ -347,21 +417,36 @@ static void cudaComputeAmoebaPmeDirectFixedEField( amoebaGpuContext amoebaGpu )
     memset( debugArray->_pSysStream[0],      0, sizeof( float )*4*paddedNumberOfAtoms*paddedNumberOfAtoms);
     debugArray->Upload();
 
-    (*gpu->psInteractionCount)[0]              = gpu->sim.workUnits;
-    gpu->psInteractionCount->Upload();
-
     // print intermediate results for the targetAtom 
 
     unsigned int targetAtom  = 0;
-#endif
 
-//#define SET_ALPHA_EWALD
-//#ifdef SET_ALPHA_EWALD
-//    amoebaGpu->gpuContext->sim.alphaEwald = 5.4459052e+00f;
-//    (void) fprintf( amoebaGpu->log, "computeCudaAmoebaPmeFixedEField: forceing alphaEwald=%15.7e\n", amoebaGpu->gpuContext->sim.alphaEwald );
-//    SetCalculateAmoebaCudaPmeFixedEFieldSim(amoebaGpu);
-//#endif
+    int maxPrint             = 3002;
+    amoebaGpu->psE_Field->Download();
+    (void) fprintf( amoebaGpu->log, "Recip EFields In\n" );
+    for( int ii = 0; ii < gpu->natoms; ii++ ){
+        (void) fprintf( amoebaGpu->log, "%5d ", ii); 
 
+        int indexOffset     = ii*3;
+
+        // E_Field
+
+        int isNan  = isNanOrInfinity( amoebaGpu->psE_Field->_pSysStream[0][indexOffset] );
+            isNan += isNanOrInfinity( amoebaGpu->psE_Field->_pSysStream[0][indexOffset+1] );
+            isNan += isNanOrInfinity( amoebaGpu->psE_Field->_pSysStream[0][indexOffset+2] );
+
+        (void) fprintf( amoebaGpu->log,"E[%16.9e %16.9e %16.9e] %s\n",
+                        amoebaGpu->psE_Field->_pSysStream[0][indexOffset],
+                        amoebaGpu->psE_Field->_pSysStream[0][indexOffset+1],
+                        amoebaGpu->psE_Field->_pSysStream[0][indexOffset+2], (isNan ? "XXX" :"") );
+   
+        if( ii == maxPrint && (gpu->natoms - maxPrint) > ii ){
+            ii = gpu->natoms - maxPrint;
+        }
+    }
+    (void) fflush( amoebaGpu->log );
+    (void) fprintf( amoebaGpu->log, "Recip EFields End\n" );
+#endif
 
     kClearFields_3( amoebaGpu, 2 );
 
@@ -413,12 +498,10 @@ static void cudaComputeAmoebaPmeDirectFixedEField( amoebaGpuContext amoebaGpu )
                         sizeof(FixedFieldParticle), sizeof(FixedFieldParticle)*amoebaGpu->nonbondThreadsPerBlock, amoebaGpu->energyOutputBuffers, 
                         (*gpu->psInteractionCount)[0], gpu->sim.workUnits, gpu->bOutputBufferPerWarp );
         (void) fflush( amoebaGpu->log );
+/*
+        (void) fprintf( amoebaGpu->log, "Out WorkArray_3_[1,2]  paddedNumberOfAtoms=%d\n",  amoebaGpu->paddedNumberOfAtoms, amoebaGpu->outputBuffers );
         amoebaGpu->psWorkArray_3_1->Download();
         amoebaGpu->psWorkArray_3_2->Download();
-        amoebaGpu->psE_Field->Download();
-        amoebaGpu->psE_FieldPolar->Download();
-        (void) fprintf( amoebaGpu->log, "Out WorkArray_3_[1,2]  paddedNumberOfAtoms=%d\n",  amoebaGpu->paddedNumberOfAtoms, amoebaGpu->outputBuffers );
-        int maxPrint        = 32;
         for( int ii = 0; ii < amoebaGpu->paddedNumberOfAtoms; ii++ ){
            (void) fprintf( amoebaGpu->log, "%5d ", ii); 
 
@@ -444,8 +527,9 @@ static void cudaComputeAmoebaPmeDirectFixedEField( amoebaGpuContext amoebaGpu )
            }
         }
         (void) fflush( amoebaGpu->log );
-        (void) fprintf( amoebaGpu->log, "EFields End\n" );
-
+*/
+        amoebaGpu->psE_Field->Download();
+        amoebaGpu->psE_FieldPolar->Download();
         for( int ii = 0; ii < gpu->natoms; ii++ ){
            (void) fprintf( amoebaGpu->log, "%5d ", ii); 
 
@@ -520,6 +604,14 @@ static void cudaComputeAmoebaPmeDirectFixedEField( amoebaGpuContext amoebaGpu )
 
 void cudaComputeAmoebaPmeFixedEField( amoebaGpuContext amoebaGpu )
 {
-    cudaComputeAmoebaPmeDirectFixedEField( amoebaGpu );
+
+// zero field?
+
+unsigned int offset = 3*amoebaGpu->paddedNumberOfAtoms*sizeof( float );
+memset( amoebaGpu->psE_Field->_pSysStream[0], 0, offset );
+amoebaGpu->psE_Field->Upload();
+
+
     kCalculateAmoebaPMEFixedMultipoleField( amoebaGpu );
+    cudaComputeAmoebaPmeDirectFixedEField( amoebaGpu );
 }

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

@@ -1985,21 +1985,33 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
 
     // load in parameters
  
-    int numberOfMultipoles            = atoi( tokens[1].c_str() );
+    unsigned int tokenIndex           = 1;
+    int numberOfMultipoles            = atoi( tokens[tokenIndex++].c_str() );
     int usePme                        = 0;
+    int bsOrder                       = 0;
+    std::vector<int> grid(3,0);
     double aewald                     = 0.0;
     double cutoffDistance             = 0.0;
-    double box[3]                     = { 10.0, 10.0, 10.0 };
+    std::vector<double> box(3,10.0);
 
     // usePme, aewald, cutoffDistance added w/ Version 1
 
     if( version > 0 ){
-        usePme                        = atoi( tokens[2].c_str() );
-        aewald                        = atof( tokens[3].c_str() );
-        cutoffDistance                = atof( tokens[4].c_str() );
-        box[0]                        = atof( tokens[5].c_str() );
-        box[1]                        = atof( tokens[6].c_str() );
-        box[2]                        = atof( tokens[7].c_str() );
+        usePme                        = atoi( tokens[tokenIndex++].c_str() );
+        aewald                        = atof( tokens[tokenIndex++].c_str() );
+        cutoffDistance                = atof( tokens[tokenIndex++].c_str() );
+
+        box[0]                        = atof( tokens[tokenIndex++].c_str() );
+        box[1]                        = atof( tokens[tokenIndex++].c_str() );
+        box[2]                        = atof( tokens[tokenIndex++].c_str() );
+
+        //double electric               = atof( tokens[tokenIndex++].c_str() );
+                                        tokenIndex++;
+
+        bsOrder                       = atoi( tokens[tokenIndex++].c_str() );
+        grid[0]                       = atoi( tokens[tokenIndex++].c_str() );
+        grid[1]                       = atoi( tokens[tokenIndex++].c_str() );
+        grid[2]                       = atoi( tokens[tokenIndex++].c_str() );
     }
   
     if( usePme ){
@@ -2009,6 +2021,8 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
     }
     multipoleForce->setCutoffDistance( cutoffDistance );
     multipoleForce->setAEwald( aewald );
+    multipoleForce->setPmeBSplineOrder( bsOrder );
+    multipoleForce->setPmeGridDimensions( grid );
     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 ){
         (void) fprintf( log, "%s number of MultipoleParameter terms=%d usePme=%d aewald=%15.7e cutoffDistance=%12.4f\n",
@@ -3605,15 +3619,19 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS
                 readVec3( filePtr, tokens, forces[AMOEBA_GK_FORCE], &lineCount, field, log );
             } else if( field == "AmoebaGkAndCavityForce" ){
                 readVec3( filePtr, tokens, forces[AMOEBA_GK_CAVITY_FORCE], &lineCount, field, log );
-            } else if( field == "AmoebaGk_A_ForceAndTorque"            || 
-                       field == "AmoebaGk_A_Force"                     ||
-                       field == "AmoebaSurfaceParameters"              ||
-                       field == "AmoebaGk_A_DrB"                       ||
-                       field == "AmoebaDBorn"                          ||
-                       field == "AmoebaBorn1Force"                     ||
-                       field == "AmoebaBornForce"                      ||
-                       field == "AmoebaGkEdiffForceAndTorque"          ||
-                       field == "AmoebaGkEdiffForce"                   ){
+            } else if( field == "AmoebaGk_A_ForceAndTorque"                     || 
+                       field == "AmoebaGk_A_Force"                              ||
+                       field == "AmoebaSurfaceParameters"                       ||
+                       field == "AmoebaGk_A_DrB"                                ||
+                       field == "AmoebaDBorn"                                   ||
+                       field == "AmoebaBorn1Force"                              ||
+                       field == "AmoebaBornForce"                               ||
+                       field == "AmoebaGkEdiffForceAndTorque"                   ||
+                       field == "AmoebaGkEdiffForce"                            ||
+                       field == "PmeDirectForceAndTorqueOutLoop"                ||
+                       field == "PmeDirectForceIncludingMappedTorqueOutLoop"    ||
+                       field == "PmeDirectForceTorqueInLoop"                    
+                    ){
                 std::vector< std::vector<double> > vectorOfDoubleVectors;
                 readVectorOfDoubleVectors( filePtr, tokens, vectorOfDoubleVectors, &lineCount, field, log );
                 supplementary[field] = vectorOfDoubleVectors;