updateParametersInContext() can modify parameter offsets (#4732)

* updateParametersInContext() can modify parameter offsets

* Reordering respects parameter offsets

* Implemented for CUDA and HIP

openmmapi/include/openmm/NonbondedForce.h

@@ -572,11 +572,13 @@ public:
      * Simply call setParticleParameters() and setExceptionParameters() to modify this object's parameters, then call
      * updateParametersInContext() to copy them over to the Context.
      *
-     * This method has several limitations.  The only information it updates is the parameters of particles and exceptions.
-     * 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, only the chargeProd, sigma, and epsilon values of an exception
-     * can be changed; the pair of particles involved in the exception cannot change.  Finally, this method cannot be used
-     * to add new particles or exceptions, only to change the parameters of existing ones.
+     * This method has several limitations.  The only information it updates is the parameters of particles and exceptions,
+     * as well as parameter offsets for particles and exceptions.  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,
+     * only the chargeProd, sigma, and epsilon values of an exception can be changed; the pair of particles involved in the
+     * exception cannot change.  Likewise, it can update charge, sigma and epsilon for a parameter offset, but not the
+     * identify of the particle or exception the offset is applied to, or which global parameter it is based on.  Finally,
+     * this method cannot be used to add new particles or exceptions, only to change the parameters of existing ones.
      */
     void updateParametersInContext(Context& context);
     /**

platforms/cpu/include/CpuKernels.h

@@ -278,6 +278,7 @@ private:
     int numParticles, num14, chargePosqIndex, ljPosqIndex;
     std::vector<std::vector<int> > bonded14IndexArray;
     std::vector<std::vector<double> > bonded14ParamArray;
+    std::map<int, int> nb14Index;
     double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, ewaldSelfEnergy, dispersionCoefficient;
     int kmax[3], gridSize[3], dispersionGridSize[3];
     bool useSwitchingFunction, exceptionsArePeriodic, useOptimizedPme, hasInitializedPme, hasInitializedDispersionPme, hasParticleOffsets, hasExceptionOffsets;

platforms/cpu/src/CpuKernels.cpp

@@ -506,7 +506,6 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
     numParticles = force.getNumParticles();
     exclusions.resize(numParticles);
     vector<int> nb14s;
-    map<int, int> nb14Index;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
@@ -753,6 +752,32 @@ void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
         bonded14IndexArray[i][0] = particle1;
         bonded14IndexArray[i][1] = particle2;
     }
+    particleParamOffsets.clear();
+    exceptionParamOffsets.clear();
+    particleParamOffsets.resize(force.getNumParticles());
+    exceptionParamOffsets.resize(force.getNumExceptions());
+    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+        string param;
+        int particle;
+        double charge, sigma, epsilon;
+        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of a particle parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        particleParamOffsets[particle].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
+    }
+    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+        string param;
+        int exception;
+        double charge, sigma, epsilon;
+        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of an exception parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        exceptionParamOffsets[nb14Index[exception]].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
+    }
     computeParameters(context, false);
     
     // Recompute the coefficient for the dispersion correction.

platforms/cuda/include/CudaKernels.h

@@ -211,6 +211,7 @@ private:
     std::vector<std::pair<int, int> > exceptionAtoms;
     std::vector<std::string> paramNames;
     std::vector<double> paramValues;
+    std::map<int, int> exceptionIndex;
     double ewaldSelfEnergy, dispersionCoefficient, alpha, dispersionAlpha;
     int interpolateForceThreads;
     int gridSizeX, gridSizeY, gridSizeZ;

platforms/cuda/src/CudaKernels.cpp

@@ -89,8 +89,35 @@ double CudaCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bo
 class CudaCalcNonbondedForceKernel::ForceInfo : public CudaForceInfo {
 public:
     ForceInfo(const NonbondedForce& force) : force(force) {
+        particleOffset.resize(force.getNumParticles(), -1);
+        exceptionOffset.resize(force.getNumExceptions(), -1);
+        for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+            string parameter;
+            int particleIndex;
+            double chargeScale, sigmaScale, epsilonScale;
+            force.getParticleParameterOffset(i, parameter, particleIndex, chargeScale, sigmaScale, epsilonScale);
+            particleOffset[particleIndex] = i;
+        }
+        for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+            string parameter;
+            int exceptionIndex;
+            double chargeProdScale, sigmaScale, epsilonScale;
+            force.getExceptionParameterOffset(i, parameter, exceptionIndex, chargeProdScale, sigmaScale, epsilonScale);
+            exceptionOffset[exceptionIndex] = i;
+        }
     }
     bool areParticlesIdentical(int particle1, int particle2) {
+        if (particleOffset[particle1] != -1 || particleOffset[particle2] != -1) {
+            if (particleOffset[particle1] == -1 || particleOffset[particle2] == -1)
+                return false;
+            string parameter1, parameter2;
+            int particleIndex1, particleIndex2;
+            double chargeScale1, sigmaScale1, epsilonScale1, chargeScale2, sigmaScale2, epsilonScale2;
+            force.getParticleParameterOffset(particleOffset[particle1], parameter1, particleIndex1, chargeScale1, sigmaScale1, epsilonScale1);
+            force.getParticleParameterOffset(particleOffset[particle2], parameter2, particleIndex2, chargeScale2, sigmaScale2, epsilonScale2);
+            if (parameter1 != parameter2 || chargeScale1 != chargeScale2 || sigmaScale1 != sigmaScale2 || epsilonScale1 != epsilonScale2)
+                return false;
+        }
         double charge1, charge2, sigma1, sigma2, epsilon1, epsilon2;
         force.getParticleParameters(particle1, charge1, sigma1, epsilon1);
         force.getParticleParameters(particle2, charge2, sigma2, epsilon2);
@@ -108,6 +135,17 @@ public:
         particles[1] = particle2;
     }
     bool areGroupsIdentical(int group1, int group2) {
+        if (exceptionOffset[group1] != -1 || exceptionOffset[group2] != -1) {
+            if (exceptionOffset[group1] == -1 || exceptionOffset[group2] == -1)
+                return false;
+            string parameter1, parameter2;
+            int exceptionIndex1, exceptionIndex2;
+            double chargeProdScale1, sigmaScale1, epsilonScale1, chargeProdScale2, sigmaScale2, epsilonScale2;
+            force.getExceptionParameterOffset(exceptionOffset[group1], parameter1, exceptionIndex1, chargeProdScale1, sigmaScale1, epsilonScale1);
+            force.getExceptionParameterOffset(exceptionOffset[group2], parameter2, exceptionIndex2, chargeProdScale2, sigmaScale2, epsilonScale2);
+            if (parameter1 != parameter2 || chargeProdScale1 != chargeProdScale2 || sigmaScale1 != sigmaScale2 || epsilonScale1 != epsilonScale2)
+                return false;
+        }
         int particle1, particle2;
         double chargeProd1, chargeProd2, sigma1, sigma2, epsilon1, epsilon2;
         force.getExceptionParameters(group1, particle1, particle2, chargeProd1, sigma1, epsilon1);
@@ -116,6 +154,7 @@ public:
     }
 private:
     const NonbondedForce& force;
+    vector<int> particleOffset, exceptionOffset;
 };
 
 class CudaCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
@@ -253,7 +292,6 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
     }
     vector<pair<int, int> > exclusions;
     vector<int> exceptions;
-    map<int, int> exceptionIndex;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
@@ -1126,12 +1164,58 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
         }
         baseExceptionParams.upload(baseExceptionParamsVec);
     }
-    
+
+    // Record parameter offsets.
+
+    vector<vector<mm_float4> > particleOffsetVec(force.getNumParticles());
+    vector<vector<mm_float4> > exceptionOffsetVec(numExceptions);
+    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+        string param;
+        int particle;
+        double charge, sigma, epsilon;
+        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of a particle parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        particleOffsetVec[particle].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+    }
+    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+        string param;
+        int exception;
+        double charge, sigma, epsilon;
+        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
+        int index = exceptionIndex[exception];
+        if (index < startIndex || index >= endIndex)
+            continue;
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of an exception parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        exceptionOffsetVec[index-startIndex].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+    }
+    if (max(force.getNumParticleParameterOffsets(), 1) != particleParamOffsets.getSize())
+        throw OpenMMException("updateParametersInContext: The number of particle parameter offsets has changed");
+    vector<mm_float4> p, e;
+    for (int i = 0; i < particleOffsetVec.size(); i++)
+        for (int j = 0; j < particleOffsetVec[i].size(); j++)
+            p.push_back(particleOffsetVec[i][j]);
+    for (int i = 0; i < exceptionOffsetVec.size(); i++)
+        for (int j = 0; j < exceptionOffsetVec[i].size(); j++)
+            e.push_back(exceptionOffsetVec[i][j]);
+    if (force.getNumParticleParameterOffsets() > 0)
+        particleParamOffsets.upload(p);
+    if (max((int) e.size(), 1) != exceptionParamOffsets.getSize())
+        throw OpenMMException("updateParametersInContext: The number of exception parameter offsets has changed");
+    if (e.size() > 0)
+        exceptionParamOffsets.upload(e);
+
     // Compute other values.
-    
+
     if (force.getUseDispersionCorrection() && cu.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
-    cu.invalidateMolecules(info, firstParticle <= lastParticle, firstException <= lastException);
+    cu.invalidateMolecules(info, firstParticle <= lastParticle || force.getNumParticleParameterOffsets() > 0,
+                           firstException <= lastException || force.getNumExceptionParameterOffsets() > 0);
     recomputeParams = true;
 }
 

platforms/hip/include/HipKernels.h

@@ -207,6 +207,7 @@ private:
     std::vector<std::pair<int, int> > exceptionAtoms;
     std::vector<std::string> paramNames;
     std::vector<double> paramValues;
+    std::map<int, int> exceptionIndex;
     double ewaldSelfEnergy, dispersionCoefficient, alpha, dispersionAlpha;
     int interpolateForceThreads;
     int gridSizeX, gridSizeY, gridSizeZ;

platforms/hip/src/HipKernels.cpp

@@ -91,8 +91,35 @@ double HipCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, boo
 class HipCalcNonbondedForceKernel::ForceInfo : public HipForceInfo {
 public:
     ForceInfo(const NonbondedForce& force) : force(force) {
+        particleOffset.resize(force.getNumParticles(), -1);
+        exceptionOffset.resize(force.getNumExceptions(), -1);
+        for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+            string parameter;
+            int particleIndex;
+            double chargeScale, sigmaScale, epsilonScale;
+            force.getParticleParameterOffset(i, parameter, particleIndex, chargeScale, sigmaScale, epsilonScale);
+            particleOffset[particleIndex] = i;
+        }
+        for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+            string parameter;
+            int exceptionIndex;
+            double chargeProdScale, sigmaScale, epsilonScale;
+            force.getExceptionParameterOffset(i, parameter, exceptionIndex, chargeProdScale, sigmaScale, epsilonScale);
+            exceptionOffset[exceptionIndex] = i;
+        }
     }
     bool areParticlesIdentical(int particle1, int particle2) {
+        if (particleOffset[particle1] != -1 || particleOffset[particle2] != -1) {
+            if (particleOffset[particle1] == -1 || particleOffset[particle2] == -1)
+                return false;
+            string parameter1, parameter2;
+            int particleIndex1, particleIndex2;
+            double chargeScale1, sigmaScale1, epsilonScale1, chargeScale2, sigmaScale2, epsilonScale2;
+            force.getParticleParameterOffset(particleOffset[particle1], parameter1, particleIndex1, chargeScale1, sigmaScale1, epsilonScale1);
+            force.getParticleParameterOffset(particleOffset[particle2], parameter2, particleIndex2, chargeScale2, sigmaScale2, epsilonScale2);
+            if (parameter1 != parameter2 || chargeScale1 != chargeScale2 || sigmaScale1 != sigmaScale2 || epsilonScale1 != epsilonScale2)
+                return false;
+        }
         double charge1, charge2, sigma1, sigma2, epsilon1, epsilon2;
         force.getParticleParameters(particle1, charge1, sigma1, epsilon1);
         force.getParticleParameters(particle2, charge2, sigma2, epsilon2);
@@ -110,6 +137,17 @@ public:
         particles[1] = particle2;
     }
     bool areGroupsIdentical(int group1, int group2) {
+        if (exceptionOffset[group1] != -1 || exceptionOffset[group2] != -1) {
+            if (exceptionOffset[group1] == -1 || exceptionOffset[group2] == -1)
+                return false;
+            string parameter1, parameter2;
+            int exceptionIndex1, exceptionIndex2;
+            double chargeProdScale1, sigmaScale1, epsilonScale1, chargeProdScale2, sigmaScale2, epsilonScale2;
+            force.getExceptionParameterOffset(exceptionOffset[group1], parameter1, exceptionIndex1, chargeProdScale1, sigmaScale1, epsilonScale1);
+            force.getExceptionParameterOffset(exceptionOffset[group2], parameter2, exceptionIndex2, chargeProdScale2, sigmaScale2, epsilonScale2);
+            if (parameter1 != parameter2 || chargeProdScale1 != chargeProdScale2 || sigmaScale1 != sigmaScale2 || epsilonScale1 != epsilonScale2)
+                return false;
+        }
         int particle1, particle2;
         double chargeProd1, chargeProd2, sigma1, sigma2, epsilon1, epsilon2;
         force.getExceptionParameters(group1, particle1, particle2, chargeProd1, sigma1, epsilon1);
@@ -118,6 +156,7 @@ public:
     }
 private:
     const NonbondedForce& force;
+    vector<int> particleOffset, exceptionOffset;
 };
 
 class HipCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
@@ -247,7 +286,6 @@ void HipCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
     }
     vector<pair<int, int> > exclusions;
     vector<int> exceptions;
-    map<int, int> exceptionIndex;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
@@ -1035,11 +1073,57 @@ void HipCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
         baseExceptionParams.upload(baseExceptionParamsVec);
     }
 
+    // Record parameter offsets.
+
+    vector<vector<mm_float4> > particleOffsetVec(force.getNumParticles());
+    vector<vector<mm_float4> > exceptionOffsetVec(numExceptions);
+    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+        string param;
+        int particle;
+        double charge, sigma, epsilon;
+        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of a particle parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        particleOffsetVec[particle].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+    }
+    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+        string param;
+        int exception;
+        double charge, sigma, epsilon;
+        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
+        int index = exceptionIndex[exception];
+        if (index < startIndex || index >= endIndex)
+            continue;
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of an exception parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        exceptionOffsetVec[index-startIndex].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+    }
+    if (max(force.getNumParticleParameterOffsets(), 1) != particleParamOffsets.getSize())
+        throw OpenMMException("updateParametersInContext: The number of particle parameter offsets has changed");
+    vector<mm_float4> p, e;
+    for (int i = 0; i < particleOffsetVec.size(); i++)
+        for (int j = 0; j < particleOffsetVec[i].size(); j++)
+            p.push_back(particleOffsetVec[i][j]);
+    for (int i = 0; i < exceptionOffsetVec.size(); i++)
+        for (int j = 0; j < exceptionOffsetVec[i].size(); j++)
+            e.push_back(exceptionOffsetVec[i][j]);
+    if (force.getNumParticleParameterOffsets() > 0)
+        particleParamOffsets.upload(p);
+    if (max((int) e.size(), 1) != exceptionParamOffsets.getSize())
+        throw OpenMMException("updateParametersInContext: The number of exception parameter offsets has changed");
+    if (e.size() > 0)
+        exceptionParamOffsets.upload(e);
+
     // Compute other values.
 
     if (force.getUseDispersionCorrection() && cu.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
-    cu.invalidateMolecules(info, firstParticle <= lastParticle, firstException <= lastException);
+    cu.invalidateMolecules(info, firstParticle <= lastParticle || force.getNumParticleParameterOffsets() > 0,
+                           firstException <= lastException || force.getNumExceptionParameterOffsets() > 0);
     recomputeParams = true;
 }
 

platforms/opencl/include/OpenCLKernels.h

@@ -214,6 +214,7 @@ private:
     std::vector<std::pair<int, int> > exceptionAtoms;
     std::vector<std::string> paramNames;
     std::vector<double> paramValues;
+    std::map<int, int> exceptionIndex;
     double ewaldSelfEnergy, dispersionCoefficient, alpha, dispersionAlpha;
     int gridSizeX, gridSizeY, gridSizeZ;
     int dispersionGridSizeX, dispersionGridSizeY, dispersionGridSizeZ;

platforms/opencl/src/OpenCLKernels.cpp

@@ -104,8 +104,35 @@ double OpenCLCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context,
 class OpenCLCalcNonbondedForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
     ForceInfo(int requiredBuffers, const NonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+        particleOffset.resize(force.getNumParticles(), -1);
+        exceptionOffset.resize(force.getNumExceptions(), -1);
+        for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+            string parameter;
+            int particleIndex;
+            double chargeScale, sigmaScale, epsilonScale;
+            force.getParticleParameterOffset(i, parameter, particleIndex, chargeScale, sigmaScale, epsilonScale);
+            particleOffset[particleIndex] = i;
+        }
+        for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+            string parameter;
+            int exceptionIndex;
+            double chargeProdScale, sigmaScale, epsilonScale;
+            force.getExceptionParameterOffset(i, parameter, exceptionIndex, chargeProdScale, sigmaScale, epsilonScale);
+            exceptionOffset[exceptionIndex] = i;
+        }
     }
     bool areParticlesIdentical(int particle1, int particle2) {
+        if (particleOffset[particle1] != -1 || particleOffset[particle2] != -1) {
+            if (particleOffset[particle1] == -1 || particleOffset[particle2] == -1)
+                return false;
+            string parameter1, parameter2;
+            int particleIndex1, particleIndex2;
+            double chargeScale1, sigmaScale1, epsilonScale1, chargeScale2, sigmaScale2, epsilonScale2;
+            force.getParticleParameterOffset(particleOffset[particle1], parameter1, particleIndex1, chargeScale1, sigmaScale1, epsilonScale1);
+            force.getParticleParameterOffset(particleOffset[particle2], parameter2, particleIndex2, chargeScale2, sigmaScale2, epsilonScale2);
+            if (parameter1 != parameter2 || chargeScale1 != chargeScale2 || sigmaScale1 != sigmaScale2 || epsilonScale1 != epsilonScale2)
+                return false;
+        }
         double charge1, charge2, sigma1, sigma2, epsilon1, epsilon2;
         force.getParticleParameters(particle1, charge1, sigma1, epsilon1);
         force.getParticleParameters(particle2, charge2, sigma2, epsilon2);
@@ -123,6 +150,17 @@ public:
         particles[1] = particle2;
     }
     bool areGroupsIdentical(int group1, int group2) {
+        if (exceptionOffset[group1] != -1 || exceptionOffset[group2] != -1) {
+            if (exceptionOffset[group1] == -1 || exceptionOffset[group2] == -1)
+                return false;
+            string parameter1, parameter2;
+            int exceptionIndex1, exceptionIndex2;
+            double chargeProdScale1, sigmaScale1, epsilonScale1, chargeProdScale2, sigmaScale2, epsilonScale2;
+            force.getExceptionParameterOffset(exceptionOffset[group1], parameter1, exceptionIndex1, chargeProdScale1, sigmaScale1, epsilonScale1);
+            force.getExceptionParameterOffset(exceptionOffset[group2], parameter2, exceptionIndex2, chargeProdScale2, sigmaScale2, epsilonScale2);
+            if (parameter1 != parameter2 || chargeProdScale1 != chargeProdScale2 || sigmaScale1 != sigmaScale2 || epsilonScale1 != epsilonScale2)
+                return false;
+        }
         int particle1, particle2;
         double chargeProd1, chargeProd2, sigma1, sigma2, epsilon1, epsilon2;
         force.getExceptionParameters(group1, particle1, particle2, chargeProd1, sigma1, epsilon1);
@@ -131,6 +169,7 @@ public:
     }
 private:
     const NonbondedForce& force;
+    vector<int> particleOffset, exceptionOffset;
 };
 
 class OpenCLCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
@@ -257,7 +296,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     }
     vector<pair<int, int> > exclusions;
     vector<int> exceptions;
-    map<int, int> exceptionIndex;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
@@ -1146,12 +1184,58 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
         }
         baseExceptionParams.upload(baseExceptionParamsVec);
     }
-    
+
+    // Record parameter offsets.
+
+    vector<vector<mm_float4> > particleOffsetVec(force.getNumParticles());
+    vector<vector<mm_float4> > exceptionOffsetVec(numExceptions);
+    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+        string param;
+        int particle;
+        double charge, sigma, epsilon;
+        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of a particle parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        particleOffsetVec[particle].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+    }
+    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+        string param;
+        int exception;
+        double charge, sigma, epsilon;
+        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
+        int index = exceptionIndex[exception];
+        if (index < startIndex || index >= endIndex)
+            continue;
+        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
+        if (paramPos == paramNames.end())
+            throw OpenMMException("updateParametersInContext: The parameter of an exception parameter offset has changed");
+        int paramIndex = paramPos-paramNames.begin();
+        exceptionOffsetVec[index-startIndex].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
+    }
+    if (max(force.getNumParticleParameterOffsets(), 1) != particleParamOffsets.getSize())
+        throw OpenMMException("updateParametersInContext: The number of particle parameter offsets has changed");
+    vector<mm_float4> p, e;
+    for (int i = 0; i < particleOffsetVec.size(); i++)
+        for (int j = 0; j < particleOffsetVec[i].size(); j++)
+            p.push_back(particleOffsetVec[i][j]);
+    for (int i = 0; i < exceptionOffsetVec.size(); i++)
+        for (int j = 0; j < exceptionOffsetVec[i].size(); j++)
+            e.push_back(exceptionOffsetVec[i][j]);
+    if (force.getNumParticleParameterOffsets() > 0)
+        particleParamOffsets.upload(p);
+    if (max((int) e.size(), 1) != exceptionParamOffsets.getSize())
+        throw OpenMMException("updateParametersInContext: The number of exception parameter offsets has changed");
+    if (e.size() > 0)
+        exceptionParamOffsets.upload(e);
+
     // Compute other values.
-    
+
     if (force.getUseDispersionCorrection() && cl.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
-    cl.invalidateMolecules(info, firstParticle <= lastParticle, firstException <= lastException);
+    cl.invalidateMolecules(info, firstParticle <= lastParticle || force.getNumParticleParameterOffsets() > 0,
+                           firstException <= lastException || force.getNumExceptionParameterOffsets() > 0);
     recomputeParams = true;
 }
 

platforms/reference/include/ReferenceKernels.h

@@ -662,6 +662,7 @@ private:
     std::vector<std::vector<double> > particleParamArray, bonded14ParamArray;
     std::vector<std::array<double, 3> > baseParticleParams, baseExceptionParams;
     std::map<std::pair<std::string, int>, std::array<double, 3> > particleParamOffsets, exceptionParamOffsets;
+    std::map<int, int> nb14Index;
     double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, dispersionCoefficient;
     int kmax[3], gridSize[3], dispersionGridSize[3];
     bool useSwitchingFunction, exceptionsArePeriodic;

platforms/reference/src/ReferenceKernels.cpp

@@ -903,7 +903,6 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
     numParticles = force.getNumParticles();
     exclusions.resize(numParticles);
     vector<int> nb14s;
-    map<int, int> nb14Index;
     for (int i = 0; i < force.getNumExceptions(); i++) {
         int particle1, particle2;
         double chargeProd, sigma, epsilon;
@@ -1038,6 +1037,18 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
 void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force, int firstParticle, int lastParticle, int firstException, int lastException) {
     if (force.getNumParticles() != numParticles)
         throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    if (force.getNumParticleParameterOffsets() != particleParamOffsets.size())
+        throw OpenMMException("updateParametersInContext: The number of particle parameter offsets has changed");
+    if (force.getNumExceptionParameterOffsets() != exceptionParamOffsets.size())
+        throw OpenMMException("updateParametersInContext: The number of exception parameter offsets has changed");
+    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+        string param;
+        int particle;
+        double charge, sigma, epsilon;
+        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
+        if (particleParamOffsets.find(make_pair(param, particle)) == particleParamOffsets.end())
+            throw OpenMMException("updateParametersInContext: The parameter or particle index of a particle parameter offset has changed");
+    }
 
     // Identify which exceptions are 1-4 interactions.
 
@@ -1048,6 +1059,8 @@ void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& con
         double charge, sigma, epsilon;
         force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
         exceptionsWithOffsets.insert(exception);
+        if (exceptionParamOffsets.find(make_pair(param, nb14Index[exception])) == exceptionParamOffsets.end())
+            throw OpenMMException("updateParametersInContext: The parameter or exception index of an exception parameter offset has changed");
     }
     vector<int> nb14s;
     for (int i = 0; i < force.getNumExceptions(); i++) {
@@ -1070,6 +1083,22 @@ void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& con
         bonded14IndexArray[i][0] = particle1;
         bonded14IndexArray[i][1] = particle2;
     }
+    particleParamOffsets.clear();
+    exceptionParamOffsets.clear();
+    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
+        string param;
+        int particle;
+        double charge, sigma, epsilon;
+        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
+        particleParamOffsets[make_pair(param, particle)] = {charge, sigma, epsilon};
+    }
+    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
+        string param;
+        int exception;
+        double charge, sigma, epsilon;
+        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
+        exceptionParamOffsets[make_pair(param, nb14Index[exception])] = {charge, sigma, epsilon};
+    }
     
     // Recompute the coefficient for the dispersion correction.
 

tests/TestNonbondedForce.h

@@ -901,6 +901,39 @@ void testParameterOffsets() {
             energy += ONE_4PI_EPS0*pairChargeProd[i][j]/dist + 4.0*pairEpsilon[i][j]*(pow(x, 12.0)-pow(x, 6.0));
         }
     ASSERT_EQUAL_TOL(energy, context.getState(State::Energy).getPotentialEnergy(), 1e-5);
+
+    // Update the offsets and see if the energy is still correct.
+
+    force->setParticleParameterOffset(0, "p1", 0, 2.0, 0.6, 0.6);
+    force->setParticleParameterOffset(1, "p2", 1, 1.1, 1.1, 2.1);
+    force->setExceptionParameterOffset(0, "p1", 1, 0.4, 0.4, 1.4);
+    force->updateParametersInContext(context);
+    particleCharge = {0.0+2.0*0.5, 1.0+1.1*1.5, -1.0, 0.5};
+    particleSigma = {1.0+0.6*0.5, 0.5+1.1*1.5, 2.0, 2.0};
+    particleEpsilon = {0.5+0.6*0.5, 0.6+2.1*1.5, 0.7, 0.8};
+    for (int i = 0; i < 4; i++)
+        for (int j = i+1; j < 4; j++) {
+            pairChargeProd[i][j] = particleCharge[i]*particleCharge[j];
+            pairSigma[i][j] = 0.5*(particleSigma[i]+particleSigma[j]);
+            pairEpsilon[i][j] = sqrt(particleEpsilon[i]*particleEpsilon[j]);
+        }
+    pairChargeProd[0][3] = 0.0;
+    pairSigma[0][3] = 1.0;
+    pairEpsilon[0][3] = 0.0;
+    pairChargeProd[2][3] = 0.5+0.4*0.5;
+    pairSigma[2][3] = 1.0+0.4*0.5;
+    pairEpsilon[2][3] = 1.5+1.4*0.5;
+    pairChargeProd[0][1] = 1.0;
+    pairSigma[0][1] = 1.5;
+    pairEpsilon[0][1] = 1.0;
+    energy = 0.0;
+    for (int i = 0; i < 4; i++)
+        for (int j = i+1; j < 4; j++) {
+            double dist = j-i;
+            double x = pairSigma[i][j]/dist;
+            energy += ONE_4PI_EPS0*pairChargeProd[i][j]/dist + 4.0*pairEpsilon[i][j]*(pow(x, 12.0)-pow(x, 6.0));
+        }
+    ASSERT_EQUAL_TOL(energy, context.getState(State::Energy).getPotentialEnergy(), 1e-5);
 }
 
 void testEwaldExceptions() {