Extened Free energy plugin to allow cutoffs; cleaned up code and added tests

plugins/freeEnergy/platforms/reference/src/ReferenceFreeEnergyKernels.cpp

@@ -246,11 +246,12 @@ void ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel::initialize(const System&
     // 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));
+        const NonbondedSoftcoreForce* nonbonded = dynamic_cast<const NonbondedSoftcoreForce*>(&system.getForce(i));
         if (nonbonded != NULL) {
-            if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
+            if (nonbonded->getNonbondedMethod() != NonbondedSoftcoreForce::NoCutoff){
                 obcParameters->setUseCutoff(static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
-            if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
+            }
+            if (nonbonded->getNonbondedMethod() == NonbondedSoftcoreForce::CutoffPeriodic) {
                 Vec3 boxVectors[3];
                 system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
                 RealOpenMM periodicBoxSize[3];
@@ -317,8 +318,10 @@ void ReferenceFreeEnergyCalcGBVISoftcoreForceKernel::initialize(const System& sy
     gBVIParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
     gBVIParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
 
-    if (force.getNonbondedMethod() != GBVISoftcoreForce::NoCutoff)
+    if (force.getNonbondedMethod() != GBVISoftcoreForce::NoCutoff){
         gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
+     }
+
     if (force.getNonbondedMethod() == GBVISoftcoreForce::CutoffPeriodic) {
         Vec3 boxVectors[3];
         system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);

plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.cpp

@@ -161,11 +161,11 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtom
                                                                 RealOpenMM** atomParameters, int** exclusions,
                                                                 RealOpenMM* fixedParameters, vector<RealVec>& forces,
                                                                 RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
-
+   std::vector<RealOpenMM> de( numberOfAtoms, 0.0 );
    if (cutoff) {
        for (int i = 0; i < (int) neighborList->size(); i++) {
            OpenMM::AtomPair pair = (*neighborList)[i];
-           calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy);
+           calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy, de);
        }
    } else {
 
@@ -185,7 +185,7 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtom
           for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
 
              if( exclusionIndices[jj] != ii ){
-                 calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, totalEnergy);
+                 calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, totalEnergy, de);
              }
           }
        }
@@ -206,9 +206,10 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtom
 
      --------------------------------------------------------------------------------------- */
 
-void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates,
-                        RealOpenMM** atomParameters, vector<RealVec>& forces,
-                        RealOpenMM* totalEnergy ) const {
+void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj,
+                                                               vector<RealVec>& atomCoordinates,
+                                                               RealOpenMM** atomParameters, vector<RealVec>& forces,
+                                                               RealOpenMM* totalEnergy, std::vector<RealOpenMM>& de ) const {
 
     // ---------------------------------------------------------------------------------------
 
@@ -246,21 +247,28 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, v
 
     RealOpenMM energy          = zero;
     RealOpenMM dEdR;
+    RealOpenMM dEdRx;
+    RealOpenMM dEdCol;
 
     if( minSoftCoreLJLambda < one ){
        calculateOneSoftCoreLJIxn( deltaR[0][ReferenceForce::RIndex], sig, eps, minSoftCoreLJLambda, &dEdR, &energy );
     } else {
        calculateOneLJIxn( inverseR, sig, eps, &dEdR, &energy );
     }
+    dEdRx = dEdR;
 
     // Coulomb
 
-    if (cutoff)
+    if (cutoff){
        dEdR       += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
-    else
+       dEdCol      = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
+    } else {
        dEdR       += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
+       dEdCol      = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
+    }
 
-    dEdR          *= inverseR*inverseR;
+    dEdCol *= inverseR*inverseR;
+    dEdR   *= inverseR*inverseR;
 
     // accumulate forces
 

plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.h

@@ -67,7 +67,7 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
           
       void calculateOneIxn( int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* totalEnergy ) const;
+                            RealOpenMM* totalEnergy, std::vector<RealOpenMM>& de ) const;
 
    public:
 
@@ -166,7 +166,7 @@ private:
            --------------------------------------------------------------------------------------- */
       
       void calculateOneLJIxn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
-                              RealOpenMM* dEdR, RealOpenMM* energy ) const;
+                              RealOpenMM* dEdR, RealOpenMM* energy) const;
       
         /**---------------------------------------------------------------------------------------
       

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

@@ -206,7 +206,8 @@ void CpuGBVISoftcore::computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadiu
 
    --------------------------------------------------------------------------------------- */
 
-void CpuGBVISoftcore::computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMMVector& bornRadii ){
+void CpuGBVISoftcore::computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates,
+                                        RealOpenMMVector& bornRadii ){
 
     // ---------------------------------------------------------------------------------------
 
@@ -254,6 +255,7 @@ void CpuGBVISoftcore::computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordi
                    continue;
 
                sum  += bornRadiusScaleFactors[atomJ]*CpuGBVISoftcore::getVolume( r, radiusI, scaledRadii[atomJ] );
+
             }
         }
 
@@ -267,8 +269,11 @@ void CpuGBVISoftcore::computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordi
             computeBornRadiiUsingQuinticSpline( atomicRadius3, sum, gbviParameters, 
                                                 bornRadii[atomI], &switchDeriviativeValue );
             switchDeriviative[atomI] = switchDeriviativeValue;
+
         }
+
     }
+
 }
 
 /**---------------------------------------------------------------------------------------
@@ -558,12 +563,7 @@ void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, ve
     // set energy/forces to zero
 
     RealOpenMMVector bornForces( numberOfAtoms, 0.0 );
-    std::vector<RealVec> forces( numberOfAtoms );
-    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
-        forces[atomI][0] = zero;
-        forces[atomI][1] = zero;
-        forces[atomI][2] = zero;
-    }
+    std::vector<RealVec> forces( numberOfAtoms, RealVec(0.0, 0.0, 0.0) );
 
     // ---------------------------------------------------------------------------------------
 
@@ -602,7 +602,6 @@ void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, ve
             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];
@@ -621,10 +620,10 @@ void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, ve
 
            }
            bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
+
        }
     }
 
-    // ---------------------------------------------------------------------------------------
 
     // second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
 
@@ -634,6 +633,9 @@ void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, ve
     const RealOpenMMVector& scaledRadii                 = gbviParameters->getScaledRadii();
     RealOpenMMVector& switchDeriviative                 = getSwitchDeriviative();
     const RealOpenMMVector& bornRadiusScaleFactors      = gbviParameters->getBornRadiusScaleFactors();
+
+    // ---------------------------------------------------------------------------------------
+
     for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
  
         RealOpenMM R        = atomicRadii[atomI];
@@ -706,6 +708,8 @@ void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, ve
 
     }
 
+    //printGbvi( atomCoordinates, partialCharges, bornRadii,bornForces, forces, "GBVI softcore post loop2", stderr );
+
     // apply prefactor tau = (1/diel_solute - 1/diel_solvent)
 
     RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());  
@@ -717,6 +721,68 @@ void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, ve
 
 }
 
+/**---------------------------------------------------------------------------------------
+
+    Print GB/VI parameters, radii, forces, ...
+
+    @param atomCoordinates     atomic coordinates
+    @param partialCharges      partial charges
+    @param bornRadii           Born radii (may be empty)
+    @param bornForces          Born forces (may be empty)
+    @param forces              forces (may be empty)
+    @param idString            id string (who is calling)
+    @param log                 log file
+
+    --------------------------------------------------------------------------------------- */
+
+void CpuGBVISoftcore::printGbvi( const std::vector<OpenMM::RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges,
+                                 const RealOpenMMVector& bornRadii,
+                                 const RealOpenMMVector& bornForces,
+                                 const std::vector<OpenMM::RealVec>& forces,
+                                 const std::string& idString, FILE* log ){
+
+    // ---------------------------------------------------------------------------------------
+
+    const GBVISoftcoreParameters* gbviParameters    = getGBVISoftcoreParameters();
+    const int numberOfAtoms                         = gbviParameters->getNumberOfAtoms();
+    const RealOpenMMVector& atomicRadii             = gbviParameters->getAtomicRadii();
+    const RealOpenMMVector& gammaParameters         = gbviParameters->getGammaParameters();
+    const RealOpenMMVector& bornRadiusScaleFactors  = gbviParameters->getBornRadiusScaleFactors();
+    const RealOpenMM preFactor                      = 2.0*gbviParameters->getElectricConstant();
+
+    // ---------------------------------------------------------------------------------------
+
+    const RealOpenMMVector& scaledRadii       = gbviParameters->getScaledRadii();
+    const RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
+    RealOpenMM tau                            = static_cast<RealOpenMM>(gbviParameters->getTau());  
+
+    (void) fprintf( log, "Reference Gbvi     %s atoms=%d\n", idString.c_str(), numberOfAtoms );
+    (void) fprintf( log, "    tau            %15.7e\n", tau ); 
+    (void) fprintf( log, "    scaleMethod    %d (QuinticEnum=%d)\n", 
+                    gbviParameters->getBornRadiusScalingSoftcoreMethod(), GBVISoftcoreParameters::QuinticSpline );
+    (void) fprintf( log, "    preFactor      %15.7e)\n", preFactor );
+ 
+    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+        (void) fprintf( log, "%6d r=%15.7e rSc=%15.7e swd=%15.7e lmda=%4.2f tau*gam=%15.7e q=%15.7e", atomI,
+                        atomicRadii[atomI],    
+                        scaledRadii[atomI],    
+                        switchDeriviative[atomI],    
+                        bornRadiusScaleFactors[atomI],
+                        tau*gammaParameters[atomI],    
+                        partialCharges[atomI] );
+        if( bornRadii.size() > atomI ){
+             (void) fprintf( log, " bR=%15.7e", bornRadii[atomI] );
+        }
+        if( bornForces.size() > atomI ){
+             (void) fprintf( log, " tau*bF=%15.7e", tau*bornForces[atomI] );    
+        }
+        (void) fprintf( log, "\n" );
+    }   
+
+    return;
+
+}
+
 /**---------------------------------------------------------------------------------------
 
    Use double precision 

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

@@ -303,6 +303,26 @@ class CpuGBVISoftcore {
                                               GBVISoftcoreParameters* gbviParameters, 
                                               RealOpenMM& bornRadius, RealOpenMM* switchDeriviative );
       
+    /**---------------------------------------------------------------------------------------
+    
+        Print GB/VI parameters, radii, forces, ...
+    
+        @param atomCoordinates     atomic coordinates
+        @param partialCharges      partial charges
+        @param bornRadii           Born radii (may be empty)
+        @param bornForces          Born forces (may be empty)
+        @param forces              forces (may be empty)
+        @param idString            id string (who is calling)
+        @param log                 log file
+    
+        --------------------------------------------------------------------------------------- */
+    
+    void printGbvi( const std::vector<OpenMM::RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges,
+                    const RealOpenMMVector& bornRadii,
+                    const RealOpenMMVector& bornForces,
+                    const std::vector<OpenMM::RealVec>& forces,
+                    const std::string& idString, FILE* log );
+    
 };
 
 // ---------------------------------------------------------------------------------------

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

@@ -209,11 +209,12 @@ void CpuObcSoftcore::computeBornRadii( const vector<RealVec>& atomCoordinates,
                     // this case (atom i completely inside atom j) is not considered in the original paper
                     // Jay Ponder and the authors of Tinker recognized this and
                     // worked out the details
-     
                     if( offsetRadiusI < (scaledRadiusJ - r) ){
                         term += two*( radiusIInverse - l_ij);
                     }
+
                     sum += nonPolarScaleFactors[atomJ]*term;
+
                 }
             }
         }
@@ -341,10 +342,7 @@ RealOpenMM CpuObcSoftcore::computeBornEnergyForces( vector<RealVec>& atomCoordin
     computeBornRadii( atomCoordinates,  bornRadii );
 
     RealOpenMM obcEnergy                 = zero;
-    RealOpenMMVector bornForces( numberOfAtoms );
-    for( int ii = 0; ii < numberOfAtoms; ii++ ){
-       bornForces[ii] = zero;
-    }
+    RealOpenMMVector bornForces( numberOfAtoms, 0.0 );
 
     // ---------------------------------------------------------------------------------------
 

plugins/freeEnergy/platforms/reference/tests/TestReferenceGBVISoftcoreForce.cpp

@@ -106,7 +106,7 @@ void testSingleParticle() {
     system.addParticle(2.0);
     LangevinIntegrator integrator(0, 0.1, 0.01);
 
-    GBVISoftcoreForce* forceField = new GBVISoftcoreForce;
+    GBVISoftcoreForce* forceField = new GBVISoftcoreForce();
 
     double charge         = 1.0;
     double radius         = 0.15;
@@ -130,7 +130,7 @@ void testSingleParticle() {
     double tau            = (1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric());
 
     double bornEnergy     = (-charge*charge/(8*PI_M*eps0))*tau/bornRadius;
-    double nonpolarEnergy = -0.1*gamma*tau*std::pow( radius/bornRadius, 3.0);
+    double nonpolarEnergy = -gamma*tau*std::pow( radius/bornRadius, 3.0);
 
     double expectedE      = (bornEnergy+nonpolarEnergy); 
     double obtainedE      = state.getPotentialEnergy(); 

plugins/freeEnergy/serialization/src/FreeEnergyGBVISoftcoreProxy.cpp

@@ -50,13 +50,6 @@ void GBVISoftcoreForceProxy::serialize(const void* object, SerializationNode& no
     node.setDoubleProperty("solventDielectric", force.getSolventDielectric());
     node.setDoubleProperty("quinticLwrLmtFctr", force.getQuinticLowerLimitFactor());
     node.setDoubleProperty("quinticUpprBrLmt", force.getQuinticUpperBornRadiusLimit());
-/*
-    double alpha, beta, gamma;
-    force.getTanhParameters( alpha, beta, gamma );
-    node.setDoubleProperty("alphaT", alpha);
-    node.setDoubleProperty("betaT", beta);
-    node.setDoubleProperty("gammaT", gamma);
-*/
     SerializationNode& particles = node.createChildNode("Particles");
     for (int i = 0; i < force.getNumParticles(); i++) {
         double charge, radius, gamma, bornRadiusScaleFactor;
@@ -84,13 +77,6 @@ void* GBVISoftcoreForceProxy::deserialize(const SerializationNode& node) const {
         force->setSolventDielectric(node.getDoubleProperty("solventDielectric"));
         force->setQuinticLowerLimitFactor(node.getDoubleProperty("quinticLwrLmtFctr"));
         force->setQuinticUpperBornRadiusLimit(node.getDoubleProperty("quinticUpprBrLmt"));
-/*
-        double alpha, beta, gamma;
-        alpha  = node.getDoubleProperty("alphaT");
-        beta   = node.getDoubleProperty("betaT");
-        gamma  = node.getDoubleProperty("gammaT");
-        force->setTanhParameters( alpha, beta, gamma );
-*/
         const SerializationNode& particles = node.getChildNode("Particles");
         for (int i = 0; i < (int) particles.getChildren().size(); i++) {
             const SerializationNode& particle = particles.getChildren()[i];