Removed SASAForce -- not currently used

Added checks for nonbonded methods
Removed unused functionality from GK force
Moved non-accessor methods from Force classes to ForceImpl classes

plugins/amoeba/openmmapi/src/AmoebaSASAForceImpl.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/AmoebaSASAForceImpl.h"
-#include "openmm/internal/ContextImpl.h"
-#include "amoebaKernels.h"
-
-using namespace OpenMM;
-using std::vector;
-
-AmoebaSASAForceImpl::AmoebaSASAForceImpl(AmoebaSASAForce& owner) : owner(owner) {
-}
-
-void AmoebaSASAForceImpl::initialize(ContextImpl& context) {
-    kernel = context.getPlatform().createKernel(CalcAmoebaSASAForceKernel::Name(), context);
-    if (owner.getNumParticles() != context.getSystem().getNumParticles())
-        throw OpenMMException("AmoebaSASAForce must have exactly as many particles as the System it belongs to.");
-    dynamic_cast<CalcAmoebaSASAForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
-}
-
-double AmoebaSASAForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    return dynamic_cast<CalcAmoebaSASAForceKernel&>(kernel.getImpl()).execute(context, includeForces, includeEnergy);
-}
-
-std::vector<std::string> AmoebaSASAForceImpl::getKernelNames() {
-    std::vector<std::string> names;
-    names.push_back(CalcAmoebaSASAForceKernel::Name());
-    return names;
-}

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -46,6 +46,11 @@ AmoebaVdwForceImpl::~AmoebaVdwForceImpl() {
 }
 
 void AmoebaVdwForceImpl::initialize(ContextImpl& context) {
+    System& system = context.getSystem();
+
+    if (owner.getNumParticles() != system.getNumParticles())
+        throw OpenMMException("AmoebaVdwForce must have exactly as many particles as the System it belongs to.");
+
     kernel = context.getPlatform().createKernel(CalcAmoebaVdwForceKernel::Name(), context);
     dynamic_cast<CalcAmoebaVdwForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }

plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForce.cpp

@@ -63,79 +63,6 @@ void AmoebaWcaDispersionForce::setParticleParameters(int particleIndex, double r
     parameters[particleIndex].epsilon         = epsilon;
 }
 
-void AmoebaWcaDispersionForce::getMaximumDispersionEnergy(int particleIndex, double& maxDispersionEnergy ) const {
-
-    const double pi  = 3.1415926535897;
-
-    // from last loop in subroutine knp in ksolv.f
-
-    double rdisp, epsi;
-    getParticleParameters( particleIndex, rdisp, epsi );
-    if( epsi <= 0.0 || rdisp <= 0.0 ){
-        maxDispersionEnergy = 0.0;
-        return;
-    }
-    double rmini     = rdisp;
-    rdisp           += getDispoff();
-    double epso      = getEpso();
-    double emixo     = std::sqrt(epso) + std::sqrt(epsi);
-           emixo     = 4.0*epso*epsi/(emixo*emixo);
-
-    double rmino     = getRmino();
-    double rmino2    = rmino*rmino;
-    double rmini2    = rmini*rmini;
-    double rmixo     = 2.0*(rmino2*rmino + rmini2*rmini) / (rmino2 + rmini2);
-
-    double rmixo3    = rmixo*rmixo*rmixo;
-    double rmixo7    = rmixo*rmixo3*rmixo3;
-    double ao        = emixo*rmixo7;
-
-    double epsh      = getEpsh();
-    double emixh     = std::sqrt(epsh) + std::sqrt(epsi);
-           emixh     = 4.0*epsh*epsi/(emixh*emixh);
-    double rminh     = getRminh();
-    double rminh2    = rminh*rminh;
-    double rmixh     = rminh*rminh + rmini2;
-           rmixh     = 2.0 * (rminh2*rminh + rmini2*rmini) / (rminh2 + rmini2);
-    double rmixh3    = rmixh*rmixh*rmixh;
-    double rmixh7    = rmixh3*rmixh3*rmixh;
-    double ah        = emixh*rmixh7;
-
-    double rdisp3    = rdisp*rdisp*rdisp;
-    double rdisp7    = rdisp*rdisp3*rdisp3;
-    double rdisp11   = rdisp7*rdisp3*rdisp;
-
-    double cdisp;
-    if( rdisp < rmixh) {
-        cdisp = -4.0*pi*emixh*(rmixh3-rdisp3)/3.0 - emixh*18.0/11.0*rmixh3*pi;
-    } else {
-        cdisp = 2.0*pi*(2.0*rmixh7-11.0*rdisp7)*ah/ (11.0*rdisp11);
-    }
-    cdisp *= 2.0;
-    if (rdisp < rmixo ) {
-        cdisp -= 4.0*pi*emixo*(rmixo3-rdisp3)/3.0;
-        cdisp -= emixo*18.0/11.0*rmixo3*pi;
-    } else { 
-        cdisp += 2.0*pi*(2.0*rmixo7-11.0*rdisp7) * ao/(11.0*rdisp11);
-    }
-    maxDispersionEnergy = getSlevy()*getAwater()*cdisp;
-//    (void) fprintf( stderr,"Wca %5d %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e\n",
-//                    particleIndex, rdisp,rmini,epsi, emixh,rmixh,emixo,rmixo,cdisp );
-
-    return;
-}
-
-double AmoebaWcaDispersionForce::getTotalMaximumDispersionEnergy( void ) const {
-
-    double totalMaximumDispersionEnergy = 0.0;
-    for( int ii = 0; ii < getNumParticles(); ii++ ){
-        double maximumDispersionEnergy;
-        getMaximumDispersionEnergy( ii, maximumDispersionEnergy );
-        totalMaximumDispersionEnergy += maximumDispersionEnergy;
-    }
-    return totalMaximumDispersionEnergy;
-}
-
 double AmoebaWcaDispersionForce::getEpso( void ) const {
     return epso;
 }

plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForceImpl.cpp

@@ -32,6 +32,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "internal/AmoebaWcaDispersionForceImpl.h"
 #include "amoebaKernels.h"
+#include <cmath>
 
 using namespace OpenMM;
 
@@ -46,6 +47,11 @@ AmoebaWcaDispersionForceImpl::~AmoebaWcaDispersionForceImpl() {
 }
 
 void AmoebaWcaDispersionForceImpl::initialize(ContextImpl& context) {
+ 
+    System& system = context.getSystem();
+    if (owner.getNumParticles() != system.getNumParticles())
+        throw OpenMMException("AmoebaWcaDispersionForce must have exactly as many particles as the System it belongs to.");
+
     kernel = context.getPlatform().createKernel(CalcAmoebaWcaDispersionForceKernel::Name(), context);
     dynamic_cast<CalcAmoebaWcaDispersionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
@@ -53,6 +59,79 @@ void AmoebaWcaDispersionForceImpl::initialize(ContextImpl& context) {
 double AmoebaWcaDispersionForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return dynamic_cast<CalcAmoebaWcaDispersionForceKernel&>(kernel.getImpl()).execute(context, includeForces, includeEnergy);
 }
+void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( const AmoebaWcaDispersionForce& force, int particleIndex, double& maxDispersionEnergy ) {
+
+    const double pi  = 3.1415926535897;
+
+    // from last loop in subroutine knp in ksolv.f
+
+    double rdisp, epsi;
+    force.getParticleParameters( particleIndex, rdisp, epsi );
+    if( epsi <= 0.0 || rdisp <= 0.0 ){
+        maxDispersionEnergy = 0.0;
+        return;
+    }
+    double rmini     = rdisp;
+    rdisp           += force.getDispoff();
+    double epso      = force.getEpso();
+    double emixo     = std::sqrt(epso) + std::sqrt(epsi);
+           emixo     = 4.0*epso*epsi/(emixo*emixo);
+
+    double rmino     = force.getRmino();
+    double rmino2    = rmino*rmino;
+    double rmini2    = rmini*rmini;
+    double rmixo     = 2.0*(rmino2*rmino + rmini2*rmini) / (rmino2 + rmini2);
+
+    double rmixo3    = rmixo*rmixo*rmixo;
+    double rmixo7    = rmixo*rmixo3*rmixo3;
+    double ao        = emixo*rmixo7;
+
+    double epsh      = force.getEpsh();
+    double emixh     = std::sqrt(epsh) + std::sqrt(epsi);
+           emixh     = 4.0*epsh*epsi/(emixh*emixh);
+    double rminh     = force.getRminh();
+    double rminh2    = rminh*rminh;
+    double rmixh     = rminh*rminh + rmini2;
+           rmixh     = 2.0 * (rminh2*rminh + rmini2*rmini) / (rminh2 + rmini2);
+    double rmixh3    = rmixh*rmixh*rmixh;
+    double rmixh7    = rmixh3*rmixh3*rmixh;
+    double ah        = emixh*rmixh7;
+
+    double rdisp3    = rdisp*rdisp*rdisp;
+    double rdisp7    = rdisp*rdisp3*rdisp3;
+    double rdisp11   = rdisp7*rdisp3*rdisp;
+
+    double cdisp;
+    if( rdisp < rmixh) {
+        cdisp = -4.0*pi*emixh*(rmixh3-rdisp3)/3.0 - emixh*18.0/11.0*rmixh3*pi;
+    } else {
+        cdisp = 2.0*pi*(2.0*rmixh7-11.0*rdisp7)*ah/ (11.0*rdisp11);
+    }
+    cdisp *= 2.0;
+    if (rdisp < rmixo ) {
+        cdisp -= 4.0*pi*emixo*(rmixo3-rdisp3)/3.0;
+        cdisp -= emixo*18.0/11.0*rmixo3*pi;
+    } else {
+        cdisp += 2.0*pi*(2.0*rmixo7-11.0*rdisp7) * ao/(11.0*rdisp11);
+    }
+    maxDispersionEnergy = force.getSlevy()*force.getAwater()*cdisp;
+
+//    (void) fprintf( stderr,"Wca %5d %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e\n",
+//                    particleIndex, rdisp,rmini,epsi, emixh,rmixh,emixo,rmixo,cdisp );
+
+    return;
+}
+
+double AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy( const AmoebaWcaDispersionForce& force ){
+
+    double totalMaximumDispersionEnergy = 0.0;
+    for( int ii = 0; ii < force.getNumParticles(); ii++ ){
+        double maximumDispersionEnergy;
+        getMaximumDispersionEnergy( force, ii, maximumDispersionEnergy );
+        totalMaximumDispersionEnergy += maximumDispersionEnergy;
+    }
+    return totalMaximumDispersionEnergy;
+}
 
 std::vector<std::string> AmoebaWcaDispersionForceImpl::getKernelNames() {
     std::vector<std::string> names;

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernelFactory.cpp

@@ -49,7 +49,6 @@ extern "C" void OPENMMCUDA_EXPORT registerKernelFactories() {
              platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaSASAForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory);
              platform.registerKernelFactory(CalcAmoebaWcaDispersionForceKernel::Name(), factory);
         }
@@ -124,9 +123,6 @@ KernelImpl* AmoebaCudaKernelFactory::createKernelImpl(std::string name, const Pl
     if (name == CalcAmoebaGeneralizedKirkwoodForceKernel::Name())
         return new CudaCalcAmoebaGeneralizedKirkwoodForceKernel(name, platform, *amoebaCudaData, context.getSystem());
 
-    if (name == CalcAmoebaSASAForceKernel::Name())
-        return new CudaCalcAmoebaSASAForceKernel(name, platform, *amoebaCudaData, context.getSystem());
-
     if (name == CalcAmoebaVdwForceKernel::Name())
         return new CudaCalcAmoebaVdwForceKernel(name, platform, *amoebaCudaData, context.getSystem());
 

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -25,12 +25,13 @@
  * -------------------------------------------------------------------------- */
 
 #include "AmoebaCudaKernels.h"
-#include "openmm/LangevinIntegrator.h"
-#include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
 #include "kernels/amoebaGpuTypes.h"
 #include "kernels/cudaKernels.h"
 #include "kernels/amoebaCudaKernels.h"
+#include "internal/AmoebaMultipoleForceImpl.h"
+#include "internal/AmoebaWcaDispersionForceImpl.h"
+
 #include <cmath>
 #ifdef _MSC_VER
 #include <windows.h>
@@ -509,7 +510,7 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
     polarizationCovalentList.push_back( AmoebaMultipoleForce::PolarizationCovalent14 );
 
     std::vector<int> covalentDegree;
-    force.getCovalentDegree( covalentDegree );
+    AmoebaMultipoleForceImpl::getCovalentDegree( force, covalentDegree );
     int dipoleIndex      = 0;
     int quadrupoleIndex  = 0;
     int maxCovalentRange = 0;
@@ -556,13 +557,13 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
         multipoleAtomCovalentInfo[i] = covalentLists;
 
         int minCovalentIndex, maxCovalentIndex;
-        force.getCovalentRange( i, covalentList, &minCovalentIndex, &maxCovalentIndex );
+        AmoebaMultipoleForceImpl::getCovalentRange( force, i, covalentList, &minCovalentIndex, &maxCovalentIndex );
         minCovalentIndices[i] = minCovalentIndex;
         if( maxCovalentRange < (maxCovalentIndex - minCovalentIndex) ){
             maxCovalentRange = maxCovalentIndex - minCovalentIndex;
         }
 
-        force.getCovalentRange( i, polarizationCovalentList, &minCovalentIndex, &maxCovalentIndex );
+        AmoebaMultipoleForceImpl::getCovalentRange( force, i, polarizationCovalentList, &minCovalentIndex, &maxCovalentIndex );
         minCovalentPolarizationIndices[i] = minCovalentIndex;
         if( maxCovalentRange < (maxCovalentIndex - minCovalentIndex) ){
             maxCovalentRange = maxCovalentIndex - minCovalentIndex;
@@ -632,43 +633,6 @@ double CudaCalcAmoebaGeneralizedKirkwoodForceKernel::execute(ContextImpl& contex
     return 0.0;
 }
 
-/* -------------------------------------------------------------------------- *
- *                       AmoebaSASA                                           *
- * -------------------------------------------------------------------------- */
-
-CudaCalcAmoebaSASAForceKernel::CudaCalcAmoebaSASAForceKernel(std::string name, const Platform& platform, AmoebaCudaData& data, System& system) :
-      CalcAmoebaSASAForceKernel(name, platform), data(data), system(system) {
-    data.incrementKernelCount();
-}
-
-CudaCalcAmoebaSASAForceKernel::~CudaCalcAmoebaSASAForceKernel() {
-    data.decrementKernelCount();
-}
-
-void CudaCalcAmoebaSASAForceKernel::initialize(const System& system, const AmoebaSASAForce& force) {
-
-/*
-    //data.hasAmoebaSASA = true;
-    int numParticles = system.getNumParticles();
-    std::vector<float> radii(numParticles);
-    std::vector<float> weights(numParticles);
-    for( int ii = 0; ii < numParticles; ii++ ){
-        double particleRadius, particleWeight;
-        force.getParticleParameters(ii, particleRadius, particleWeight);
-        radii[ii]    = static_cast<float>( particleRadius );
-        weights[ii]  = static_cast<float>( particleWeight);
-    }   
-
-fprintf( stderr, "\nIn CudaCalcAmoebaSASAForceKernel::initialize %d\n", numParticles );
-fflush( stderr );
-    gpuSetAmoebaSASAParameters( data.amoebaGpu, static_cast<float>(force.getProbeRadius() ), radii, weights);
-*/
-}
-
-double CudaCalcAmoebaSASAForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    return 0.0;
-}
-
 static void computeAmoebaVdwForce( AmoebaCudaData& data ) {
 
     amoebaGpuContext gpu = data.getAmoebaGpu();
@@ -771,7 +735,7 @@ void CudaCalcAmoebaWcaDispersionForceKernel::initialize(const System& system, co
         radii[ii]         = static_cast<float>( radius );
         epsilons[ii]      = static_cast<float>( epsilon );
     }   
-    float totalMaximumDispersionEnergy =  static_cast<float>( force.getTotalMaximumDispersionEnergy( ) );
+    float totalMaximumDispersionEnergy =  static_cast<float>( AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy( force ) );
     gpuSetAmoebaWcaDispersionParameters( data.getAmoebaGpu(), radii, epsilons, totalMaximumDispersionEnergy,
                                           static_cast<float>( force.getEpso( )),
                                           static_cast<float>( force.getEpsh( )),

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h

@@ -326,34 +326,6 @@ private:
     System& system;
 };
 
-/**
- * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
- */
-class CudaCalcAmoebaSASAForceKernel : public CalcAmoebaSASAForceKernel {
-public:
-    CudaCalcAmoebaSASAForceKernel(std::string name, const Platform& platform, AmoebaCudaData& data, System& system);
-    ~CudaCalcAmoebaSASAForceKernel();
-    /**
-     * Initialize the kernel.
-     * 
-     * @param system     the System this kernel will be applied to
-     * @param force      the AmoebaMultipoleForce this kernel will be used for
-     */
-    void initialize(const System& system, const AmoebaSASAForce& force);
-    /**
-     * Execute the kernel to calculate the forces and/or energy.
-     *
-     * @param context        the context in which to execute this kernel
-     * @param includeForces  true if forces should be calculated
-     * @param includeEnergy  true if the energy should be calculated
-     * @return the potential energy due to the force
-     */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-private:
-    AmoebaCudaData& data;
-    System& system;
-};
-
 /**
  * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
  */

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

@@ -65,12 +65,6 @@ extern void SetCalculateAmoebaCudaWcaDispersionSim(amoebaGpuContext gpu);
 extern void GetCalculateAmoebaCudaWcaDispersionSim(amoebaGpuContext gpu);
 extern void kCalculateAmoebaWcaDispersionForces(amoebaGpuContext amoebaGpu );
 
-// SASA
-
-extern void SetCalculateAmoebaCudaSASAForcesSim(amoebaGpuContext gpu);
-extern void GetCalculateAmoebaCudaSASAForcesSim(amoebaGpuContext gpu);
-extern void kCalculateAmoebaSASAForces(amoebaGpuContext amoebaGpu );
-
 // fixed electric field
 
 extern void SetCalculateAmoebaCudaFixedEFieldSim(amoebaGpuContext gpu);

plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp

@@ -2736,7 +2736,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
             wcaDispersionForce->setParticleParameters( ii, radius, epsilon );
 
             if( ii < static_cast<int>(maxDispersionEnergyVector.size()) ){
-                wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy );
+                AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( *wcaDispersionForce, ii, maxDispersionEnergy );
                 double tinkerValue  = maxDispersionEnergyVector[ii];
                        tinkerValue *= CalToJoule;
                 double delta        = fabs( maxDispersionEnergy - tinkerValue );
@@ -2769,7 +2769,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
             wcaDispersionForce->getParticleParameters( ii, radius, epsilon );
             (void) fprintf( log, "%8d %10.4f %10.4f", ii, radius, epsilon );
             if( ii < maxDispersionEnergyVector.size() ){
-                wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy );
+                AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( *wcaDispersionForce, ii, maxDispersionEnergy );
                 if( useOpenMMUnits )maxDispersionEnergy /= CalToJoule;
                 double delta = fabs( maxDispersionEnergy - maxDispersionEnergyVector[ii] );
                 const char* error  = (delta > 1.0e-05) ? "XXX" : "";
@@ -2791,7 +2791,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
         int errors = 0;
         for( unsigned int ii = 0; ii < arraySize && ii < maxDispersionEnergyVector.size(); ii++ ){
             double maxDispersionEnergy;
-            wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy );
+            AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( *wcaDispersionForce, ii, maxDispersionEnergy );
             if( useOpenMMUnits )maxDispersionEnergy /= CalToJoule;
             double delta = fabs( maxDispersionEnergy - maxDispersionEnergyVector[ii] );
             if( delta > 1.0e-05 ){
@@ -2810,133 +2810,6 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
     return wcaDispersionForce->getNumParticles();
 }
 
-/**---------------------------------------------------------------------------------------
-
-    Read Amoeba surface parameters
-
-    @param filePtr              file pointer to parameter file
-    @param forceFlag            flag signalling whether force is to be added to system
-                                if force == 0 || forceFlag & AMOEBA_HARMONIC_ANGLE_FORCE, then included
-    @param tokens               array of strings from first line of parameter file for this block of parameters
-    @param system               System reference
-    @param lineCount            used to track line entries read from parameter file
-    @param log                  log file pointer -- may be NULL
-
-    @return number of multipole parameters
-
-    --------------------------------------------------------------------------------------- */
-
-static int readAmoebaSurfaceParameters( FILE* filePtr, MapStringInt& forceMap, const StringVector& tokens,
-                                         System& system, int* lineCount,
-                                         MapStringVectorOfVectors& supplementary, FILE* log ){
-
-// ---------------------------------------------------------------------------------------
-
-     static const std::string methodName      = "readAmoebaSurfaceParameters";
-    
-// ---------------------------------------------------------------------------------------
-
-     // validate number of tokens
- 
-     if( tokens.size() < 1 ){
-         char buffer[1024];
-         (void) sprintf( buffer, "%s no surface entries???\n", methodName.c_str() );
-         throwException(__FILE__, __LINE__, buffer );
-         exit(-1);
-     }
- 
-     // create force
- 
-     AmoebaSASAForce* sasaForce = new AmoebaSASAForce();
-
-//globalSasaForce = sasaForce;
-
-     MapStringIntI forceActive  = forceMap.find( AMOEBA_SASA_FORCE );
-     if( forceActive != forceMap.end() && (*forceActive).second ){
-         system.addForce( sasaForce );
-         if( log ){
-             (void) fprintf( log, "Amoeba SASA force is being included.\n" );
-         }    
-     } else if( log ){
-         (void) fprintf( log, "Amoeba SASA force is not being included.\n" );
-     }
-
-     // load in parameters
-  
-     int numberOfParticles            = atoi( tokens[1].c_str() );
-     if( log ){
-         (void) fprintf( log, "%s number of sasaForce terms=%d\n", methodName.c_str(), numberOfParticles );
-     }
-     for( int ii = 0; ii < numberOfParticles; ii++ ){
-         StringVector lineTokens;
-         int isNotEof = readLine( filePtr, lineTokens, lineCount, log );
-         if( lineTokens.size() > 2 ){
-             int tokenIndex       = 0;
-             int index            = atoi( lineTokens[tokenIndex++].c_str() );
-             double radius        = atof( lineTokens[tokenIndex++].c_str() );
-             double weight        = atof( lineTokens[tokenIndex++].c_str() );
-             sasaForce->addParticle( radius, weight );
-         } else {
-             (void) fprintf( log, "%s AmoebaSASAForce tokens incomplete at line=%d\n", methodName.c_str(), *lineCount );
-             exit(-1);
-         }
-     }
- 
-     int isNotEof                   = 1;
-     int hits                       = 0;
-     while( hits < 1 ){
-        StringVector tokens;
-        isNotEof = readLine( filePtr, tokens, lineCount, log );
-        if( isNotEof && tokens.size() > 0 ){
- 
-           std::string field       = tokens[0];
-           if( field == "AmoebaSurfaceProbe" ){
-              sasaForce->setProbeRadius( atof( tokens[1].c_str() ) );
-              hits++;
-           } else {
-                 char buffer[1024];
-                 (void) sprintf( buffer, "%s read past SASA block at line=%d\n", methodName.c_str(), *lineCount );
-                 throwException(__FILE__, __LINE__, buffer );
-                 exit(-1);
-           }
-        } else {
-           char buffer[1024];
-           (void) sprintf( buffer, "%s invalid token count at line=%d?\n", methodName.c_str(), *lineCount );
-           throwException(__FILE__, __LINE__, buffer );
-           exit(-1);
-        }
-     }
-
-     // diagnostics
- 
-     if( log ){
-
-         //static const unsigned int maxPrint   = MAX_PRINT;
-         static const unsigned int maxPrint   = 15;
-         unsigned int arraySize               = static_cast<unsigned int>(sasaForce->getNumParticles());
-         //(void) fprintf( log, "%s: %u sample of AmoebaSASAForce parameters in %s units; probe radius=%10.4f\n",
-                         //methodName.c_str(), arraySize, (useOpenMMUnits ? "OpenMM" : "Amoeba"),
-         (void) fprintf( log, "%s: %u sample of AmoebaSASAForce parameters; probe radius=%10.4f\n",
-                         methodName.c_str(), arraySize,
-                         sasaForce->getProbeRadius() );
-  
-         for( unsigned int ii = 0; ii < arraySize;  ii++ ){
-             double radius, weight;
-             sasaForce->getParticleParameters( ii, radius, weight );
-             (void) fprintf( log, "%8d %10.4f %10.4f\n", ii, radius, weight );
-
-             // skip to end
-    
-             if( ii == maxPrint && (arraySize - maxPrint) > ii ){
-                 ii = arraySize - maxPrint - 1;
-             } 
-         }
-         (void) fflush( log );
-     }
- 
-     return sasaForce->getNumParticles();
-}
-
 /**---------------------------------------------------------------------------------------
 
     Read Constraints
@@ -3354,18 +3227,6 @@ static void getStringForceMap( System& system, MapStringForce& forceMap, FILE* l
             }
         }
     
-        // SASA Force
-
-        if( !hit ){
-
-            try {
-               AmoebaSASAForce& sasaForce        = dynamic_cast<AmoebaSASAForce&>(force);
-               forceMap[AMOEBA_SASA_FORCE]       = &force;
-               hit++;
-            } catch( std::bad_cast ){
-            }
-        }
-    
         // stretch bend force
 
         if( !hit ){
@@ -3557,6 +3418,10 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS
                     (void) fprintf( log, "CMMotionRemover added w/ frequency=%d at line=%d\n", frequency, lineCount );
                 }
    
+            // not used any longer -- was used in SASA force
+
+            } else if( field == "AmoebaSurfaceProbe" ){
+           
             // All forces/energy
   
             } else if( field == ALL_FORCES ){
@@ -3675,6 +3540,7 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS
                 readVec3( filePtr, tokens, forces[AMOEBA_GK_CAVITY_FORCE], &lineCount, field, log );
             } else if( field == "AmoebaGk_A_ForceAndTorque"            || 
                        field == "AmoebaGk_A_Force"                     ||
+                       field == "AmoebaSurfaceParameters"              ||
                        field == "AmoebaGk_A_DrB"                       ||
                        field == "AmoebaDBorn"                          ||
                        field == "AmoebaBorn1Force"                     ||
@@ -3702,11 +3568,6 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS
                     potentialEnergy[AMOEBA_GK_CAVITY_FORCE] = value; 
                 }
  
-            // Amoeba SASA
- 
-            } else if( field == "AmoebaSurfaceParameters" ){
-                readAmoebaSurfaceParameters( filePtr, forceMap, tokens, system, &lineCount, supplementary, log );
- 
             // Amoeba Vdw
  
             } else if( field == "AmoebaVdw14_7SigEpsTable"  || field == "AmoebaVdw14_7Reduction" ){

plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.h

@@ -47,8 +47,8 @@
 #include "AmoebaGeneralizedKirkwoodForce.h"
 #include "AmoebaVdwForce.h"
 #include "AmoebaWcaDispersionForce.h"
-#include "AmoebaSASAForce.h"
 #include "internal/windowsExport.h"
+#include "internal/AmoebaWcaDispersionForceImpl.h"
 
 #include <ctime>
 #include <vector>