Implemented updateParametersInContext() for DrudeForce

plugins/drude/openmmapi/include/openmm/DrudeForce.h

@@ -125,7 +125,7 @@ public:
      */
     void setParticleParameters(int index, int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34);
     /**
-     * Add an interaction to the list of screenedPairs.
+     * Add an interaction to the list of screened pairs.
      *
      * @param particle1  the index within this Force of the first particle involved in the interaction
      * @param particle2  the index within this Force of the second particle involved in the interaction
@@ -152,10 +152,10 @@ public:
      */
     void setScreenedPairParameters(int index, int particle1, int particle2, double thole);
     /**
-     * Update the particle and screenedPair parameters in a Context to match those stored in this Force object.  This method
+     * Update the particle and screened pair 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() and setScreenedPairParameters() to modify this object's parameters, then call
-     * updateParametersInState() to copy them over to the Context.
+     * updateParametersInContext() to copy them over to the Context.
      * 
      * This method has several limitations.  It can be used to modify the numeric parameters associated with a particle or
      * screened pair (polarizability, thole, etc.), but not the identities of the particles they involve.  It also cannot

plugins/drude/platforms/cuda/src/CudaDrudeKernels.cpp

@@ -179,7 +179,59 @@ double CudaCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForce
 }
 
 void CudaCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) {
+    int numContexts = cu.getPlatformData().contexts.size();
+    
+    // Set the particle parameters.
     
+    int startParticleIndex = cu.getContextIndex()*force.getNumParticles()/numContexts;
+    int endParticleIndex = (cu.getContextIndex()+1)*force.getNumParticles()/numContexts;
+    int numParticles = endParticleIndex-startParticleIndex;
+    if (numParticles > 0) {
+        if (particleParams == NULL || numParticles != particleParams->getSize())
+            throw OpenMMException("updateParametersInContext: The number of Drude particles has changed");
+        vector<float4> paramVector(numParticles);
+        for (int i = 0; i < numParticles; i++) {
+            int p, p1, p2, p3, p4;
+            double charge, polarizability, aniso12, aniso34;
+            force.getParticleParameters(startParticleIndex+i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
+            double a1 = (p2 == -1 ? 1 : aniso12);
+            double a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34);
+            double a3 = 3-a1-a2;
+            double k3 = charge*charge/(polarizability*a3);
+            double k1 = charge*charge/(polarizability*a1) - k3;
+            double k2 = charge*charge/(polarizability*a2) - k3;
+            if (p2 == -1)
+                k1 = 0;
+            if (p3 == -1 || p4 == -1)
+                k2 = 0;
+            paramVector[i] = make_float4((float) k1, (float) k2, (float) k3, 0.0f);
+        }
+        particleParams->upload(paramVector);
+    }
+    
+    // Set the pair parameters.
+    
+    int startPairIndex = cu.getContextIndex()*force.getNumScreenedPairs()/numContexts;
+    int endPairIndex = (cu.getContextIndex()+1)*force.getNumScreenedPairs()/numContexts;
+    int numPairs = endPairIndex-startPairIndex;
+    if (numPairs > 0) {
+        if (pairParams == NULL || numPairs != pairParams->getSize())
+            throw OpenMMException("updateParametersInContext: The number of screened pairs has changed");
+        vector<float2> paramVector(numPairs);
+        for (int i = 0; i < numPairs; i++) {
+            int drude1, drude2;
+            double thole;
+            force.getScreenedPairParameters(startPairIndex+i, drude1, drude2, thole);
+            int p, p1, p2, p3, p4;
+            double charge1, charge2, polarizability1, polarizability2, aniso12, aniso34;
+            force.getParticleParameters(drude1, p, p1, p2, p3, p4, charge1, polarizability1, aniso12, aniso34);
+            force.getParticleParameters(drude2, p, p1, p2, p3, p4, charge2, polarizability2, aniso12, aniso34);
+            double screeningScale = thole/pow(polarizability1*polarizability2, 1.0/6.0);
+            double energyScale = ONE_4PI_EPS0*charge1*charge2;
+            paramVector[i] = make_float2((float) screeningScale, (float) energyScale);
+        }
+        pairParams->upload(paramVector);
+    }
 }
 
 CudaIntegrateDrudeLangevinStepKernel::~CudaIntegrateDrudeLangevinStepKernel() {

plugins/drude/platforms/cuda/tests/TestCudaDrudeForce.cpp

@@ -156,6 +156,43 @@ void testThole() {
     validateForce(system, positions, energySpring1+energySpring2+energyDipole);
 }
 
+void testChangingParameters() {
+    const double k = 1.5;
+    const double charge = 0.1;
+    const double alpha = charge*charge/k;
+    Platform& platform = Platform::getPlatformByName("CUDA");
+    
+    // Create the system.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    DrudeForce* drude = new DrudeForce();
+    drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
+    system.addForce(drude);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    
+    // Check the energy.
+    
+    VerletIntegrator integ(1.0);
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*k*3*3, state.getPotentialEnergy(), 1e-5);
+    
+    // Modify the parameters.
+    
+    const double k2 = 2.2;
+    const double charge2 = 0.3;
+    const double alpha2 = charge2*charge2/k2;
+    drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
+    drude->updateParametersInContext(context);
+    state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*k2*3*3, state.getPotentialEnergy(), 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerDrudeCudaKernelFactories();
@@ -164,6 +201,7 @@ int main(int argc, char* argv[]) {
         testSingleParticle();
         testAnisotropicParticle();
         testThole();
+        testChangingParameters();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/drude/platforms/opencl/src/OpenCLDrudeKernels.cpp

@@ -184,7 +184,59 @@ double OpenCLCalcDrudeForceKernel::execute(ContextImpl& context, bool includeFor
 }
 
 void OpenCLCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) {
+    int numContexts = cl.getPlatformData().contexts.size();
+    
+    // Set the particle parameters.
     
+    int startParticleIndex = cl.getContextIndex()*force.getNumParticles()/numContexts;
+    int endParticleIndex = (cl.getContextIndex()+1)*force.getNumParticles()/numContexts;
+    int numParticles = endParticleIndex-startParticleIndex;
+    if (numParticles > 0) {
+        if (particleParams == NULL || numParticles != particleParams->getSize())
+            throw OpenMMException("updateParametersInContext: The number of Drude particles has changed");
+        vector<mm_float4> paramVector(numParticles);
+        for (int i = 0; i < numParticles; i++) {
+            int p, p1, p2, p3, p4;
+            double charge, polarizability, aniso12, aniso34;
+            force.getParticleParameters(startParticleIndex+i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
+            double a1 = (p2 == -1 ? 1 : aniso12);
+            double a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34);
+            double a3 = 3-a1-a2;
+            double k3 = charge*charge/(polarizability*a3);
+            double k1 = charge*charge/(polarizability*a1) - k3;
+            double k2 = charge*charge/(polarizability*a2) - k3;
+            if (p2 == -1)
+                k1 = 0;
+            if (p3 == -1 || p4 == -1)
+                k2 = 0;
+            paramVector[i] = mm_float4((float) k1, (float) k2, (float) k3, 0.0f);
+        }
+        particleParams->upload(paramVector);
+    }
+    
+    // Set the pair parameters.
+    
+    int startPairIndex = cl.getContextIndex()*force.getNumScreenedPairs()/numContexts;
+    int endPairIndex = (cl.getContextIndex()+1)*force.getNumScreenedPairs()/numContexts;
+    int numPairs = endPairIndex-startPairIndex;
+    if (numPairs > 0) {
+        if (pairParams == NULL || numPairs != pairParams->getSize())
+            throw OpenMMException("updateParametersInContext: The number of screened pairs has changed");
+        vector<mm_float2> paramVector(numPairs);
+        for (int i = 0; i < numPairs; i++) {
+            int drude1, drude2;
+            double thole;
+            force.getScreenedPairParameters(startPairIndex+i, drude1, drude2, thole);
+            int p, p1, p2, p3, p4;
+            double charge1, charge2, polarizability1, polarizability2, aniso12, aniso34;
+            force.getParticleParameters(drude1, p, p1, p2, p3, p4, charge1, polarizability1, aniso12, aniso34);
+            force.getParticleParameters(drude2, p, p1, p2, p3, p4, charge2, polarizability2, aniso12, aniso34);
+            double screeningScale = thole/pow(polarizability1*polarizability2, 1.0/6.0);
+            double energyScale = ONE_4PI_EPS0*charge1*charge2;
+            paramVector[i] = mm_float2((float) screeningScale, (float) energyScale);
+        }
+        pairParams->upload(paramVector);
+    }
 }
 
 OpenCLIntegrateDrudeLangevinStepKernel::~OpenCLIntegrateDrudeLangevinStepKernel() {

plugins/drude/platforms/opencl/tests/TestOpenCLDrudeForce.cpp

@@ -156,6 +156,43 @@ void testThole() {
     validateForce(system, positions, energySpring1+energySpring2+energyDipole);
 }
 
+void testChangingParameters() {
+    const double k = 1.5;
+    const double charge = 0.1;
+    const double alpha = charge*charge/k;
+    Platform& platform = Platform::getPlatformByName("OpenCL");
+    
+    // Create the system.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    DrudeForce* drude = new DrudeForce();
+    drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
+    system.addForce(drude);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    
+    // Check the energy.
+    
+    VerletIntegrator integ(1.0);
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*k*3*3, state.getPotentialEnergy(), 1e-5);
+    
+    // Modify the parameters.
+    
+    const double k2 = 2.2;
+    const double charge2 = 0.3;
+    const double alpha2 = charge2*charge2/k2;
+    drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
+    drude->updateParametersInContext(context);
+    state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*k2*3*3, state.getPotentialEnergy(), 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerDrudeOpenCLKernelFactories();
@@ -164,6 +201,7 @@ int main(int argc, char* argv[]) {
         testSingleParticle();
         testAnisotropicParticle();
         testThole();
+        testChangingParameters();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/drude/platforms/reference/tests/TestReferenceDrudeForce.cpp

@@ -154,11 +154,49 @@ void testThole() {
     validateForce(system, positions, energySpring1+energySpring2+energyDipole);
 }
 
+void testChangingParameters() {
+    const double k = 1.5;
+    const double charge = 0.1;
+    const double alpha = charge*charge/k;
+    Platform& platform = Platform::getPlatformByName("Reference");
+    
+    // Create the system.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    DrudeForce* drude = new DrudeForce();
+    drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
+    system.addForce(drude);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    
+    // Check the energy.
+    
+    VerletIntegrator integ(1.0);
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*k*3*3, state.getPotentialEnergy(), 1e-5);
+    
+    // Modify the parameters.
+    
+    const double k2 = 2.2;
+    const double charge2 = 0.3;
+    const double alpha2 = charge2*charge2/k2;
+    drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
+    drude->updateParametersInContext(context);
+    state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*k2*3*3, state.getPotentialEnergy(), 1e-5);
+}
+
 int main() {
     try {
         testSingleParticle();
         testAnisotropicParticle();
         testThole();
+        testChangingParameters();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;