Reference code for GB/VI

olla/include/kernels.h

@@ -36,6 +36,7 @@
 #include "BrownianIntegrator.h"
 #include "CMMotionRemover.h"
 #include "GBSAOBCForce.h"
+#include "GBVIForce.h"
 #include "HarmonicAngleForce.h"
 #include "HarmonicBondForce.h"
 #include "KernelImpl.h"
@@ -277,6 +278,39 @@ public:
     virtual double executeEnergy(OpenMMContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by GBVIForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcGBVIForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcGBVIForces";
+    }
+    CalcGBVIForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * 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 scaled radii
+     */
+    virtual void initialize(const System& system, const GBVIForce& force, const std::vector<double>& scaledRadii) = 0;
+    /**
+     * Execute the kernel to calculate the forces.
+     * 
+     * @param context    the context in which to execute this kernel
+     */
+    virtual void executeForces(OpenMMContextImpl& context) = 0;
+    /**
+     * 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
+     */
+    virtual double executeEnergy(OpenMMContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/GBVIForce.h

+#ifndef OPENMM_GBVIFORCEFIELD_H_
+#define OPENMM_GBVIFORCEFIELD_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class implements an implicit solvation force using the GB/VI model.
+ * <p>
+ * If the System also contains a NonbondedForce, this force will use the cutoffs
+ * and periodic boundary conditions specified in it.
+ */
+
+class OPENMM_EXPORT GBVIForce : public Force {
+public:
+    /*
+     * Create a GBVIForce.
+     * 
+     * @param numParticles    the number of particles in the system
+     */
+    GBVIForce(int numParticles);
+    /**
+     * Get the number of particles in the system.
+     */
+    int getNumParticles() const {
+        return particles.size();
+    }
+    /**
+     * Get the force field parameters for a particle.
+     * 
+     * @param index          the index of the particle for which to get parameters
+     * @param charge         the charge of the particle, measured in units of the proton charge
+     * @param radius         the GBSA radius of the particle, measured in nm
+     * @param gamma          the gamma parameter
+     */
+    void getParticleParameters(int index, double& charge, double& radius, double& gamma) const;
+    /**
+     * Set the force field parameters for a particle.
+     * 
+     * @param index          the index of the particle for which to set parameters
+     * @param charge         the charge of the particle, measured in units of the proton charge
+     * @param radius         the GB/VI radius of the particle, measured in nm
+     * @param gamma          the gamma parameter
+     */
+    void setParticleParameters(int index, double charge, double radius, double gamma);
+    /**
+     * Get the dielectric constant for the solvent.
+     */
+    double getSolventDielectric() const {
+        return solventDielectric;
+    }
+    /**
+     * Set the dielectric constant for the solvent.
+     */
+    void setSolventDielectric(double dielectric) {
+        solventDielectric = dielectric;
+    }
+    /**
+     * Get the dielectric constant for the solute.
+     */
+    double getSoluteDielectric() const {
+        return soluteDielectric;
+    }
+    /**
+     * Set the dielectric constant for the solute.
+     */
+    void setSoluteDielectric(double dielectric) {
+        soluteDielectric = dielectric;
+    }
+protected:
+    ForceImpl* createImpl();
+private:
+    class ParticleInfo;
+    double solventDielectric, soluteDielectric;
+
+// 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;
+
+#if defined(_MSC_VER)
+#pragma warning(pop)
+#endif
+
+};
+
+class GBVIForce::ParticleInfo {
+public:
+    double charge, radius, gamma;
+    ParticleInfo() {
+        charge = radius = gamma = 0.0;
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_GBVIFORCEFIELD_H_*/

openmmapi/include/internal/GBVIForceImpl.h

+#ifndef OPENMM_GBVIFORCEFIELDIMPL_H_
+#define OPENMM_GBVIFORCEFIELDIMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "ForceImpl.h"
+#include "GBVIForce.h"
+#include "Kernel.h"
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of GBVIForce.
+ */
+
+class GBVIForceImpl : public ForceImpl {
+public:
+    GBVIForceImpl(GBVIForce& owner);
+    void initialize(OpenMMContextImpl& context);
+    GBVIForce& getOwner() {
+        return owner;
+    }
+    
+    // calculate scaled radii (Eq. 5 of Labute paper [JCC 29 1693-1698 2008])
+
+    void findScaledRadii( int numberOfParticles, const std::vector<std::vector<int> >& bondIndices,
+                          const std::vector<double> & bondLengths, std::vector<double> & scaledRadii) const;
+
+    void updateContextState(OpenMMContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    void calcForces(OpenMMContextImpl& context, Stream& forces);
+    double calcEnergy(OpenMMContextImpl& context);
+    std::map<std::string, double> getDefaultParameters() {
+        return std::map<std::string, double>(); // This force field doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+private:
+    GBVIForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_GBVIFORCEFIELDIMPL_H_*/

openmmapi/src/GBVIForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include "OpenMMException.h"
+#include "GBVIForce.h"
+#include "internal/GBVIForceImpl.h"
+
+using namespace OpenMM;
+
+GBVIForce::GBVIForce(int numParticles) : particles(numParticles), solventDielectric(78.3), soluteDielectric(1.0) {
+}
+
+void GBVIForce::getParticleParameters(int index, double& charge, double& radius, double& gamma) const {
+    charge = particles[index].charge;
+    radius = particles[index].radius;
+    gamma  = particles[index].gamma;
+}
+
+void GBVIForce::setParticleParameters(int index, double charge, double radius, double gamma) {
+    particles[index].charge = charge;
+    particles[index].radius = radius;
+    particles[index].gamma  = gamma;
+}
+
+ForceImpl* GBVIForce::createImpl() {
+    return new GBVIForceImpl(*this);
+}

openmmapi/src/GBVIForceImpl.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "internal/GBVIForceImpl.h"
+#include "internal/OpenMMContextImpl.h"
+#include "OpenMMException.h"
+#include "kernels.h"
+#include <vector>
+#include <math.h>
+
+using namespace OpenMM;
+using std::vector;
+
+GBVIForceImpl::GBVIForceImpl(GBVIForce& owner) : owner(owner) {
+}
+
+void GBVIForceImpl::initialize(OpenMMContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcGBVIForceKernel::Name(), context);
+
+    // load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
+    
+    System& system            = context.getSystem();
+    int     numberOfParticles = owner.getNumParticles();
+
+    vector<vector<int> > bondIndices;
+    vector<double> bondLengths;
+
+    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);
+            }
+            break;
+        }
+    }   
+
+    // Also treat constrained distances as bonds.
+
+    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);
+    }   
+
+    vector<double> scaledRadii;
+    scaledRadii.resize(numberOfParticles);
+    findScaledRadii( numberOfParticles, bondIndices, bondLengths, scaledRadii);
+
+    dynamic_cast<CalcGBVIForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner, scaledRadii);
+}
+
+#define GBVIDebug 1
+
+void GBVIForceImpl::findScaledRadii( int numberOfParticles, const std::vector<std::vector<int> >& bondIndices,
+                                     const std::vector<double> & bondLengths, std::vector<double> & scaledRadii) const {
+
+    // load 1-2 indicies for each atom 
+
+    std::vector<std::vector<int> > bonded12(numberOfParticles);
+
+    for (int i = 0; i < (int) bondIndices.size(); ++i) {
+        bonded12[bondIndices[i][0]].push_back(i);
+        bonded12[bondIndices[i][1]].push_back(i);
+    }
+
+    int errors = 0;
+
+    // compute scaled radii (Eq. 5 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+    for (int j = 0; j < (int) bonded12.size(); ++j){
+
+        double charge;
+        double gamma;
+        double radiusJ;
+        double scaledRadiusJ;
+     
+        owner.getParticleParameters(j, charge, radiusJ, gamma); 
+
+        if(  bonded12[j].size() == 0 ){
+            (void) fprintf( stderr, "Warning GBVIForceImpl::findScaledRadii atom %d has no covalent bonds; using atomic radius=%.3f.\n", j, radiusJ );
+            scaledRadiusJ = radiusJ;
+//             errors++;
+        } else {
+
+            double rJ2    = radiusJ*radiusJ;
+    
+            // loop over bonded neighbors of atom j, applying Eq. 5 in Labute
+
+            scaledRadiusJ = 0.0;
+            for (int i = 0; i < (int) bonded12[j].size(); ++i){
+    
+               int index            = bonded12[j][i];
+               int bondedAtomIndex  = (j == bondIndices[index][0]) ? bondIndices[index][1] : bondIndices[index][0];
+              
+               double radiusI;
+               owner.getParticleParameters(bondedAtomIndex, charge, radiusI, gamma); 
+               double rI2           = radiusI*radiusI;
+    
+               double a_ij          = (radiusI - bondLengths[index]);
+                      a_ij         *= a_ij;
+                      a_ij          = (rJ2 - a_ij)/(2.0*bondLengths[index]);
+    
+               double a_ji          = radiusJ - bondLengths[index];
+                      a_ji         *= a_ji;
+                      a_ji          = (rI2 - a_ji)/(2.0*bondLengths[index]);
+    
+               scaledRadiusJ       += a_ij*a_ij*(3.0*radiusI - a_ij) + a_ji*a_ji*( 3.0*radiusJ - a_ji );
+            }
+    
+            scaledRadiusJ  = (radiusJ*radiusJ*radiusJ) - 0.125*scaledRadiusJ; 
+            if( scaledRadiusJ > 0.0 ){
+                scaledRadiusJ  = 0.95*pow( scaledRadiusJ, (1.0/3.0) );
+            } else {
+                scaledRadiusJ  = 0.0;
+            }
+        }
+        scaledRadii[j] = scaledRadiusJ;
+
+    }
+
+    // abort if errors
+
+    if( errors ){
+        throw OpenMMException("GBVIForceImpl::findScaledRadii errors -- aborting");
+    }
+
+#if GBVIDebug
+    (void) fprintf( stderr, "                  R              q          gamma   scaled radii no. bnds\n" );
+    double totalQ = 0.0;
+    for( int i = 0; i < (int) scaledRadii.size(); i++ ){
+
+        double charge;
+        double gamma;
+        double radiusI;
+     
+        owner.getParticleParameters(i, charge, radiusI, gamma); 
+        totalQ += charge;
+        (void) fprintf( stderr, "%4d %14.5e %14.5e %14.5e %14.5e %d\n", i, radiusI, charge, gamma, scaledRadii[i], (int) bonded12[i].size() );
+    }
+    (void) fprintf( stderr, "Total charge=%e\n", totalQ );
+    (void) fflush( stderr );
+#endif
+
+#undef GBVIDebug
+
+}
+
+void GBVIForceImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {
+    dynamic_cast<CalcGBVIForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double GBVIForceImpl::calcEnergy(OpenMMContextImpl& context) {
+    return dynamic_cast<CalcGBVIForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+std::vector<std::string> GBVIForceImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(CalcGBVIForceKernel::Name());
+    return names;
+}

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -52,6 +52,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcRBTorsionForceKernel(name, platform);
     else if (name == CalcGBSAOBCForceKernel::Name())
         return new ReferenceCalcGBSAOBCForceKernel(name, platform);
+    else if (name == CalcGBVIForceKernel::Name())
+        return new ReferenceCalcGBVIForceKernel(name, platform);
     else if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform);
     else if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -32,6 +32,7 @@
 #include "ReferenceKernels.h"
 #include "ReferenceFloatStreamImpl.h"
 #include "gbsa/CpuObc.h"
+#include "gbsa/CpuGBVI.h"
 #include "SimTKReference/ReferenceAndersenThermostat.h"
 #include "SimTKReference/ReferenceAngleBondIxn.h"
 #include "SimTKReference/ReferenceBondForce.h"
@@ -442,7 +443,7 @@ void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBS
         const NonbondedForce* nonbonded = dynamic_cast<const NonbondedForce*>(&system.getForce(i));
         if (nonbonded != NULL) {
             if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
-                obcParameters->setUseCutoff(nonbonded->getCutoffDistance());
+                obcParameters->setUseCutoff(static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
             if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
                 Vec3 boxVectors[3];
                 nonbonded->getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
@@ -473,6 +474,79 @@ double ReferenceCalcGBSAOBCForceKernel::executeEnergy(OpenMMContextImpl& context
     return obc->getEnergy();
 }
 
+ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
+    if (gbvi) {
+        delete gbvi;
+    }
+}
+
+void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii ) {
+
+    int numParticles = system.getNumParticles();
+
+    charges.resize(numParticles);
+    vector<RealOpenMM> atomicRadii(numParticles);
+    vector<RealOpenMM> scaledRadii(numParticles);
+    vector<RealOpenMM> gammas(numParticles);
+
+    for (int i = 0; i < numParticles; ++i) {
+        double charge, radius, gamma;
+        force.getParticleParameters(i, charge, radius, gamma);
+        charges[i]       = static_cast<RealOpenMM>(charge);
+        atomicRadii[i]   = static_cast<RealOpenMM>(radius);
+        gammas[i]        = static_cast<RealOpenMM>(gamma);
+        scaledRadii[i]   = static_cast<RealOpenMM>(inputScaledRadii[i]);
+    }
+
+    GBVIParameters * gBVIParameters = new GBVIParameters(numParticles);
+    gBVIParameters->setAtomicRadii(atomicRadii);
+    gBVIParameters->setGammaParameters(gammas);
+    gBVIParameters->setScaledRadii(scaledRadii);
+    gBVIParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
+    gBVIParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
+
+    // If there is a NonbondedForce in this system, use it to initialize cutoffs and periodic boundary conditions.
+
+    for (int i = 0; i < system.getNumForces(); i++) {
+        const NonbondedForce* nonbonded = dynamic_cast<const NonbondedForce*>(&system.getForce(i));
+        if (nonbonded != NULL) {
+            if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
+                gBVIParameters->setUseCutoff( static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
+            if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
+                Vec3 boxVectors[3];
+                nonbonded->getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+                RealOpenMM periodicBoxSize[3];
+                periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+                periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+                periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
+                gBVIParameters->setPeriodic(periodicBoxSize);
+            }
+            break;
+        }
+    }
+    gbvi = new CpuGBVI(gBVIParameters);
+
+}
+
+void ReferenceCalcGBVIForceKernel::executeForces(OpenMMContextImpl& context) {
+
+    RealOpenMM** posData   = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
+    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
+    RealOpenMM* bornRadii  = new RealOpenMM[context.getSystem().getNumParticles()];
+    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
+    gbvi->computeBornForces(bornRadii, posData, &charges[0], forceData);
+    delete[] bornRadii;
+}
+
+double ReferenceCalcGBVIForceKernel::executeEnergy(OpenMMContextImpl& context) {
+    RealOpenMM** posData  = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
+    RealOpenMM* bornRadii = new RealOpenMM[context.getSystem().getNumParticles()];
+    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
+    RealOpenMM energy     = gbvi->computeBornEnergy(bornRadii ,posData, &charges[0]);
+    delete[] bornRadii;
+    return static_cast<double>(energy);
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferenceKernels.h

@@ -37,6 +37,7 @@
 #include "SimTKReference/ReferenceNeighborList.h"
 
 class CpuObc;
+class CpuGBVI;
 class ReferenceAndersenThermostat;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
@@ -276,6 +277,40 @@ private:
     std::vector<RealOpenMM> charges;
 };
 
+/**
+ * This kernel is invoked by GBVIForce to calculate the forces acting on the system.
+ */
+class ReferenceCalcGBVIForceKernel : public CalcGBVIForceKernel {
+public:
+    ReferenceCalcGBVIForceKernel(std::string name, const Platform& platform) : CalcGBVIForceKernel(name, platform) {
+    }
+    ~ReferenceCalcGBVIForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system       the System this kernel will be applied to
+     * @param force        the GBVIForce this kernel will be used for
+     * @param scaled radii 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(OpenMMContextImpl& 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(OpenMMContextImpl& context);
+private:
+    CpuGBVI * gbvi;
+    std::vector<RealOpenMM> charges;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -45,6 +45,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
+    registerKernelFactory(CalcGBVIForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/reference/src/gbsa/CpuGBVI.cpp

+
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "CpuGBVI.h"
+#include "../SimTKReference/ReferenceForce.h"
+#include <math.h>
+
+static const double CAL_TO_JOULE = 0.4184;
+
+/**---------------------------------------------------------------------------------------
+
+   CpuGBVI constructor
+
+   gbviParameters      gbviParameters object
+   
+   --------------------------------------------------------------------------------------- */
+
+CpuGBVI::CpuGBVI( ImplicitSolventParameters* gbviParameters ) : CpuImplicitSolvent( gbviParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::CpuGBVI";
+
+   // ---------------------------------------------------------------------------------------
+
+   _initializeGBVIDataMembers( );
+
+   _gbviParameters = static_cast<GBVIParameters*> (gbviParameters);
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   CpuGBVI destructor
+
+   --------------------------------------------------------------------------------------- */
+
+CpuGBVI::~CpuGBVI( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::~CpuGBVI";
+
+   // ---------------------------------------------------------------------------------------
+
+   //if( _gbviParameters != NULL ){
+     // delete _gbviParameters;
+   //}
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Initialize data members
+
+   --------------------------------------------------------------------------------------- */
+
+void CpuGBVI::_initializeGBVIDataMembers( void ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::initializeDataMembers";
+
+   // ---------------------------------------------------------------------------------------
+
+   _gbviParameters  = NULL;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get GBVIParameters reference
+
+   @return GBVIParameters reference
+
+   --------------------------------------------------------------------------------------- */
+
+GBVIParameters* CpuGBVI::getGBVIParameters( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::getGBVIParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _gbviParameters;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set GBVIParameters reference
+
+   @param GBVIParameters reference
+
+   @return SimTKOpenMMCommon::DefaultReturn;
+
+   --------------------------------------------------------------------------------------- */
+
+int CpuGBVI::setGBVIParameters(  GBVIParameters* gbviParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::setGBVIParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   _gbviParameters = gbviParameters;
+   return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   @param atomCoordinates     atomic coordinates
+   @param bornRadii           output array of Born radii
+   @param chain               not used here
+
+   @return SimTKOpenMMCommon::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+#define GBVIDebug 0
+
+int CpuGBVI::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii, RealOpenMM* chain ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const RealOpenMM zero          = (RealOpenMM) 0.0;
+   static const RealOpenMM one           = (RealOpenMM) 1.0;
+   static const RealOpenMM minusThree    = (RealOpenMM) -3.0;
+   static const RealOpenMM oneEighth     = (RealOpenMM) 0.125;
+   static const RealOpenMM minusOneThird = (RealOpenMM) (-1.0/3.0);
+   static const RealOpenMM three         = (RealOpenMM) 3.0;
+
+   static const char* methodName         = "CpuGBVI::computeBornRadii";
+
+   // ---------------------------------------------------------------------------------------
+
+   GBVIParameters* gbviParameters        = getGBVIParameters();
+   int numberOfAtoms                     = gbviParameters->getNumberOfAtoms();
+   RealOpenMM* atomicRadii               = gbviParameters->getAtomicRadii();
+   const RealOpenMM* scaledRadii         = gbviParameters->getScaledRadii();
+
+   // ---------------------------------------------------------------------------------------
+
+#if GBVIDebug
+   FILE* logFile                         = stderr;
+   (void) fprintf( logFile, "\n%s\n", methodName );
+#endif
+
+   // calculate Born radii
+
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+     
+      RealOpenMM radiusI         = atomicRadii[atomI];
+      RealOpenMM sum             = zero;
+
+      // sum over volumes
+
+      for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
+
+         if( atomJ != atomI ){
+
+            RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+            if (_gbviParameters->getPeriodic())
+                ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
+            else
+                ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+
+            RealOpenMM r               = deltaR[ReferenceForce::RIndex];
+
+            if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
+                continue;
+
+            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 );
+#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] );
+#endif
+
+   }
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+#undef GBVIDebug
+
+/**---------------------------------------------------------------------------------------
+
+   Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return volume
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM CpuGBVI::getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::getVolume";
+
+   static const RealOpenMM zero         = (RealOpenMM)  0.0;
+   static const RealOpenMM minusThree   = (RealOpenMM) -3.0;
+
+   RealOpenMM              diff         = (S - R);
+   if( FABS( diff ) < r ){
+
+      RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
+
+      return (CpuGBVI::getL( r, (r + S),    S ) -
+              CpuGBVI::getL( r, lowerBound, S ));
+
+   } else if( r <= diff ){
+
+      return CpuGBVI::getL( r, (r + S), S ) -
+             CpuGBVI::getL( r, (r - S), S ) + 
+             POW( R, minusThree );
+
+   } else {
+      return zero;
+   }
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get L (used in analytical solution for volume integrals) 
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM CpuGBVI::getL( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::getL";
+
+   static const RealOpenMM one           = (RealOpenMM) 1.0;
+   static const RealOpenMM threeHalves   = (RealOpenMM) 1.5;
+   static const RealOpenMM third         = (RealOpenMM) (1.0/3.0);
+   static const RealOpenMM fourth        = (RealOpenMM) 0.25;
+   static const RealOpenMM eighth        = (RealOpenMM) 0.125;
+
+   // ---------------------------------------------------------------------------------------
+
+   RealOpenMM rInv   = one/r;
+
+   RealOpenMM xInv   = one/x;
+   RealOpenMM xInv2  = xInv*xInv;
+   RealOpenMM xInv3  = xInv2*xInv;
+
+   RealOpenMM diff2  = (r + S)*(r - S);
+
+   return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get partial derivative of L wrt r
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM CpuGBVI::dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::dL_dr";
+
+   static const RealOpenMM one           = (RealOpenMM) 1.0;
+   static const RealOpenMM threeHalves   = (RealOpenMM) 1.5;
+   static const RealOpenMM threeEights   = (RealOpenMM) 0.375;
+   static const RealOpenMM third         = (RealOpenMM) (1.0/3.0);
+   static const RealOpenMM fourth        = (RealOpenMM) 0.25;
+   static const RealOpenMM eighth        = (RealOpenMM) 0.125;
+
+   // ---------------------------------------------------------------------------------------
+
+   RealOpenMM rInv   = one/r;
+   RealOpenMM rInv2  = rInv*rInv;
+
+   RealOpenMM xInv   = one/x;
+   RealOpenMM xInv2  = xInv*xInv;
+   RealOpenMM xInv3  = xInv2*xInv;
+
+   RealOpenMM diff2  = (r + S)*(r - S);
+
+   return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get partial derivative of L wrt x
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM CpuGBVI::dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::dL_dx";
+
+   static const RealOpenMM one           = (RealOpenMM)  1.0;
+   static const RealOpenMM half          = (RealOpenMM)  0.5;
+   static const RealOpenMM threeHalvesM  = (RealOpenMM) -1.5;
+   static const RealOpenMM third         = (RealOpenMM)  (1.0/3.0);
+
+   // ---------------------------------------------------------------------------------------
+
+   RealOpenMM rInv   = one/r;
+
+   RealOpenMM xInv   = one/x;
+   RealOpenMM xInv2  = xInv*xInv;
+   RealOpenMM xInv3  = xInv2*xInv;
+
+   RealOpenMM diff   = (r + S)*(r - S);
+
+   return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Sgb function
+
+   @param t                   r*r*G_i*G_j
+
+   @return Sgb (top of p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM CpuGBVI::Sgb( RealOpenMM t ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::Sgb";
+
+   static const RealOpenMM zero    = (RealOpenMM) 0.0;
+   static const RealOpenMM one     = (RealOpenMM) 1.0;
+   static const RealOpenMM fourth  = (RealOpenMM) 0.25;
+
+   // ---------------------------------------------------------------------------------------
+
+   return ( (t != zero) ? one/SQRT( (one + (fourth*EXP( -t ))/t) ) : zero);
+}
+
+#define GBVIDebug 0
+
+/**---------------------------------------------------------------------------------------
+
+   Get GB/VI energy
+
+   @param bornRadii           Born radii
+   @param atomCoordinates     atomic coordinates
+   @param partialCharges      partial charges
+
+   @return energy
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM CpuGBVI::computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
+                                       const RealOpenMM* partialCharges ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "CpuGBVI::computeBornEnergy";
+
+   static const RealOpenMM zero    = (RealOpenMM) 0.0;
+   static const RealOpenMM one     = (RealOpenMM) 1.0;
+   static const RealOpenMM two     = (RealOpenMM) 2.0;
+   static const RealOpenMM three   = (RealOpenMM) 3.0;
+   static const RealOpenMM four    = (RealOpenMM) 4.0;
+   static const RealOpenMM half    = (RealOpenMM) 0.5;
+   static const RealOpenMM fourth  = (RealOpenMM) 0.25;
+   static const RealOpenMM eighth  = (RealOpenMM) 0.125;
+
+   // ---------------------------------------------------------------------------------------
+
+   const GBVIParameters* gbviParameters = getGBVIParameters();
+   const RealOpenMM preFactor           = gbviParameters->getElectricConstant();
+   const int numberOfAtoms              = gbviParameters->getNumberOfAtoms();
+   const RealOpenMM* atomicRadii        = gbviParameters->getAtomicRadii();
+   const RealOpenMM* gammaParameters    = gbviParameters->getGammaParameters();
+
+   if( bornRadii == NULL ){
+      bornRadii   = getBornRadii();
+   }
+
+#if GBVIDebug
+   FILE* logFile                        = stderr;
+   (void) fprintf( logFile, "\n%s\n", methodName );
+   (void) fflush( logFile );
+#endif
+
+   // ---------------------------------------------------------------------------------------
+
+   // Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
+
+   RealOpenMM energy                 = zero;
+
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+ 
+      RealOpenMM partialChargeI   = preFactor*partialCharges[atomI];
+      RealOpenMM partialChargeI2  = two*partialChargeI;
+
+      // self-energy term
+
+      energy                      += partialChargeI*partialCharges[atomI]/bornRadii[atomI];
+
+      // cavity term
+
+      RealOpenMM ratio             = (atomicRadii[atomI]/bornRadii[atomI]);
+      energy                      -= gammaParameters[atomI]*ratio*ratio*ratio;
+
+/*
+RealOpenMM e1 = partialChargeI*partialCharges[atomI]/bornRadii[atomI];
+RealOpenMM e2 = gammaParameters[atomI]*ratio*ratio*ratio;
+(void) fprintf( stderr, "E %d self=%.4e gamma=%.4e e=%.4e\n", atomI, e1, e2, energy );
+*/
+      for( int atomJ = atomI + 1; atomJ < numberOfAtoms; atomJ++ ){
+
+         RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+         if (_gbviParameters->getPeriodic())
+             ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
+         else
+             ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+         if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
+             continue;
+
+         RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
+         RealOpenMM t                  = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);         
+         energy                       += partialChargeI2*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 );
+*/
+      }
+   }
+
+#if GBVIDebug
+   (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] );
+   }
+   (void) fflush( logFile );
+#endif
+
+   RealOpenMM conversion = CAL_TO_JOULE*gbviParameters->getTau();  
+   return (conversion*energy);
+ 
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get GB/VI forces
+
+   @param bornRadii           Born radii
+   @param atomCoordinates     atomic coordinates
+   @param partialCharges      partial charges
+   @param forces              forces
+
+   @return SimTKOpenMMCommon::DefaultReturn;
+
+   --------------------------------------------------------------------------------------- */
+
+
+int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
+                                const RealOpenMM* partialCharges, RealOpenMM** inputForces ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName        = "CpuGBVI::computeBornEnergyForces";
+
+   static const RealOpenMM zero         = (RealOpenMM) 0.0;
+   static const RealOpenMM one          = (RealOpenMM) 1.0;
+   static const RealOpenMM two          = (RealOpenMM) 2.0;
+   static const RealOpenMM three        = (RealOpenMM) 3.0;
+   static const RealOpenMM four         = (RealOpenMM) 4.0;
+   static const RealOpenMM half         = (RealOpenMM) 0.5;
+   static const RealOpenMM oneThird     = (RealOpenMM) (1.0/3.0);
+   static const RealOpenMM fourth       = (RealOpenMM) 0.25;
+   static const RealOpenMM eighth       = (RealOpenMM) 0.125;
+
+   // ---------------------------------------------------------------------------------------
+
+#if GBVIDebug
+   FILE* logFile                        = stderr;
+   (void) fprintf( logFile, "\n%s\n", methodName );
+   (void) fflush( logFile );
+#endif
+
+   const GBVIParameters* gbviParameters = getGBVIParameters();
+   const int numberOfAtoms              = gbviParameters->getNumberOfAtoms();
+   const RealOpenMM* atomicRadii        = gbviParameters->getAtomicRadii();
+   const RealOpenMM* gammaParameters    = gbviParameters->getGammaParameters();
+
+   if( bornRadii == NULL ){
+      bornRadii   = getBornRadii();
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // constants
+
+   const RealOpenMM preFactor           = two*gbviParameters->getElectricConstant();
+
+   // ---------------------------------------------------------------------------------------
+
+   // set energy/forces to zero
+
+   const unsigned int arraySzInBytes    = sizeof( RealOpenMM )*numberOfAtoms;
+
+   RealOpenMM** forces  = (RealOpenMM**) malloc( sizeof( RealOpenMM* )*numberOfAtoms );
+   RealOpenMM*  block   = (RealOpenMM*)  malloc( sizeof( RealOpenMM )*numberOfAtoms*3 );
+	memset( block, 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+	RealOpenMM* blockPtr = block;
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      forces[ii] = blockPtr;
+		blockPtr  += 3;
+   }
+
+   RealOpenMM* bornForces = getBornForce();
+   memset( bornForces, 0, arraySzInBytes );
+
+   // ---------------------------------------------------------------------------------------
+
+   // first main loop
+
+   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)
+
+      RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
+      for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
+
+         RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+         if (_gbviParameters->getPeriodic())
+             ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
+         else
+             ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+         if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
+             continue;
+
+         RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
+
+         RealOpenMM deltaX             = deltaR[ReferenceForce::XIndex];
+         RealOpenMM deltaY             = deltaR[ReferenceForce::YIndex];
+         RealOpenMM deltaZ             = deltaR[ReferenceForce::ZIndex];
+
+         RealOpenMM alpha2_ij          = bornRadii[atomI]*bornRadii[atomJ];
+         RealOpenMM D_ij               = r2/(four*alpha2_ij);
+
+         RealOpenMM expTerm            = EXP( -D_ij );
+         RealOpenMM denominator2       = r2 + alpha2_ij*expTerm; 
+         RealOpenMM denominator        = SQRT( denominator2 ); 
+         
+         RealOpenMM Gpol               = (partialChargeI*partialCharges[atomJ])/denominator; 
+         RealOpenMM dGpol_dr           = -Gpol*( one - fourth*expTerm )/denominator2;  
+
+         RealOpenMM dGpol_dalpha2_ij   = -half*Gpol*expTerm*( one + D_ij )/denominator2;
+
+         if( atomI != atomJ ){
+
+             bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
+
+             deltaX            *= dGpol_dr;
+             deltaY            *= dGpol_dr;
+             deltaZ            *= dGpol_dr;
+
+             forces[atomI][0]  += deltaX;
+             forces[atomI][1]  += deltaY;
+             forces[atomI][2]  += deltaZ;
+
+             forces[atomJ][0]  -= deltaX;
+             forces[atomJ][1]  -= deltaY;
+             forces[atomJ][2]  -= deltaZ;
+
+         }
+
+         // 3 FLOP
+
+         bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
+
+      }
+   }
+
+#if GBVIDebug
+   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] );
+   }
+   (void) fflush( logFile );
+#endif
+
+   // ---------------------------------------------------------------------------------------
+
+   // second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
+
+   // dGpol/dBornRadius) = bornForces[]
+   // dBornRadius/dr     = (1/3)*(bR**4)*(dV/dr)
+
+   const RealOpenMM* scaledRadii         = gbviParameters->getScaledRadii();
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+ 
+      RealOpenMM R        = atomicRadii[atomI];
+
+      RealOpenMM b2       = bornRadii[atomI]*bornRadii[atomI];
+      bornForces[atomI]  *= oneThird*b2*b2;
+
+      for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
+
+         if( atomJ != atomI ){
+
+            RealOpenMM deltaX             = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
+            RealOpenMM deltaY             = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
+            RealOpenMM deltaZ             = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
+    
+            RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+            if (_gbviParameters->getPeriodic())
+                ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
+            else
+                ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+            if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
+                continue;
+   
+            RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
+                       deltaX             = deltaR[ReferenceForce::XIndex];
+                       deltaY             = deltaR[ReferenceForce::YIndex];
+                       deltaZ             = deltaR[ReferenceForce::ZIndex];
+
+            RealOpenMM r                  = SQRT( r2 );
+ 
+            RealOpenMM S                  = scaledRadii[atomJ];
+            RealOpenMM diff               = (S - R);
+
+            RealOpenMM de                 = zero;
+
+            // find dRb/dr, where Rb is the Born radius
+
+            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 ) );   
+            }
+
+if( 0 ){
+   RealOpenMM delta = (RealOpenMM) 1.0e-02;
+   (void) fprintf( stderr, "\n" );
+   for( int kk = 0; kk < 5; kk++ ){
+      RealOpenMM V1    = CpuGBVI::getVolume( r, R, S );
+      RealOpenMM V2    = CpuGBVI::getVolume( r+delta, R, S );
+      RealOpenMM df    = (V2-V1)/delta;
+      (void) fprintf( stderr, "df %d %d [%14.6e %14.6e] V[%14.6e %14.6e] %.2e\n", atomI, atomJ, de, df, V2, V1, delta );
+      delta *= (RealOpenMM) 0.1;
+   }
+
+   double deltaD = 1.0e-02;
+   double ded = CpuGBVI::dL_drD( (double) r, r+S, S ) + CpuGBVI::dL_dxD( r, r+S, S ) - ( CpuGBVI::dL_drD( r, (r-S), S ) + CpuGBVI::dL_dxD( r, (r-S), S ) );
+   for( int kk = 0; kk < 5; kk++ ){
+      double V1    = CpuGBVI::getVolumeD( r, R, S );
+      double V2    = CpuGBVI::getVolumeD( r+deltaD, R, S );
+      double df    = (V2-V1)/deltaD;
+      (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;
+   }
+}
+
+
+
+
+            // de = (dG/dRb)(dRb/dr)
+
+            de                      *= bornForces[atomI]/r;
+
+            deltaX                  *= de;
+            deltaY                  *= de;
+            deltaZ                  *= de;
+   
+            forces[atomI][0]        += deltaX;
+            forces[atomI][1]        += deltaY;
+            forces[atomI][2]        += deltaZ;
+  
+            forces[atomJ][0]        -= deltaX;
+            forces[atomJ][1]        -= deltaY;
+            forces[atomJ][2]        -= deltaZ;
+
+         }
+      }
+
+   }
+
+#if GBVIDebug
+   (void) fprintf( logFile, "Atom      BornRadii      BornForce                                         Forces\n" );
+   double forceSum[3] = { 0.0, 0.0, 0.0 };
+   for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+      forceSum[0] += forces[atomI][0];
+      forceSum[1] += forces[atomI][1];
+      forceSum[2] += forces[atomI][2];
+      (void) fprintf( logFile, "%4d %14.6e %14.6e [%14.6e %14.6e %14.6e]\n", atomI, bornRadii[atomI], bornForces[atomI],  forces[atomI][0], forces[atomI][1], forces[atomI][2] );
+   }   
+   (void) fprintf( logFile, "F sum=[%14.6e %14.6e %14.6e]\n", forceSum[0], forceSum[1], forceSum[2] );
+   (void) fflush( logFile );
+#endif
+
+#undef GBVIDebug
+
+   // convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
+
+   RealOpenMM conversion = 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];
+   }
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Get string w/ state 
+   
+   @param title               title (optional)
+      
+   @return string containing state
+      
+   --------------------------------------------------------------------------------------- */
+
+std::string CpuGBVI::getStateString( const char* title ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuImplicitSolvent::getStateString";
+
+   // ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   message << CpuImplicitSolvent::getStateString( title );
+
+   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 CpuGBVI::writeBornEnergyForces( RealOpenMM** atomCoordinates,
+                                   const RealOpenMM* partialCharges, RealOpenMM** forces,
+                                   const std::string& resultsFileName ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nCpuGBVI::writeBornEnergyForces";
+
+   // ---------------------------------------------------------------------------------------
+/*
+   ImplicitSolventParameters* implicitSolventParameters = getImplicitSolventParameters();
+   const GBVIParameters* gbviParameters                   = static_cast<const GBVIParameters*>(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 SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   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 CpuGBVI::writeForceLoop1( int numberOfAtoms, RealOpenMM** forces, const RealOpenMM* bornForce,
+                             const std::string& outputFileName ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nCpuGBVI::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 CpuGBVI::writeForceLoop( int numberOfAtoms, const IntVector& chunkSizes,
+                            const RealOpenMMPtrPtrVector& realRealOpenMMVector, 
+                            const RealOpenMMPtrVector& realVector,
+                            const std::string& outputFileName ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nCpuGBVI::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 GBVIParameters 
+                              entry is used
+   @param atomCoordinates     atomic coordinates
+   @param partialCharges      partial charges
+   @param forces              forces
+
+   @return SimTKOpenMMCommon::DefaultReturn;
+
+   The array bornRadii is also updated and the obcEnergy
+
+   --------------------------------------------------------------------------------------- */
+
+int CpuGBVI::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
+                                      const RealOpenMM* partialCharges, RealOpenMM** forces ){
+ 
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuGBVI::computeBornEnergyForcesPrint";
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Use double precision 
+
+   Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return volume
+
+   --------------------------------------------------------------------------------------- */
+
+double CpuGBVI::getVolumeD( double r, double R, double S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::getVolume";
+
+   static const double zero         = 0.0;
+   static const double minusThree   = -3.0;
+
+   double              diff    = (S - R);
+   if( fabs( diff ) < r ){
+
+      double lowerBound = (R > (r - S)) ? R : (r - S);
+
+      return (CpuGBVI::getLD( r, (r + S),    S ) -
+              CpuGBVI::getLD( r, lowerBound, S ));
+
+   } else if( r < diff ){
+
+      return CpuGBVI::getLD( r, (r + S), S ) -
+             CpuGBVI::getLD( r, (r - S), S ) + 
+             pow( R, minusThree );
+
+   } else {
+      return zero;
+   }
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Use double precision 
+
+   Get L (used in analytical solution for volume integrals) 
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+double CpuGBVI::getLD( double r, double x, double S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::getL";
+
+   static const double one           =  1.0;
+   static const double threeHalves   =  1.5;
+   static const double third         =  1.0/3.0;
+   static const double fourth        =  0.25;
+   static const double eighth        =  0.125;
+
+   // ---------------------------------------------------------------------------------------
+
+   double rInv   = one/r;
+
+   double xInv   = one/x;
+   double xInv2  = xInv*xInv;
+   double xInv3  = xInv2*xInv;
+
+   double diff2  = (r + S)*(r - S);
+
+   return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Use double precision 
+
+   Get partial derivative of L wrt r
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+double CpuGBVI::dL_drD( double r, double x, double S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::dL_dr";
+
+   static const double one           =  1.0;
+   static const double threeHalves   =  1.5;
+   static const double threeEights   =  0.375;
+   static const double third         =  1.0/3.0;
+   static const double fourth        =  0.25;
+   static const double eighth        =  0.125;
+
+   // ---------------------------------------------------------------------------------------
+
+   double rInv   = one/r;
+   double rInv2  = rInv*rInv;
+
+   double xInv   = one/x;
+   double xInv2  = xInv*xInv;
+   double xInv3  = xInv2*xInv;
+
+   double diff2  = (r + S)*(r - S);
+
+   return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Use double precision 
+
+   Get partial derivative of L wrt x
+
+   @param r                   distance between atoms i & j
+   @param R                   atomic radius
+   @param S                   scaled atomic radius
+
+   @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+   --------------------------------------------------------------------------------------- */
+
+double CpuGBVI::dL_dxD( double r, double x, double S ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "CpuGBVI::dL_dx";
+
+   static const double one           =   1.0;
+   static const double half          =   0.5;
+   static const double threeHalvesM  =  -1.5;
+   static const double third         =   1.0/3.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   double rInv   = one/r;
+
+   double xInv   = one/x;
+   double xInv2  = xInv*xInv;
+   double xInv3  = xInv2*xInv;
+
+   double diff   = (r + S)*(r - S);
+
+   return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
+}

platforms/reference/src/gbsa/CpuGBVI.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __CpuGBVI_H__
+#define __CpuGBVI_H__
+
+#include "GBVIParameters.h"
+#include "CpuImplicitSolvent.h"
+
+// ---------------------------------------------------------------------------------------
+
+class CpuGBVI : public CpuImplicitSolvent {
+
+   private:
+
+      // GB/VI parameters
+
+      GBVIParameters* _gbviParameters;
+
+      // initialize data members (more than
+      // one constructor, so centralize intialization here)
+
+      void _initializeGBVIDataMembers( void );
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+      
+         @param implicitSolventParameters    ImplicitSolventParameters reference
+      
+         @return CpuImplicitSolvent object
+      
+         --------------------------------------------------------------------------------------- */
+
+       CpuGBVI( ImplicitSolventParameters* gbviParameters );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~CpuGBVI( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Return GBVIParameters
+      
+         @return GBVIParameters
+      
+         --------------------------------------------------------------------------------------- */
+
+      GBVIParameters* getGBVIParameters( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set ImplicitSolventParameters
+      
+         @param ImplicitSolventParameters
+      
+         @return SimTKOpenMMCommon::DefaultReturn 
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setGBVIParameters( GBVIParameters* gbviParameters );
+ 
+      /**---------------------------------------------------------------------------------------
+      
+         Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
+
+         @param atomCoordinates   atomic coordinates
+         @param bornRadii         output array of Born radii
+         @param gbviChain         not used
+      
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii,
+                            RealOpenMM* gbviChain = 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 
+      
+         title             title (optional)
+      
+         @return state string
+      
+         --------------------------------------------------------------------------------------- */
+      
+      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])
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return volume
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get L (analytical solution for volume integrals) 
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM getL( RealOpenMM r, RealOpenMM x, RealOpenMM S );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get partial derivative of L wrt r
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get partial derivative of L wrt x
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Sgb function
+      
+         @param t                   r*r*G_i*G_j
+      
+         @return Sgb (p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static RealOpenMM Sgb( RealOpenMM t );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get GB/VI energy
+      
+         @param bornRadii           Born radii
+         @param atomCoordinates     atomic coordinates
+         @param partialCharges      partial charges
+      
+         @return energy
+      
+         --------------------------------------------------------------------------------------- */
+      
+      RealOpenMM computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
+                                    const RealOpenMM* partialCharges );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get GB/VI forces
+      
+         @param bornRadii           Born radii
+         @param atomCoordinates     atomic coordinates
+         @param partialCharges      partial charges
+         @param forces              output forces
+      
+         @return SimTKOpenMMCommon::DefaultReturn;
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
+                             const RealOpenMM* partialCharges, RealOpenMM** inputForces );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return volume
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static double getVolumeD( double r, double R, double S );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get L (analytical solution for volume integrals) 
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static double getLD( double r, double x, double S );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get partial derivative of L wrt r
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static double dL_drD( double r, double x, double S );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get partial derivative of L wrt x
+      
+         @param r                   distance between atoms i & j
+         @param R                   atomic radius
+         @param S                   scaled atomic radius
+      
+         @return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
+      
+         --------------------------------------------------------------------------------------- */
+      
+      static double dL_dxD( double r, double x, double S );
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __CpuGBVI_H__

platforms/reference/src/gbsa/GBVIParameters.cpp

+
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <math.h>
+#include <sstream>
+
+#include "GBVIParameters.h"
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+
+// #define UseGromacsMalloc 1
+
+#ifdef UseGromacsMalloc
+extern "C" {
+#include "smalloc.h" 
+}
+#endif
+
+const std::string GBVIParameters::ParameterFileName = std::string( "params.agb" );
+
+/**---------------------------------------------------------------------------------------
+
+   GBVIParameters:
+
+		Calculates for each atom
+
+			(1) the van der Waal radii
+         (2) volume
+         (3) fixed terms in Obc equation gPol
+         (4) list of atoms that should be excluded in calculating
+				 force -- nonbonded atoms (1-2, and 1-3 atoms)
+
+	Implementation:
+
+		Slightly different sequence of calls when running on CPU vs GPU.
+		Difference arise because the CPU-side data arrays for the Brook
+		streams are allocated by the BrookStreamWrapper objects. These
+		arrays are then used by GBVIParameters when initializing the
+		the values (vdwRadii, volume, ...) to be used in the calculation.
+
+		Cpu:
+			 GBVIParameters* gb_VIParameters = new GBVIParameters( numberOfAtoms, log );
+          gb_VIParameters->initializeParameters( top );
+
+		Gpu:
+
+			gb_VIParameters   = new GBVIParameters( gpu->natoms, log );
+			
+			// set arrays for cpu using stream data field; 
+			// initializeParameters() only allocates space for arrays if they are not set (==NULL)
+			// also set flag so that GBVIParameters destructor does not free arrays 
+			
+			gb_VIParameters->setVdwRadii(  getBrookStreamWrapperAtIndex( GpuObc::gb_VIVdwRadii  )->getData() );
+			gb_VIParameters->setVolume(    getBrookStreamWrapperAtIndex( GpuObc::gb_VIVolume    )->getData() );
+			gb_VIParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuObc::gb_VIGpolFixed )->getData() );
+			gb_VIParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuObc::gb_VIBornRadii )->getData() );
+			
+			gb_VIParameters->setFreeArrays( false );
+			
+			gb_VIParameters->initializeParameters( top );
+ 
+
+   Issues:
+
+		Tinker's atom radii are used. 
+      The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
+      only tested for generic proteins
+		see mapGmxAtomNameToTinkerAtomNumber()
+
+   --------------------------------------------------------------------------------------- */
+
+
+/**---------------------------------------------------------------------------------------
+
+   GBVIParameters constructor (Simbios) 
+
+   @param numberOfAtoms       number of atoms
+
+   --------------------------------------------------------------------------------------- */
+
+GBVIParameters::GBVIParameters( int numberOfAtoms ) : ImplicitSolventParameters( numberOfAtoms ), cutoff(false), periodic(false) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::GBVIParameters";
+   
+   // ---------------------------------------------------------------------------------------
+
+   _ownScaledRadii       = 0;
+   _scaledRadii          = NULL;
+   _ownGammaParameters   = 0;
+   _gammaParameters      = NULL;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   GBVIParameters destructor (Simbios) 
+
+   --------------------------------------------------------------------------------------- */
+
+GBVIParameters::~GBVIParameters( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::~GBVIParameters";
+   
+   // ---------------------------------------------------------------------------------------
+
+   // in GPU runs, arrays may be 'owned' by BrookStreamWrapper -- hence they should not
+   // be freed here, i.e., _freeArrays should be 'false'   
+
+#ifdef UseGromacsMalloc
+
+/*
+   if( _freeArrays ){
+
+      if( _vdwRadii != NULL ){
+         save_free( "_vdwRadii", __FILE__, __LINE__, _vdwRadii );
+      }
+   
+   } */
+
+#else
+
+   if( _ownScaledRadii ){
+      delete[] _scaledRadii;
+   }
+   delete _gammaParameters;
+/*
+   if( getFreeArrays() ){
+
+   } */
+
+#endif
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get AtomicRadii array
+
+   @return array of atomic radii
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM* GBVIParameters::getAtomicRadii( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
+
+   // ---------------------------------------------------------------------------------------
+
+   RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
+
+   return atomicRadii;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set AtomicRadii array
+
+   @param atomicRadii array of atomic radii
+
+   @return SimTKOpenMMCommon::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::setAtomicRadii:";
+
+   // ---------------------------------------------------------------------------------------
+
+   return ImplicitSolventParameters::setAtomicRadii( atomicRadii );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set AtomicRadii array
+
+   @param atomicRadii vector of atomic radii
+
+   @return SimTKOpenMMCommon::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\nGBVIParameters::setAtomicRadii:";
+
+   // ---------------------------------------------------------------------------------------
+
+   return ImplicitSolventParameters::setAtomicRadii( atomicRadii );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Return scaled radii
+   If not previously set, allocate space
+
+   @return array 
+
+   --------------------------------------------------------------------------------------- */
+
+const RealOpenMM* GBVIParameters::getScaledRadii( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::getScaledRadii";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _scaledRadii == NULL ){
+      GBVIParameters* localThis  = const_cast<GBVIParameters* const>(this);
+      localThis->_scaledRadii    = new RealOpenMM[getNumberOfAtoms()];
+      localThis->_ownScaledRadii = true;
+      memset( _scaledRadii, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
+   }   
+   return _scaledRadii;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set flag indicating whether scale factors array should be deleted
+
+   @param ownScaledRadii flag indicating whether scale factors 
+                                 array should be deleted
+
+   @return SimTKOpenMMCommon::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setOwnScaledRadii( int ownScaledRadii ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::setOwnScaleFactors";
+
+   // ---------------------------------------------------------------------------------------
+
+   _ownScaledRadii = ownScaledRadii;
+
+   return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set scaled radii
+
+   @param scaledRadii  scaledRadii
+
+   @return SimTKOpenMMCommon::DefaultReturn always
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setScaledRadii( RealOpenMM* scaledRadii ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuObc::setScaledRadii";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _ownScaledRadii && _scaledRadii != scaledRadii ){
+      delete[] _scaledRadii;
+      _ownScaledRadii = false;
+   }
+
+   _scaledRadii = scaledRadii;
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+#if RealOpenMMType == 2
+
+/**---------------------------------------------------------------------------------------
+
+   Set scaled radii
+
+   @param scaledRadii  scaledRadii
+
+   @return SimTKOpenMMCommon::DefaultReturn always
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setScaledRadii( float* scaledRadii ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::setScaledRadii";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _scaledRadii == NULL ){
+      _scaledRadii    = new RealOpenMM[getNumberOfAtoms()];
+      _ownScaledRadii = true;
+   }   
+   for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
+      _scaledRadii[ii] = (RealOpenMM) scaledRadii[ii];
+   }
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Set scaled radii
+
+   @param scaledRadii  scaledRadii
+
+   @return SimTKOpenMMCommon::DefaultReturn always
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setScaledRadii( const RealOpenMMVector& scaledRadii ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuObc::setScaledRadii";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _ownScaledRadii && _scaledRadii != NULL ){
+      delete[] _scaledRadii;
+   }
+   _ownScaledRadii = true;
+   _scaledRadii    = new RealOpenMM[getNumberOfAtoms()];
+   for( int ii = 0; ii < (int) scaledRadii.size(); ii++ ){
+      _scaledRadii[ii] = scaledRadii[ii];
+   }
+
+   return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Return gamma parameters
+   If not previously set, allocate space
+
+   @return array 
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM* GBVIParameters::getGammaParameters( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::getGammaParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _gammaParameters == NULL ){
+      GBVIParameters* localThis = const_cast<GBVIParameters* const>(this);
+      localThis->_gammaParameters    = new RealOpenMM[getNumberOfAtoms()];
+      localThis->_ownGammaParameters = true;
+      memset( _gammaParameters, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
+   }   
+   return _gammaParameters;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set flag indicating whether scale factors array should be deleted
+
+   @param ownGammaParameters   flag indicating whether gamma parameter
+                               array should be deleted
+
+   @return SimTKOpenMMCommon::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setOwnGammaParameters( int ownGammaParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::setOwnScaleFactors";
+
+   // ---------------------------------------------------------------------------------------
+
+   _ownGammaParameters = ownGammaParameters;
+
+   return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set gamma parameters
+
+   @param gammas  gamma parameters
+
+   @return SimTKOpenMMCommon::DefaultReturn always
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setGammaParameters( RealOpenMM* gammas ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuObc::setGammas";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _ownGammaParameters && _gammaParameters != gammas ){
+      delete[] _gammaParameters;
+      _ownGammaParameters = false;
+   }
+
+   _gammaParameters = gammas;
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+#if RealOpenMMType == 2
+
+/**---------------------------------------------------------------------------------------
+
+   Set gamma parameters
+
+   @param gammas  gammas
+
+   @return SimTKOpenMMCommon::DefaultReturn always
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setGammaParameters( float* gammas ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::setGammas";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _gammaParameters == NULL ){
+      _gammaParameters    = new RealOpenMM[getNumberOfAtoms()];
+      _ownGammaParameters = true;
+   }   
+   for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
+      _gammaParameters[ii] = (RealOpenMM) gammas[ii];
+   }
+
+   return SimTKOpenMMCommon::DefaultReturn;
+
+}
+
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Set gamma parameters
+
+   @param gammas  gammas
+
+   @return SimTKOpenMMCommon::DefaultReturn always
+
+   --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setGammaParameters( const RealOpenMMVector& gammas ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nCpuObc::setGammas";
+
+   // ---------------------------------------------------------------------------------------
+
+   if( _ownGammaParameters && _gammaParameters != NULL ){
+      delete[] _gammaParameters;
+   }
+   _ownGammaParameters = true;
+
+   _gammaParameters    = new RealOpenMM[getNumberOfAtoms()];
+   for( int ii = 0; ii < (int) gammas.size(); ii++ ){
+      _gammaParameters[ii] = gammas[ii];
+   }
+
+   return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Get string w/ state
+   
+   @param title               title (optional)
+      
+   @return string
+      
+   --------------------------------------------------------------------------------------- */
+
+std::string GBVIParameters::getStateString( const char* title ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::getStateString";
+
+   // ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   message << ImplicitSolventParameters::getStateString( title );
+
+   std::string tab           = getStringTab();
+
+   return message.str();
+
+}
+
+/**---------------------------------------------------------------------------------------
+            
+   Return zero value if all parameters set; else return nonzero
+         
+   @return ready status
+            
+   --------------------------------------------------------------------------------------- */
+ 
+int GBVIParameters::isNotReady( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\nGBVIParameters::isNotReady";
+
+   // ---------------------------------------------------------------------------------------
+
+   int isReady = ImplicitSolventParameters::isNotReady();
+
+   int errors  = 0;
+   std::stringstream message;
+   message << methodName;
+
+   const RealOpenMM* scaledRadii = getScaledRadii();
+   if( scaledRadii == NULL || scaledRadii[0] <= 0.0 ){
+      errors++;
+      message << "\n   scaledRadii are not set";
+   }
+
+   const RealOpenMM* gamma = getGammaParameters();
+   if( gamma == NULL ){
+      errors++;
+      message << "\n   gamma parameters are not set";
+   }
+
+   // check scale factors are in correct units
+
+   RealOpenMM average, stdDev, maxValue, minValue;
+   int minIndex, maxIndex;
+   SimTKOpenMMUtilities::getArrayStatistics( getNumberOfAtoms(), scaledRadii, &average,
+                                             &stdDev, &minValue, &minIndex,
+                                             &maxValue, &maxIndex );
+
+   if( average < 0.3 || average > 2.0 || minValue < 0.1 ){
+      errors++;
+      message << "\n   scale factors for atomic radii appear not to be set correctly -- radii should be in nm";
+      message << "\n   average radius=" << average << " min radius=" << minValue << " at atom index=" << minIndex;
+   }   
+
+   if( errors ){
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   errors += isReady;
+
+   return errors;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setUseCutoff( RealOpenMM distance ) {
+
+    cutoff = true;
+    cutoffDistance = distance;
+    return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get whether to use a cutoff.
+
+     --------------------------------------------------------------------------------------- */
+
+bool GBVIParameters::getUseCutoff() {
+    return cutoff;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get the cutoff distance.
+
+     --------------------------------------------------------------------------------------- */
+
+RealOpenMM GBVIParameters::getCutoffDistance() {
+    return cutoffDistance;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Set the force to use periodic boundary conditions.  This requires that a cutoff has
+     also been set, and the smallest side of the periodic box is at least twice the cutoff
+     distance.
+
+     @param boxSize             the X, Y, and Z widths of the periodic box
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+int GBVIParameters::setPeriodic( RealOpenMM* boxSize ) {
+
+    assert(cutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+    return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get whether to use periodic boundary conditions.
+
+     --------------------------------------------------------------------------------------- */
+
+bool GBVIParameters::getPeriodic() {
+    return periodic;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get the periodic box dimension
+
+     --------------------------------------------------------------------------------------- */
+
+const RealOpenMM* GBVIParameters::getPeriodicBox() {
+    return periodicBoxSize;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get tau prefactor
+
+   @return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM GBVIParameters::getTau( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nGBVIParameters::getTau:";
+
+   static const RealOpenMM zero = 0.0;
+   static const RealOpenMM one  = 1.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   RealOpenMM tau;
+   if( getSoluteDielectric() != zero && getSolventDielectric() != zero ){
+      tau = (one/getSoluteDielectric()) - (one/getSolventDielectric());
+   } else {
+      tau = zero;
+   }   
+
+   return tau;
+}

platforms/reference/src/gbsa/GBVIParameters.h

+
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __GBVIParameters_H__
+#define __GBVIParameters_H__
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "ImplicitSolventParameters.h"
+
+// ---------------------------------------------------------------------------------------
+
+class GBVIParameters : public ImplicitSolventParameters {
+
+   public:
+
+       static const std::string ParameterFileName;
+
+   private:
+
+      // scaled radii
+
+      int _ownScaledRadii;
+      RealOpenMM* _scaledRadii;
+
+      // gamma parameters
+      int _ownGammaParameters;
+      RealOpenMM* _gammaParameters;
+
+      // cutoff and periodic boundary conditions
+      
+      bool cutoff;
+      bool periodic;
+      RealOpenMM periodicBoxSize[3];
+      RealOpenMM cutoffDistance;
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         GBVIParameters constructor (Simbios) 
+      
+         @param numberOfAtoms       number of atoms
+      
+         --------------------------------------------------------------------------------------- */
+      
+       GBVIParameters( int numberOfAtoms );
+
+      /**---------------------------------------------------------------------------------------
+      
+         GBVIParameters destructor (Simbios) 
+      
+         --------------------------------------------------------------------------------------- */
+      
+       ~GBVIParameters( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Return scaled radii
+      
+         @return array
+      
+         --------------------------------------------------------------------------------------- */
+      
+      const RealOpenMM* getScaledRadii( void ) const;
+        
+      /**---------------------------------------------------------------------------------------
+      
+         Return scaled radii
+      
+         @return array
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setScaledRadii( RealOpenMM* scaledRadii );
+#if RealOpenMMType == 2
+      int setScaledRadii( float* scaledRadii );
+#endif
+      int setScaledRadii( const RealOpenMMVector& scaledRadii );
+        
+      /**---------------------------------------------------------------------------------------
+      
+         Set flag indicating whether scaled radii array should be deleted
+      
+         @param ownScaledRadiusFactors flag indicating whether scaled radii
+                                       array should be deleted
+
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setOwnScaledRadii( int ownScaledRadii );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get AtomicRadii array w/ dielectric offset applied
+      
+         @return array of atom volumes
+      
+         --------------------------------------------------------------------------------------- */
+
+      RealOpenMM* getAtomicRadii( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set AtomicRadii array
+      
+         @param atomicRadii array of atomic radii
+      
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setAtomicRadii( RealOpenMM* atomicRadii );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set AtomicRadii array
+      
+         @param atomicRadii vector of atomic radii
+      
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setAtomicRadii( const RealOpenMMVector& atomicRadii );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set flag indicating whether gamma parameter array should be deleted
+      
+         @param ownGammaParameters   flag indicating whether gamma parameter
+                                     array should be deleted
+
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setOwnGammaParameters( int ownGammaParameters );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get GammaParameters array
+      
+         @return array of gamma values
+      
+         --------------------------------------------------------------------------------------- */
+
+      RealOpenMM* getGammaParameters( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set GammaParameters array
+      
+         @param gammaParameters    array of gamma parameters
+      
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setGammaParameters( RealOpenMM* gammaParameters );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set GammaParameters array
+      
+         @param gammaParameters   array of gamma parameters
+      
+         @return SimTKOpenMMCommon::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+
+      int setGammaParameters( const RealOpenMMVector& gammaParameters );
+
+      /**---------------------------------------------------------------------------------------
+            
+         Get string w/ state
+         
+         @param title               title (optional)
+            
+         @return string
+            
+         --------------------------------------------------------------------------------------- */
+      
+      std::string getStateString( const char* title ) const;
+
+      /**---------------------------------------------------------------------------------------
+                  
+         Return zero value if all parameters set; else return nonzero
+               
+         @return ready status
+                  
+         --------------------------------------------------------------------------------------- */
+       
+      int isNotReady( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use a cutoff.
+
+         @param distance            the cutoff distance
+
+         @return SimTKOpenMMCommon::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int setUseCutoff( RealOpenMM distance );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get whether to use a cutoff.
+
+         --------------------------------------------------------------------------------------- */
+
+      bool getUseCutoff();
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the cutoff distance.
+
+         --------------------------------------------------------------------------------------- */
+
+      RealOpenMM getCutoffDistance();
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use periodic boundary conditions.  This requires that a cutoff has
+         already been set, and the smallest side of the periodic box is at least twice the cutoff
+         distance.
+
+         @param boxSize             the X, Y, and Z widths of the periodic box
+
+         @return SimTKOpenMMCommon::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int setPeriodic( RealOpenMM* boxSize );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get whether to use periodic boundary conditions.
+
+         --------------------------------------------------------------------------------------- */
+
+      bool getPeriodic();
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the periodic box dimension
+
+         --------------------------------------------------------------------------------------- */
+
+      const RealOpenMM* getPeriodicBox();
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get tau prefactor
+      
+         @return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric
+      
+         --------------------------------------------------------------------------------------- */
+      
+      RealOpenMM getTau( void ) const;
+      
+
+};
+   
+#endif // __GBVIParameters_H__

platforms/reference/tests/TestReferenceGBVIForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the reference implementation of GBVIForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "HarmonicBondForce.h"
+#include "GBVIForce.h"
+#include "GBSAOBCForce.h"
+#include "System.h"
+#include "LangevinIntegrator.h"
+#include "NonbondedForce.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testSingleParticle() {
+    ReferencePlatform platform;
+    System system(1, 0);
+    system.setParticleMass(0, 2.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+
+    GBVIForce* forceField = new GBVIForce(1);
+
+    double charge         = 0.0;
+    double radius         = 0.15;
+    double gamma          = 1.0;
+    forceField->setParticleParameters(0, charge, radius, gamma);
+    system.addForce(forceField);
+
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(1);
+    positions[0] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+
+    double bornRadius     = radius; 
+    double eps0           = EPSILON0;
+    double tau            = (1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric());
+
+    double bornEnergy     = (-charge*charge/(8*PI_M*eps0))*tau/bornRadius;
+    double nonpolarEnergy = -CAL2JOULE*gamma*tau*std::pow( radius/bornRadius, 3.0);
+
+    double expectedE      = (bornEnergy+nonpolarEnergy); 
+    double obtainedE      = state.getPotentialEnergy(); 
+    double diff           = fabs( obtainedE - expectedE );
+    (void) fprintf( stderr, "testSingleParticle expected=%14.6e obtained=%14.6e diff=%14.6e breakdown:[%14.6e %14.6e]\n",
+                    expectedE, obtainedE, diff, bornEnergy, nonpolarEnergy );
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+}
+
+/*
+void testCutoffAndPeriodic() {
+    ReferencePlatform platform;
+    System system(2, 0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* gbsa = new GBSAOBCForce(2);
+    NonbondedForce* nonbonded = new NonbondedForce(2, 0);
+    gbsa->setParticleParameters(0, -1, 0.15, 1);
+    nonbonded->setParticleParameters(0, -1, 1, 0);
+    gbsa->setParticleParameters(1, 1, 0.15, 1);
+    nonbonded->setParticleParameters(1, 1, 1, 0);
+    const double cutoffDistance = 3.0;
+    const double boxSize = 10.0;
+    nonbonded->setCutoffDistance(cutoffDistance);
+    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+
+    // Calculate the forces for both cutoff and periodic with two different atom positions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state2 = context.getState(State::Forces);
+    positions[1][0]+= boxSize;
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state4 = context.getState(State::Forces);
+
+    // All forces should be identical, exception state3 which should be zero.
+
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
+}
+*/
+
+void testEnergyEthane() {
+
+    ReferencePlatform platform;
+    const int numParticles = 8;
+    System system(numParticles, 0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+
+    //void HarmonicBondForce::getBondParameters(int index, int& particle1, int& particle2, double& length, double& k)
+    double C_HBondDistance   = 0.1097;
+    double C_CBondDistance   = 0.1504;
+    HarmonicBondForce* bonds = new HarmonicBondForce(7);
+    bonds->setBondParameters(0, 0, 1, C_HBondDistance, 0.0); 
+    bonds->setBondParameters(1, 2, 1, C_HBondDistance, 0.0); 
+    bonds->setBondParameters(2, 3, 1, C_HBondDistance, 0.0); 
+
+    bonds->setBondParameters(3, 1, 4, C_CBondDistance, 0.0); 
+
+    bonds->setBondParameters(4, 5, 4, C_HBondDistance, 0.0); 
+    bonds->setBondParameters(5, 6, 4, C_HBondDistance, 0.0); 
+    bonds->setBondParameters(6, 7, 4, C_HBondDistance, 0.0); 
+
+    system.addForce(bonds);
+
+    double C_radius, C_gamma, C_charge, H_radius, H_gamma, H_charge;
+
+    int AM1_BCC = 1;
+    H_charge    = -0.053;
+    C_charge    = -3.0*H_charge;
+    if( AM1_BCC ){
+       C_radius =  0.180;
+       C_gamma  = -0.2863;
+       H_radius =  0.125;
+       H_gamma  =  0.2437;
+    } else {
+       C_radius =  0.215;
+       C_gamma  = -1.1087;
+       H_radius =  0.150;
+       H_gamma  =  0.1237;
+    }
+
+    int VI = 1;
+    if( VI ){
+       (void) fprintf( stderr, "Applying GB/VI\n" );
+       GBVIForce* forceField = new GBVIForce(numParticles);
+       for( int i = 0; i < numParticles; i++ ){
+          forceField->setParticleParameters(i, H_charge, H_radius, H_gamma);
+       }
+       forceField->setParticleParameters(1, C_charge, C_radius, C_gamma);
+       forceField->setParticleParameters(4, C_charge, C_radius, C_gamma);
+       system.addForce(forceField);
+    } else {
+       (void) fprintf( stderr, "Applying GBSA OBC\n" );
+       GBSAOBCForce* forceField = new GBSAOBCForce(numParticles);
+       double H_scale           =  0.85;
+       double C_scale           =  0.72;
+       for( int i = 0; i < numParticles; i++ ){
+          forceField->setParticleParameters(i, H_charge, H_radius, H_scale );
+       }
+       forceField->setParticleParameters(1, C_charge, C_radius, C_scale);
+       forceField->setParticleParameters(4, C_charge, C_radius, C_scale);
+       system.addForce(forceField);
+    }
+
+    OpenMMContext context(system, integrator, platform);
+    
+/*    
+    0.5480    1.7661    0.0000 H   0  0  0  0  0  0           0  0  0
+    0.7286    0.8978    0.6468 C   0  0  0  0  0  0           0  0  0
+    0.4974    0.0000    0.0588 H   0  0  0  0  0  0           0  0  0
+    0.0000    0.9459    1.4666 H   0  0  0  0  0  0           0  0  0
+    2.1421    0.8746    1.1615 C   0  0  0  0  0  0           0  0  0
+    2.3239    0.0050    1.8065 H   0  0  0  0  0  0           0  0  0
+    2.8705    0.8295    0.3416 H   0  0  0  0  0  0           0  0  0
+    2.3722    1.7711    1.7518 H   0  0  0  0  0  0           0  0  0
+*/
+
+    vector<Vec3> positions(numParticles);
+    positions[0] = Vec3(0.5480,    1.7661,    0.0000);
+    positions[1] = Vec3(0.7286,    0.8978,    0.6468);
+    positions[2] = Vec3(0.4974,    0.0000,    0.0588);
+    positions[3] = Vec3(0.0000,    0.9459,    1.4666);
+    positions[4] = Vec3(2.1421,    0.8746,    1.1615);
+    positions[5] = Vec3(2.3239,    0.0050,    1.8065);
+    positions[6] = Vec3(2.8705,    0.8295,    0.3416);
+    positions[7] = Vec3(2.3722,    1.7711,    1.7518);
+    context.setPositions(positions);
+
+    State state = context.getState(State::Forces | State::Energy);
+    (void) fprintf( stderr, "Energy %.4e\n", state.getPotentialEnergy() );
+    
+    // Take a small step in the direction of the energy gradient.
+    
+    double norm        = 0.0;
+    double forceSum[3] = { 0.0, 0.0, 0.0 };
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f  = state.getForces()[i];
+        (void) fprintf( stderr, "F %d [%14.6e %14.6e %14.6e]\n", i, f[0], f[1], f[2] );
+        norm   += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+        forceSum[0] += f[0];
+        forceSum[1] += f[1];
+        forceSum[2] += f[2];
+    }
+    norm               = std::sqrt(norm);
+    (void) fprintf( stderr, "Fsum [%14.6e %14.6e %14.6e] norm=%14.6e\n", forceSum[0], forceSum[1], forceSum[2], norm );
+
+    const double delta = 1e-4;
+    double step = delta/norm;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = state.getForces()[i];
+        positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+    }
+    context.setPositions(positions);
+    
+    State state2 = context.getState(State::Energy);
+
+    (void) fprintf( stderr, "Energies %.8e %.8e\n", state.getPotentialEnergy(), state2.getPotentialEnergy() );
+
+    // See whether the potential energy changed by the expected amount.
+    
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 0.01)
+}
+
+void testEnergyTwoParticle() {
+
+    ReferencePlatform platform;
+    const int numParticles = 2;
+    System system(numParticles, 0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+
+    //void HarmonicBondForce::getBondParameters(int index, int& particle1, int& particle2, double& length, double& k)
+    double C_HBondDistance   = 3.0;
+    HarmonicBondForce* bonds = new HarmonicBondForce(1);
+    bonds->setBondParameters(0, 0, 1, C_HBondDistance, 0.0); 
+    system.addForce(bonds);
+
+    double C_radius, C_gamma, C_charge, H_radius, H_gamma, H_charge;
+/*
+    H_charge    = -1.0;
+    C_charge    =  1.0;
+
+    H_gamma     =  1.0;
+    C_gamma     =  1.0;
+
+    H_radius    =  1.0;
+    C_radius    =  1.0;
+*/ 
+    H_charge    = -0.5;
+    C_charge    =  0.5;
+
+    H_gamma     =  0.5;
+    C_gamma     =  0.5;
+
+    H_radius    =  1.5;
+    C_radius    =  1.5;
+ 
+    int VI = 1;
+    if( VI ){
+       (void) fprintf( stderr, "Applying GB/VI\n" );
+       GBVIForce* forceField = new GBVIForce(numParticles);
+       forceField->setParticleParameters(0, H_charge, H_radius, H_gamma);
+       forceField->setParticleParameters(1, C_charge, C_radius, C_gamma);
+       system.addForce(forceField);
+    } else {
+       (void) fprintf( stderr, "Applying GBSA OBC\n" );
+       GBSAOBCForce* forceField = new GBSAOBCForce(numParticles);
+       forceField->setParticleParameters(0, H_charge, H_radius, 1.0);
+       forceField->setParticleParameters(1, C_charge, C_radius, 1.0);
+       system.addForce(forceField);
+    }
+
+    OpenMMContext context(system, integrator, platform);
+    
+    vector<Vec3> positions(numParticles);
+    positions[0] = Vec3(         0.0000,    0.0000,    0.0000);
+    positions[1] = Vec3(C_HBondDistance,    0.0000,    0.0000);
+    context.setPositions(positions);
+
+    State state = context.getState(State::Forces | State::Energy);
+    
+    // Take a small step in the direction of the energy gradient.
+    
+    double norm        = 0.0;
+    double forceSum[3] = { 0.0, 0.0, 0.0 };
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f  = state.getForces()[i];
+        (void) fprintf( stderr, "F %d [%14.6e %14.6e %14.6e]\n", i, f[0], f[1], f[2] );
+        norm   += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+        forceSum[0] += f[0];
+        forceSum[1] += f[1];
+        forceSum[2] += f[2];
+    }
+    norm               = std::sqrt(norm);
+    (void) fprintf( stderr, "Fsum [%14.6e %14.6e %14.6e] norm=%14.6e\n", forceSum[0], forceSum[1], forceSum[2], norm );
+
+    const double delta = 1e-4;
+    double step = delta/norm;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = state.getForces()[i];
+        positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+    }
+    context.setPositions(positions);
+    
+    State state2 = context.getState(State::Energy);
+
+    double diff = fabs( norm - (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta );
+    (void) fprintf( stderr, "Energies %14.6e %14.6e diff=%14.6e [%14.6e %14.6e]\n",
+                    state.getPotentialEnergy(), state2.getPotentialEnergy(), diff, norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta );
+
+    // See whether the potential energy changed by the expected amount.
+    
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 0.01)
+}
+
+int main() {
+    try {
+/*
+        testSingleParticle();
+        testCutoffAndPeriodic();
+        testForce();
+        testEnergyEthane();
+        testEnergy1();
+*/
+//        testSingleParticle();
+        testEnergyTwoParticle();
+        testEnergyEthane();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}