Fixed problem w/ cavity force

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

@@ -1899,10 +1899,26 @@ void gpuSetAmoebaObcParameters( amoebaGpuContext amoebaGpu, float innerDielectri
 
     gpu->sim.preFactor               = -amoebaGpu->amoebaSim.electric*((1.0f/innerDielectric)-(1.0f/solventDielectric));
 
+    if( amoebaGpu->log ){
+        (void) fprintf( amoebaGpu->log,"gpuSetAmoebaObcParameters: cavity=%d dielectricOffset=%15.7e probeRadius=%15.7e surfaceAreaFactor=%15.7e\n", 
+                        includeCavityTerm, dielectricOffset, probeRadius, surfaceAreaFactor );
+        (void) fprintf( amoebaGpu->log,"                           gkc=%12.3f solventDielectric=%15.7e innerDielectric=%15.7e sim.preFactor=%15.7e\n", 
+                        amoebaGpu->amoebaSim.gkc, amoebaGpu->amoebaSim.dwater, amoebaGpu->amoebaSim.dielec, gpu->sim.preFactor );
+        (void) fprintf( amoebaGpu->log,"                           fc=%15.7e fd=%15.7e fq=%15.7e\n",
+                        amoebaGpu->amoebaSim.fc, amoebaGpu->amoebaSim.fq, amoebaGpu->amoebaSim.fq );
+        (void) fprintf( amoebaGpu->log,"\nRadius (r-off) scl*(r-off) scl\n" );
+        for (unsigned int i = 0; i < amoebaGpu->paddedNumberOfAtoms && i < 10; i++) 
+        {
+            (void) fprintf( amoebaGpu->log,"%6d %15.7e %15.7e %15.7e %15.7e\n", i,
+                            radius[i] , (*gpu->psObcData)[i].x, (*gpu->psObcData)[i].y, scale[i] );
+        }
+    }
+
     gpuRotationToLabFrameAllocate( amoebaGpu );
     gpuFixedEFieldAllocate( amoebaGpu );
     gpuElectrostaticAllocate( amoebaGpu );
     gpuKirkwoodAllocate( amoebaGpu );
+
 }
 
 static int encodeCell( unsigned int x, unsigned int y ){
@@ -3764,6 +3780,41 @@ void cudaLoadCudaFloatArray( int numberOfParticles, int entriesPerParticle, CUDA
     }
 }
 
+/**---------------------------------------------------------------------------------------
+
+   Load contents of arrays into vector
+
+   @param numberOfParticles    number of particles
+   @param entriesPerParticle   entries/particles array
+   @param array                cuda array
+   @param outputVector         output vector
+
+   --------------------------------------------------------------------------------------- */
+
+void cudaLoadCudaFloat2Array( int numberOfParticles, int entriesPerParticle, CUDAStream<float2>* array, VectorOfDoubleVectors& outputVector ) 
+{
+    // ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName = "cudaLoadCudaFloat2Array";
+
+    // ---------------------------------------------------------------------------------------
+
+    array->Download();
+    int runningIndex  = 0;
+    
+    outputVector.resize( numberOfParticles ); 
+ 
+    for( int ii = 0; ii < numberOfParticles; ii++ ){ 
+        if( entriesPerParticle > 0 ){
+            outputVector[ii].push_back( array->_pSysStream[0][runningIndex].x );
+        }
+        if( entriesPerParticle > 1 ){
+            outputVector[ii].push_back( array->_pSysStream[0][runningIndex].y );
+        }
+        runningIndex++;
+    }
+}
+
 /**---------------------------------------------------------------------------------------
 
    Load contents of arrays into vector
@@ -3779,7 +3830,7 @@ void cudaLoadCudaFloat4Array( int numberOfParticles, int entriesPerParticle, CUD
 {
     // ---------------------------------------------------------------------------------------
 
-    // static const std::string methodName = "cudaLoadCudaFloatArray";
+    // static const std::string methodName = "cudaLoadCudaFloat4Array";
 
     // ---------------------------------------------------------------------------------------
 
@@ -4089,6 +4140,7 @@ void trackMutualInducedIterations( amoebaGpuContext amoebaGpu, int iteration){
 // ---------------------------------------------------------------------------------------
 
     if( amoebaGpu->log == NULL || currentStep > 20000 )return;
+    //if( amoebaGpu->log == NULL )return;
 
     gpuContext gpu                       = amoebaGpu->gpuContext;
     currentStep++;
@@ -4117,7 +4169,7 @@ void trackMutualInducedIterations( amoebaGpuContext amoebaGpu, int iteration){
             fileId.push_back( currentStep );
         }    
         if( (currentStep % 20) == 0 || fileId[0] > 20 ){
-             (void) fprintf( amoebaGpu->log, "step=%d fileId=%d\n", currentStep, fileId[0] );
+             (void) fprintf( amoebaGpu->log, "step=%d fileId=%d iterations=%d\n", currentStep, fileId[0], iteration );
         }
         (void) fflush( amoebaGpu->log );
         VectorOfDoubleVectors outputVector;
@@ -4129,7 +4181,7 @@ void trackMutualInducedIterations( amoebaGpuContext amoebaGpu, int iteration){
             cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psInducedDipoleS,     outputVector );
             cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psInducedDipolePolarS,outputVector );
 */
-        cudaWriteVectorOfDoubleVectorsToFile( "CudaMIT", fileId, outputVector );
+        cudaWriteVectorOfDoubleVectorsToFile( "CudaMI", fileId, outputVector );
         int nansPresent = isNanOrInfinity( amoebaGpu->mutualInducedCurrentEpsilon );
         if( nansPresent == 0 ){ 
             for( int ii = 0; ii < gpu->natoms && nansPresent == 0; ii++ ){

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

@@ -140,6 +140,7 @@ extern void cudaWriteFloat1AndFloat1ArraysToFile( int numberOfAtoms, char* fname
 extern void readFile( std::string fileName, StringVectorVector& fileContents );
  
 extern void cudaLoadCudaFloatArray( int numberOfParticles, int entriesPerParticle, CUDAStream<float>* array, VectorOfDoubleVectors& outputVector );
+extern void cudaLoadCudaFloat2Array( int numberOfParticles, int entriesPerParticle, CUDAStream<float2>* array, VectorOfDoubleVectors& outputVector );
 extern void cudaLoadCudaFloat4Array( int numberOfParticles, int entriesPerParticle, CUDAStream<float4>* array, VectorOfDoubleVectors& outputVector );
 extern void cudaWriteVectorOfDoubleVectorsToFile( char* fname, std::vector<int>& fileId, VectorOfDoubleVectors& outputVector );
 

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

@@ -4,6 +4,7 @@
 
 #include "amoebaGpuTypes.h"
 #include "amoebaCudaKernels.h"
+#include "cudaKernels.h"
 #include "kCalculateAmoebaCudaUtilities.h"
 #include "kCalculateAmoebaCudaKirkwoodParticle.h"
 extern void kCalculateObcGbsaForces2(gpuContext gpu);
@@ -1693,7 +1694,7 @@ __launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
 __launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
 #endif
  void kReduceToBornForcePrefactor_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* fieldIn1, float* fieldIn2,
-                                                    float* fieldOut )
+                                          float* fieldOut )
 {
     unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
 
@@ -1750,7 +1751,7 @@ __launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
 __launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
 #endif
 void kReduceToBornForcePrefactorAndSASA_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* fieldIn1, float* fieldIn2,
-                                                           float* fieldOut )
+                                                float* fieldOut )
 {
     unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
 
@@ -1802,13 +1803,13 @@ void kReduceToBornForcePrefactorAndSASA_kernel( unsigned int fieldComponents, un
               ratio6   = ratio6*ratio6;
         float saTerm   = cSim.surfaceAreaFactor * r * r * ratio6;
 
-        totalForce     += saTerm / bornRadius;
-        totalForce     *= bornRadius * bornRadius * obcChain;
+        totalForce    += saTerm / bornRadius;
+        totalForce    *= bornRadius * bornRadius * obcChain;
 
-        energy         += saTerm;
+        fieldOut[pos]  = totalForce;
 
-        fieldOut[pos]   = totalForce*bornRadius*bornRadius*obcChain;
-        pos            += gridDim.x * blockDim.x;
+        energy        += saTerm;
+        pos           += gridDim.x * blockDim.x;
     }
 
     cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy / -6.0f;
@@ -1830,11 +1831,48 @@ static void kReduceToBornForcePrefactor( amoebaGpuContext amoebaGpu )
 {
 
     if( amoebaGpu->includeObcCavityTerm ){
+
         kReduceToBornForcePrefactorAndSASA_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
                                                     amoebaGpu->paddedNumberOfAtoms, amoebaGpu->outputBuffers,
                                                     amoebaGpu->psWorkArray_1_1->_pDevStream[0],
                                                     amoebaGpu->psWorkArray_1_2->_pDevStream[0],
                                                     amoebaGpu->gpuContext->psBornForce->_pDevStream[0] );
+#ifdef AMOEBA_DEBUG
+    if( amoebaGpu->log ){
+
+        // kClearEnergy() should be called prior to kReduceToBornForcePrefactorAndSASA_kernel
+
+        double energy = kReduceEnergy( amoebaGpu->gpuContext );
+        amoebaGpu->gpuContext->psBornForce->Download();
+        amoebaGpu->gpuContext->psObcData->Download();
+        amoebaGpu->gpuContext->psBornRadii->Download();
+        (void) fprintf( amoebaGpu->log, "Born force w/ cavity energy=%15.7e.\n", energy ); (void) fflush( amoebaGpu->log );
+
+        for( int ii = 0; ii < amoebaGpu->gpuContext->natoms; ii++ ){
+           (void) fprintf( amoebaGpu->log, "%5d ", ii);
+           (void) fprintf( amoebaGpu->log,"bF %16.9e obc=%16.9e bR=%16.9e\n",
+                           amoebaGpu->gpuContext->psBornForce->_pSysStream[0][ii],
+                           amoebaGpu->gpuContext->psObcData->_pSysStream[0][ii].x,
+                           amoebaGpu->gpuContext->psBornRadii->_pSysStream[0][ii] );
+        }
+        (void) fflush( amoebaGpu->log );
+        if( 1 ){
+            std::vector<int> fileId;
+            //fileId.push_back( 0 );
+            VectorOfDoubleVectors outputVector;
+            cudaLoadCudaFloat4Array( amoebaGpu->gpuContext->natoms,  3, amoebaGpu->gpuContext->psPosq4,       outputVector );
+            cudaLoadCudaFloatArray(  amoebaGpu->gpuContext->natoms,  1, amoebaGpu->gpuContext->psBornRadii,   outputVector );
+            cudaLoadCudaFloat2Array( amoebaGpu->gpuContext->natoms,  2, amoebaGpu->gpuContext->psObcData,     outputVector );
+            cudaLoadCudaFloatArray(  amoebaGpu->gpuContext->natoms,  1, amoebaGpu->gpuContext->psBornForce,   outputVector );
+            cudaWriteVectorOfDoubleVectorsToFile( "CudaBornForce", fileId, outputVector );
+            (void) fprintf( amoebaGpu->log, "kReduceToBornForcePrefactor: exiting.\n" );
+            (void) fprintf( stderr, "kReduceToBornForcePrefactor: exiting.\n" ); (void) fflush( stderr );
+            exit(0);
+        }
+
+    }
+#endif
+
     } else {
         kReduceToBornForcePrefactor_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
                                              amoebaGpu->paddedNumberOfAtoms, amoebaGpu->outputBuffers,