Decouple multipole & vdw cutoff

plugins/amoeba/openmmapi/include/AmoebaVdwForce.h

@@ -157,6 +157,20 @@ public:
      */
     double getCutoff( void ) const;
 
+    /**
+     * Set flag for using neighbor list for vdw ixn
+     * 
+     * @param neighboristFlag neighbor list flag
+     */
+    void setUseNeighborList( int neighborListFlag );
+
+    /**
+     * Get neighbor list flag for vdw ixn
+     * 
+     * @return neighbor list flag
+     */
+    int getUseNeighborList( void ) const;
+
     /**
      * Set flag for employing periodic boundary conditions
      * 
@@ -177,7 +191,7 @@ private:
 
     class VdwInfo;
     int usePBC;
-    double cutoff;
+    int useNeighborList; double cutoff;
     std::string sigmaCombiningRule;
     std::string epsilonCombiningRule;
     std::vector< std::vector<int> > exclusions;

plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp

@@ -38,7 +38,7 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-AmoebaVdwForce::AmoebaVdwForce() : usePBC(0), cutoff(1.0e+10) {
+AmoebaVdwForce::AmoebaVdwForce() : usePBC(0), cutoff(1.0e+10), useNeighborList(0) {
 }
 
 int AmoebaVdwForce::addParticle(int ivIndex, int classIndex, double sigma, double epsilon, double reductionFactor ) {
@@ -112,6 +112,14 @@ double AmoebaVdwForce::getCutoff( void ) const {
     return cutoff;
 }
 
+void AmoebaVdwForce::setUseNeighborList( int useNeighborListFlag ){
+    useNeighborList = useNeighborListFlag;
+}
+
+int AmoebaVdwForce::getUseNeighborList( void ) const {
+    return useNeighborList;
+}
+
 void AmoebaVdwForce::setPBC( int pbcFlag ){
     usePBC = pbcFlag;
 }

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

@@ -42,7 +42,8 @@ AmoebaCudaData::AmoebaCudaData( CudaPlatform::PlatformData& data ) : cudaPlatfor
     log                           = NULL;
     contextImpl                   = NULL;
     gpuInitialized                = false;
-    applyCutoff                   = 0;
+    applyMultipoleCutoff          = 0;
+    useVdwNeighborList            = 0;
     multipoleForceCount           = 0;
 }   
 
@@ -130,12 +131,20 @@ int AmoebaCudaData::getMultipoleForceCount( void ) const {
     return multipoleForceCount;
 }
 
-void AmoebaCudaData::setApplyCutoff( int inputApplyCutoff ) {
-    applyCutoff = inputApplyCutoff;
+void AmoebaCudaData::setApplyMultipoleCutoff( int inputApplyMultipoleCutoff ) {
+    applyMultipoleCutoff = inputApplyMultipoleCutoff;
 }
 
-int AmoebaCudaData::getApplyCutoff( void ) const {
-    return applyCutoff;
+int AmoebaCudaData::getApplyMultipoleCutoff( void ) const {
+    return applyMultipoleCutoff;
+}
+
+void AmoebaCudaData::setUseVdwNeighborList( int inputUseVdwNeighborList ) {
+    useVdwNeighborList = inputUseVdwNeighborList;
+}
+
+int AmoebaCudaData::getUseVdwNeighborList( void ) const {
+    return useVdwNeighborList;
 }
 
 }

plugins/amoeba/platforms/cuda/src/AmoebaCudaData.h

@@ -154,18 +154,32 @@ public:
     void incrementMultipoleForceCount( void ); 
 
     /**
-     * Get multipole force count 
+     * Get multipole cutoff
      * 
-     * @return multipole force count
+     * @return multipole cutoff
      */
-    int getApplyCutoff( ) const; 
+    int getApplyMultipoleCutoff( ) const; 
 
     /**
-     * Get multipole force count 
+     * Set multipole cutoff
      * 
-     * @return multipole force count
+     * @return multipole cutoff
+     */
+    void setApplyMultipoleCutoff( int applyMultipoleCutoff ); 
+
+    /**
+     * Get vdw cutoff
+     * 
+     * @return vdw cutoff
+     */
+    int getUseVdwNeighborList( ) const; 
+
+    /**
+     * Set vdw cutoff
+     * 
+     * @return vdw cutoff
      */
-    void setApplyCutoff( int applyCutoff ); 
+    void setUseVdwNeighborList( int useVdwNeighborList ); 
 
     CudaPlatform::PlatformData& cudaPlatformData;
 
@@ -174,7 +188,8 @@ private:
     amoebaGpuContext amoebaGpu;
     bool hasAmoebaBonds, hasAmoebaGeneralizedKirkwood, hasAmoebaMultipole;
     int multipoleForceCount;
-    int applyCutoff;
+    int applyMultipoleCutoff;
+    int useVdwNeighborList;
     KernelImpl* localForceKernel;
     unsigned int kernelCount;
     void* contextImpl;

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

@@ -983,7 +983,7 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
             zsize = pmeGridDimension[2];
         }
         gpuSetAmoebaPMEParameters(data.getAmoebaGpu(), (float) alpha, xsize, ysize, zsize);
-        data.setApplyCutoff( 1 );
+        data.setApplyMultipoleCutoff( 1 );
         data.cudaPlatformData.nonbondedMethod = PARTICLE_MESH_EWALD;
         amoebaGpuContext amoebaGpu            = data.getAmoebaGpu();
         gpuContext gpu                        = amoebaGpu->gpuContext;
@@ -1063,7 +1063,7 @@ static void computeAmoebaVdwForce( AmoebaCudaData& data ) {
 
     // Vdw14_7F
 
-    kCalculateAmoebaVdw14_7Forces(gpu, data.getApplyCutoff());
+    kCalculateAmoebaVdw14_7Forces(gpu, data.getUseVdwNeighborList());
 }
 
 /* -------------------------------------------------------------------------- *
@@ -1129,6 +1129,7 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
                                force.getSigmaCombiningRule(), force.getEpsilonCombiningRule(),
                                allExclusions, force.getPBC(), static_cast<float>(force.getCutoff()) );
     data.getAmoebaGpu()->gpuContext->forces.push_back(new ForceInfo(force));
+    data.setUseVdwNeighborList( force.getUseNeighborList() );
 }
 
 double CudaCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

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

@@ -139,7 +139,6 @@ struct cudaAmoebaGmxSimulation {
     float scalingDistanceCutoff;                    // scaling cutoff
     float2*         pDampingFactorAndThole;         // Thole & damping factors
 
-    float* pRotationMatrix; 
     int4*  pMultipoleParticlesIdsAndAxisType; 
     int*   pMultipoleAxisOffset; 
     float* pMolecularDipole; 

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

@@ -36,6 +36,13 @@ void GetCalculateAmoebaCudaUtilitiesSim(amoebaGpuContext amoebaGpu)
 #undef METHOD_NAME
 #undef USE_PERIODIC
 
+#undef METHOD_NAME
+#define USE_PERIODIC
+#define METHOD_NAME(a, b) a##Periodic##b
+#include "kFindInteractingBlocksVdw.h"
+#undef METHOD_NAME
+#undef USE_PERIODIC
+
 __global__ 
 #if (__CUDA_ARCH__ >= 200)
 __launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)

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

@@ -12,4 +12,8 @@ extern __global__ void kFindBlocksWithInteractionsPeriodic_kernel();
 extern __global__ void kFindInteractionsWithinBlocksPeriodic_kernel(unsigned int*);
 
 
+extern __global__ void kFindBlocksWithInteractionsVdwPeriodic_kernel();
+extern __global__ void kFindInteractionsWithinBlocksVdwPeriodic_kernel(unsigned int*);
+
+
 #endif

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

@@ -36,8 +36,8 @@ void GetCalculateAmoebaCudaVdw14_7Sim(amoebaGpuContext amoebaGpu)
     RTERROR(status, "GetCalculateAmoebaCudaVdw14_7Sim: cudaMemcpyFromSymbol: SetSim copy from cAmoebaSim failed");
 }
 
-//#define AMOEBA_DEBUG
-#undef AMOEBA_DEBUG
+#define AMOEBA_DEBUG
+//#undef AMOEBA_DEBUG
 
 __device__ void zeroVdw14_7SharedForce( struct Vdw14_7Particle* sA ) 
 {
@@ -533,15 +533,15 @@ void kCalculateAmoebaVdw14_7Forces( amoebaGpuContext amoebaGpu, int applyCutoff
     if( applyCutoff ){
 
         kFindBlockBoundsPeriodic_kernel<<<(gpu->psGridBoundingBox->_length+63)/64, 64>>>();
-        LAUNCHERROR("kFindBlockBoundsPeriodic");
-        kFindBlocksWithInteractionsPeriodic_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
-        LAUNCHERROR("kFindBlocksWithInteractionsPeriodic");
+        LAUNCHERROR("kFindBlockBoundsVdwPeriodic");
+        kFindBlocksWithInteractionsVdwPeriodic_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
+        LAUNCHERROR("kFindBlocksWithInteractionsVdwPeriodic");
         compactStream(gpu->compactPlan, gpu->sim.pInteractingWorkUnit, amoebaGpu->amoebaSim.pVdwWorkUnit, gpu->sim.pInteractionFlag, gpu->sim.workUnits, gpu->sim.pInteractionCount);
-        kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
+        kFindInteractionsWithinBlocksVdwPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                     sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
-        LAUNCHERROR("kFindInteractionsWithinBlocksPeriodic");
+        LAUNCHERROR("kFindInteractionsWithinBlocksVdwPeriodic");
 
-if( 0 ){  
+if( 1 ){  
     gpu->psInteractionCount->Download();
     gpu->psInteractingWorkUnit->Download();
     gpu->psInteractionFlag->Download();

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

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2009 Stanford University and the Authors.           *
+ * Authors: Scott Le Grand, Peter Eastman                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This file contains the kernels for identifying interacting blocks.  It is included
+ * several times in kCalculateCDLJForces.cu with different #defines to generate
+ * different versions of the kernels.
+ */
+
+/**
+ * Find a bounding box for the atoms in each block.
+ */
+/*
+__global__ void METHOD_NAME(kFindBlockBounds, _kernel)()
+{
+    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
+    unsigned int base = pos << GRIDBITS;
+    if (base < cSim.atoms)
+    {
+        float4 apos = cSim.pPosq[base];
+#ifdef USE_PERIODIC
+        apos.x -= floor(apos.x*cSim.invPeriodicBoxSizeX)*cSim.periodicBoxSizeX;
+        apos.y -= floor(apos.y*cSim.invPeriodicBoxSizeY)*cSim.periodicBoxSizeY;
+        apos.z -= floor(apos.z*cSim.invPeriodicBoxSizeZ)*cSim.periodicBoxSizeZ;
+        float4 firstPoint = apos;
+#endif
+        float minx = apos.x;
+        float maxx = apos.x;
+        float miny = apos.y;
+        float maxy = apos.y;
+        float minz = apos.z;
+        float maxz = apos.z;
+        for (unsigned int i = 1; i < GRID; i++)
+        {
+            apos = cSim.pPosq[base+i];
+#ifdef USE_PERIODIC
+            apos.x -= floor((apos.x-firstPoint.x)*cSim.invPeriodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
+            apos.y -= floor((apos.y-firstPoint.y)*cSim.invPeriodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
+            apos.z -= floor((apos.z-firstPoint.z)*cSim.invPeriodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
+#endif
+            minx = min(minx, apos.x);
+            maxx = max(maxx, apos.x);
+            miny = min(miny, apos.y);
+            maxy = max(maxy, apos.y);
+            minz = min(minz, apos.z);
+            maxz = max(maxz, apos.z);
+        }
+        cSim.pGridBoundingBox[pos] = make_float4(0.5f*(maxx-minx), 0.5f*(maxy-miny), 0.5f*(maxz-minz), 0);
+        cSim.pGridCenter[pos] = make_float4(0.5f*(maxx+minx), 0.5f*(maxy+miny), 0.5f*(maxz+minz), 0);
+    }
+}
+*/
+
+/**
+ * Compare the bounding boxes for each pair of blocks.  If they are sufficiently far apart,
+ * mark them as non-interacting.
+ */
+__global__ void METHOD_NAME(kFindBlocksWithInteractionsVdw, _kernel)()
+{
+    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
+    while (pos < cSim.workUnits)
+    {
+        // Extract cell coordinates from appropriate work unit
+
+        unsigned int x = cSim.pWorkUnit[pos];
+        unsigned int y = ((x >> 2) & 0x7fff);
+        x = (x >> 17);
+
+        // Find the distance between the bounding boxes of the two cells.
+
+        float4 centera = cSim.pGridCenter[x];
+        float4 centerb = cSim.pGridCenter[y];
+        float dx = centera.x-centerb.x;
+        float dy = centera.y-centerb.y;
+        float dz = centera.z-centerb.z;
+#ifdef USE_PERIODIC
+        dx -= floor(dx*cSim.invPeriodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
+        dy -= floor(dy*cSim.invPeriodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
+        dz -= floor(dz*cSim.invPeriodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
+#endif
+        float4 boxSizea = cSim.pGridBoundingBox[x];
+        float4 boxSizeb = cSim.pGridBoundingBox[y];
+        dx = max(0.0f, abs(dx)-boxSizea.x-boxSizeb.x);
+        dy = max(0.0f, abs(dy)-boxSizea.y-boxSizeb.y);
+        dz = max(0.0f, abs(dz)-boxSizea.z-boxSizeb.z);
+        cSim.pInteractionFlag[pos] = (dx*dx+dy*dy+dz*dz > cAmoebaSim.vdwCutoff2 ? 0 : 1);
+        pos += gridDim.x*blockDim.x;
+    }
+}
+
+/**
+ * Compare each atom in one block to the bounding box of another block, and set
+ * flags for which ones are interacting.
+ */
+__global__ void METHOD_NAME(kFindInteractionsWithinBlocksVdw, _kernel)(unsigned int* workUnit)
+{
+    extern __shared__ unsigned int flags[];
+    unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
+    unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
+    unsigned int numWorkUnits = cSim.pInteractionCount[0];
+    unsigned int pos = warp*numWorkUnits/totalWarps;
+    unsigned int end = (warp+1)*numWorkUnits/totalWarps;
+    unsigned int index = threadIdx.x & (GRID - 1);
+
+    unsigned int lasty = 0xFFFFFFFF;
+    float4 apos;
+    while (pos < end)
+    {
+        // Extract cell coordinates from appropriate work unit
+        unsigned int x = workUnit[pos];
+        unsigned int y = ((x >> 2) & 0x7fff);
+        bool bExclusionFlag = (x & 0x1);
+        x = (x >> 17);
+        if (x == y || bExclusionFlag)
+        {
+            // Assume this block will be dense.
+
+            if (index == 0)
+                cSim.pInteractionFlag[pos] = 0xFFFFFFFF;
+        }
+        else
+        {
+            // Load the bounding box for x and the atom positions for y.
+
+            float4 center = cSim.pGridCenter[x];
+            float4 boxSize = cSim.pGridBoundingBox[x];
+            if (y != lasty)
+            {
+                apos = cSim.pPosq[(y<<GRIDBITS)+index];
+            }
+
+            // Find the distance of the atom from the bounding box.
+
+            float dx = apos.x-center.x;
+            float dy = apos.y-center.y;
+            float dz = apos.z-center.z;
+#ifdef USE_PERIODIC
+            dx -= floor(dx*cSim.invPeriodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
+            dy -= floor(dy*cSim.invPeriodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
+            dz -= floor(dz*cSim.invPeriodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
+#endif
+            dx = max(0.0f, abs(dx)-boxSize.x);
+            dy = max(0.0f, abs(dy)-boxSize.y);
+            dz = max(0.0f, abs(dz)-boxSize.z);
+            flags[threadIdx.x] = (dx*dx+dy*dy+dz*dz > cAmoebaSim.vdwCutoff2 ? 0 : 1 << index);
+
+            // Sum the flags.
+
+            if (index % 2 == 0)
+                flags[threadIdx.x] += flags[threadIdx.x+1];
+            if (index % 4 == 0)
+                flags[threadIdx.x] += flags[threadIdx.x+2];
+            if (index % 8 == 0)
+                flags[threadIdx.x] += flags[threadIdx.x+4];
+            if (index % 16 == 0)
+                flags[threadIdx.x] += flags[threadIdx.x+8];
+            if (index == 0)
+            {
+                unsigned int allFlags = flags[threadIdx.x] + flags[threadIdx.x+16];
+
+                // Count how many flags are set, and based on that decide whether to compute all interactions
+                // or only a fraction of them.
+
+                unsigned int bits = (allFlags&0x55555555) + ((allFlags>>1)&0x55555555);
+                bits = (bits&0x33333333) + ((bits>>2)&0x33333333);
+                bits = (bits&0x0F0F0F0F) + ((bits>>4)&0x0F0F0F0F);
+                bits = (bits&0x00FF00FF) + ((bits>>8)&0x00FF00FF);
+                bits = (bits&0x0000FFFF) + ((bits>>16)&0x0000FFFF);
+                cSim.pInteractionFlag[pos] = (bits > 12 ? 0xFFFFFFFF : allFlags);
+            }
+            lasty = y;
+        }
+        pos++;
+    }
+}