dded GB/VI to Cuda platform

Free energy plugin added
Plugin will not run w/ Obc or GB/VI forces unless line 2004 of gpu.cpp (gpu->sim.totalNonbondOutputBuffers  = 2*gpu->sim.nonbondOutputBuffers;) is commented in -- working on removing this constraint
Also unit tests for GB/VI currently fail 

CMakeLists.txt

@@ -215,6 +215,8 @@ IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
     SET(STATIC_TARGET ${STATIC_TARGET}_d)
 ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
 
+# used by plugin
+SET(OPENMM_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 
 # These are all the places to search for header files which are
 # to be part of the API.
@@ -303,6 +305,16 @@ IF(OPENMM_BUILD_OPENCL_LIB)
     ADD_SUBDIRECTORY(platforms/opencl)
 ENDIF(OPENMM_BUILD_OPENCL_LIB)
 
+# FreeEnergy plugin
+
+SET(OPENMM_BUILD_FREE_ENERGY_PLUGIN OFF CACHE BOOL "Build FreeEnergy plugin for Nvidia GPUs")
+SET(OPENMM_BUILD_FREE_ENERGY_PATH)
+IF(OPENMM_BUILD_FREE_ENERGY_PLUGIN)
+   SET(OPENMM_BUILD_FREE_ENERGY_PATH ${CMAKE_CURRENT_SOURCE_DIR}/plugins/freeEnergy)
+#MESSAGE("\nIn0 OpenMM: OPENMM_BUILD_FREE_ENERGY_PATH=" ${OPENMM_BUILD_FREE_ENERGY_PATH} )
+   ADD_SUBDIRECTORY(plugins/freeEnergy)
+ENDIF(OPENMM_BUILD_FREE_ENERGY_PLUGIN)
+
 INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_TARGET})
 INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_TARGET})
 FILE(GLOB CORE_HEADERS     include/*.h          */include/*.h)

olla/include/openmm/kernels.h

@@ -378,7 +378,7 @@ public:
 class CalcGBSAOBCForceKernel : public KernelImpl {
 public:
     static std::string Name() {
-        return "CalcGBSAOBCForces";
+        return "CalcGBSAOBCForce";
     }
     CalcGBSAOBCForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
     }
@@ -410,7 +410,7 @@ public:
 class CalcGBVIForceKernel : public KernelImpl {
 public:
     static std::string Name() {
-        return "CalcGBVIForces";
+        return "CalcGBVIForce";
     }
     CalcGBVIForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
     }

openmmapi/include/openmm/GBVIForce.h

@@ -50,7 +50,7 @@ namespace OpenMM {
 
 class OPENMM_EXPORT GBVIForce : public Force {
 public:
-    /**
+    /** 
      * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
      */
     enum NonbondedMethod {
@@ -68,7 +68,7 @@ public:
          * each other particle.  Interactions beyond the cutoff distance are ignored.
          */
         CutoffPeriodic = 2,
-    };
+    };  
     /*
      * Create a GBVIForce.
      */
@@ -107,6 +107,41 @@ public:
      * @param gamma          the gamma parameter
      */
     void setParticleParameters(int index, double charge, double radius, double gamma);
+    /**
+     * Add a bond 
+     *
+     * @param particle1 the index of the first particle 
+     * @param particle2 the index of the second particle
+     * @param distance  the distance between the two particles, measured in nm
+     * @return the index of the bond that was added
+     */
+    int addBond(int particle1, int particle2, double distance);
+
+    /** 
+     * Get the parameters defining a bond
+     * 
+     * @param index     the index of the bond for which to get parameters
+     * @param particle1 the index of the first particle involved in the bond
+     * @param particle2 the index of the second particle involved in the bond
+     * @param distance  the distance between the two particles, measured in nm
+     */
+    void getBondParameters(int index, int& particle1, int& particle2, double& distance) const;
+    /**
+     * Set 1-2 bonds
+     * 
+     * @param index          index of the bond for which to set parameters
+     * @param particle1      index of first atom in bond
+     * @param particle2      index of second atom in bond
+     * @param bondLength     bond length
+     */
+    void setBondParameters( int index, int particle1, int particle2, double bondLength);
+    /** 
+     * Get number of bonds
+     * 
+     * @return number of bonds
+     */
+    int getNumBonds( void ) const;
+
     /**
      * Get the dielectric constant for the solvent.
      */
@@ -131,20 +166,20 @@ public:
     void setSoluteDielectric(double dielectric) {
         soluteDielectric = dielectric;
     }
-    /**
+    /** 
      * Get the method used for handling long range nonbonded interactions.
      */
     NonbondedMethod getNonbondedMethod() const;
-    /**
+    /** 
      * Set the method used for handling long range nonbonded interactions.
      */
     void setNonbondedMethod(NonbondedMethod method);
-    /**
+    /** 
      * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
      * is NoCutoff, this value will have no effect.
      */
     double getCutoffDistance() const;
-    /**
+    /** 
      * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
      * is NoCutoff, this value will have no effect.
      */
@@ -155,7 +190,23 @@ private:
     class ParticleInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance, solventDielectric, soluteDielectric;
+    class BondInfo;
+
+// Retarded visual studio compiler complains about being unable to 
+// export private stl class members.
+// This stanza explains that it should temporarily shut up.
+#if defined(_MSC_VER)
+#pragma warning(push)
+#pragma warning(disable:4251)
+#endif
+
     std::vector<ParticleInfo> particles;
+    std::vector<BondInfo> bonds;
+
+#if defined(_MSC_VER)
+#pragma warning(pop)
+#endif
+
 };
 
 class GBVIForce::ParticleInfo {
@@ -169,6 +220,20 @@ public:
     }
 };
 
+class GBVIForce::BondInfo {
+public:
+    int particle1, particle2;
+    double bondLength;
+    BondInfo() {
+        bondLength     = 0.0;
+        particle1      = -1;
+        particle2      = -1;
+    }
+    BondInfo(int atomIndex1, int atomIndex2, double bondLength) :
+             particle1(atomIndex1), particle2(atomIndex2), bondLength(bondLength) {
+    }
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_GBVIFORCEFIELD_H_*/

openmmapi/include/openmm/internal/GBVIForceImpl.h

@@ -56,6 +56,9 @@ public:
     void findScaledRadii( int numberOfParticles, const std::vector<std::vector<int> >& bondIndices,
                           const std::vector<double> & bondLengths, std::vector<double> & scaledRadii) const;
 
+    // if bond info not set, then use bond forces/constraints
+    int getBondsFromForces(ContextImpl& context);
+
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }

openmmapi/src/GBVIForce.cpp

@@ -33,6 +33,7 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/GBVIForce.h"
 #include "openmm/internal/GBVIForceImpl.h"
+#include <sstream>
 
 using namespace OpenMM;
 
@@ -72,6 +73,27 @@ void GBVIForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
+int GBVIForce::addBond(int particle1, int particle2, double bondLength) {
+    bonds.push_back(BondInfo(particle1, particle2, bondLength));
+    return bonds.size()-1;
+}
+
+void GBVIForce::setBondParameters( int index, int particle1, int particle2, double bondLength) {
+    bonds[index].particle1  = particle1;
+    bonds[index].particle2  = particle2;
+    bonds[index].bondLength = bondLength;
+}
+
+int GBVIForce::getNumBonds( void ) const {
+   return (int) bonds.size();
+}
+
+void GBVIForce::getBondParameters(int index, int& bondIndex1, int& bondIndex2, double& bondLength) const {
+    bondIndex1 = bonds[index].particle1;
+    bondIndex2 = bonds[index].particle2;
+    bondLength = bonds[index].bondLength;
+}
+
 ForceImpl* GBVIForce::createImpl() {
     return new GBVIForceImpl(*this);
 }

openmmapi/src/GBVIForceImpl.cpp

@@ -36,6 +36,7 @@
 #include <vector>
 #include <cmath>
 #include <cstdio>
+#include <sstream>
 
 using namespace OpenMM;
 using std::vector;
@@ -47,52 +48,94 @@ void GBVIForceImpl::initialize(ContextImpl& context) {
     kernel = context.getPlatform().createKernel(CalcGBVIForceKernel::Name(), context);
     if (owner.getNumParticles() != context.getSystem().getNumParticles())
         throw OpenMMException("GBVIForce must have exactly as many particles as the System it belongs to.");
-
-    // load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
     
     System& system            = context.getSystem();
     int     numberOfParticles = owner.getNumParticles();
+    int numberOfBonds         = owner.getNumBonds();
+
+    // load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
+    // numberOfBonds < 1, indicating they were not set by the user
+
+    if( numberOfBonds < 1 ){
+        (void) fprintf( stderr, "Warning: no covalent bonds set for GB/VI force!\n" ); 
+//        getBondsFromForces( context );
+//        numberOfBonds = owner.getNumBonds();
+    }
+
+    std::vector< std::vector<int> > bondIndices;
+    bondIndices.resize( numberOfBonds );
+
+    std::vector<double> bondLengths;
+    bondLengths.resize( numberOfBonds );
+
+    for (int i = 0; i < numberOfBonds; i++) {
+        int particle1, particle2;
+        double bondLength;
+        owner.getBondParameters(i, particle1, particle2, bondLength);
+        if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
+            std::stringstream msg;
+            msg << "GBVISoftcoreForce: Illegal particle index: ";
+            msg << particle1;
+            throw OpenMMException(msg.str());
+        }
+        if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
+            std::stringstream msg;
+            msg << "GBVISoftcoreForce: Illegal particle index: ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        if (bondLength < 0 ) { 
+            std::stringstream msg;
+            msg << "GBVISoftcoreForce: negative bondlength: ";
+            msg << bondLength;
+            throw OpenMMException(msg.str());
+        }
+        bondIndices[i].push_back( particle1 );  
+        bondIndices[i].push_back( particle2 );
+        bondLengths[i] = bondLength;
+    }   
+
+    vector<double> scaledRadii;
+    scaledRadii.resize(numberOfParticles);
+    findScaledRadii( numberOfParticles, bondIndices, bondLengths, scaledRadii);
+
+    dynamic_cast<CalcGBVIForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner, scaledRadii);
+}
 
-    vector<vector<int> > bondIndices;
-    vector<double> bondLengths;
+int GBVIForceImpl::getBondsFromForces(ContextImpl& context) {
 
+    // load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
+    
+    System& system            = context.getSystem();
     for (int i = 0; i < system.getNumForces(); i++) {
         if (dynamic_cast<HarmonicBondForce*>(&system.getForce(i)) != NULL) {
             const HarmonicBondForce& force = dynamic_cast<const HarmonicBondForce&>(system.getForce(i));
-            bondIndices.resize(force.getNumBonds());
             for (int j = 0; j < force.getNumBonds(); ++j) {
                 int particle1, particle2;
                 double length, k;
                 force.getBondParameters(j, particle1, particle2, length, k); 
-                bondIndices[j].push_back(particle1);
-                bondIndices[j].push_back(particle2);
-                bondLengths.push_back(length);
+                owner.addBond( particle1, particle2, length );
             }
             break;
         }
     }   
 
-    // Also treat constrained distances as bonds.
+    // Also treat constrained distances as bonds if mass of one particle is < (2 + epsilon) (~2=deuterium) 
 
-    int numBonds = bondIndices.size();
-    bondIndices.resize(numBonds+system.getNumConstraints());
     for (int j = 0; j < system.getNumConstraints(); j++) {
         int particle1, particle2;
         double distance;
         system.getConstraintParameters(j, particle1, particle2, distance);
-        bondIndices[j+numBonds].push_back(particle1);
-        bondIndices[j+numBonds].push_back(particle2);
-        bondLengths.push_back(distance);
+        double mass1 = system.getParticleMass( particle1 );
+        double mass2 = system.getParticleMass( particle2 );
+        if( mass1 < 2.1 || mass2 < 2.1 ){
+            owner.addBond( particle1, particle2, distance );
+        }
     }   
 
-    vector<double> scaledRadii;
-    scaledRadii.resize(numberOfParticles);
-    findScaledRadii( numberOfParticles, bondIndices, bondLengths, scaledRadii);
-
-    dynamic_cast<CalcGBVIForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner, scaledRadii);
 }
 
-#define GBVIDebug 1
+#define GBVIDebug 0
 
 void GBVIForceImpl::findScaledRadii( int numberOfParticles, const std::vector<std::vector<int> >& bondIndices,
                                      const std::vector<double> & bondLengths, std::vector<double> & scaledRadii) const {

platforms/cuda/src/CudaKernelFactory.cpp

@@ -51,6 +51,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomNonbondedForceKernel(name, platform, data, context.getSystem());
     if (name == CalcGBSAOBCForceKernel::Name())
         return new CudaCalcGBSAOBCForceKernel(name, platform, data);
+    if (name == CalcGBVIForceKernel::Name())
+        return new CudaCalcGBVIForceKernel(name, platform, data);
     if (name == IntegrateVerletStepKernel::Name())
         return new CudaIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -51,11 +51,16 @@ void CudaCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context)
 
 void CudaCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
     _gpuContext* gpu = data.gpu;
-    if (gpu->bIncludeGBSA) {
+
+    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI) {
         gpu->bRecalculateBornRadii = true;
         kCalculateCDLJObcGbsaForces1(gpu);
         kReduceObcGbsaBornForces(gpu);
-        kCalculateObcGbsaForces2(gpu);
+        if (gpu->bIncludeGBSA ) { 
+           kCalculateObcGbsaForces2(gpu);
+        } else {
+           kCalculateGBVIForces2(gpu);
+        }
     }
     else if (data.hasNonbonded)
         kCalculateCDLJForces(gpu);
@@ -75,19 +80,21 @@ void CudaCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context)
 
 double CudaCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
     _gpuContext* gpu = data.gpu;
-    if (gpu->bIncludeGBSA) {
+    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI) {
         gpu->bRecalculateBornRadii = true;
         kCalculateCDLJObcGbsaForces1(gpu);
         kReduceObcGbsaBornForces(gpu);
-        kCalculateObcGbsaForces2(gpu);
+        if (gpu->bIncludeGBSA ) {
+           kCalculateObcGbsaForces2(gpu);
+        } else {
+           kCalculateGBVIForces2(gpu);
+        }
     }
     else if (data.hasNonbonded)
         kCalculateCDLJForces(gpu);
     if (data.hasCustomNonbonded)
         kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
     kCalculateLocalForces(gpu);
-    if (gpu->bIncludeGBSA)
-        kReduceBornSumAndForces(gpu);
     return kReduceEnergy(gpu)+data.ewaldSelfEnergy;
 }
 
@@ -533,6 +540,38 @@ void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCF
 void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
 }
 
+CudaCalcGBVIForceKernel::~CudaCalcGBVIForceKernel() {
+}
+
+void CudaCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii) {
+
+    int numParticles = system.getNumParticles();
+    _gpuContext* gpu = data.gpu;
+
+    vector<int> particle(numParticles);
+    vector<float> radius(numParticles);
+    vector<float> scaledRadii(numParticles);
+    vector<float> gammas(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        double charge, particleRadius, gamma, bornRadiusScaleFactor;
+        force.getParticleParameters(i, charge, particleRadius, gamma );
+        particle[i]                  = i;
+        radius[i]                    = (float) particleRadius;
+        gammas[i]                    = (float) gamma;
+        scaledRadii[i]               = (float) inputScaledRadii[i];
+    }
+    gpuSetGBVIParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), particle,
+                         radius, gammas, scaledRadii );
+}
+
+void CudaCalcGBVIForceKernel::executeForces(ContextImpl& context) {
+}
+
+double CudaCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
+    return 0.0;
+}
+
 static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
 
     // Initialize any terms that haven't already been handled by a Force.
@@ -643,7 +682,6 @@ void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langev
     double stepSize = integrator.getStepSize();
     if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
         // Initialize the GPU parameters.
-        
         double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
         gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
         gpuSetConstants(gpu);

platforms/cuda/src/CudaKernels.h

@@ -389,6 +389,39 @@ private:
     CudaPlatform::PlatformData& data;
 };
 
+/**
+ * This kernel is invoked by GBVIForce to calculate the forces acting on the system.
+ */
+class CudaCalcGBVIForceKernel : public CalcGBVIForceKernel {
+public:
+    CudaCalcGBVIForceKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : CalcGBVIForceKernel(name, platform), data(data) {
+    }
+    ~CudaCalcGBVIForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system      the System this kernel will be applied to
+     * @param force       the GBVIForce this kernel will be used for
+     * @param scaledRadii the scaled radii (Eq. 5 of Labute paper)
+     */
+    void initialize(const System& system, const GBVIForce& force, const std::vector<double> & scaledRadii);
+    /**
+     * Execute the kernel to calculate the forces.
+     * 
+     * @param context    the context in which to execute this kernel
+     */
+    void executeForces(ContextImpl& context);
+    /**
+     * Execute the kernel to calculate the energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the GBVIForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    CudaPlatform::PlatformData& data;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/cuda/src/kernels/cudaKernels.h

@@ -30,16 +30,22 @@
 extern void kClearForces(gpuContext gpu);
 extern void kClearEnergy(gpuContext gpu);
 extern void kClearBornForces(gpuContext gpu);
+extern void kClearObcGbsaBornSum(gpuContext gpu);
 extern void kCalculateObcGbsaBornSum(gpuContext gpu);
 extern void kReduceObcGbsaBornSum(gpuContext gpu);
+extern void kCalculateGBVIBornSum(gpuContext gpu);
+extern void kReduceGBVIBornSum(gpuContext gpu);
+extern void kClearGBVIBornSum( gpuContext gpu );
 extern void kGenerateRandoms(gpuContext gpu);
 
 // Main loop
 extern void kCalculateCDLJObcGbsaForces1(gpuContext gpu);
+extern void kCalculateCDLJGBVIForces1(gpuContext gpu);
 extern void kCalculateCDLJForces(gpuContext gpu);
 extern void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid);
 extern void kReduceObcGbsaBornForces(gpuContext gpu);
 extern void kCalculateObcGbsaForces2(gpuContext gpu);
+extern void kCalculateGBVIForces2(gpuContext gpu);
 extern void kCalculateLocalForces(gpuContext gpu);
 extern void kCalculateAndersenThermostat(gpuContext gpu);
 extern void kReduceBornSumAndForces(gpuContext gpu);
@@ -73,8 +79,12 @@ extern void SetCalculateLocalForcesSim(gpuContext gpu);
 extern void GetCalculateLocalForcesSim(gpuContext gpu);
 extern void SetCalculateObcGbsaBornSumSim(gpuContext gpu);
 extern void GetCalculateObcGbsaBornSumSim(gpuContext gpu);
+extern void SetCalculateGBVIBornSumSim(gpuContext gpu);
+extern void GetCalculateGBVIBornSumSim(gpuContext gpu);
 extern void SetCalculateObcGbsaForces2Sim(gpuContext gpu);
 extern void GetCalculateObcGbsaForces2Sim(gpuContext gpu);
+extern void SetCalculateGBVIForces2Sim(gpuContext gpu);
+extern void GetCalculateGBVIForces2Sim(gpuContext gpu);
 extern void SetCalculateAndersenThermostatSim(gpuContext gpu);
 extern void GetCalculateAndersenThermostatSim(gpuContext gpu);
 extern void SetCalculatePMESim(gpuContext gpu);

platforms/cuda/src/kernels/gpu.cpp

@@ -524,7 +524,8 @@ void gpuSetLJ14Parameters(gpuContext gpu, float epsfac, float fudge, const vecto
     psLJ14Parameter->Upload();
 }
 
-static void setExclusions(gpuContext gpu, const vector<vector<int> >& exclusions) {
+extern "C"
+void setExclusions(gpuContext gpu, const vector<vector<int> >& exclusions) {
     if (gpu->exclusions.size() > 0) {
         bool ok = (exclusions.size() == gpu->exclusions.size());
         for (int i = 0; i < exclusions.size() && ok; i++) {
@@ -569,7 +570,8 @@ void gpuSetCoulombParameters(gpuContext gpu, float epsfac, const vector<int>& at
     }
 
     // Dummy out extra atom data
-    for (unsigned int i = coulombs; i < gpu->sim.paddedNumberOfAtoms; i++)
+
+    for (unsigned int i = gpu->natoms; i < gpu->sim.paddedNumberOfAtoms; i++)
     {
         (*gpu->psPosq4)[i].x       = 100000.0f + i * 10.0f;
         (*gpu->psPosq4)[i].y       = 100000.0f + i * 10.0f;
@@ -578,7 +580,6 @@ void gpuSetCoulombParameters(gpuContext gpu, float epsfac, const vector<int>& at
         (*gpu->psSigEps2)[i].x     = 0.0f;
         (*gpu->psSigEps2)[i].y     = 0.0f;
     }
-
     gpu->psPosq4->Upload();
     gpu->psSigEps2->Upload();
 }
@@ -916,6 +917,53 @@ void gpuSetObcParameters(gpuContext gpu, float innerDielectric, float solventDie
     gpu->sim.preFactor = 2.0f*electricConstant*((1.0f/innerDielectric)-(1.0f/solventDielectric))*gpu->sim.forceConversionFactor;
 }
 
+extern "C"
+void gpuSetGBVIParameters(gpuContext gpu, float innerDielectric, float solventDielectric, const vector<int>& atom, const vector<float>& radius, 
+                          const vector<float>& gamma, const vector<float>& scaledRadii )
+{
+    unsigned int atoms = atom.size();
+    gpu->bIncludeGBVI  = true;
+    double tau         = ((1.0f/innerDielectric)-(1.0f/solventDielectric)); 
+    for (unsigned int i = 0; i < atoms; i++)
+    {
+            (*gpu->psGBVIData)[i].x = radius[i];
+            (*gpu->psGBVIData)[i].y = scaledRadii[i];
+            (*gpu->psGBVIData)[i].z = tau*gamma[i];
+            (*gpu->psGBVIData)[i].w = 1.0f;
+
+(*gpu->psObcData)[i].x  = radius[i];
+(*gpu->psObcData)[i].y  = 0.9f*radius[i];
+
+#undef DUMP_PARAMETERS
+#define DUMP_PARAMETERS 0
+#if (DUMP_PARAMETERS == 1)
+        (void) fprintf( stderr,"GBVI param: %5u R=%14.7e scaledR=%14.7e gamma*tau=%14.7e bornRadiusScaleFactor=%14.7e\n",
+                        i, (*gpu->psGBVIData)[i].x, (*gpu->psGBVIData)[i].y,
+                        (*gpu->psGBVIData)[i].z, (*gpu->psGBVIData)[i].w ); 
+#endif
+    }
+//(void) fprintf( stderr, "gpuSetGBVIParameters: setting Obc parameters!!!! should be removed.\n" );
+    // Dummy out extra atom data
+    for (unsigned int i = atoms; i < gpu->sim.paddedNumberOfAtoms; i++)
+    {
+        (*gpu->psBornRadii)[i]      = 0.2f;
+        (*gpu->psGBVIData)[i].x     = 0.01f;
+        (*gpu->psGBVIData)[i].y     = 0.01f;
+        (*gpu->psGBVIData)[i].z     = 0.01f;
+        (*gpu->psGBVIData)[i].w     = 1.00f;
+    }
+
+    gpu->psBornRadii->Upload();
+    gpu->psGBVIData->Upload();
+gpu->psObcData->Upload();
+    gpu->sim.preFactor = 2.0f*electricConstant*((1.0f/innerDielectric)-(1.0f/solventDielectric))*gpu->sim.forceConversionFactor;
+
+#if (DUMP_PARAMETERS == 1)
+(void) fprintf( stderr, "gpuSetGBVIParameters: preFactor=%14.6e elecCnstnt=%.4f frcCnvrsnFctr=%.4f tau=%.4f.\n",
+                gpu->sim.preFactor, 2.0f*electricConstant, gpu->sim.forceConversionFactor, ((1.0f/innerDielectric)-(1.0f/solventDielectric)) );
+#endif
+}
+
 static void markShakeClusterInvalid(ShakeCluster& cluster, map<int, ShakeCluster>& allClusters, vector<bool>& invalidForShake)
 {
     cluster.valid = false;
@@ -1435,6 +1483,8 @@ int gpuAllocateInitialBuffers(gpuContext gpu)
     gpu->sim.pAttr                      = gpu->psSigEps2->_pDevStream[0];
     gpu->psObcData                      = new CUDAStream<float2>(gpu->sim.paddedNumberOfAtoms, 1, "ObcData");
     gpu->sim.pObcData                   = gpu->psObcData->_pDevStream[0];
+    gpu->psGBVIData                     = new CUDAStream<float4>(gpu->sim.paddedNumberOfAtoms, 1, "GBVIData");
+    gpu->sim.pGBVIData                  = gpu->psGBVIData->_pDevStream[0];
     gpu->psStepSize                     = new CUDAStream<float2>(1, 1, "StepSize");
     gpu->sim.pStepSize                  = gpu->psStepSize->_pDevStream[0];
     (*gpu->psStepSize)[0] = make_float2(0.0f, 0.0f);
@@ -1671,6 +1721,7 @@ void* gpuInit(int numAtoms, unsigned int device, bool useBlockingSync)
     gpu->bRemoveCM                  = false;
     gpu->bRecalculateBornRadii      = true;
     gpu->bIncludeGBSA               = false;
+    gpu->bIncludeGBVI               = false;
     gpuInitializeRandoms(gpu);
 
     // To be determined later
@@ -1883,6 +1934,7 @@ void gpuShutDown(gpuContext gpu)
     if (gpu->psTabulatedErfc != NULL)
         delete gpu->psTabulatedErfc;
     delete gpu->psObcData;
+    delete gpu->psGBVIData;
     delete gpu->psObcChain;
     delete gpu->psBornForce;
     delete gpu->psBornRadii;
@@ -1957,10 +2009,10 @@ int gpuBuildOutputBuffers(gpuContext gpu)
         gpu->bOutputBufferPerWarp           = false;
         gpu->sim.nonbondOutputBuffers       = gpu->sim.paddedNumberOfAtoms / GRID;
     }
-    gpu->sim.totalNonbondOutputBuffers  = (gpu->bIncludeGBSA ? 2 * gpu->sim.nonbondOutputBuffers : gpu->sim.nonbondOutputBuffers);
+    gpu->sim.totalNonbondOutputBuffers  = ( (gpu->bIncludeGBSA || gpu->bIncludeGBVI) ? 2 * gpu->sim.nonbondOutputBuffers : gpu->sim.nonbondOutputBuffers);
+//gpu->sim.totalNonbondOutputBuffers  = 2*gpu->sim.nonbondOutputBuffers;
     gpu->sim.outputBuffers              = gpu->sim.totalNonbondOutputBuffers;
 
-
     unsigned int outputBuffers = gpu->sim.totalNonbondOutputBuffers;
     for (unsigned int i = 0; i < gpu->sim.paddedNumberOfAtoms; i++)
     {
@@ -2234,7 +2286,9 @@ int gpuSetConstants(gpuContext gpu)
     SetCalculateCustomNonbondedForcesSim(gpu);
     SetCalculateLocalForcesSim(gpu);
     SetCalculateObcGbsaBornSumSim(gpu);
+    SetCalculateGBVIBornSumSim(gpu);
     SetCalculateObcGbsaForces2Sim(gpu);
+    SetCalculateGBVIForces2Sim(gpu);
     SetCalculateAndersenThermostatSim(gpu);
     SetCalculatePMESim(gpu);
     SetForcesSim(gpu);

platforms/cuda/src/kernels/gputypes.h

@@ -62,7 +62,6 @@ enum SM_VERSION
  * to gromacs functions*/
 struct _gpuContext {
     
-    
     //Cache this here so that it doesn't
     //have to be repeatedly passed around
     int natoms;
@@ -87,6 +86,7 @@ struct _gpuContext {
     bool bRecalculateBornRadii;
     bool bOutputBufferPerWarp;
     bool bIncludeGBSA;
+    bool bIncludeGBVI;
     bool tabulatedFunctionsChanged;
     unsigned long seed;
     SM_VERSION sm_version;
@@ -114,6 +114,7 @@ struct _gpuContext {
     CUDAStream<int>* psPmeAtomRange;           // The range of sorted atoms at each grid point
     CUDAStream<float2>* psPmeAtomGridIndex;    // The grid point each atom is at
     CUDAStream<float2>* psObcData;
+    CUDAStream<float4>* psGBVIData;
     CUDAStream<float>* psObcChain;
     CUDAStream<float>* psBornForce;
     CUDAStream<float>* psBornRadii;
@@ -223,6 +224,10 @@ void gpuSetPeriodicBoxSize(gpuContext gpu, float xsize, float ysize, float zsize
 extern "C"
 void gpuSetObcParameters(gpuContext gpu, float innerDielectric, float solventDielectric, const std::vector<float>& radius, const std::vector<float>& scale, const std::vector<float>& charge);
 
+extern "C" 
+void gpuSetGBVIParameters(gpuContext gpu, float innerDielectric, float solventDielectric, const std::vector<int>& atom, const std::vector<float>& radius,
+                          const std::vector<float>& gammas, const std::vector<float>& scaledRadii);
+
 extern "C"
 void gpuSetConstraintParameters(gpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<float>& distance,
         const std::vector<float>& invMass1, const std::vector<float>& invMass2, float constraintTolerance);
@@ -275,4 +280,7 @@ int gpuSetConstants(gpuContext gpu);
 extern "C"
 void gpuReorderAtoms(gpuContext gpu);
 
+extern "C"
+void setExclusions(gpuContext gpu, const std::vector<std::vector<int> >& exclusions);
+
 #endif //__GPUTYPES_H__

platforms/cuda/src/kernels/kCalculateCDLJObcGbsaForces1.cu

@@ -55,6 +55,7 @@ void SetCalculateCDLJObcGbsaForces1Sim(gpuContext gpu)
     cudaError_t status;
     status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));     
     RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
+
 }
 
 void GetCalculateCDLJObcGbsaForces1Sim(gpuContext gpu)
@@ -142,15 +143,22 @@ void kCalculateCDLJObcGbsaForces1(gpuContext gpu)
         case NO_CUTOFF:
             if (gpu->bRecalculateBornRadii)
             {
-                kCalculateObcGbsaBornSum(gpu);
-                kReduceObcGbsaBornSum(gpu);
+                if( gpu->bIncludeGBVI ){
+                   kCalculateGBVIBornSum(gpu);
+                   kReduceGBVIBornSum(gpu);
+                } else {
+                   kCalculateObcGbsaBornSum(gpu);
+                   kReduceObcGbsaBornSum(gpu);
+                }
+
             }
             if (gpu->bOutputBufferPerWarp)
-                kCalculateCDLJObcGbsaN2ByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
-                        sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
-            else
-                kCalculateCDLJObcGbsaN2Forces1_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
-                        sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
+                   kCalculateCDLJObcGbsaN2ByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
+                           sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
+               else
+                   kCalculateCDLJObcGbsaN2Forces1_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
+                           sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
+   
             LAUNCHERROR("kCalculateCDLJObcGbsaN2Forces1");
             break;
         case CUTOFF:

platforms/cuda/src/kernels/kCalculateObcGbsaBornSum.cu

@@ -45,7 +45,6 @@ struct Atom {
     float r;
     float sr;
     float sum;
-    float junk;
 };
 
 static __constant__ cudaGmxSimulation cSim;
@@ -150,14 +149,21 @@ void kReduceObcGbsaBornSum(gpuContext gpu)
     LAUNCHERROR("kReduceObcGbsaBornSum");
 }
 
+void kClearObcGbsaBornSum(gpuContext gpu)
+{
+  //  printf("kClearObcGbsaBornSum\n");
+    kClearObcGbsaBornSum_kernel<<<gpu->sim.blocks, 384>>>();
+}
+
 void kCalculateObcGbsaBornSum(gpuContext gpu)
 {
   //  printf("kCalculateObcgbsaBornSum\n");
-    kClearObcGbsaBornSum_kernel<<<gpu->sim.blocks, 384>>>();
+    kClearObcGbsaBornSum(gpu);
     LAUNCHERROR("kClearBornSum");
     switch (gpu->sim.nonbondedMethod)
     {
         case NO_CUTOFF:
+
             if (gpu->bOutputBufferPerWarp)
                 kCalculateObcGbsaN2ByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                         sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);

platforms/cuda/src/kernels/kForces.cu

@@ -43,14 +43,14 @@ void SetForcesSim(gpuContext gpu)
 {
     cudaError_t status;
     status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));     
-    RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
+    RTERROR(status, "cudaMemcpyToSymbol: SetForcesSim copy to cSim failed");
 }
 
 void GetForcesSim(gpuContext gpu)
 {
     cudaError_t status;
     status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));     
-    RTERROR(status, "cudaMemcpyFromSymbol: SetSim copy from cSim failed");
+    RTERROR(status, "cudaMemcpyFromSymbol: GetForcesSim copy from cSim failed");
 }
 
 __global__ void kClearForces_kernel()
@@ -243,7 +243,7 @@ __global__ void kReduceForces_kernel()
         {
             totalForce += *pFt;
         }
-        
+
         pFt = (float*)cSim.pForce4 + pos;
         *pFt = totalForce;
         pos += gridDim.x * blockDim.x;
@@ -259,11 +259,13 @@ void kReduceForces(gpuContext gpu)
 
 double kReduceEnergy(gpuContext gpu)
 {
-//    printf("kReduceEnergy\n");
+    //printf("kReduceEnergy\n");
     gpu->psEnergy->Download();
-    double sum = 0.0f;
-    for (int i = 0; i < gpu->sim.energyOutputBuffers; i++)
+    double sum = 0.0;
+    for (int i = 0; i < gpu->sim.energyOutputBuffers; i++){
         sum += (*gpu->psEnergy)[i];
+    }
+
     return sum;
 }
 
@@ -308,29 +310,87 @@ __global__ void kReduceObcGbsaBornForces_kernel()
         }
         float r            = (obcData.x + cSim.dielectricOffset + cSim.probeRadius);
         float ratio6       = pow((obcData.x + cSim.dielectricOffset) / bornRadius, 6.0f);
-        //float saTerm       = cSim.surfaceAreaFactor * r * r * ratio6;
         float saTerm       = cSim.surfaceAreaFactor * r * r * ratio6;
-        totalForce        += saTerm / bornRadius; // 1.102 == Temp mysterious fudge factor, FIX FIX FIX
 
-        /* E */
-        energy += saTerm;
+        totalForce        += saTerm / bornRadius;
+        totalForce        *= bornRadius * bornRadius * obcChain;
+
+        energy            += saTerm;
+
+        pFt                = cSim.pBornForce + pos;
+        *pFt               = totalForce;
+
+        pos               += gridDim.x * blockDim.x;
+    }
+
+    // correct for surface area factor of -6
+    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy / -6.0f;
+}
+
+__global__ void kReduceGBVIBornForces_kernel()
+{
+    unsigned int pos = (blockIdx.x * blockDim.x + threadIdx.x);
+    float energy = 0.0f;
+    while (pos < cSim.atoms)
+    {
+        float bornRadius  = cSim.pBornRadii[pos];
+        float4 gbviData   = cSim.pGBVIData[pos];
+        float totalForce  = 0.0f;
+        float* pFt        = cSim.pBornForce + pos;
+
+        int i = cSim.nonbondOutputBuffers;
+        while (i >= 4)
+        {
+            float f1    = *pFt;
+            pFt        += cSim.stride;
+            float f2    = *pFt;
+            pFt        += cSim.stride;
+            float f3    = *pFt;
+            pFt        += cSim.stride;
+            float f4    = *pFt;
+            pFt        += cSim.stride;
+            totalForce += f1 + f2 + f3 + f4;
+            i -= 4;
+        }
+        if (i >= 2)
+        {
+            float f1    = *pFt;
+            pFt        += cSim.stride;
+            float f2    = *pFt;
+            pFt        += cSim.stride;
+            totalForce += f1 + f2;
+            i -= 2;
+        }
+        if (i > 0)
+        {
+            totalForce += *pFt;
+        }
 
-        totalForce *= bornRadius * bornRadius * obcChain;
+        float ratio         = (gbviData.x/bornRadius);
+        float ratio3        = ratio*ratio*ratio;
+        energy             -= gbviData.z*ratio3;
+        totalForce         += (3.0f*gbviData.z*ratio3)/bornRadius; // 'cavity' term
+        float br2           = bornRadius*bornRadius;
+        totalForce         *= (1.0f/3.0f)*br2*br2;
 
         pFt = cSim.pBornForce + pos;
         *pFt = totalForce;
         pos += gridDim.x * blockDim.x;
     }
-    /* E */
-    // correct for surface area factor of -6
-    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy / -6.0f;
+    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy;
 }
 
 void kReduceObcGbsaBornForces(gpuContext gpu)
 {
     //printf("kReduceObcGbsaBornForces\n");
-    kReduceObcGbsaBornForces_kernel<<<gpu->sim.blocks, gpu->sim.bf_reduce_threads_per_block>>>();
-    LAUNCHERROR("kReduceObcGbsaBornForces");
+    if( gpu->bIncludeGBSA ){
+       kReduceObcGbsaBornForces_kernel<<<gpu->sim.blocks, gpu->sim.bf_reduce_threads_per_block>>>();
+       LAUNCHERROR("kReduceObcGbsaBornForces");
+    } else if( gpu->bIncludeGBVI ){
+       kReduceGBVIBornForces_kernel<<<gpu->sim.blocks, gpu->sim.bf_reduce_threads_per_block>>>();
+       LAUNCHERROR("kReduceGBVIBornForces");
+    }   
+
 }
 
 

platforms/reference/src/SimTKReference/ReferenceBondIxn.cpp

@@ -155,7 +155,7 @@ RealOpenMM ReferenceBondIxn::getNormedDotProduct( RealOpenMM* vector1, RealOpenM
 // for angles near pi, double is required due to the 'steepness' of acos()
 // in this regime.
   
-#define USE_DOUBLE_FOR_NORMED_DOT_PRODUCT
+//#define USE_DOUBLE_FOR_NORMED_DOT_PRODUCT
 
 #if defined USE_DOUBLE_FOR_NORMED_DOT_PRODUCT
    double v1D[3];

platforms/reference/src/SimTKReference/ReferenceRbDihedralBond.cpp

@@ -109,7 +109,7 @@ int ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
 
    // debug flag
 
-   static const int debug              = 0;
+   static int debug                    = 0;
 
    static const int LastAtomIndex      = 4;
 
@@ -244,9 +244,12 @@ int ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
       }
       message << std::endl;
 
-      message << " k="     << parameters[0];
-      message << " a="     << parameters[1];
-      message << " m="     << parameters[2];
+      message << std::endl << numberOfParameters << " Parameters: [";
+      for( int ii = 0; ii < numberOfParameters; ii++ ){
+         message << parameters[ii] << " ";
+      }
+      message << " ]" << std::endl;
+
       message << " ang="   << dihederalAngle;
       message << " dotD="  << cosPhi;
       message << " sign="  << signOfAngle;

platforms/reference/src/SimTKUtilities/TwoDimArraySimTk.h

@@ -27,7 +27,8 @@
 
 // ---------------------------------------------------------------------------------------
 
-#include "UtilitiesSimTk.h"
+//#include "SimTKOpenMMUtilities.h"
+#include "SimTKOpenMMUtilities.h"
 
 /**---------------------------------------------------------------------------------------
 
@@ -116,11 +117,11 @@ TwoDimArraySimTk<T>::TwoDimArraySimTk( int rowSize, int columnSize ){
    _columnSize        = columnSize;
 
 #ifdef useXMalloc 
-   _2Darray           = (T**) UtilitiesSimTk::Xmalloc( "TwoDimArraySimTk", __FILE__, __LINE__, rowSize*sizeof( T* ) );
-   _2DMemoryBlock     = (T*)  UtilitiesSimTk::Xmalloc( "block",   __FILE__, __LINE__, columnSize*rowSize*sizeof( T ) );
+   _2Darray           = (T**) SimTKOpenMMUtilities::Xmalloc( "TwoDimArraySimTk", __FILE__, __LINE__, rowSize*sizeof( T* ) );
+   _2DMemoryBlock     = (T*)  SimTKOpenMMUtilities::Xmalloc( "block",   __FILE__, __LINE__, columnSize*rowSize*sizeof( T ) );
 #else
-   _2Darray           = (T**) UtilitiesSimTk::Xmalloc( "TwoDimArraySimTk", __FILE__, __LINE__, rowSize*sizeof( T* ) );
-   _2DMemoryBlock     = (T*)  UtilitiesSimTk::Xmalloc( "block",   __FILE__, __LINE__, columnSize*rowSize*sizeof( T ) );
+   _2Darray           = (T**) SimTKOpenMMUtilities::Xmalloc( "TwoDimArraySimTk", __FILE__, __LINE__, rowSize*sizeof( T* ) );
+   _2DMemoryBlock     = (T*)  SimTKOpenMMUtilities::Xmalloc( "block",   __FILE__, __LINE__, columnSize*rowSize*sizeof( T ) );
 #endif
 
    T* blockPtr        = _2DMemoryBlock;
@@ -157,10 +158,10 @@ template <typename T> TwoDimArraySimTk<T>::~TwoDimArraySimTk( ){
    // ---------------------------------------------------------------------------------------
 
    if( _2DMemoryBlock ){
-      UtilitiesSimTk::Xfree(  "2DMemoryBlock", __FILE__, __LINE__, _2DMemoryBlock );
+      SimTKOpenMMUtilities::Xfree(  "2DMemoryBlock", __FILE__, __LINE__, _2DMemoryBlock );
    }
    if( _2Darray ){
-      UtilitiesSimTk::Xfree( "2Darray", __FILE__, __LINE__, _2Darray );
+      SimTKOpenMMUtilities::Xfree( "2Darray", __FILE__, __LINE__, _2Darray );
    }
       
 }

platforms/reference/src/gbsa/CpuGBVI.cpp

@@ -163,14 +163,14 @@ int CpuGBVI::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRad
 
    // ---------------------------------------------------------------------------------------
 
-   GBVIParameters* gbviParameters        = getGBVIParameters();
-   int numberOfAtoms                     = gbviParameters->getNumberOfAtoms();
-   RealOpenMM* atomicRadii               = gbviParameters->getAtomicRadii();
-   const RealOpenMM* scaledRadii         = gbviParameters->getScaledRadii();
+   GBVIParameters* gbviParameters                   = getGBVIParameters();
+   int numberOfAtoms                                = gbviParameters->getNumberOfAtoms();
+   RealOpenMM* atomicRadii                          = gbviParameters->getAtomicRadii();
+   const RealOpenMM* scaledRadii                    = gbviParameters->getScaledRadii();
 
    // ---------------------------------------------------------------------------------------
 
-#if GBVIDebug
+#if( GBVIDebug == 1 )
    FILE* logFile                         = stderr;
    (void) fprintf( logFile, "\n%s\n", methodName );
 #endif
@@ -201,18 +201,24 @@ int CpuGBVI::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRad
 
             sum  += CpuGBVI::getVolume( r, radiusI, scaledRadii[atomJ] );
 
-#if GBVIDebug
-            (void) fprintf( logFile, "%d r=%.5f add for atom J=%d scR=%.4e %.6e sum=%14.6e\n",
-                            atomI, r, atomJ, scaledRadii[atomJ], getVolume( r, radiusI, scaledRadii[atomJ] ), sum );
+#if( GBVIDebug == 1 )
+if( atomI == 1568 || atomJ == 1568 ){
+            (void) fprintf( logFile, "%d addJ=%d scR=%14.6e %14.6e sum=%14.6e rI=%14.6e r=%14.6e S-R=%14.6e\n",
+                            atomI, atomJ, scaledRadii[atomJ], getVolume( r, radiusI, scaledRadii[atomJ] ), sum,
+                            radiusI, r, (scaledRadii[atomJ]-radiusI) );
+}
 #endif
          }
       }
 
+#if( GBVIDebug == 1 )
+      (void) fprintf( logFile, "%d Born radius sum=%14.6e %14.6e %14.6e ", atomI, sum, POW( radiusI, minusThree ), (POW( radiusI, minusThree ) - sum) );
+#endif
       sum              = POW( radiusI, minusThree ) - sum;
       bornRadii[atomI] = POW( sum, minusOneThird );
 
-#if GBVIDebug
-      (void) fprintf( logFile, "%d Born radius %14.6e\n", atomI, bornRadii[atomI] );
+#if( GBVIDebug == 1 )
+      (void) fprintf( logFile, "br=%14.6e\n", atomI, bornRadii[atomI] );
 #endif
 
    }
@@ -442,7 +448,7 @@ RealOpenMM CpuGBVI::computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM**
       bornRadii   = getBornRadii();
    }
 
-#if GBVIDebug
+#if( GBVIDebug == 1 )
    FILE* logFile                        = stderr;
    (void) fprintf( logFile, "\n%s\n", methodName );
    (void) fflush( logFile );
@@ -451,22 +457,24 @@ RealOpenMM CpuGBVI::computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM**
    // ---------------------------------------------------------------------------------------
 
    // Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
+   // to minimze roundoff error sum cavityEnergy separately since in general much
+   // smaller than other contributions
 
    RealOpenMM energy                 = zero;
+   RealOpenMM cavityEnergy           = zero;
 
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
  
-      RealOpenMM partialChargeI   = preFactor*partialCharges[atomI];
-      RealOpenMM partialChargeI2  = two*partialChargeI;
+      RealOpenMM partialChargeI   = partialCharges[atomI];
 
       // self-energy term
 
-      energy                      += partialChargeI*partialCharges[atomI]/bornRadii[atomI];
+      RealOpenMM  atomIEnergy     = half*partialChargeI/bornRadii[atomI];
 
       // cavity term
 
-      RealOpenMM ratio             = (atomicRadii[atomI]/bornRadii[atomI]);
-      energy                      -= gammaParameters[atomI]*ratio*ratio*ratio;
+      RealOpenMM ratio            = (atomicRadii[atomI]/bornRadii[atomI]);
+      cavityEnergy               += gammaParameters[atomI]*ratio*ratio*ratio;
 
 /*
 RealOpenMM e1 = partialChargeI*partialCharges[atomI]/bornRadii[atomI];
@@ -485,15 +493,19 @@ RealOpenMM e2 = gammaParameters[atomI]*ratio*ratio*ratio;
 
          RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
          RealOpenMM t                  = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);         
-         energy                       += partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
+         atomIEnergy                  += partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
 /*
 RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
 (void) fprintf( stderr, "E %d %d e3=%.4e r2=%4e t=%.3e sgb=%.4e e=%.5e\n", atomI, atomJ, e3, r2, t, Sgb( t ), energy );
 */
       }
+
+      energy += two*partialChargeI*atomIEnergy;
    }
+   energy *= CAL_TO_JOULE*preFactor;
+   energy -= cavityEnergy;
 
-#if GBVIDebug
+#if( GBVIDebug == 1 )
    (void) fprintf( logFile, "ElectricConstant=%.4e Tau=%.4e e=%.5e eOut=%.5e\n", preFactor, gbviParameters->getTau(), energy, gbviParameters->getTau()*energy );
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
       (void) fprintf( logFile, "bR %d bR=%16.8e\n", atomI, bornRadii[atomI] );
@@ -501,11 +513,15 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference
    (void) fflush( logFile );
 #endif
 
-   RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau());  
+   RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());  
    return (conversion*energy);
  
 }
 
+#undef GBVIDebug
+
+#define GBVIDebug 0
+
 /**---------------------------------------------------------------------------------------
 
    Get GB/VI forces
@@ -539,7 +555,7 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo
 
    // ---------------------------------------------------------------------------------------
 
-#if GBVIDebug
+#if( GBVIDebug == 1 )
    FILE* logFile                        = stderr;
    (void) fprintf( logFile, "\n%s\n", methodName );
    (void) fflush( logFile );
@@ -584,11 +600,6 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo
 
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
  
-      // partial of cavity term wrt Born radius
-
-      RealOpenMM          ratio = (atomicRadii[atomI]/bornRadii[atomI]);
-      bornForces[atomI]        += (three*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI]; 
-
       // partial of polar term wrt Born radius
       // and (dGpol/dr)(dr/dx)
 
@@ -641,16 +652,50 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo
 
          // 3 FLOP
 
+#if 0
+if( atomI == 0 ){
+   (void) fprintf( logFile, "bFCalc: %6d %6d %14.6e %14.6e   %14.6e %14.6e\n", atomI, atomJ, dGpol_dalpha2_ij, bornRadii[atomJ], bornForces[atomI], bornRadii[atomI] );
+}
+#endif
          bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
 
       }
    }
 
-#if GBVIDebug
+#if( GBVIDebug == 1 )
+{
+   double stupidFactor                   = three/CAL_TO_JOULE;
+   RealOpenMM conversion                 = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau());  
+   int maxPrint                          = 10;
+   const RealOpenMM* scaledRadii         = gbviParameters->getScaledRadii();
+
+   (void) fprintf( logFile, "Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, CAL_TO_JOULE, gbviParameters->getTau() );
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-       (void) fprintf( logFile, "F1 %d bF=%14.6e [%14.6e %14.6e %14.6e]\n", atomI, bornForces[atomI], forces[atomI][0], forces[atomI][1],  forces[atomI][2] );
+      RealOpenMM R        = atomicRadii[atomI];
+
+      // partial of cavity term wrt Born radius
+     
+      RealOpenMM  ratio   = (atomicRadii[atomI]/bornRadii[atomI]);
+      bornForces[atomI]  += (stupidFactor*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI]; 
+
+      RealOpenMM b2       = bornRadii[atomI]*bornRadii[atomI];
+      double xx           = bornForces[atomI]*oneThird*b2*b2;
+
+      // xx*conversion should agree w/ values pulled out of kReduceGBVIBornForces_kernel in kForces.cu
+
+      (void) fprintf( logFile, "F1 %6d r/sclR[%14.6e %14.6e] bR=%14.6e bF=%14.6e %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], bornForces[atomI], xx*conversion,
+//                      forces[atomI][0], forces[atomI][1], forces[atomI][2],
+                      conversion*forces[atomI][0], conversion*forces[atomI][1],  conversion*forces[atomI][2],
+                      atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
+      if( atomI == maxPrint ){
+         atomI = numberOfAtoms - maxPrint;
+         if( atomI < maxPrint )atomI = maxPrint;
+      }
    }
    (void) fflush( logFile );
+}
 #endif
 
    // ---------------------------------------------------------------------------------------
@@ -660,11 +705,28 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo
    // dGpol/dBornRadius) = bornForces[]
    // dBornRadius/dr     = (1/3)*(bR**4)*(dV/dr)
 
+#if 0
+   (void) fprintf( logFile, "Clearing forces before loop2 periodic=%d cutoff=%d cutoffR=%14.7e\n",
+                   _gbviParameters->getPeriodic(), _gbviParameters->getUseCutoff(),  _gbviParameters->getCutoffDistance() );
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+      forces[atomI][0]  = zero;
+      forces[atomI][1]  = zero;
+      forces[atomI][2]  = zero;
+   } 
+   (void) fflush( logFile );
+#endif
+
    const RealOpenMM* scaledRadii         = gbviParameters->getScaledRadii();
+   double stupidFactor                   = three/CAL_TO_JOULE;
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
  
       RealOpenMM R        = atomicRadii[atomI];
 
+      // partial of cavity term wrt Born radius
+     
+      RealOpenMM  ratio   = (atomicRadii[atomI]/bornRadii[atomI]);
+      bornForces[atomI]  += (stupidFactor*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI]; 
+
       RealOpenMM b2       = bornRadii[atomI]*bornRadii[atomI];
       bornForces[atomI]  *= oneThird*b2*b2;
 
@@ -698,22 +760,18 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo
 
             // find dRb/dr, where Rb is the Born radius
 
+            de = CpuGBVI::dL_dr( r, r+S, S ) + CpuGBVI::dL_dx( r, r+S, S );   
             if( FABS( diff ) < r ){
-               de = CpuGBVI::dL_dr( r, r+S, S ) + CpuGBVI::dL_dx( r, r+S, S );   
                if( R > (r - S) ){
-//(void) fprintf( stderr, "F2 %d %d block 0\n", atomI, atomJ );
                   de -= CpuGBVI::dL_dr( r, R, S );  
                } else {
-//(void) fprintf( stderr, "F2 %d %d block 1\n", atomI, atomJ );
                   de -= ( CpuGBVI::dL_dr( r, (r-S), S ) + CpuGBVI::dL_dx( r, (r-S), S ) );
                }
             } else if( r < (S - R) ){
-//(void) fprintf( stderr, "F3 %d %d block 2\n", atomI, atomJ );
-               de =   CpuGBVI::dL_dr( r, r+S, S ) + CpuGBVI::dL_dx( r, r+S, S ) -
-                    ( CpuGBVI::dL_dr( r, r-S, S ) + CpuGBVI::dL_dx( r, r-S, S ) );   
+               de -= ( CpuGBVI::dL_dr( r, r-S, S ) + CpuGBVI::dL_dx( r, r-S, S ) );   
             }
 
-if( 0 ){
+#if 0
    RealOpenMM delta = (RealOpenMM) 1.0e-02;
    (void) fprintf( stderr, "\n" );
    for( int kk = 0; kk < 5; kk++ ){
@@ -733,12 +791,10 @@ if( 0 ){
       (void) fprintf( stderr, "df %d %d [%14.6e %14.6e] V[%14.6e %14.6e] %.2e\n", atomI, atomJ, ded, df, V2, V1, deltaD );
       deltaD *= 0.1;
    }
-}
-
-
+#endif
 
 
-            // de = (dG/dRb)(dRb/dr)
+             // de = (dG/dRb)(dRb/dr)
 
             de                      *= bornForces[atomI]/r;
 
@@ -759,7 +815,7 @@ if( 0 ){
 
    }
 
-#if GBVIDebug
+#if( GBVIDebug == 1 )
    (void) fprintf( logFile, "Atom      BornRadii      BornForce                                         Forces\n" );
    double forceSum[3] = { 0.0, 0.0, 0.0 };
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
@@ -772,17 +828,33 @@ if( 0 ){
    (void) fflush( logFile );
 #endif
 
-#undef GBVIDebug
-
    // convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
 
    RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau());  
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-      inputForces[atomI][0] = conversion*forces[atomI][0];
-      inputForces[atomI][1] = conversion*forces[atomI][1];
-      inputForces[atomI][2] = conversion*forces[atomI][2];
+      inputForces[atomI][0] += conversion*forces[atomI][0];
+      inputForces[atomI][1] += conversion*forces[atomI][1];
+      inputForces[atomI][2] += conversion*forces[atomI][2];
    }
 
+#if( GBVIDebug == 1 )
+{
+   (void) fprintf( logFile, "\nPost conversion\n" );
+   (void) fprintf( logFile, "Atom      BornRadii      BornForce                                         Forces\n" );
+   int maxPrint = 10;
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+      (void) fprintf( logFile, "%4d %14.6e %14.6e [%14.6e %14.6e %14.6e]\n", atomI, bornRadii[atomI], conversion*bornForces[atomI],
+                      inputForces[atomI][0], inputForces[atomI][1], inputForces[atomI][2] );
+      if( atomI == maxPrint ){
+         atomI = numberOfAtoms - maxPrint;
+         if( atomI < maxPrint )atomI = numberOfAtoms;
+      }
+   }
+   (void) fflush( logFile );
+}
+#endif
+#undef GBVIDebug
+
    delete[] forces;
    delete[] block;