DrudeForce supports periodic boundary conditions (#4523)

* DrudeForce supports periodic boundary conditions

* Fixed uninitialized memory

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

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -134,7 +134,7 @@ public:
      */
     int addScreenedPair(int particle1, int particle2, double thole);
     /**
-     * Get the force field parameters for screened pair.
+     * Get the force field parameters for a screened pair.
      *
      * @param index           the index of the pair for which to get parameters
      * @param[out] particle1  the index within this Force of the first particle involved in the interaction
@@ -143,7 +143,7 @@ public:
      */
     void getScreenedPairParameters(int index, int& particle1, int& particle2, double& thole) const;
     /**
-     * Set the force field parameters for screened pair.
+     * Set the force field parameters for a screened pair.
      *
      * @param index      the index of the pair for which to get parameters
      * @param particle1  the index within this Force of the first particle involved in the interaction
@@ -162,15 +162,19 @@ public:
      * be used to add new particles or screenedPairs, only to change the parameters of existing ones.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Set whether this force should apply periodic boundary conditions when calculating displacements.
+     * Usually this is not appropriate for bonded forces, but there are situations when it can be useful.
+     * 
+     * Periodic boundary conditions are only applied to screened pairs.  They are never used for the
+     * force between a Drude particle and its parent particle, regardless of this setting.
+     */
+    void setUsesPeriodicBoundaryConditions(bool periodic);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.
-     *
-     * @returns true if nonbondedMethod uses PBC and false otherwise
      */
-    bool usesPeriodicBoundaryConditions() const {
-        return false;
-    }
+    bool usesPeriodicBoundaryConditions() const;
 protected:
     ForceImpl* createImpl() const;
 private:
@@ -178,6 +182,7 @@ private:
     class ScreenedPairInfo;
     std::vector<ParticleInfo> particles;
     std::vector<ScreenedPairInfo> screenedPairs;
+    bool usePeriodic;
 };
 
 /**

plugins/drude/openmmapi/src/DrudeForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -38,7 +38,7 @@
 using namespace OpenMM;
 using namespace std;
 
-DrudeForce::DrudeForce() {
+DrudeForce::DrudeForce() : usePeriodic(false) {
 }
 
 int DrudeForce::addParticle(int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34) {
@@ -105,3 +105,11 @@ ForceImpl* DrudeForce::createImpl() const {
 void DrudeForce::updateParametersInContext(Context& context) {
     dynamic_cast<DrudeForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
 }
+
+void DrudeForce::setUsesPeriodicBoundaryConditions(bool periodic) {
+    usePeriodic = periodic;
+}
+
+bool DrudeForce::usesPeriodicBoundaryConditions() const {
+    return usePeriodic;
+}

plugins/drude/platforms/common/src/CommonDrudeKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -158,6 +158,7 @@ void CommonCalcDrudeForceKernel::initialize(const System& system, const DrudeFor
         }
         pairParams.upload(paramVector);
         map<string, string> replacements;
+        replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
         replacements["PARAMS"] = cc.getBondedUtilities().addArgument(pairParams, "float2");
         cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonDrudeKernelSources::drudePairForce, replacements), force.getForceGroup());
     }

plugins/drude/platforms/common/src/kernels/drudePairForce.cc

@@ -7,6 +7,9 @@ real3 force4 = make_real3(0);
 // First pair.
 
 real3 delta = make_real3(pos1.x-pos3.x, pos1.y-pos3.y, pos1.z-pos3.z);
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(delta)
+#endif
 real rInv = RSQRT(dot(delta, delta));
 real r = RECIP(rInv);
 real u = drudeParams.x*r;
@@ -20,6 +23,9 @@ force3 -= f;
 // Second pair.
 
 delta = make_real3(pos1.x-pos4.x, pos1.y-pos4.y, pos1.z-pos4.z);
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(delta)
+#endif
 rInv = RSQRT(dot(delta, delta));
 r = RECIP(rInv);
 u = drudeParams.x*r;
@@ -33,6 +39,9 @@ force4 -= f;
 // Third pair.
 
 delta = make_real3(pos2.x-pos3.x, pos2.y-pos3.y, pos2.z-pos3.z);
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(delta)
+#endif
 rInv = RSQRT(dot(delta, delta));
 r = RECIP(rInv);
 u = drudeParams.x*r;
@@ -46,6 +55,9 @@ force3 -= f;
 // Fourth pair.
 
 delta = make_real3(pos2.x-pos4.x, pos2.y-pos4.y, pos2.z-pos4.z);
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(delta)
+#endif
 rInv = RSQRT(dot(delta, delta));
 r = RECIP(rInv);
 u = drudeParams.x*r;

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -36,6 +36,7 @@
 #include "SimTKOpenMMUtilities.h"
 #include "ReferenceConstraints.h"
 #include "ReferenceVirtualSites.h"
+#include "ReferenceForce.h"
 #include <set>
 
 using namespace OpenMM;
@@ -56,6 +57,11 @@ static vector<Vec3>& extractForces(ContextImpl& context) {
     return *data->forces;
 }
 
+static Vec3* extractBoxVectors(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return data->periodicBoxVectors;
+}
+
 static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
     return *data->constraints;
@@ -103,6 +109,7 @@ void ReferenceCalcDrudeForceKernel::initialize(const System& system, const Drude
     pairThole.resize(numPairs);
     for (int i = 0; i < numPairs; i++)
         force.getScreenedPairParameters(i, pair1[i], pair2[i], pairThole[i]);
+    periodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -170,6 +177,7 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
     // Compute the screened interaction between bonded dipoles.
     
     int numPairs = pair1.size();
+    Vec3* boxVectors = extractBoxVectors(context);
     for (int i = 0; i < numPairs; i++) {
         int dipole1 = pair1[i];
         int dipole2 = pair2[i];
@@ -181,8 +189,13 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
                 int p1 = dipole1Particles[j];
                 int p2 = dipole2Particles[k];
                 double chargeProduct = charge[dipole1]*charge[dipole2]*(j == k ? 1 : -1);
-                Vec3 delta = pos[p1]-pos[p2];
-                double r = sqrt(delta.dot(delta));
+                double deltaR[ReferenceForce::LastDeltaRIndex];
+                if (periodic)
+                    ReferenceForce::getDeltaRPeriodic(pos[p2], pos[p1], boxVectors, deltaR);
+                else
+                    ReferenceForce::getDeltaR(pos[p2], pos[p1], deltaR);
+                Vec3 delta(deltaR[ReferenceForce::XIndex], deltaR[ReferenceForce::YIndex], deltaR[ReferenceForce::ZIndex]);
+                double r = deltaR[ReferenceForce::RIndex];
                 double u = r*uscale;
                 double screening = 1.0 - (1.0+0.5*u)*exp(-u);
                 energy += ONE_4PI_EPS0*chargeProduct*screening/r;

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -77,6 +77,7 @@ private:
     std::vector<double> charge, polarizability, aniso12, aniso34;
     std::vector<int> pair1, pair2;
     std::vector<double> pairThole;
+    bool periodic;
 };
 
 /**

plugins/drude/tests/TestDrudeForce.h

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -127,26 +127,29 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
     return 1.0-(1.0+u/2)*exp(-u);
 }
 
-void testThole() {
+void testThole(bool periodic) {
     const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
     const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double thole = 2.5;
+    const double box = 2.0;
     System system;
     system.addParticle(1.0);
     system.addParticle(1.0);
     system.addParticle(1.0);
     system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(box, 0, 0), Vec3(0, box, 0), Vec3(0, 0, box));
     DrudeForce* drude = new DrudeForce();
     drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
     drude->addParticle(3, 2, -1, -1, -1, charge, alpha, 1, 1);
     drude->addScreenedPair(0, 1, thole);
+    drude->setUsesPeriodicBoundaryConditions(periodic);
     system.addForce(drude);
     vector<Vec3> positions(4);
     positions[0] = Vec3(0, 0, 0);
     positions[1] = Vec3(0, -0.5, 0);
-    positions[2] = Vec3(1, 0, 0);
-    positions[3] = Vec3(1, 0, 0.3);
+    positions[2] = Vec3(1.1, 0, 0);
+    positions[3] = Vec3(1.1, 0, 0.3);
     double energySpring1 = 0.5*k*0.5*0.5;
     double energySpring2 = 0.5*k*0.3*0.3;
     double energyDipole = 0.0;
@@ -154,6 +157,11 @@ void testThole() {
     for (int i = 0; i < 2; i++)
         for (int j = 2; j < 4; j++) {
             Vec3 delta = positions[i]-positions[j];
+            if (periodic) {
+                delta[0] -= box*floor(delta[0]/box+0.5);
+                delta[1] -= box*floor(delta[1]/box+0.5);
+                delta[2] -= box*floor(delta[2]/box+0.5);
+            }
             double r = sqrt(delta.dot(delta));
             energyDipole += ONE_4PI_EPS0*q[i]*q[j]*computeScreening(r, thole, alpha, alpha)/r;
         }
@@ -204,7 +212,8 @@ int main(int argc, char* argv[]) {
         setupKernels(argc, argv);
         testSingleParticle();
         testAnisotropicParticle();
-        testThole();
+        testThole(false);
+        testThole(true);
         testChangingParameters();
         runPlatformTests();
     }