Fix for particles w/o multipole axis (e.g. ions)

-This line, and th -se below, will be ignored--

M    amoebaCudaGpu.cpp
M    kCalculateAmoebaCudaMapTorques.cu
M    kCalculateAmoebaCudaRotateFrame.cu

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

@@ -1646,8 +1646,8 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
         amoebaGpu->amoebaSim.dielec               = 1.0f;
     }
 
-    static const int maxAxisType = 5;
-    int axisTypeCount[maxAxisType+1] = { 0, 0, 0, 0, 0, 0 };
+    static const int maxAxisType = 6;
+    int axisTypeCount[maxAxisType+1] = { 0, 0, 0, 0, 0, 0, 0 };
     for( int ii = 0; ii < static_cast<int>(charges.size()); ii++ ){
 
         // axis type & multipole particles ids
@@ -1666,19 +1666,23 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
         // and need test system
 
         int axisParticleIndex                                                     = multipoleParticleZ[ii];
-        if( maxIndices[axisParticleIndex] < ii ){
-            maxIndices[axisParticleIndex] = ii;
-        }
-        if( amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] > ii ){
-            amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] = ii;
+        if( axisParticleIndex > -1 ){
+            if( maxIndices[axisParticleIndex] < ii ){
+                maxIndices[axisParticleIndex] = ii;
+            }
+            if( amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] > ii ){
+                amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] = ii;
+            }
         }
 
         axisParticleIndex                                                         = multipoleParticleX[ii];
-        if( maxIndices[axisParticleIndex] < ii ){
-            maxIndices[axisParticleIndex] = ii;
-        }
-        if( amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] > ii ){
-            amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] = ii;
+        if( axisParticleIndex > -1 ){
+            if( maxIndices[axisParticleIndex] < ii ){
+                maxIndices[axisParticleIndex] = ii;
+            }
+            if( amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] > ii ){
+                amoebaGpu->psMultipoleAxisOffset->_pSysData[axisParticleIndex] = ii;
+            }
         }
 
         axisParticleIndex                                                         = multipoleParticleY[ii];
@@ -1692,7 +1696,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
         }
 
         if( 0 && amoebaGpu->log )
-            fprintf( amoebaGpu->log, "Z1 %4d [%4d %4d] %4d %4d %4d %4d   %d %d\n", ii,
+            fprintf( amoebaGpu->log, "Z1 %4d %d [%4d %4d] %4d %4d %4d %4d   %d %d\n", ii,  axisType[ii],
                      multipoleParticleX[ii], multipoleParticleY[ii], multipoleParticleZ[ii],
                      amoebaGpu->psMultipoleAxisOffset->_pSysData[multipoleParticleZ[ii]],
                      maxIndices[multipoleParticleZ[ii]],
@@ -1931,7 +1935,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
     }
 
     if( amoebaGpu->log ){
-        std::string axisLabel[maxAxisType+1] = {  "ZThenX", "Bisector", "ZBisect", "ThreeFold", "ZOnly", "Unknown"};
+        std::string axisLabel[maxAxisType+1] = {  "ZThenX", "Bisector", "ZBisect", "ThreeFold", "ZOnly", "NoAxisType", "Unknown"};
         for( unsigned int kk = 0; kk < (maxAxisType+1); kk++ ){
             (void) fprintf( amoebaGpu->log, "%2u %10s atom count=%d\n", kk, axisLabel[kk].c_str(), axisTypeCount[kk] );
         }
@@ -1969,6 +1973,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
 
         int diff = maxIndices[ii] - amoebaGpu->psMultipoleAxisOffset->_pSysData[ii];
         if( diff > amoebaGpu->maxMapTorqueDifference ){
+            //fprintf( stderr, "maxMapTorqueDifference %d maxInd=%d off=%d diff=%d\n", ii, maxIndices[ii], amoebaGpu->psMultipoleAxisOffset->_pSysData[ii], diff );
             amoebaGpu->maxMapTorqueDifference = diff;
         }
     }

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

@@ -101,6 +101,18 @@ void amoebaMapTorqueToForce_kernel( float* torque, int maxDiff, float* tempElecF
     int axisAtom                  = multiPoleAtoms[threadId].z;
     int axisType                  = multiPoleAtoms[threadId].w;
 
+    // NoAxisType
+
+    if( axisType == 5 )
+    { 
+       int min                       = multiPoleAtomOffset[threadId];
+       int offset                    = 3*(threadId*(maxDiff + 1)-min); 
+       tempElecForce[offset  ]       = 0.0f;
+       tempElecForce[offset+1]       = 0.0f;
+       tempElecForce[offset+2]       = 0.0f;
+       return;
+    }
+
     vector[U][0]                  = atomCoord[threadId].x - atomCoord[axisAtom].x;
     vector[U][1]                  = atomCoord[threadId].y - atomCoord[axisAtom].y;
     vector[U][2]                  = atomCoord[threadId].z - atomCoord[axisAtom].z;
@@ -689,7 +701,7 @@ void cudaComputeAmoebaMapTorques( amoebaGpuContext amoebaGpu, CUDAStream<float>*
     // check that BLOCK_SIZE >= amoebaGpu->maxMapTorqueDifference
 
     if( amoebaGpu->maxMapTorqueDifference > BLOCK_SIZE ){
-        (void) fprintf( amoebaGpu->log, "block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
+        (void) fprintf( amoebaGpu->log, "cudaComputeAmoebaMapTorques: block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
                         BLOCK_SIZE, amoebaGpu->maxMapTorqueDifference );
         exit(-1);  
     }
@@ -859,7 +871,7 @@ void cudaComputeAmoebaMapTorquesAndAddTotalForce( amoebaGpuContext amoebaGpu,
     // check that BLOCK_SIZE >= amoebaGpu->maxMapTorqueDifference
 
     if( amoebaGpu->maxMapTorqueDifference > BLOCK_SIZE ){
-        (void) fprintf( amoebaGpu->log, "block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
+        (void) fprintf( amoebaGpu->log, "cudaComputeAmoebaMapTorquesAndAddTotalForce: block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
                         BLOCK_SIZE, amoebaGpu->maxMapTorqueDifference );
         exit(-1);  
     }
@@ -1061,7 +1073,7 @@ void cudaComputeAmoebaMapTorquesAndAddTotalForce2( amoebaGpuContext amoebaGpu,
     // check that BLOCK_SIZE >= amoebaGpu->maxMapTorqueDifference
 
     if( amoebaGpu->maxMapTorqueDifference > BLOCK_SIZE ){
-        (void) fprintf( amoebaGpu->log, "block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
+        (void) fprintf( amoebaGpu->log, "cudaComputeAmoebaMapTorquesAndAddTotalForce2: block size (%d) in amoebaMapTorqueReduce_kernel is too small ( > %d)! -- aborting.\n",
                         BLOCK_SIZE, amoebaGpu->maxMapTorqueDifference );
         exit(-1);
     }

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

@@ -53,10 +53,12 @@ __device__ static float normVector3( float* vector )
 
 #undef AMOEBA_DEBUG
 
-// ZThenX    == 0
-// Bisector  == 1
-// ZBisect   == 2
-// ThreeFold == 3
+// ZThenX     == 0
+// Bisector   == 1
+// ZBisect    == 2
+// ThreeFold  == 3
+// ZOnly      == 4
+// NoAxisType == 5
 
 __global__
 #if (__CUDA_ARCH__ >= 200)
@@ -293,7 +295,7 @@ void kCudaComputeLabFrameMoments_kernel( void )
             sum                 = normVector3( vectorZ );
         }
  
-    } else if( axisType == 4 ){ 
+    } else if( axisType >= 4 ){ 
 
         vectorX[0]             = 0.1f;
         vectorX[1]             = 0.1f;