bornRadiusScaleFactors were not included in second GB/VI force loop

removed some cruft

plugins/freeEnergy/CMakeLists.txt

@@ -172,6 +172,7 @@ SET(OPENMM_FREE_ENERGY_BUILD_VERSION 0)
 ADD_DEFINITIONS(-DOPENMM_FREE_ENERGY_LIBRARY_NAME=${OPENMM_FREE_ENERGY_LIBRARY_NAME}
                 -DOPENMM_FREE_ENERGY_MAJOR_VERSION=${OPENMM_FREE_ENERGY_MAJOR_VERSION}
                 -DOPENMM_FREE_ENERGY_MINOR_VERSION=${OPENMM_FREE_ENERGY_MINOR_VERSION}
+                -DINCLUDE_FREE_ENERGY_PLUGIN
 		          -DOPENMM_FREE_ENERGY_BUILD_VERSION=${OPENMM_FREE_ENERGY_BUILD_VERSION})
 
 # CMake quotes automatically when building Visual Studio projects but we need

plugins/freeEnergy/platforms/cuda/src/CudaFreeEnergyKernels.cpp

@@ -386,8 +386,8 @@ void CudaFreeEnergyCalcNonbondedSoftcoreForceKernel::executeForces(ContextImpl&
 // ---------------------------------------------------------------------------------------
 
     _gpuContext* gpu = data.gpu;
-//    static int call  = 0;
-//    call++;
+    static int call  = 0;
+    call++;
 
     // write array, ... address's to board
 
@@ -425,7 +425,9 @@ void CudaFreeEnergyCalcNonbondedSoftcoreForceKernel::executeForces(ContextImpl&
   
     // local LJ-14 forces
 
+//kPrintForces( gpu, "Pre  kCalculateLocalSoftcoreForces ", call );
     kCalculateLocalSoftcoreForces(gpu);
+//kPrintForces( gpu, "Post kCalculateLocalSoftcoreForces ", call );
 //kPrintForces(gpu, "Post kCalculateLocalSoftcoreForces", call );
 //kReduceForces(gpu);
 }
@@ -681,14 +683,11 @@ void CudaFreeEnergyCalcGBVISoftcoreForceKernel::executeForces(ContextImpl& conte
     if( setSim == 0 ){
        setSim++;
        SetCalculateGBVISoftcoreBornSumGpuSim( gpu );
+       SetCalculateObcGbsaSoftcoreBornSumSim( gpu );
        SetCalculateCDLJObcGbsaSoftcoreGpu1Sim( gpu );
        SetCalculateGBVISoftcoreForces2Sim( gpu );
     }
 
-    // required?
-
-    gpu->bRecalculateBornRadii = true; // fixed
-
     // calculate Born radii and first loop of GB/VI forces
 
     if( debug && log ){
@@ -700,16 +699,28 @@ void CudaFreeEnergyCalcGBVISoftcoreForceKernel::executeForces(ContextImpl& conte
         }
     }
 
+    // In kCalculateObcGbsaSoftcoreBornSum: SetCalculateObcGbsaSoftcoreBornSumSim
     kClearSoftcoreBornForces(gpu);
 
+    // In kCalculateGBVISoftcoreBornSum: SetCalculateGBVISoftcoreBornSumGpuSim
     kCalculateGBVISoftcoreBornSum(gpu);
 
     if( getQuinticScaling() ){
+
+        // kCalculateGBVISoftcoreBornSum.cu
+
         kReduceGBVIBornSumQuinticScaling(gpu, gpuGBVISoftcore );
     } else {
+
+        // In kCalculateGBVISoftcoreBornSum.cu
+
         kReduceGBVISoftcoreBornSum(gpu);
     }
 
+    // In kCalculateCDLJObcGbsaSoftcoreForces1.cu
+    //    SetCalculateCDLJObcGbsaSoftcoreGpu1Sim
+    //    SetCalculateCDLJObcGbsaSoftcoreSupplementary1Sim (called in GpuNonbondedSoftcore.cpp)
+
     kCalculateCDLJObcGbsaSoftcoreForces1(gpu);
 
     if( debug && log ){
@@ -721,13 +732,16 @@ void CudaFreeEnergyCalcGBVISoftcoreForceKernel::executeForces(ContextImpl& conte
     // compute Born forces
 
     if( getQuinticScaling() ){
+
+        // In kCalculateGBVISoftcoreBornSum.cu
+
         kReduceGBVIBornForcesQuinticScaling(gpu);
     } else {
-        gpu->bIncludeGBVI = true;
+
+        // In kCalculateGBVISoftcoreBornSum.cu
         kReduceGBVISoftcoreBornForces(gpu);
-        gpu->bIncludeGBVI = false;
-    }
 
+    }
     if( debug && log ){
         (void) fprintf( log, "\n%s: calling kCalculateGBVIForces2\n", methodName.c_str() );
         (void) fflush( log );
@@ -735,7 +749,10 @@ void CudaFreeEnergyCalcGBVISoftcoreForceKernel::executeForces(ContextImpl& conte
 
     // second loop of GB/VI forces
 
+    // In kCalculateGBVISoftcoreForces2.cu (SetCalculateGBVISoftcoreForces2Sim)
+
     kCalculateGBVISoftcoreForces2(gpu);
+//kPrintForces( gpu, "Post GBVISoftcoreForces2", call );
 }
 
 double CudaFreeEnergyCalcGBVISoftcoreForceKernel::executeEnergy(ContextImpl& context) {

plugins/freeEnergy/platforms/cuda/src/kernels/GpuFreeEnergyCudaKernels.h

@@ -114,6 +114,7 @@ extern void kReduceGBVIBornForcesQuinticScaling( gpuContext gpu );
 extern void kCalculateGBVISoftcoreForces2( gpuContext gpu );
 extern void kReduceGBVISoftcoreBornForces(gpuContext gpu);
 extern void kReduceGBVISoftcoreBornSum(gpuContext gpu);
+extern void kPrintGBVISoftcore(gpuContext gpu, GpuGBVISoftcore* gpuGBVISoftcore, std::string callId, int call);
 
 extern void kClearSoftcoreBornForces(gpuContext gpu);
 

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateCDLJObcGbsaSoftcoreForces1.cu

@@ -86,8 +86,8 @@ void GetCalculateCDLJObcGbsaSoftcoreForces1Sim( gpuContext gpu )
     RTERROR(status, "cudaMemcpyFromSymbol: SetSim copy from cSim failed");
 
 }
-texture<float, 1, cudaReadModeElementType> tabulatedErfcRef;
 
+#if 0
 __device__ float fastErfc(float r)
 {
     float normalized = cSim.tabulatedErfcScale*r;
@@ -96,6 +96,7 @@ __device__ float fastErfc(float r)
     float fract1 = 1.0f-fract2;
     return fract1*tex1Dfetch(tabulatedErfcRef, index) + fract2*tex1Dfetch(tabulatedErfcRef, index+1);
 }
+#endif
 
 // Include versions of the kernel for N^2 calculations.
 

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreAux.h

@@ -29,8 +29,8 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#ifndef __GpuGBVIAUX_H__
-#define __GpuGBVIAUX_H__
+#ifndef __Gpu_GBVI_SOFTCORE_AUX_H__
+#define __Gpu_GBVI_SOFTCORE_AUX_H__
 
 /**
  * This file contains subroutines used in evaluating quantities associated w/ the GB/VI function
@@ -130,4 +130,4 @@ __device__ float getGBVIBornForce2( float bornRadius, float R, float bornForce,
 
 }
 
-#endif // __GpuGBVIAUX_H__
+#endif // __Gpu_GBVI_SOFTCORE_AUX_H__

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreBornSum.cu

@@ -502,52 +502,47 @@ void kReduceGBVIBornForcesQuinticScaling(gpuContext gpu)
     //printf("kReduceObcGbsaBornForces\n");
     kReduceGBVIBornForcesQuinticScaling_kernel<<<gpu->sim.blocks, gpu->sim.bf_reduce_threads_per_block>>>();
     LAUNCHERROR("kReduceGBVIBornForcesQuinticScaling");
-
 }
 
-void kCalculateGBVISoftcoreBornSum(gpuContext gpu)
+void kPrintGBVISoftcore(gpuContext gpu, GpuGBVISoftcore* gpuGBVISoftcore, std::string callId, int call)
 {
-    //printf("kCalculateGBVIBornSum\n");
-    kClearGBVISoftcoreBornSum( gpu );
-    LAUNCHERROR("kClearGBVIBornSum from kCalculateGBVISoftcoreBornSum");
-
-    //size_t numWithInteractions;
-    switch (gpu->sim.nonbondedMethod)
-    {
-        case NO_CUTOFF:
-
-#define GBVI 0
-#if GBVI == 1
-int maxPrint = 31;
-gpu->psWorkUnit->Download();
-fprintf( stderr, "kCalculateGBVISoftcoreBornSum: bOutputBufferPerWarp=%u blks=%u th/blk=%u wu=%u %u shrd=%u\n", gpu->bOutputBufferPerWarp,
-                 gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block, gpu->sim.workUnits, gpu->psWorkUnit->_pSysStream[0][0],
-        sizeof(Atom)*gpu->sim.nonbond_threads_per_block );
-
+    int maxPrint = 20;
+    (void) fprintf( stderr, "kPrintGBVgSoftcore %s %d\n", callId.c_str(), call );
     gpu->psGBVIData->Download();
-               gpu->psBornSum->Download();
+    gpu->psBornRadii->Download();
+    gpu->psBornForce->Download();
     gpu->psPosq4->Download();
+    CUDAStream<float>* switchDeriviative = gpuGBVISoftcore-> getSwitchDerivative( );
+    switchDeriviative->Download();
 
-                (void) fprintf( stderr, "\nkCalculateGBVISoftcoreBornSum: pre BornSum %s Born radii & params\n",
-                               (gpu->bIncludeGBVI ? "GBVI" : "Obc") );
+    (void) fprintf( stderr, "BornSum Born radii & params\n" );
     for( int ii = 0; ii < gpu->natoms; ii++ ){
-                   (void) fprintf( stderr, "%6d bSum=%14.6e param[%14.6e %14.6e %14.6e] x[%14.6f %14.6f %14.6f %14.6f]\n",
+        (void) fprintf( stderr, "%6d prm[%14.6e %14.6e %14.6e] bR=%14.6e bF=%14.6e swDrv=%14.6e x[%14.6f %14.6f %14.6f %14.6f]\n",
                         ii,
-                                   gpu->psBornSum->_pSysStream[0][ii],
                         gpu->psGBVIData->_pSysStream[0][ii].x,
                         gpu->psGBVIData->_pSysStream[0][ii].y,
                         gpu->psGBVIData->_pSysStream[0][ii].z,
+                        gpu->psBornRadii->_pSysStream[0][ii],
+                        gpu->psBornForce->_pSysStream[0][ii],
+                        switchDeriviative->_pSysStream[0][ii],
                         gpu->psPosq4->_pSysStream[0][ii].x, gpu->psPosq4->_pSysStream[0][ii].y,
-                                   gpu->psPosq4->_pSysStream[0][ii].z, gpu->psPosq4->_pSysStream[0][ii].w
-                                 );
+                        gpu->psPosq4->_pSysStream[0][ii].z, gpu->psPosq4->_pSysStream[0][ii].w );
         if( (ii == maxPrint) && ( ii < (gpu->natoms - maxPrint)) ){
             ii = gpu->natoms - maxPrint;
         }
     }
+}
 
-#endif
-#undef GBVI
+void kCalculateGBVISoftcoreBornSum(gpuContext gpu)
+{
+    //printf("kCalculateGBVIBornSum\n");
+    kClearGBVISoftcoreBornSum( gpu );
+    LAUNCHERROR("kClearGBVIBornSum from kCalculateGBVISoftcoreBornSum");
 
+    //size_t numWithInteractions;
+    switch (gpu->sim.nonbondedMethod)
+    {
+        case NO_CUTOFF:
 
             if (gpu->bOutputBufferPerWarp){
                 kCalculateGBVISoftcoreN2ByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
@@ -576,5 +571,5 @@ fprintf( stderr, "kCalculateGBVISoftcoreBornSum: bOutputBufferPerWarp=%u blks=%u
             break;
 #endif
     }
-    LAUNCHERROR("kCalculateGBVIBornSum");
+    LAUNCHERROR("kCalculateGBVISoftcoreBornSum");
 }

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreForces2.cu

@@ -51,6 +51,7 @@ struct Atom {
     float fy;
     float fz;
     float fb;
+    float bornRadiusScaleFactor;
 };
 
 
@@ -61,6 +62,7 @@ void SetCalculateGBVISoftcoreForces2Sim(gpuContext gpu)
     cudaError_t status;
     status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));     
     RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
+    //(void) fprintf( stderr, "SetCalculateGBVISoftcoreForces2Sim called.\n" );
 }
 
 void GetCalculateGBVISoftcoreForces2Sim(gpuContext gpu)

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreForces2.h

@@ -29,7 +29,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "kCalculateGBVIAux.h"
+#include "kCalculateGBVISoftcoreAux.h"
 
 /**
  * This file contains the kernel for evalauating the second stage of GBSA.  It is included
@@ -77,6 +77,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
             sA[threadIdx.x].z                       = apos.z;
             sA[threadIdx.x].r                       = ar.x;
             sA[threadIdx.x].sr                      = ar.y;
+            sA[threadIdx.x].bornRadiusScaleFactor   = ar.w;
             sA[threadIdx.x].fb                      = fb;
 
             for (unsigned int j = (tgx+1)&(GRID-1); j != tgx; j = (j+1)&(GRID-1))
@@ -93,7 +94,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
                 float r                 = sqrt(r2);
 
                 // Atom I Born forces and sum
-                float dE                = getGBVI_dE2( r, ar.x, psA[j].sr, fb );
+                float dE                = psA[j].bornRadiusScaleFactor*getGBVI_dE2( r, ar.x, psA[j].sr, fb );
                
 #if defined USE_PERIODIC
                 if (i >= cSim.atoms || x+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
@@ -152,6 +153,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
                 sA[threadIdx.x].z                       = temp.z;
                 sA[threadIdx.x].r                       = temp1.x;
                 sA[threadIdx.x].sr                      = temp1.y;
+                sA[threadIdx.x].bornRadiusScaleFactor   = temp1.w;
             }
 #ifdef USE_CUTOFF
             unsigned int flags = cSim.pInteractionFlag[pos];
@@ -177,7 +179,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
                     float r2                = dx * dx + dy * dy + dz * dz;
                     float r                 = sqrt(r2);
 
-                    float dE                = getGBVI_dE2( r, ar.x, psA[tj].sr, fb );
+                    float dE                = psA[tj].bornRadiusScaleFactor*getGBVI_dE2( r, ar.x, psA[tj].sr, fb );
 
 #if defined USE_PERIODIC
                     if (i >= cSim.atoms || y+tj >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
@@ -203,7 +205,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
                     psA[tj].fz             += d;
 
                     // Atom J Born sum term
-                    dE                      = getGBVI_dE2( r, psA[tj].r, ar.y, psA[tj].fb );
+                    dE                      = ar.w*getGBVI_dE2( r, psA[tj].r, ar.y, psA[tj].fb );
 
 #ifdef USE_PERIODIC
                     if (i >= cSim.atoms || y+tj >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
@@ -250,7 +252,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
                         float r                 = sqrt(r2);
 
                         // Interleaved Atom I and J Born Forces and sum components
-                        float dE                = getGBVI_dE2( r, ar.x, psA[j].sr, fb );
+                        float dE                = psA[j].bornRadiusScaleFactor*getGBVI_dE2( r, ar.x, psA[j].sr, fb );
 
 #if defined USE_PERIODIC
                         if (i >= cSim.atoms || y+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
@@ -276,7 +278,7 @@ __global__ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int
                         tempBuffer[threadIdx.x].z = d;
 
                         // Atom J Born sum term
-                        dE                      = getGBVI_dE2( r, psA[j].r, ar.y, psA[j].fb );
+                        dE                      = ar.w*getGBVI_dE2( r, psA[j].r, ar.y, psA[j].fb );
 
 #ifdef USE_PERIODIC
                         if (i >= cSim.atoms || y+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateNonbondedSoftcore.cu

@@ -427,7 +427,7 @@ fprintf( stderr, "kCalculateCDLJSoftcoreForces: bOutputBufferPerWarp=%u blks=%u
 void kPrintForces(gpuContext gpu, std::string idString, int call )
 {
  //   printf("kReduceForces\n");
-#define GBVI_DEBUG 0
+#define GBVI_DEBUG 4
 #if ( GBVI_DEBUG == 4 )
 
                 gpu->psBornRadii->Download();
@@ -436,9 +436,9 @@ void kPrintForces(gpuContext gpu, std::string idString, int call )
                 gpu->psBornForce->Download();
                 gpu->psForce4->Download();
                 gpu->psPosq4->Download();
-                int maxPrint = 0; 
+                int maxPrint = 30; 
 int   nanHit       = 0;
-int   targetIndex  = 852;
+int   targetIndex  = -852;
 float maxForce     = 3.0e+04;
 float maxPosition  = 2.0e+02;
                 for( int ii = 0; ii < gpu->natoms; ii++ ){
@@ -471,8 +471,8 @@ if( isnan( gpu->psBornForce->_pSysStream[0][ii] ) ||
    nanHit = 1;
 }
 
-                //if( hit || ii < maxPrint || ii >= (gpu->natoms - maxPrint) )
-                if( hit ){
+                if( hit || ii < maxPrint || ii >= (gpu->natoms - maxPrint) ){
+                //if( hit ){
                     static int firstHit = 1;
                     if( firstHit ){
                        firstHit = 0;

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreBornSum.cu

@@ -97,7 +97,7 @@ __global__ void kClearSoftcoreBornForces_kernel()
 
 void kClearSoftcoreBornForces(gpuContext gpu)
 {
-  //  printf("kClearBornForces\n");
+  //  printf("kClearSoftcoreBornForces\n");
     kClearSoftcoreBornForces_kernel<<<gpu->sim.blocks, 384>>>();
     LAUNCHERROR("kClearSoftcoreBornForces");
 }

plugins/freeEnergy/platforms/cuda/tests/CMakeLists.txt

 #
 # Testing
 #
-
 ENABLE_TESTING()
-
 INCLUDE(${CMAKE_SOURCE_DIR}/platforms/cuda/cuda-cmake/FindCuda.cmake)
 INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
 INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/include)
@@ -48,9 +46,9 @@ ENDFOREACH(TEST_PROG ${TEST_PROGS})
 
 # TestCudaUsingParameterFile customized w/ command-line argument (input file name used in test) 
 
-#ADD_EXECUTABLE(TestCudaUsingParameterFile TstCudaUsingParameterFile.cpp)
-#TARGET_LINK_LIBRARIES(TestCudaUsingParameterFile ${SHARED_TARGET})
-#ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/lambdaSdObcParameters.txt")
+ADD_EXECUTABLE(TestFreeEnergyCudaUsingParameterFile TstFreeEnergyCudaUsingParameterFile.cpp)
+TARGET_LINK_LIBRARIES(TestFreeEnergyCudaUsingParameterFile ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET})
+ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/lambdaSdObcParameters.txt")
 #ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfNoPbcParameters.txt")
 #
 #SET(TEST_ROOT TestCudaUsingParameterFile)

plugins/freeEnergy/platforms/cuda/tests/TstFreeEnergyCudaUsingParameterFile.cpp

@@ -5054,12 +5054,15 @@ void testReferenceCudaForces( std::string parameterFileName, MapStringInt& force
 
    // Run several steps and see if relative force difference is within tolerance
    
+   int includeEnergy = 1;
    for( int step = 0; step < numberOfSteps; step++ ){
 
       // pull info out of contexts
 
       int types                                       = State::Positions | State::Velocities | State::Forces;
-      //int types                                       = State::Positions | State::Velocities | State::Forces | State::Energy;
+      if( includeEnergy ){
+         types                                       |= State::Energy;
+      }
 
       State cudaState                                 =      cudaContext->getState( types );
       State referenceState                            = referenceContext->getState( types );
@@ -5068,20 +5071,29 @@ void testReferenceCudaForces( std::string parameterFileName, MapStringInt& force
       std::vector<Vec3> referenceVelocities           = referenceState.getVelocities();
       std::vector<Vec3> referenceForces               = referenceState.getForces();
 
-//      double referenceKineticEnergy                   = referenceState.getKineticEnergy();
-//      double referencePotentialEnergy                 = referenceState.getPotentialEnergy();
-
-double referenceKineticEnergy                   = 0.0;
-double referencePotentialEnergy                 = 0.0;
+      double referenceKineticEnergy;
+      double referencePotentialEnergy;
+      if( includeEnergy ){
+         referenceKineticEnergy                       = referenceState.getKineticEnergy();
+         referencePotentialEnergy                     = referenceState.getPotentialEnergy();
+      } else {
+         referenceKineticEnergy                       = 0.0;
+         referencePotentialEnergy                     = 0.0;
+      }
 
       std::vector<Vec3> cudaCoordinates               = cudaState.getPositions();
       std::vector<Vec3> cudaVelocities                = cudaState.getVelocities();
       std::vector<Vec3> cudaForces                    = cudaState.getForces();
 
-//      double cudaKineticEnergy                        = cudaState.getKineticEnergy();
-//      double cudaPotentialEnergy                      = cudaState.getPotentialEnergy();
-double cudaKineticEnergy                        = 0.0;
-double cudaPotentialEnergy                      = 0.0;
+      double cudaKineticEnergy                        = cudaState.getKineticEnergy();
+      double cudaPotentialEnergy                      = cudaState.getPotentialEnergy();
+      if( includeEnergy ){
+         cudaKineticEnergy                            = cudaState.getKineticEnergy();
+         cudaPotentialEnergy                          = cudaState.getPotentialEnergy();
+      } else {
+         cudaKineticEnergy                            = 0.0;
+         cudaPotentialEnergy                          = 0.0;
+      }
 
       // diagnostics
 

plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.cpp

@@ -727,7 +727,7 @@ int CpuGBVISoftcore::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM*
 
    // ---------------------------------------------------------------------------------------
 
-   static const char* methodName        = "CpuGBVISoftcore::computeBornEnergyForces";
+   static const char* methodName        = "CpuGBVISoftcore::computeBornForces";
 
    static const RealOpenMM zero         = (RealOpenMM) 0.0;
    static const RealOpenMM one          = (RealOpenMM) 1.0;
@@ -852,7 +852,7 @@ if( atomI == 0 ){
 {
    double stupidFactor                   = three/CAL_TO_JOULE;
    RealOpenMM conversion                 = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau());  
-   int maxPrint                          = 200000;
+   int maxPrint                          = 20;
    const RealOpenMM* scaledRadii         = gbviParameters->getScaledRadii();
    RealOpenMM* switchDeriviative         = getSwitchDeriviative();
 
@@ -866,9 +866,9 @@ if( atomI == 0 ){
       double xx           = switchDeriviative[atomI]*bF*oneThird*b2*b2;
 
       // xx*conversion should agree w/ values pulled out of kReduceGBVISoftcoreBornForces_kernel in kForces.cu
-      (void) fprintf( logFile, "F1 %6d r/sclR[%14.6e %14.6e] bR=%14.6e bF=%14.6e sw=%14.6e %14.6e f[%14.6e %14.6e %14.6e](cnvrtd)"
+      (void) fprintf( logFile, "F1 %6d r/sclR[%14.6e %14.6e] bR=%14.6e bF=%14.6e sw=%14.6e f[%14.6e %14.6e %14.6e](cnvrtd)"
                                " x[%14.6e %14.6e %14.6e]\n",
-                      atomI, atomicRadii[atomI], scaledRadii[atomI], bornRadii[atomI], bF, switchDeriviative[atomI], xx*conversion,
+                      atomI, atomicRadii[atomI], scaledRadii[atomI], bornRadii[atomI], xx*conversion, switchDeriviative[atomI],
                       conversion*forces[atomI][0], conversion*forces[atomI][1],  conversion*forces[atomI][2],
                       atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
       if( atomI == maxPrint ){
@@ -877,6 +877,15 @@ if( atomI == 0 ){
       }
    }
    (void) fflush( logFile );
+   int clearForces = 0;
+   if( clearForces ){
+      (void) fprintf( logFile, "Forces cleared after loop 1\n" );
+      for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+          forces[atomI][0] = 0.0f;
+          forces[atomI][1] = 0.0f;
+          forces[atomI][2] = 0.0f;
+      }
+   }
 }
 #endif
 
@@ -935,6 +944,7 @@ if( atomI == 0 ){
    const RealOpenMM* scaledRadii                    = gbviParameters->getScaledRadii();
    RealOpenMM* switchDeriviative                    = getSwitchDeriviative();
    double stupidFactor                              = three/CAL_TO_JOULE;
+   const RealOpenMM* bornRadiusScaleFactors         = gbviParameters->getBornRadiusScaleFactors();
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
  
       RealOpenMM R        = atomicRadii[atomI];
@@ -973,11 +983,9 @@ if( atomI == 0 ){
             RealOpenMM S                  = scaledRadii[atomJ];
             RealOpenMM diff               = (S - R);
 
-            RealOpenMM de                 = zero;
-
             // find dRb/dr, where Rb is the Born radius
 
-            de = CpuGBVISoftcore::dL_dr( r, r+S, S ) + CpuGBVISoftcore::dL_dx( r, r+S, S );   
+            RealOpenMM de = CpuGBVISoftcore::dL_dr( r, r+S, S ) + CpuGBVISoftcore::dL_dx( r, r+S, S );   
             if( FABS( diff ) < r ){
                if( R > (r - S) ){
                   de -= CpuGBVISoftcore::dL_dr( r, R, S );  
@@ -988,9 +996,8 @@ if( atomI == 0 ){
                de -= ( CpuGBVISoftcore::dL_dr( r, r-S, S ) + CpuGBVISoftcore::dL_dx( r, r-S, S ) );   
             }
 
+
 #if 0
-   RealOpenMM delta = (RealOpenMM) 1.0e-02;
-   (void) fprintf( stderr, "\n" );
    for( int kk = 0; kk < 5; kk++ ){
       RealOpenMM V1    = CpuGBVISoftcore::getVolume( r, R, S );
       RealOpenMM V2    = CpuGBVISoftcore::getVolume( r+delta, R, S );
@@ -1010,10 +1017,9 @@ if( atomI == 0 ){
    }
 #endif
 
-
              // de = (dG/dRb)(dRb/dr)
 
-            de                      *= bornForces[atomI]/r;
+            de                      *= bornRadiusScaleFactors[atomJ]*bornForces[atomI]/r;
 
             deltaX                  *= de;
             deltaY                  *= de;
@@ -1027,12 +1033,23 @@ if( atomI == 0 ){
             forces[atomJ][1]        -= deltaY;
             forces[atomJ][2]        -= deltaZ;
 
+#if 0
+   if( atomI == 2613 ){
+   (void) fprintf( stderr, "AtomJ %5d r=%14.7e de=%14.7e bfI=%14.7e finalDe=%14.7e [%14.7e %14.7e %14.7e]\n",
+                   atomJ, r, de, bornForces[atomI], (de*bornForces[atomI]/r), 
+                   forces[atomI][0], forces[atomI][1], forces[atomI][2]  );
+   } else if( atomJ == 2613 ){
+   (void) fprintf( stderr, "AtomI %5d r=%14.7e de=%14.7e bfI=%14.7e finalDe=%14.7e [%14.7e %14.7e %14.7e]\n",
+                   atomI, r, de, bornForces[atomI], (de*bornForces[atomI]/r), 
+                   forces[atomJ][0], forces[atomJ][1], forces[atomJ][2]  );
+   }
+#endif
          }
       }
 
    }
 
-#if( GBVISoftcoreDebug == 1 )
+#if( GBVISoftcoreDebug == 9 )
 {
    (void) fprintf( logFile, "\nPre conversion\n" );
    (void) fprintf( logFile, "Atom        ScaledRadii    BornRadii      BornForce      SwitchDrv                                         Forces\n" );
@@ -1065,10 +1082,13 @@ if( atomI == 0 ){
 {
    (void) fprintf( logFile, "\nPost conversion\n" );
    (void) fprintf( logFile, "Atom        BornRadii      BornForce      SwitchDrv                                         Forces\n" );
-   int maxPrint = 1000000;
+   int maxPrint = 20;
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-      (void) fprintf( logFile, "%6d %14.6e %14.6e %14.6e [%14.6e %14.6e %14.6e] %s\n", atomI, bornRadii[atomI], conversion*bornForces[atomI],
-                      switchDeriviative[atomI], inputForces[atomI][0], inputForces[atomI][1], inputForces[atomI][2], (fabs( switchDeriviative[atomI] - 1.0 ) > 1.0e-05 ? "SWWWWW" : "") );
+      (void) fprintf( logFile, "%6d %14.6e %14.6e %14.6e 2[%14.6e %14.6e %14.6e] ttlF[%14.6e %14.6e %14.6e] %s\n",
+                      atomI, bornRadii[atomI], conversion*bornForces[atomI], switchDeriviative[atomI],
+                      conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*forces[atomI][2], 
+                      inputForces[atomI][0], inputForces[atomI][1], inputForces[atomI][2], 
+                      (fabs( switchDeriviative[atomI] - 1.0 ) > 1.0e-05 ? "SWWWWW" : "") );
       if( atomI == maxPrint ){
          atomI = numberOfAtoms - maxPrint;
          if( atomI < maxPrint )atomI = numberOfAtoms;
@@ -1120,246 +1140,6 @@ std::string CpuGBVISoftcore::getStateString( const char* title ) const {
    return message.str();
 }
 
-/**---------------------------------------------------------------------------------------
-
-   Write Born energy and forces (Simbios)
-
-   @param atomCoordinates     atomic coordinates
-   @param partialCharges      partial charges
-   @param forces              forces
-   @param resultsFileName     output file name
-
-   @return SimTKOpenMMCommon::DefaultReturn unless
-           file cannot be opened
-           in which case return SimTKOpenMMCommon::ErrorReturn
-
-   --------------------------------------------------------------------------------------- */
-
-int CpuGBVISoftcore::writeBornEnergyForces( RealOpenMM** atomCoordinates,
-                                   const RealOpenMM* partialCharges, RealOpenMM** forces,
-                                   const std::string& resultsFileName ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   static const char* methodName  = "\nCpuGBVISoftcore::writeBornEnergyForces";
-
-   // ---------------------------------------------------------------------------------------
-/*
-   ImplicitSolventParameters* implicitSolventParameters = getImplicitSolventParameters();
-   const GBVISoftcoreParameters* gbviParameters                   = static_cast<const GBVISoftcoreParameters*>(implicitSolventParameters);
-   
-
-   int numberOfAtoms                    = gbviParameters->getNumberOfAtoms();
-   const RealOpenMM* atomicRadii        = gbviParameters->getAtomicRadii();
-   const RealOpenMM* bornRadii          = getBornRadiiConst();
-   const RealOpenMM* scaledRadii        = gbviParameters->getScaledRadiusFactors();
-   const RealOpenMM* gbviChain           = getObcChainConst();
-   const RealOpenMM  energy             = getEnergy();
-
-   // ---------------------------------------------------------------------------------------
-
-   // open file -- return if unsuccessful
-
-   FILE* implicitSolventResultsFile = NULL;
-#ifdef WIN32
-   fopen_s( &implicitSolventResultsFile, resultsFileName.c_str(), "w" );
-#else
-   implicitSolventResultsFile = fopen( resultsFileName.c_str(), "w" );
-#endif
-
-   // diganostics
-
-   std::stringstream message;
-   message << methodName;
-   if( implicitSolventResultsFile != NULL ){
-      std::stringstream message;
-      message << methodName;
-      message << " Opened file=<" << resultsFileName << ">.";
-      SimTKOpenMMLog::printMessage( message );
-   } else {
-      std::stringstream message;
-      message << methodName;
-      message << "  could not open file=<" << resultsFileName << "> -- abort output.";
-      SimTKOpenMMLog::printMessage( message );
-      return SimTKOpenMMCommon::ErrorReturn;
-   }
-
-   // header
-
-   (void) fprintf( implicitSolventResultsFile, "# %d atoms E=%.7e   format: coords(3) bornRadii(input) q atomicRadii scaleFactors forces gbviChain\n",
-                   numberOfAtoms, energy );
-
-   RealOpenMM forceConversion  = (RealOpenMM) 1.0;
-   RealOpenMM lengthConversion = (RealOpenMM) 1.0;
-
-   // output
-
-   if( forces != NULL && atomCoordinates != NULL && partialCharges != NULL && atomicRadii != NULL ){
-      for( int ii = 0; ii < numberOfAtoms; ii++ ){
-            (void) fprintf( implicitSolventResultsFile, "%.7e %.7e %.7e %.7e %.5f %.5f %.5f %.7e %.7e %.7e %.7e\n",
-                            lengthConversion*atomCoordinates[ii][0],
-                            lengthConversion*atomCoordinates[ii][1], 
-                            lengthConversion*atomCoordinates[ii][2],
-                           (bornRadii != NULL ? lengthConversion*bornRadii[ii] : 0.0),
-                            partialCharges[ii], lengthConversion*atomicRadii[ii], scaledRadii[ii],
-                            forceConversion*forces[ii][0],
-                            forceConversion*forces[ii][1],
-                            forceConversion*forces[ii][2],
-                            forceConversion*gbviChain[ii]
-                          );
-      }
-   }
-   (void) fclose( implicitSolventResultsFile );
-
-*/
-   return 0;
-
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Write  results from first loop
-
-   @param numberOfAtoms       number of atoms
-   @param forces              forces
-   @param bornForce           Born force prefactor
-   @param outputFileName      output file name
-
-   @return SimTKOpenMMCommon::DefaultReturn unless
-           file cannot be opened
-           in which case return SimTKOpenMMCommon::ErrorReturn
-
-   --------------------------------------------------------------------------------------- */
-
-int CpuGBVISoftcore::writeForceLoop1( int numberOfAtoms, RealOpenMM** forces, const RealOpenMM* bornForce,
-                             const std::string& outputFileName ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName  = "\nCpuGBVISoftcore::writeForceLoop1";
-
-   // ---------------------------------------------------------------------------------------
-
-   int chunkSize;
-   if( bornForce ){
-      chunkSize = 3;
-   } else {
-      chunkSize = 4;
-   }
-
-   StringVector lineVector;
-   std::stringstream header;
-   lineVector.push_back( "# bornF F" );
-   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-      std::stringstream line;
-      line << (atomI+1) << " ";
-      SimTKOpenMMUtilities::formatRealStringStream( line, forces[atomI], chunkSize );
-      if( bornForce ){
-         line << " " << bornForce[atomI];
-      }
-      lineVector.push_back( line.str() );
-   }
-   return SimTKOpenMMUtilities::writeFile( lineVector, outputFileName );
-
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Write results
-
-   @param numberOfAtoms        number of atoms
-   @param chunkSizes           vector of chunk sizes for realRealOpenMMVector
-   @param realRealOpenMMVector vector of RealOpenMM**
-   @param realVector           vector of RealOpenMM*
-   @param outputFileName       output file name
-
-   @return SimTKOpenMMCommon::DefaultReturn unless
-           file cannot be opened
-           in which case return SimTKOpenMMCommon::ErrorReturn
-
-   --------------------------------------------------------------------------------------- */
-
-int CpuGBVISoftcore::writeForceLoop( int numberOfAtoms, const IntVector& chunkSizes,
-                            const RealOpenMMPtrPtrVector& realRealOpenMMVector, 
-                            const RealOpenMMPtrVector& realVector,
-                            const std::string& outputFileName ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName  = "\nCpuGBVISoftcore::writeForceLoop";
-
-   static const int maxChunks = 10;
-   int chunks[maxChunks];
-
-   // ---------------------------------------------------------------------------------------
-
-   for( int ii = 0; ii < (int) chunkSizes.size(); ii++ ){
-      chunks[ii] = chunkSizes[ii];
-   }
-   for( int ii = (int) chunkSizes.size(); ii < maxChunks; ii++ ){
-      chunks[ii] = 3;
-   }
-
-   StringVector lineVector;
-   std::stringstream header;
-   // lineVector.push_back( "# " );
-
-   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-
-      std::stringstream line;
-		char buffer[128];
-
-		(void) sprintf( buffer, "%4d ", atomI );
-		line << buffer;
-
-      int index = 0;
-      for( RealOpenMMPtrPtrVectorCI ii = realRealOpenMMVector.begin(); ii != realRealOpenMMVector.end(); ii++ ){
-         RealOpenMM** forces = *ii;
-			(void) sprintf( buffer, "%11.5f %11.5f %11.5f ", forces[atomI][0], forces[atomI][1], forces[atomI][2] );
-			line << buffer;
-//         SimTKOpenMMUtilities::formatRealStringStream( line, forces[atomI], chunks[index++] );
-//         line << " ";
-      }
-
-      for( RealOpenMMPtrVectorCI ii = realVector.begin(); ii != realVector.end(); ii++ ){
-         RealOpenMM* array = *ii;
-			(void) sprintf( buffer, "%11.5f ", array[atomI] );
-         line << buffer;
-      }
-
-      lineVector.push_back( line.str() );
-   }
-   return SimTKOpenMMUtilities::writeFile( lineVector, outputFileName );
-
-}
-
-/**---------------------------------------------------------------------------------------
-
-   Get Obc Born energy and forces -- used debugging
-
-   @param bornRadii           Born radii -- optional; if NULL, then GBVISoftcoreParameters 
-                              entry is used
-   @param atomCoordinates     atomic coordinates
-   @param partialCharges      partial charges
-   @param forces              forces
-
-   @return 0;
-
-   The array bornRadii is also updated and the obcEnergy
-
-   --------------------------------------------------------------------------------------- */
-
-int CpuGBVISoftcore::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
-                                      const RealOpenMM* partialCharges, RealOpenMM** forces ){
- 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nCpuGBVISoftcore::computeBornEnergyForcesPrint";
-
-   return 0;
-
-}
-
 /**---------------------------------------------------------------------------------------
 
    Use double precision 

plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.h

@@ -118,25 +118,6 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
       int computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii,
                             RealOpenMM* switchDeriviative = NULL );
       
-      /**---------------------------------------------------------------------------------------
-      
-         Get Born energy and forces (not used)
-      
-         @param bornRadii         Born radii
-         @param atomCoordinates   atomic coordinates
-         @param partialCharges    partial charges
-         @param forces            forces
-      
-         @return force array
-      
-         --------------------------------------------------------------------------------------- */
-      
-      int computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
-                                   const RealOpenMM* partialCharges, RealOpenMM** forces );
-      
-      int computeBornEnergyForcesPrint( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
-                                        const RealOpenMM* partialCharges, RealOpenMM** forces );
-      
       /**---------------------------------------------------------------------------------------
       
          Get state 
@@ -149,61 +130,6 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
       
       std::string getStateString( const char* title ) const;
 
-      /**---------------------------------------------------------------------------------------
-      
-         Write Born energy and forces (Simbios)
-      
-         @param atomCoordinates   atomic coordinates
-         @param partialCharges    partial atom charges
-         @param forces            force array
-         @param resultsFileName   output file name
-      
-         @return SimTKOpenMMCommon::DefaultReturn if file opened; else return SimTKOpenMMCommon::ErrorReturn
-      
-         --------------------------------------------------------------------------------------- */
-          
-      int writeBornEnergyForces( RealOpenMM** atomCoordinates,
-                                 const RealOpenMM* partialCharges, RealOpenMM** forces,
-                                 const std::string& resultsFileName ) const;
-
-      /**---------------------------------------------------------------------------------------
-      
-         Write  results from first loop
-      
-         @param atomCoordinates     atomic coordinates
-         @param RealOpenMM forces         forces
-         @param outputFileName      output file name
-      
-         @return SimTKOpenMMCommon::DefaultReturn unless
-                 file cannot be opened
-                 in which case return SimTKOpenMMCommon::ErrorReturn
-      
-         --------------------------------------------------------------------------------------- */
-      
-      static int writeForceLoop1( int numberOfAtoms, RealOpenMM** forces, const RealOpenMM* bornForce,
-                                  const std::string& outputFileName );
-      
-      /**---------------------------------------------------------------------------------------
-      
-         Write results
-      
-         @param numberOfAtoms       number of atoms
-         @param chunkSizes          vector of chunk sizes for realRealOpenMMVector
-         @param realRealOpenMMVector      vector of RealOpenMM**
-         @param realVector          vector of RealOpenMM*
-         @param outputFileName      output file name
-      
-         @return SimTKOpenMMCommon::DefaultReturn unless
-                 file cannot be opened
-                 in which case return SimTKOpenMMCommon::ErrorReturn
-      
-         --------------------------------------------------------------------------------------- */
-      
-      static int writeForceLoop( int numberOfAtoms, const IntVector& chunkSizes,
-                                 const RealOpenMMPtrPtrVector& realRealOpenMMVector, 
-                                 const RealOpenMMPtrVector& realVector,
-                                 const std::string& outputFileName );
-      
       /**---------------------------------------------------------------------------------------
       
          Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])