Implemented updateParametersInContext() for GBSAOBCForce

olla/include/openmm/kernels.h

@@ -575,6 +575,13 @@ public:
      * @return the potential energy due to the force
      */
     virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GBSAOBCForce to copy the parameters from
+     */
+    virtual void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) = 0;
 };
 
 /**

openmmapi/include/openmm/GBSAOBCForce.h

@@ -153,6 +153,18 @@ public:
      * @param distance    the cutoff distance, measured in nm
      */
     void setCutoffDistance(double distance);
+    /**
+     * Update the particle parameters in a Context to match those stored in this Force object.  This method
+     * provides an efficient method to update certain parameters in an existing Context without needing to
+     * reinitialize it.  Simply call setParticleParameters() to modify this object's parameters, then call
+     * updateParametersInState() to copy them over to the Context.
+     * 
+     * The only information this method updates is the values of per-particle parameters.  All other aspects
+     * of the Force (the nonbonded method, the cutoff distance, etc.) are unaffected and can only be changed
+     * by reinitializing the Context.  Furthermore, this method cannot be used to add new particles, only to
+     * change the parameters of existing ones.
+     */
+    void updateParametersInContext(Context& context);
 protected:
     ForceImpl* createImpl();
 private:

openmmapi/include/openmm/internal/GBSAOBCForceImpl.h

@@ -58,6 +58,7 @@ public:
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }
     std::vector<std::string> getKernelNames();
+    void updateParametersInContext(ContextImpl& context);
 private:
     GBSAOBCForce& owner;
     Kernel kernel;

openmmapi/src/GBSAOBCForce.cpp

@@ -78,3 +78,7 @@ void GBSAOBCForce::setCutoffDistance(double distance) {
 ForceImpl* GBSAOBCForce::createImpl() {
     return new GBSAOBCForceImpl(*this);
 }
+
+void GBSAOBCForce::updateParametersInContext(Context& context) {
+    dynamic_cast<GBSAOBCForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
+}

openmmapi/src/GBSAOBCForceImpl.cpp

@@ -66,3 +66,7 @@ std::vector<std::string> GBSAOBCForceImpl::getKernelNames() {
     names.push_back(CalcGBSAOBCForceKernel::Name());
     return names;
 }
+
+void GBSAOBCForceImpl::updateParametersInContext(ContextImpl& context) {
+    kernel.getAs<CalcGBSAOBCForceKernel>().copyParametersToContext(context, owner);
+}

platforms/cuda/src/CudaKernels.cpp

@@ -1044,6 +1044,10 @@ double CudaCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeFor
 	return 0.0;
 }
 
+void CudaCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) {
+    throw OpenMMException("CudaPlatform does not support copyParametersToContext");
+}
+
 class CudaCalcGBVIForceKernel::ForceInfo : public CudaForceInfo {
 public:
     ForceInfo(const GBVIForce& force) : force(force) {

platforms/cuda/src/CudaKernels.h

@@ -581,6 +581,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GBSAOBCForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
     class ForceInfo;
     CudaPlatform::PlatformData& data;

platforms/opencl/src/OpenCLKernels.cpp

@@ -1797,6 +1797,34 @@ double OpenCLCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeF
     return 0.0;
 }
 
+void OpenCLCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) {
+    // Make sure the new parameters are acceptable.
+    
+    int numParticles = force.getNumParticles();
+    if (numParticles != cl.getNumAtoms())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    
+    // Record the per-particle parameters.
+    
+    OpenCLArray<mm_float4>& posq = cl.getPosq();
+    posq.download();
+    vector<mm_float2> paramsVector(numParticles);
+    const double dielectricOffset = 0.009;
+    for (int i = 0; i < numParticles; i++) {
+        double charge, radius, scalingFactor;
+        force.getParticleParameters(i, charge, radius, scalingFactor);
+        radius -= dielectricOffset;
+        paramsVector[i] = mm_float2((float) radius, (float) (scalingFactor*radius));
+        posq[i].w = (float) charge;
+    }
+    posq.upload();
+    params->upload(paramsVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 class OpenCLCustomGBForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomGBForceInfo(int requiredBuffers, const CustomGBForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {

platforms/opencl/src/OpenCLKernels.h

@@ -636,6 +636,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GBSAOBCForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
     double prefactor;
     bool hasCreatedKernels;

platforms/opencl/tests/TestOpenCLGBSAOBCForce.cpp

@@ -73,6 +73,17 @@ void testSingleParticle() {
     double extendedRadius = bornRadius+0.14; // probe radius
     double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+    
+    // Change the parameters and see if it is still correct.
+    
+    gbsa->setParticleParameters(0, 0.4, 0.25, 1);
+    gbsa->updateParametersInContext(context);
+    state = context.getState(State::Energy);
+    bornRadius = 0.25-0.009; // dielectric offset
+    bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
+    extendedRadius = bornRadius+0.14;
+    nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.25/bornRadius, 6.0);
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
 void testCutoffAndPeriodic() {

platforms/reference/src/ReferenceKernels.cpp

@@ -951,6 +951,27 @@ double ReferenceCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool inclu
     return obc->computeBornEnergyForces(posData, charges, forceData);
 }
 
+void ReferenceCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) {
+    int numParticles = force.getNumParticles();
+    ObcParameters* obcParameters = obc->getObcParameters();
+    if (numParticles != obcParameters->getAtomicRadii().size())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    vector<RealOpenMM> atomicRadii(numParticles);
+    vector<RealOpenMM> scaleFactors(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        double charge, radius, scalingFactor;
+        force.getParticleParameters(i, charge, radius, scalingFactor);
+        charges[i] = (RealOpenMM) charge;
+        atomicRadii[i] = (RealOpenMM) radius;
+        scaleFactors[i] = (RealOpenMM) scalingFactor;
+    }
+    obcParameters->setAtomicRadii(atomicRadii);
+    obcParameters->setScaledRadiusFactors(scaleFactors);
+}
+
 ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
     if (gbvi) {
         GBVIParameters * gBVIParameters = gbvi->getGBVIParameters();

platforms/reference/src/ReferenceKernels.h

@@ -585,6 +585,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GBSAOBCForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
     CpuObc* obc;
     std::vector<RealOpenMM> charges;

platforms/reference/tests/TestReferenceGBSAOBCForce.cpp

@@ -69,6 +69,17 @@ void testSingleParticle() {
     double extendedRadius = bornRadius+0.14; // probe radius
     double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+    
+    // Change the parameters and see if it is still correct.
+    
+    forceField->setParticleParameters(0, 0.4, 0.25, 1);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Energy);
+    bornRadius = 0.25-0.009; // dielectric offset
+    bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
+    extendedRadius = bornRadius+0.14;
+    nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.25/bornRadius, 6.0);
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
 void testCutoffAndPeriodic() {