Reference implementation for AmoebaGeneralizedKirkwoodForce

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp

@@ -61,9 +61,7 @@ extern "C" void initAmoebaReferenceKernels() {
              platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(),        factory);
              platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(),                   factory);
              platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(),             factory);
-/*
              platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
-*/
              platform.registerKernelFactory(CalcAmoebaWcaDispersionForceKernel::Name(), factory);
         }
     }
@@ -101,10 +99,8 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
     if (name == CalcAmoebaMultipoleForceKernel::Name())
         return new ReferenceCalcAmoebaMultipoleForceKernel(name, platform, context.getSystem());
 
-/*
     if (name == CalcAmoebaGeneralizedKirkwoodForceKernel::Name())
         return new ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(name, platform, context.getSystem());
-*/
 
     if (name == CalcAmoebaWcaDispersionForceKernel::Name())
         return new ReferenceCalcAmoebaWcaDispersionForceKernel(name, platform, context.getSystem());

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -35,6 +35,7 @@
 #include "AmoebaReferenceMultipoleForce.h"
 #include "AmoebaReferenceVdwForce.h"
 #include "AmoebaReferenceWcaDispersionForce.h"
+#include "AmoebaReferenceGeneralizedKirkwoodForce.h"
 #include "openmm/internal/AmoebaTorsionTorsionForceImpl.h"
 #include "openmm/internal/AmoebaWcaDispersionForceImpl.h"
 #include "ReferencePlatform.h"
@@ -42,6 +43,7 @@
 #include "openmm/AmoebaMultipoleForce.h"
 #include "openmm/internal/AmoebaMultipoleForceImpl.h"
 #include "openmm/internal/AmoebaVdwForceImpl.h"
+#include "openmm/internal/AmoebaGeneralizedKirkwoodForceImpl.h"
 
 #include <cmath>
 #ifdef _MSC_VER
@@ -453,27 +455,76 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
     if( nonbondedMethod != 0 && nonbondedMethod != 1 ){
          throw OpenMMException("AmoebaMultipoleForce nonbonded method not recognized.\n");
     }
-
-    polarizationType = static_cast<int>(force.getPolarizationType());
-    if( polarizationType != 0 && polarizationType != 1 ){ 
-         throw OpenMMException("AmoebaMultipoleForce polarization type not recognized.\n");
-    }    
+    polarizationType = force.getPolarizationType();
 
 }
 
 double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
 
-    AmoebaReferenceMultipoleForce amoebaReferenceMultipoleForce( AmoebaReferenceMultipoleForce::NoCutoff );
-    amoebaReferenceMultipoleForce.setMutualInducedDipoleTargetEpsilon( mutualInducedTargetEpsilon );
-    amoebaReferenceMultipoleForce.setMaximumMutualInducedDipoleIterations( mutualInducedMaxIterations );
+    ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel* gkKernel = NULL;
+    for (unsigned int ii = 0; ii < context.getForceImpls().size() && gkKernel == NULL; ii++) {
+        AmoebaGeneralizedKirkwoodForceImpl* gkImpl = dynamic_cast<AmoebaGeneralizedKirkwoodForceImpl*>(context.getForceImpls()[ii]);
+        if (gkImpl != NULL) {
+            gkKernel = dynamic_cast<ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel*>(&gkImpl->getKernel().getImpl());
+        }
+    }    
+
+    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = NULL;
+    if( gkKernel ){
+
+        // amoebaReferenceGeneralizedKirkwoodForce is deleted in AmoebaReferenceGeneralizedKirkwoodMultipoleForce
+        // destructor
+
+        AmoebaReferenceGeneralizedKirkwoodForce* amoebaReferenceGeneralizedKirkwoodForce = new AmoebaReferenceGeneralizedKirkwoodForce();
+        amoebaReferenceGeneralizedKirkwoodForce->setNumParticles( gkKernel->getNumParticles() );
+        amoebaReferenceGeneralizedKirkwoodForce->setSoluteDielectric( gkKernel->getSoluteDielectric() );
+        amoebaReferenceGeneralizedKirkwoodForce->setSolventDielectric( gkKernel->getSolventDielectric() );
+        amoebaReferenceGeneralizedKirkwoodForce->setDielectricOffset( gkKernel->getDielectricOffset() );
+        amoebaReferenceGeneralizedKirkwoodForce->setProbeRadius( gkKernel->getProbeRadius() );
+        amoebaReferenceGeneralizedKirkwoodForce->setSurfaceAreaFactor( gkKernel->getSurfaceAreaFactor() );
+        amoebaReferenceGeneralizedKirkwoodForce->setIncludeCavityTerm( gkKernel->getIncludeCavityTerm() );
+        amoebaReferenceGeneralizedKirkwoodForce->setDirectPolarization( gkKernel->getDirectPolarization() );
+
+        vector<RealOpenMM> parameters; 
+        gkKernel->getAtomicRadii( parameters );
+        amoebaReferenceGeneralizedKirkwoodForce->setAtomicRadii( parameters );
 
-    RealOpenMM energy      = amoebaReferenceMultipoleForce.calculateForceAndEnergy( posData, 
-                                                                                    charges, dipoles, quadrupoles, tholes,
-                                                                                    dampingFactors, polarity, axisTypes, 
-                                                                                    multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
-                                                                                    multipoleAtomCovalentInfo, polarizationType, forceData);
+        gkKernel->getScaleFactors( parameters );
+        amoebaReferenceGeneralizedKirkwoodForce->setScaleFactors( parameters );
+
+        gkKernel->getCharges( parameters );
+        amoebaReferenceGeneralizedKirkwoodForce->setCharges( parameters );
+
+        // calculate Grycuk Born radii
+
+        amoebaReferenceGeneralizedKirkwoodForce->calculateGrycukBornRadii( posData );
+
+        amoebaReferenceMultipoleForce = new AmoebaReferenceGeneralizedKirkwoodMultipoleForce( amoebaReferenceGeneralizedKirkwoodForce );
+
+    } else {
+         amoebaReferenceMultipoleForce = new AmoebaReferenceMultipoleForce( AmoebaReferenceMultipoleForce::NoCutoff );
+    }
+
+    amoebaReferenceMultipoleForce->setMutualInducedDipoleTargetEpsilon( mutualInducedTargetEpsilon );
+    amoebaReferenceMultipoleForce->setMaximumMutualInducedDipoleIterations( mutualInducedMaxIterations );
+    AmoebaReferenceMultipoleForce::PolarizationType refPolarizationType;
+    if( polarizationType == AmoebaMultipoleForce::Mutual ){
+        refPolarizationType = AmoebaReferenceMultipoleForce::Mutual;
+    } else if( polarizationType == AmoebaMultipoleForce::Direct ){
+        refPolarizationType = AmoebaReferenceMultipoleForce::Direct;
+    } else {
+        throw OpenMMException("Polarization type not recognzied." );
+    }
+
+    RealOpenMM energy = amoebaReferenceMultipoleForce->calculateForceAndEnergy( posData, charges, dipoles, quadrupoles, tholes,
+                                                                                dampingFactors, polarity, axisTypes, 
+                                                                                multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
+                                                                                multipoleAtomCovalentInfo, refPolarizationType, forceData);
+
+    delete amoebaReferenceMultipoleForce;
 
     return static_cast<double>(energy);
 }
@@ -487,48 +538,115 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextI
     return;
 }
 
-///* -------------------------------------------------------------------------- *
-// *                       AmoebaGeneralizedKirkwood                            *
-// * -------------------------------------------------------------------------- */
-//
-//ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, System& system) : 
-//           CalcAmoebaGeneralizedKirkwoodForceKernel(name, platform), system(system) {
-//    data.incrementKernelCount();
-//}
-//
-//ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel() {
-//    data.decrementKernelCount();
-//}
-//
-//void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force) {
-//
-//    data.setHasAmoebaGeneralizedKirkwood( true );
-//
-//    int numParticles = system.getNumParticles();
-//
-//    std::vector<RealOpenMM> radius(numParticles);
-//    std::vector<RealOpenMM> scale(numParticles);
-//    std::vector<RealOpenMM> charge(numParticles);
-//
-//    for( int ii = 0; ii < numParticles; ii++ ){
-//        double particleCharge, particleRadius, scalingFactor;
-//        force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor);
-//        radius[ii]  = static_cast<RealOpenMM>( particleRadius );
-//        scale[ii]   = static_cast<RealOpenMM>( scalingFactor );
-//        charge[ii]  = static_cast<RealOpenMM>( particleCharge );
-//    }   
-//    gpuSetAmoebaObcParameters( data.getAmoebaGpu(), static_cast<RealOpenMM>(force.getSoluteDielectric() ), 
-//                               static_cast<RealOpenMM>( force.getSolventDielectric() ), 
-//                               static_cast<RealOpenMM>( force.getDielectricOffset() ), radius, scale, charge,
-//                               force.getIncludeCavityTerm(),
-//                               static_cast<RealOpenMM>( force.getProbeRadius() ), 
-//                               static_cast<RealOpenMM>( force.getSurfaceAreaFactor() ) ); 
-//}
-//
-//double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-//    // handled in computeAmoebaMultipoleForce()
-//    return 0.0;
-//}
+/* -------------------------------------------------------------------------- *
+ *                       AmoebaGeneralizedKirkwood                            *
+ * -------------------------------------------------------------------------- */
+
+ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, System& system) : 
+           CalcAmoebaGeneralizedKirkwoodForceKernel(name, platform), system(system) {
+}
+
+ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel() {
+}
+
+int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getNumParticles( void ) const {
+    return numParticles;
+}
+
+int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getIncludeCavityTerm( void ) const {
+    return includeCavityTerm;
+}
+
+int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDirectPolarization( void ) const {
+    return directPolarization;
+}
+
+RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSoluteDielectric( void ) const {
+    return soluteDielectric;
+}
+
+RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSolventDielectric( void ) const {
+    return solventDielectric;
+}
+
+RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDielectricOffset( void ) const {
+    return dielectricOffset;
+}
+
+RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getProbeRadius( void ) const {
+    return probeRadius;
+}
+
+RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSurfaceAreaFactor( void ) const {
+    return surfaceAreaFactor;
+}
+
+void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getAtomicRadii( vector<RealOpenMM>& outputAtomicRadii ) const {
+    outputAtomicRadii.resize( atomicRadii.size() );
+    copy( atomicRadii.begin(), atomicRadii.end(), outputAtomicRadii.begin() );
+}
+
+void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getScaleFactors( vector<RealOpenMM>& outputScaleFactors ) const {
+    outputScaleFactors.resize( scaleFactors.size() );
+    copy( scaleFactors.begin(), scaleFactors.end(), outputScaleFactors.begin() );
+}
+
+void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getCharges( vector<RealOpenMM>& outputCharges ) const {
+    outputCharges.resize( charges.size() );
+    copy( charges.begin(), charges.end(), outputCharges.begin() );
+}
+
+void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force) {
+
+    // check that AmoebaMultipoleForce is present
+
+    const AmoebaMultipoleForce* amoebaMultipoleForce = NULL;
+    for (int ii = 0; ii < system.getNumForces() && amoebaMultipoleForce == NULL; ii++) {
+        amoebaMultipoleForce = dynamic_cast<const AmoebaMultipoleForce*>(&system.getForce(ii));
+    }
+
+    if (amoebaMultipoleForce == NULL) {
+        throw OpenMMException("AmoebaGeneralizedKirkwoodForce requires the System to also contain an AmoebaMultipoleForce.");
+    }
+
+    if (amoebaMultipoleForce->getNonbondedMethod() != AmoebaMultipoleForce::NoCutoff ) {
+        throw OpenMMException("AmoebaGeneralizedKirkwoodForce requires the AmoebaMultipoleForce use the NoCutoff nonbonded method.");
+    }
+
+    numParticles = system.getNumParticles();
+
+    for( int ii = 0; ii < numParticles; ii++ ){
+
+        double particleCharge, particleRadius, scalingFactor;
+        force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor);
+        atomicRadii.push_back( static_cast<RealOpenMM>( particleRadius ) );
+        scaleFactors.push_back( static_cast<RealOpenMM>( scalingFactor ) );
+        charges.push_back( static_cast<RealOpenMM>( particleCharge ) );
+
+        // Make sure the charge matches the one specified by the AmoebaMultipoleForce.
+
+        double charge2, thole, damping, polarity;
+        int axisType, atomX, atomY, atomZ;
+        vector<double> dipole, quadrupole;
+        amoebaMultipoleForce->getMultipoleParameters( ii, charge2, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
+        if ( particleCharge != charge2 ){
+            throw OpenMMException("AmoebaGeneralizedKirkwoodForce and AmoebaMultipoleForce must specify the same charge for every atom.");
+        }
+
+    }   
+    includeCavityTerm  = force.getIncludeCavityTerm();
+    soluteDielectric   = static_cast<RealOpenMM>( force.getSoluteDielectric() );
+    solventDielectric  = static_cast<RealOpenMM>( force.getSolventDielectric() );
+    dielectricOffset   = static_cast<RealOpenMM>( 0.009 );
+    probeRadius        = static_cast<RealOpenMM>( force.getProbeRadius() ), 
+    surfaceAreaFactor  = static_cast<RealOpenMM>( force.getSurfaceAreaFactor() ); 
+    directPolarization = amoebaMultipoleForce->getPolarizationType() == AmoebaMultipoleForce::Direct ? 1 : 0;
+}
+
+double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    // handled in AmoebaReferenceGeneralizedKirkwoodMultipoleForce, a derived class of the class AmoebaReferenceMultipoleForce
+    return 0.0;
+}
 
 ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
        CalcAmoebaVdwForceKernel(name, platform), system(system) {

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -29,6 +29,7 @@
 
 #include "openmm/System.h"
 #include "openmm/amoebaKernels.h"
+//#include "openmm/AmoebaMultipoleForce.h"
 #include "SimTKReference/ReferenceNeighborList.h"
 #include "SimTKUtilities/SimTKOpenMMRealType.h"
 
@@ -338,7 +339,7 @@ public:
 
 private:
     int numMultipoles;
-    int polarizationType;
+    AmoebaMultipoleForce::PolarizationType polarizationType;
     std::vector<RealOpenMM> charges;
     std::vector<RealOpenMM> dipoles;
     std::vector<RealOpenMM> quadrupoles;
@@ -439,6 +440,132 @@ private:
     System& system;
 };
 
+/**
+ * This kernel is invoked to calculate the Gerneralized Kirkwood forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
+public:
+    ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, System& system);
+    ~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel();
+    /**
+     * 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 AmoebaGeneralizedKirkwoodForce& 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);
+ 
+    /**
+     *  Get include-cavity term flag 
+     *
+     *  @return includeCavityTerm
+     */
+    int getIncludeCavityTerm( void ) const;
+
+    /**
+     *  Get number of particles
+     *
+     *  @return number of particles
+     */
+    int getNumParticles( void ) const;
+
+    /**
+     *  Get directPolarization flag 
+     *
+     *  @return directPolarization
+     *
+     */
+    int getDirectPolarization( void ) const;
+
+    /**
+     *  Get solute dielectric
+     *
+     *  @return soluteDielectric
+     *
+     */
+    RealOpenMM getSoluteDielectric( void ) const;
+
+    /**
+     *  Get solvent dielectric
+     *
+     *  @return solventDielectric
+     *
+     */
+    RealOpenMM getSolventDielectric( void ) const;
+
+    /**
+     *  Get dielectric offset
+     *
+     *  @return dielectricOffset
+     *
+     */
+    RealOpenMM getDielectricOffset( void ) const;
+
+    /**
+     *  Get probeRadius
+     *
+     *  @return probeRadius
+     *
+     */
+    RealOpenMM getProbeRadius( void ) const;
+
+    /**
+     *  Get surfaceAreaFactor
+     *
+     *  @return surfaceAreaFactor
+     *
+     */
+    RealOpenMM getSurfaceAreaFactor( void ) const;
+
+    /**
+     *  Get atomic radii
+     *
+     *  @param atomicRadii vector of atomic radii
+     *
+     */
+    void getAtomicRadii( std::vector<RealOpenMM>& atomicRadii ) const;
+
+    /**
+     *  Get scale factors
+     *
+     *  @param scaleFactors vector of scale factors
+     *
+     */
+    void getScaleFactors( std::vector<RealOpenMM>& scaleFactors ) const;
+
+    /**
+     *  Get charges
+     *
+     *  @param charges vector of charges
+     *
+     */
+    void getCharges( std::vector<RealOpenMM>& charges ) const;
+
+private:
+
+    int numParticles;
+    std::vector<RealOpenMM> atomicRadii;
+    std::vector<RealOpenMM> scaleFactors;
+    std::vector<RealOpenMM> charges;
+    RealOpenMM soluteDielectric;
+    RealOpenMM solventDielectric;
+    RealOpenMM dielectricOffset;
+    RealOpenMM probeRadius;
+    RealOpenMM surfaceAreaFactor;
+    int includeCavityTerm;
+    int directPolarization;
+    System& system;
+};
+
 } // namespace OpenMM
 
 #endif /*AMOEBA_OPENMM_REFERENCE_KERNELS_H*/

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceGeneralizedKirkwoodForce.cpp

+/* 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.
+ */
+
+#include "AmoebaReferenceGeneralizedKirkwoodForce.h"
+
+using std::vector;
+using OpenMM::RealVec;
+
+AmoebaReferenceGeneralizedKirkwoodForce::AmoebaReferenceGeneralizedKirkwoodForce( ) : _numParticles(0),
+                                                                                      _includeCavityTerm(1),
+                                                                                      _directPolarization(0),
+                                                                                      _soluteDielectric(1.0),
+                                                                                      _solventDielectric(78.3),
+                                                                                      _dielectricOffset(0.009),
+                                                                                      _probeRadius(0.14),
+                                                                                      _surfaceAreaFactor(0.0054) {
+
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setNumParticles( int numParticles ){
+    _numParticles = numParticles;
+}
+
+int AmoebaReferenceGeneralizedKirkwoodForce::getNumParticles( void ) const {
+    return _numParticles;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setIncludeCavityTerm( int includeCavityTerm ){
+    _includeCavityTerm = includeCavityTerm;
+}
+
+int AmoebaReferenceGeneralizedKirkwoodForce::getIncludeCavityTerm( void ) const {
+    return _includeCavityTerm;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setDirectPolarization( int directPolarization ){
+    _directPolarization = directPolarization;
+}
+
+int AmoebaReferenceGeneralizedKirkwoodForce::getDirectPolarization( void ) const {
+    return _directPolarization;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setSoluteDielectric( RealOpenMM soluteDielectric ){
+    _soluteDielectric  = soluteDielectric;
+}
+
+RealOpenMM AmoebaReferenceGeneralizedKirkwoodForce::getSoluteDielectric( void ) const {
+    return _soluteDielectric;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setSolventDielectric( RealOpenMM solventDielectric ){
+    _solventDielectric  = solventDielectric;
+}
+
+RealOpenMM AmoebaReferenceGeneralizedKirkwoodForce::getSolventDielectric( void ) const {
+    return _solventDielectric;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setDielectricOffset( RealOpenMM dielectricOffset ){
+    _dielectricOffset  = dielectricOffset;
+}
+
+RealOpenMM AmoebaReferenceGeneralizedKirkwoodForce::getDielectricOffset( void ) const {
+    return _dielectricOffset;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setProbeRadius( RealOpenMM probeRadius ){
+    _probeRadius  = probeRadius;
+}
+
+RealOpenMM AmoebaReferenceGeneralizedKirkwoodForce::getProbeRadius( void ) const {
+    return _probeRadius;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setSurfaceAreaFactor( RealOpenMM surfaceAreaFactor ){
+    _surfaceAreaFactor  = surfaceAreaFactor;
+}
+
+RealOpenMM AmoebaReferenceGeneralizedKirkwoodForce::getSurfaceAreaFactor( void ) const {
+    return _surfaceAreaFactor;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setAtomicRadii( const vector<RealOpenMM>& atomicRadii ){
+    _atomicRadii.resize( atomicRadii.size() );
+    copy( atomicRadii.begin(), atomicRadii.end(), _atomicRadii.begin() );
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::getAtomicRadii( vector<RealOpenMM>& atomicRadii ) const {
+    atomicRadii.resize( _atomicRadii.size() );
+    copy( _atomicRadii.begin(), _atomicRadii.end(), atomicRadii.begin() );
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setScaleFactors( const vector<RealOpenMM>& scaleFactors ){
+    _scaleFactors.resize( scaleFactors.size() );
+    copy( scaleFactors.begin(), scaleFactors.end(), _scaleFactors.begin() );
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::getScaleFactors( vector<RealOpenMM>& scaleFactors ) const {
+    scaleFactors.resize( _scaleFactors.size() );
+    copy( _scaleFactors.begin(), _scaleFactors.end(), scaleFactors.begin() );
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setCharges( const vector<RealOpenMM>& charges ){
+    _charges.resize( charges.size() );
+    copy( charges.begin(), charges.end(), _charges.begin() );
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::getGrycukBornRadii( vector<RealOpenMM>& bornRadii ) const {
+    bornRadii.resize( _bornRadii.size() );
+    copy( _bornRadii.begin(), _bornRadii.end(), bornRadii.begin() );
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::calculateGrycukBornRadii( const vector<RealVec>& particlePositions ) {
+
+    const RealOpenMM zero      = 0.0;
+    const RealOpenMM one       = 1.0;
+    const RealOpenMM three     = 3.0;
+    const RealOpenMM six       = 6.0;
+    const RealOpenMM eight     = 8.0;
+    const RealOpenMM sixteen   = 16.0;
+    const RealOpenMM oneThird  = 1.0/3.0;
+    const RealOpenMM bigRadius = 1000.0;
+
+    _bornRadii.resize( _numParticles );
+    for( unsigned int ii = 0; ii < _numParticles; ii++ ){
+
+        if( _atomicRadii[ii] <= zero ){
+            _bornRadii[ii] = bigRadius;
+            continue;
+        }
+
+        RealOpenMM bornSum = zero;
+        for( unsigned int jj = 0; jj < _numParticles; jj++ ){
+
+            if( ii == jj || _atomicRadii[jj] < zero )continue;
+          
+            RealOpenMM xr       = particlePositions[jj][0] - particlePositions[ii][0];
+            RealOpenMM yr       = particlePositions[jj][1] - particlePositions[ii][1];
+            RealOpenMM zr       = particlePositions[jj][2] - particlePositions[ii][2];
+
+            RealOpenMM r2       = xr*xr + yr*yr + zr*zr;
+            RealOpenMM r        = SQRT( r2 );
+
+            RealOpenMM sk       = _atomicRadii[jj]*_scaleFactors[jj];
+            RealOpenMM sk2      = sk*sk;
+
+            if( (_atomicRadii[ii] + r) < sk ){
+                RealOpenMM lik       = _atomicRadii[ii];
+                RealOpenMM uik       = sk - r;  
+                RealOpenMM lik3      = lik*lik*lik;
+                RealOpenMM uik3      = uik*uik*uik;
+                bornSum             -= (one/uik3 - one/lik3);
+            }   
+        
+            RealOpenMM uik = r + sk; 
+            RealOpenMM lik;
+            if( (_atomicRadii[ii] + r) < sk ){
+                lik = sk - r;  
+            } else if( r < (_atomicRadii[ii] + sk) ){
+                lik = _atomicRadii[ii];
+            } else {
+                lik = r - sk; 
+            }   
+        
+            RealOpenMM l2          = lik*lik; 
+            RealOpenMM l4          = l2*l2;
+            RealOpenMM lr          = lik*r;
+            RealOpenMM l4r         = l4*r;
+        
+            RealOpenMM u2          = uik*uik;
+            RealOpenMM u4          = u2*u2;
+            RealOpenMM ur          = uik*r;
+            RealOpenMM u4r         = u4*r;
+        
+            RealOpenMM term        = (three*(r2-sk2) + six*u2 - eight*ur)/u4r - (three*(r2-sk2) + six*l2 - eight*lr)/l4r;
+            bornSum               += term/sixteen;
+        
+        }
+        bornSum        = one/(_atomicRadii[ii]*_atomicRadii[ii]*_atomicRadii[ii]) - bornSum;
+        _bornRadii[ii] = (bornSum <= zero) ? bigRadius : POW( bornSum, -oneThird );
+    }
+
+    return;
+}

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceGeneralizedKirkwoodForce.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 __AmoebaReferenceGeneralizedKirkwoodForce_H__
+#define __AmoebaReferenceGeneralizedKirkwoodForce_H__
+
+#include "SimTKUtilities/RealVec.h"
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+// ---------------------------------------------------------------------------------------
+
+class AmoebaReferenceGeneralizedKirkwoodForce {
+
+public:
+
+    /**
+     *  Constructor
+     *  
+     */
+    AmoebaReferenceGeneralizedKirkwoodForce( );
+ 
+    /**
+     *  Destructor
+     *  
+     */
+    ~AmoebaReferenceGeneralizedKirkwoodForce( ){};
+ 
+    /**
+     *  Get number of particles 
+     *
+     *  @return numParticles
+     *
+     */
+    int getNumParticles( void ) const;
+
+    /**
+     *  Set numParticles
+     *
+     *  @param numParticles
+     *
+     */
+    void setNumParticles( int numParticles );
+
+    /**
+     * Get includeCavityTerm flag 
+     *
+     * @return includeCavityTerm
+     *
+     */
+    int getIncludeCavityTerm( void ) const;
+
+    /**
+     *  Set includeCavityTerm flag
+     *
+     *  @param includeCavityTerm flag indicating whether surface area term is to be included
+     *
+     */
+    void setIncludeCavityTerm( int includeCavityTerm );
+
+    /**
+     *  Get directPolarization flag 
+     *
+     *  @return directPolarization
+     *
+     */
+    int getDirectPolarization( void ) const;
+
+    /**
+     *  Set directPolarization flag
+     * 
+     *  @param directPolarization nonzero if direct as opposed to mutual polarization
+     */
+    void setDirectPolarization( int directPolarization );
+
+    /**
+     *  Get solute dielectric
+     *
+     *  @return soluteDielectric
+     */
+    RealOpenMM getSoluteDielectric( void ) const;
+
+    /**
+     *  Set solute dielectric
+     * 
+     *  @param soluteDielectric solute dielectric
+     * 
+     */
+    void setSoluteDielectric( RealOpenMM soluteDielectric );
+
+    /**
+     *  Get solvent dielectric
+     *
+     *  @return solventDielectric
+     *
+     */
+    RealOpenMM getSolventDielectric( void ) const;
+
+    /**
+     *  Set solvent dielectric 
+     *
+     *  @param solventDielectric solvent dielectric
+     *
+     */
+    void setSolventDielectric( RealOpenMM solventDielectric );
+
+    /**
+     *  Get dielectric offset
+     *
+     *  @return dielectricOffset
+     *
+     */
+    RealOpenMM getDielectricOffset( void ) const;
+
+    /**
+     * Set dielectric offset
+     *
+     * @param dielectricOffset dielectric offset
+     *
+     */
+    void setDielectricOffset( RealOpenMM dielectricOffset );
+
+    /**
+     * Get probeRadius
+     *
+     * @return probeRadius
+     *
+     */
+    RealOpenMM getProbeRadius( void ) const;
+
+    /**
+     * Set probe radius
+     *
+     * @param probeRadius probe radiue
+     *
+     */
+    void setProbeRadius( RealOpenMM probeRadius );
+
+    /**
+     * Get surfaceAreaFactor
+     *
+     * @return surfaceAreaFactor
+     *
+     */
+    RealOpenMM getSurfaceAreaFactor( void ) const;
+
+    /**
+     * Set surface area factor
+     *
+     * @param surfaceAreaFactor surface area factor
+     *
+     */
+    void setSurfaceAreaFactor( RealOpenMM surfaceAreaFactor );
+
+    /**
+     * Set atomic radii
+     *
+     * @param atomicRadii input vector of atomic radii
+     *
+     */
+    void setAtomicRadii( const vector<RealOpenMM>& atomicRadii );
+
+    /**
+     * Get atomic radii
+     *
+     * @param atomicRadii output vector of atomic radii
+     *
+     */
+    void getAtomicRadii( vector<RealOpenMM>& atomicRadii ) const;
+
+    /**
+     * Set scale factors
+     *
+     * @param scaleFactors input vector of scale factors
+     *
+     */
+    void setScaleFactors( const vector<RealOpenMM>& scaleFactors );
+
+    /**
+     * Get scale factors
+     *
+     * @param scaleFactors output vector of scale factors
+     *
+     */
+    void getScaleFactors( vector<RealOpenMM>& scaleFactors ) const;
+
+    /**
+     * Set charges
+     *
+     * @param charges input vector of charges
+     *
+     */
+    void setCharges( const vector<RealOpenMM>& charges );
+
+    /**
+     * Calculate Grycuk Born radii
+     *
+     * @param particlePositions particle positions
+     *
+     */
+    void calculateGrycukBornRadii( const vector<RealVec>& particlePositions );
+         
+    /**
+     * Get Grycik Born radii (must have called calculateGrycukBornRadii())
+     *
+     * @param bornRadii vector of Born radii
+     *
+     */
+    void getGrycukBornRadii( vector<RealOpenMM>& bornRadii ) const;     
+
+private:
+
+    int _numParticles;
+    int _includeCavityTerm;
+    int _directPolarization;
+
+    RealOpenMM _soluteDielectric;
+    RealOpenMM _solventDielectric;
+    RealOpenMM _dielectricOffset;
+    RealOpenMM _probeRadius;
+    RealOpenMM _surfaceAreaFactor;
+
+    std::vector<RealOpenMM> _atomicRadii;
+    std::vector<RealOpenMM> _scaleFactors;
+    std::vector<RealOpenMM> _charges;
+
+    std::vector<RealOpenMM> _bornRadii;
+
+};
+
+#endif // _AmoebaReferenceGeneralizedKirkwoodForce___