DPDIntegrator (#4799)

* Created DPDIntegrator class

* Reference implementation of DPDIntegrator

* Build neighbor list for DPDIntegrator

* Minor fixes

* Documentation for DPDIntegrator

* Python API for DPDIntegrator

* Preliminary OpenCL implementation of DPDIntegrator

* Enable USE_PERIODIC

* Use updated positions in DPD thermostat

* Working on neighbor list for OpenCL DPDIntegrator

* ReorderListener for particle types

* Serialization for DPDIntegrator

* CUDA implementation of DPDIntegrator

* HIP implementation of DPDIntegrator

* Fixed compile error in Python wrapper

* Fixed compile error in wrappers

* Fixed uninitialized memory in reference neighbor list

* Added DPDIntegrator to C++ API docs

* Fixed incorrect launch size

* Fixed nan in DPD random number generator

* Minor optimizations

* Improved load balancing

* Fixed an indexing error

* Neighbor list uses the maximum cutoff of any force

* Fixed HIP compilation error

* Fixed access to invalid memory

* Added test case for diffusion coefficient

* Try to debug segfaults on CI

* Debugging

* Debugging

* Debugging

* Debugging

* Debugging

* Debugging

* Possible fix

* Debugging

* Debugging

* Debugging

* Use correct block size on CPU OpenCL

* Workaround for bug in Intel's OpenCL for CPUs

* Removed an unnecessary define

* Removed debugging code

* Include Dart

* More Intel workarounds

* Workaround for error in NVIDIA OpenCL

platforms/hip/src/HipNonbondedUtilities.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) 2009-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2025 Stanford University and the Authors.      *
  * Portions copyright (c) 2020-2023 Advanced Micro Devices, Inc.              *
  * Authors: Peter Eastman, Nicholas Curtis                                    *
  * Contributors:                                                              *
@@ -62,7 +62,7 @@ public:
 };
 
 HipNonbondedUtilities::HipNonbondedUtilities(HipContext& context) : context(context), useCutoff(false), usePeriodic(false), useNeighborList(false), anyExclusions(false), usePadding(true),
-        blockSorter(NULL), pinnedCountBuffer(NULL), forceRebuildNeighborList(true), lastCutoff(0.0), groupFlags(0), canUsePairList(true), tilesAfterReorder(0) {
+        blockSorter(NULL), pinnedCountBuffer(NULL), forceRebuildNeighborList(true), groupFlags(0), canUsePairList(true), tilesAfterReorder(0) {
     // Decide how many thread blocks to use.
 
     string errorMessage = "Error initializing nonbonded utilities";
@@ -276,6 +276,7 @@ void HipNonbondedUtilities::initialize(const System& system) {
 
     // Create data structures for the neighbor list.
 
+    maxCutoff = getMaxCutoffDistance();
     if (useCutoff) {
         // Select a size for the arrays that hold the neighbor list.  We have to make a fairly
         // arbitrary guess, but if this turns out to be too small we'll increase it later.
@@ -429,7 +430,7 @@ void HipNonbondedUtilities::prepareInteractions(int forceGroups) {
     KernelSet& kernels = groupKernels[forceGroups];
     if (useCutoff && usePeriodic) {
         double4 box = context.getPeriodicBoxSize();
-        double minAllowedSize = 1.999999*kernels.cutoffDistance;
+        double minAllowedSize = 1.999999*maxCutoff;
         if (box.x < minAllowedSize || box.y < minAllowedSize || box.z < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
     }
@@ -440,15 +441,12 @@ void HipNonbondedUtilities::prepareInteractions(int forceGroups) {
 
     // Compute the neighbor list.
 
-    if (lastCutoff != kernels.cutoffDistance)
-        forceRebuildNeighborList = true;
     context.executeKernelFlat(kernels.findBlockBoundsKernel, &findBlockBoundsArgs[0], context.getPaddedNumAtoms(), context.getSIMDWidth());
     context.executeKernelFlat(kernels.computeSortKeysKernel, &computeSortKeysArgs[0], context.getNumAtomBlocks());
     blockSorter->sort(sortedBlocks);
     context.executeKernelFlat(kernels.sortBoxDataKernel, &sortBoxDataArgs[0], context.getNumAtoms(), 64);
     context.executeKernelFlat(kernels.findInteractingBlocksKernel, &findInteractingBlocksArgs[0], context.getNumAtomBlocks() * context.getSIMDWidth() * numTilesInBatch, findInteractingBlocksThreadBlockSize);
     forceRebuildNeighborList = false;
-    lastCutoff = kernels.cutoffDistance;
     context.executeKernelFlat(kernels.copyInteractionCountsKernel, &copyInteractionCountsArgs[0], 1, 1);
     hipEventRecord(downloadCountEvent, context.getCurrentStream());
 }
@@ -525,26 +523,23 @@ void HipNonbondedUtilities::setAtomBlockRange(double startFraction, double endFr
 
 void HipNonbondedUtilities::createKernelsForGroups(int groups) {
     KernelSet kernels;
-    double cutoff = 0.0;
     string source;
     for (int i = 0; i < 32; i++) {
         if ((groups&(1<<i)) != 0) {
-            cutoff = max(cutoff, groupCutoff[i]);
             source += groupKernelSource[i];
         }
     }
     kernels.hasForces = (source.size() > 0);
-    kernels.cutoffDistance = cutoff;
     kernels.source = source;
     kernels.forceKernel = kernels.energyKernel = kernels.forceEnergyKernel = NULL;
     if (useCutoff) {
-        double paddedCutoff = padCutoff(cutoff);
+        double paddedCutoff = padCutoff(maxCutoff);
         map<string, string> defines;
         defines["TILE_SIZE"] = context.intToString(HipContext::TileSize);
         defines["NUM_BLOCKS"] = context.intToString(context.getNumAtomBlocks());
         defines["NUM_ATOMS"] = context.intToString(context.getNumAtoms());
         defines["PADDED_NUM_ATOMS"] = context.intToString(context.getPaddedNumAtoms());
-        defines["PADDING"] = context.doubleToString(paddedCutoff-cutoff);
+        defines["PADDING"] = context.doubleToString(paddedCutoff-maxCutoff);
         defines["PADDED_CUTOFF"] = context.doubleToString(paddedCutoff);
         defines["PADDED_CUTOFF_SQUARED"] = context.doubleToString(paddedCutoff*paddedCutoff);
         defines["NUM_TILES_WITH_EXCLUSIONS"] = context.intToString(exclusionTiles.getSize());

platforms/hip/src/HipPlatform.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) 2008-2024 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
  * Portions copyright (c) 2020 Advanced Micro Devices, Inc.                   *
  * Authors: Peter Eastman, Nicholas Curtis                                    *
  * Contributors:                                                              *
@@ -104,6 +104,7 @@ HipPlatform::HipPlatform() {
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateDPDStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
     registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);

platforms/hip/tests/TestHipDPDIntegrator.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2025 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "HipTests.h"
+#include "TestDPDIntegrator.h"
+
+void runPlatformTests() {
+}

platforms/opencl/include/OpenCLNonbondedUtilities.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) 2009-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -334,7 +334,7 @@ private:
     std::vector<std::string> energyParameterDerivatives;
     std::map<int, double> groupCutoff;
     std::map<int, std::string> groupKernelSource;
-    double lastCutoff;
+    double maxCutoff;
     bool useCutoff, usePeriodic, deviceIsCpu, anyExclusions, usePadding, useNeighborList, forceRebuildNeighborList, useLargeBlocks, isAMD;
     int startTileIndex, startBlockIndex, numBlocks, maxExclusions, numForceThreadBlocks;
     int forceThreadBlockSize, interactingBlocksThreadBlockSize, groupFlags, numBlockSizes;
@@ -350,7 +350,6 @@ private:
 class OpenCLNonbondedUtilities::KernelSet {
 public:
     bool hasForces;
-    double cutoffDistance;
     std::string source;
     cl::Kernel forceKernel, energyKernel, forceEnergyKernel;
     cl::Kernel findBlockBoundsKernel;

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -138,6 +138,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new CommonIntegrateVariableLangevinStepKernel(name, platform, cl);
     if (name == IntegrateCustomStepKernel::Name())
         return new CommonIntegrateCustomStepKernel(name, platform, cl);
+    if (name == IntegrateDPDStepKernel::Name())
+        return new CommonIntegrateDPDStepKernel(name, platform, cl);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new CommonApplyAndersenThermostatKernel(name, platform, cl);
     if (name == IntegrateNoseHooverStepKernel::Name())

platforms/opencl/src/OpenCLNonbondedUtilities.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) 2009-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -53,7 +53,7 @@ public:
 };
 
 OpenCLNonbondedUtilities::OpenCLNonbondedUtilities(OpenCLContext& context) : context(context), useCutoff(false), usePeriodic(false), useNeighborList(false), anyExclusions(false), usePadding(true),
-        blockSorter(NULL), pinnedCountBuffer(NULL), pinnedCountMemory(NULL), forceRebuildNeighborList(true), lastCutoff(0.0), groupFlags(0), tilesAfterReorder(0) {
+        blockSorter(NULL), pinnedCountBuffer(NULL), pinnedCountMemory(NULL), forceRebuildNeighborList(true), groupFlags(0), tilesAfterReorder(0) {
     // Decide how many thread blocks and force buffers to use.
 
     deviceIsCpu = (context.getDevice().getInfo<CL_DEVICE_TYPE>() == CL_DEVICE_TYPE_CPU);
@@ -282,6 +282,7 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
 
     // Create data structures for the neighbor list.
 
+    maxCutoff = getMaxCutoffDistance();
     if (useCutoff) {
         // Select a size for the arrays that hold the neighbor list.  We have to make a fairly
         // arbitrary guess, but if this turns out to be too small we'll increase it later.
@@ -351,7 +352,7 @@ void OpenCLNonbondedUtilities::prepareInteractions(int forceGroups) {
     KernelSet& kernels = groupKernels[forceGroups];
     if (useCutoff && usePeriodic) {
         mm_float4 box = context.getPeriodicBoxSize();
-        double minAllowedSize = 1.999999*kernels.cutoffDistance;
+        double minAllowedSize = 1.999999*maxCutoff;
         if (box.x < minAllowedSize || box.y < minAllowedSize || box.z < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
     }
@@ -362,8 +363,6 @@ void OpenCLNonbondedUtilities::prepareInteractions(int forceGroups) {
 
     // Compute the neighbor list.
 
-    if (lastCutoff != kernels.cutoffDistance)
-        forceRebuildNeighborList = true;
     setPeriodicBoxArgs(context, kernels.findBlockBoundsKernel, 1);
     context.executeKernel(kernels.findBlockBoundsKernel, context.getNumAtomBlocks());
     context.executeKernel(kernels.computeSortKeysKernel, context.getNumAtomBlocks());
@@ -375,7 +374,6 @@ void OpenCLNonbondedUtilities::prepareInteractions(int forceGroups) {
     setPeriodicBoxArgs(context, kernels.findInteractingBlocksKernel, 0);
     context.executeKernel(kernels.findInteractingBlocksKernel, context.getNumAtoms(), interactingBlocksThreadBlockSize);
     forceRebuildNeighborList = false;
-    lastCutoff = kernels.cutoffDistance;
     context.getQueue().enqueueReadBuffer(interactionCount.getDeviceBuffer(), CL_FALSE, 0, sizeof(int), pinnedCountMemory, NULL, &downloadCountEvent);
     if (isAMD)
         context.getQueue().flush();
@@ -495,23 +493,20 @@ void OpenCLNonbondedUtilities::setAtomBlockRange(double startFraction, double en
 
 void OpenCLNonbondedUtilities::createKernelsForGroups(int groups) {
     KernelSet kernels;
-    double cutoff = 0.0;
     string source;
     for (int i = 0; i < 32; i++) {
         if ((groups&(1<<i)) != 0) {
-            cutoff = max(cutoff, groupCutoff[i]);
             source += groupKernelSource[i];
         }
     }
     kernels.hasForces = (source.size() > 0);
-    kernels.cutoffDistance = cutoff;
     kernels.source = source;
     if (useCutoff) {
-        double paddedCutoff = padCutoff(cutoff);
+        double paddedCutoff = padCutoff(maxCutoff);
         map<string, string> defines;
         defines["TILE_SIZE"] = context.intToString(OpenCLContext::TileSize);
         defines["NUM_ATOMS"] = context.intToString(context.getNumAtoms());
-        defines["PADDING"] = context.doubleToString(paddedCutoff-cutoff);
+        defines["PADDING"] = context.doubleToString(paddedCutoff-maxCutoff);
         defines["PADDED_CUTOFF"] = context.doubleToString(paddedCutoff);
         defines["PADDED_CUTOFF_SQUARED"] = context.doubleToString(paddedCutoff*paddedCutoff);
         defines["NUM_TILES_WITH_EXCLUSIONS"] = context.intToString(exclusionTiles.getSize());

platforms/opencl/src/OpenCLPlatform.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) 2008-2024 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -95,6 +95,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateDPDStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
     registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);

platforms/opencl/tests/TestOpenCLDPDIntegrator.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2025 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "OpenCLTests.h"
+#include "TestDPDIntegrator.h"
+
+void runPlatformTests() {
+}

platforms/reference/include/ReferenceDPDDynamics.h

+
+/* Portions copyright (c) 2006-2025 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 __ReferenceDPDDynamics_H__
+#define __ReferenceDPDDynamics_H__
+
+#include "ReferenceDynamics.h"
+#include "ReferenceNeighborList.h"
+#include "openmm/DPDIntegrator.h"
+#include "openmm/Vec3.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/windowsExport.h"
+#include <vector>
+
+namespace OpenMM {
+
+class OPENMM_EXPORT ReferenceDPDDynamics : public ReferenceDynamics {
+protected:
+      std::vector<Vec3> xPrime, oldx;
+      std::vector<double> masses, inverseMasses;
+      bool periodic;
+      Vec3 periodicBoxVectors[3];
+      NeighborList neighborList;
+      std::vector<int> particleType;
+      std::vector<std::vector<double> > frictionTable, cutoffTable;
+      double maxCutoff;
+public:
+    /**
+     * Constructor
+     * 
+     * @param system       the system to simulate
+     * @param integrator   the integrator being used
+     */
+    ReferenceDPDDynamics(const System& system, const DPDIntegrator& integrator);
+
+    /**
+     * Destructor
+     */
+    ~ReferenceDPDDynamics();
+
+    /**
+     * Set the periodic box vectors.
+     */
+    void setPeriodicBoxVectors(OpenMM::Vec3* vectors);
+    /**
+     * Get the maximum cutoff distance for any pair of types.
+     */
+    double getMaxCutoff() const;
+    /**
+     * Perform a time step, updating the positions and velocities.
+     * 
+     * @param context             the context this integrator is updating
+     * @param system              the System to be integrated
+     * @param atomCoordinates     atom coordinates
+     * @param velocities          velocities
+     * @param masses              atom masses
+     * @param tolerance           the constraint tolerance
+     */
+    void update(OpenMM::ContextImpl& context, std::vector<OpenMM::Vec3>& atomCoordinates,
+                std::vector<OpenMM::Vec3>& velocities, std::vector<double>& masses, double tolerance);
+
+    /**
+     * The first stage of the update algorithm.
+     */
+    virtual void updatePart1(int numParticles, std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces);
+
+    /**
+     * The second stage of the update algorithm.
+     */
+    virtual void updatePart2(int numParticles, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
+                             std::vector<OpenMM::Vec3>& xPrime);
+
+    /**
+     * The third stage of the update algorithm.
+     */
+    virtual void updatePart3(OpenMM::ContextImpl& context, int numParticles, std::vector<OpenMM::Vec3>& atomCoordinates,
+                             std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& xPrime);
+};
+
+} // namespace OpenMM
+
+#endif // __ReferenceDPDDynamics_H__

platforms/reference/include/ReferenceDynamics.h

@@ -136,6 +136,14 @@ class OPENMM_EXPORT ReferenceDynamics {
     
       double getTemperature() const;
 
+      /**---------------------------------------------------------------------------------------
+      
+         Set the temperature
+      
+         --------------------------------------------------------------------------------------- */
+    
+      void setTemperature(double temperature);
+
       /**---------------------------------------------------------------------------------------
       
          Update

platforms/reference/include/ReferenceKernels.h

@@ -68,6 +68,7 @@ class ReferenceVariableStochasticDynamics;
 class ReferenceVariableVerletDynamics;
 class ReferenceVerletDynamics;
 class ReferenceCustomDynamics;
+class ReferenceDPDDynamics;
 
 /**
  * This kernel is invoked at the beginning and end of force and energy computations.  It gives the
@@ -1509,6 +1510,42 @@ private:
     std::vector<std::vector<OpenMM::Vec3> > perDofValues; 
 };
 
+/**
+ * This kernel is invoked by DPDIntegrator to take one time step.
+ */
+class ReferenceIntegrateDPDStepKernel : public IntegrateDPDStepKernel {
+public:
+    ReferenceIntegrateDPDStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateDPDStepKernel(name, platform),
+        data(data), dynamics(NULL) {
+    }
+    ~ReferenceIntegrateDPDStepKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the DPDIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const DPDIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the DPDIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const DPDIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the DPDIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const DPDIntegrator& integrator);
+private:
+    ReferencePlatform::PlatformData& data;
+    ReferenceDPDDynamics* dynamics;
+    std::vector<double> masses;
+};
+
 /**
  * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
  */

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -104,6 +104,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceIntegrateVariableVerletStepKernel(name, platform, data);
     if (name == IntegrateCustomStepKernel::Name())
         return new ReferenceIntegrateCustomStepKernel(name, platform, data);
+    if (name == IntegrateDPDStepKernel::Name())
+        return new ReferenceIntegrateDPDStepKernel(name, platform, data);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new ReferenceApplyAndersenThermostatKernel(name, platform);
     if (name == ApplyMonteCarloBarostatKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -49,6 +49,7 @@
 #include "ReferenceCustomNonbondedIxn.h"
 #include "ReferenceCustomManyParticleIxn.h"
 #include "ReferenceCustomTorsionIxn.h"
+#include "ReferenceDPDDynamics.h"
 #include "ReferenceGayBerneForce.h"
 #include "ReferenceHarmonicBondIxn.h"
 #include "ReferenceLangevinMiddleDynamics.h"
@@ -2862,6 +2863,36 @@ void ReferenceIntegrateCustomStepKernel::setPerDofVariable(ContextImpl& context,
         perDofValues[variable][i] = values[i];
 }
 
+ReferenceIntegrateDPDStepKernel::~ReferenceIntegrateDPDStepKernel() {
+    if (dynamics != NULL)
+        delete dynamics;
+}
+
+void ReferenceIntegrateDPDStepKernel::initialize(const System& system, const DPDIntegrator& integrator) {
+    masses.resize(system.getNumParticles());
+    for (int i = 0; i < system.getNumParticles(); ++i)
+        masses[i] = system.getParticleMass(i);
+    dynamics = new ReferenceDPDDynamics(system, integrator);
+    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
+}
+
+void ReferenceIntegrateDPDStepKernel::execute(ContextImpl& context, const DPDIntegrator& integrator) {
+    dynamics->setTemperature(integrator.getTemperature());
+    dynamics->setDeltaT(integrator.getStepSize());
+    dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
+    dynamics->setVirtualSites(extractVirtualSites(context));
+    dynamics->setPeriodicBoxVectors(extractBoxVectors(context));
+    vector<Vec3>& posData = extractPositions(context);
+    vector<Vec3>& velData = extractVelocities(context);
+    dynamics->update(context, posData, velData, masses, integrator.getConstraintTolerance());
+    data.time += integrator.getStepSize();
+    data.stepCount++;
+}
+
+double ReferenceIntegrateDPDStepKernel::computeKineticEnergy(ContextImpl& context, const DPDIntegrator& integrator) {
+    return computeShiftedKineticEnergy(context, masses, 0.0);
+}
+
 ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
     if (thermostat)
         delete thermostat;

platforms/reference/src/ReferencePlatform.cpp

@@ -74,6 +74,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateDPDStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
     registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceDPDDynamics.cpp

+
+/* Portions copyright (c) 2025 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 "SimTKOpenMMUtilities.h"
+#include "ReferenceDPDDynamics.h"
+#include "ReferencePlatform.h"
+#include "ReferenceVirtualSites.h"
+#include "openmm/Integrator.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/DPDIntegrator.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/DPDIntegratorUtilities.h"
+#include "ReferenceForce.h"
+#include <set>
+
+using namespace std;
+using namespace OpenMM;
+
+ReferenceDPDDynamics::ReferenceDPDDynamics(const System& system, const DPDIntegrator& integrator) :
+           ReferenceDynamics(system.getNumParticles(), integrator.getStepSize(), integrator.getTemperature()) {
+    int numParticles = system.getNumParticles();
+    xPrime.resize(numParticles);
+    oldx.resize(numParticles);
+    periodic = system.usesPeriodicBoundaryConditions();
+    int numTypes;
+    DPDIntegratorUtilities::createTypeTables(integrator, numParticles, numTypes, particleType, frictionTable, cutoffTable, maxCutoff);
+}
+
+ReferenceDPDDynamics::~ReferenceDPDDynamics() {
+}
+
+void ReferenceDPDDynamics::setPeriodicBoxVectors(OpenMM::Vec3* vectors) {
+    periodicBoxVectors[0] = vectors[0];
+    periodicBoxVectors[1] = vectors[1];
+    periodicBoxVectors[2] = vectors[2];
+}
+
+double ReferenceDPDDynamics::getMaxCutoff() const {
+    return maxCutoff;
+}
+
+void ReferenceDPDDynamics::updatePart1(int numParticles, vector<Vec3>& velocities, vector<Vec3>& forces) {
+    for (int i = 0; i < numParticles; i++)
+        if (inverseMasses[i] != 0.0)
+            velocities[i] += (getDeltaT()*inverseMasses[i])*forces[i];
+}
+
+void ReferenceDPDDynamics::updatePart2(int numParticles, vector<Vec3>& atomCoordinates, vector<Vec3>& velocities,
+                                       vector<Vec3>& xPrime) {
+    const double halfdt = 0.5*getDeltaT();
+    const double kT = BOLTZ*getTemperature();
+
+    // First position update.
+
+    for (int i = 0; i < numParticles; i++) {
+        xPrime[i] = atomCoordinates[i];
+        if (inverseMasses[i] != 0.0)
+            xPrime[i] += velocities[i]*halfdt;
+    }
+
+    // Apply friction and noise to velocities.
+
+    vector<set<int> > exclusions(numParticles);
+    computeNeighborListVoxelHash(neighborList, numParticles, atomCoordinates, exclusions, periodicBoxVectors, periodic, maxCutoff, 0.0);
+    for (auto& pair : neighborList) {
+        int i = pair.first;
+        int j = pair.second;
+        if (masses[i] == 0.0 && masses[j] == 0.0)
+            continue;
+        int type1 = particleType[i];
+        int type2 = particleType[j];
+        double friction = frictionTable[type1][type2];
+        double cutoff = cutoffTable[type1][type2];
+        double deltaR[ReferenceForce::LastDeltaRIndex];
+        if (periodic)
+            ReferenceForce::getDeltaRPeriodic(xPrime[i], xPrime[j], periodicBoxVectors, deltaR);
+        else
+            ReferenceForce::getDeltaR(xPrime[i], xPrime[j], deltaR);
+        double r = deltaR[ReferenceForce::RIndex];
+        if (r >= cutoff)
+            continue;
+        double m = masses[i]*masses[j]/(masses[i]+masses[j]);
+        double omega = 1.0-(r/cutoff);
+        double vscale = exp(-getDeltaT()*2*friction*omega*omega);
+        double noisescale = sqrt(1-vscale*vscale);
+        Vec3 dir = Vec3(deltaR[ReferenceForce::XIndex], deltaR[ReferenceForce::YIndex], deltaR[ReferenceForce::ZIndex])/r;
+        Vec3 v = velocities[j]-velocities[i];
+        double dv = (1.0-vscale)*v.dot(dir) + noisescale*sqrt(kT/m)*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+        if (masses[i] != 0.0)
+            velocities[i] += (m/masses[i])*dv*dir;
+        if (masses[j] != 0.0)
+            velocities[j] -= (m/masses[j])*dv*dir;
+    }
+
+    // Second position update.
+
+    for (int i = 0; i < numParticles; i++)
+        if (inverseMasses[i] != 0.0) {
+            xPrime[i] = xPrime[i] + velocities[i]*halfdt;
+            oldx[i] = xPrime[i];
+        }
+}
+
+void ReferenceDPDDynamics::updatePart3(OpenMM::ContextImpl& context, int numParticles, vector<Vec3>& atomCoordinates,
+                                         vector<Vec3>& velocities, vector<Vec3>& xPrime) {
+    for (int i = 0; i < numParticles; i++) {
+        if (inverseMasses[i] != 0.0) {
+            velocities[i] += (xPrime[i]-oldx[i])/getDeltaT();
+            atomCoordinates[i] = xPrime[i];
+        }
+    }
+}
+
+void ReferenceDPDDynamics::update(ContextImpl& context, vector<Vec3>& atomCoordinates,
+                                    vector<Vec3>& velocities, vector<double>& masses, double tolerance) {
+    int numParticles = context.getSystem().getNumParticles();
+    ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+    if (this->masses.size() == 0) {
+        this->masses = masses;
+        inverseMasses.resize(masses.size());
+        for (int i = 0; i < masses.size(); i++) {
+            if (masses[i] == 0.0)
+                inverseMasses[i] = 0.0;
+            else
+                inverseMasses[i] = 1.0/masses[i];
+        }
+    }
+
+    // 1st update
+
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    updatePart1(numParticles, velocities, *data->forces);
+    if (referenceConstraintAlgorithm)
+        referenceConstraintAlgorithm->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
+
+    // 2nd update
+
+    updatePart2(numParticles, atomCoordinates, velocities, xPrime);
+    if (referenceConstraintAlgorithm)
+        referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
+
+    // 3rd update
+
+    updatePart3(context, numParticles, atomCoordinates, velocities, xPrime);
+    getVirtualSites().computePositions(context.getSystem(), atomCoordinates);
+    incrementTimeStep();
+}

platforms/reference/src/SimTKReference/ReferenceDynamics.cpp

@@ -131,6 +131,16 @@ double ReferenceDynamics::getTemperature() const {
    return _temperature;
 }
 
+/**---------------------------------------------------------------------------------------
+
+   Set the temperature
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceDynamics::setTemperature(double temperature) {
+    _temperature = temperature;
+}
+
 /**---------------------------------------------------------------------------------------
 
    Get ReferenceConstraint

platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp

@@ -84,7 +84,7 @@ typedef std::vector< VoxelItem > Voxel;
 class VoxelHash
 {
 public:
-    VoxelHash(double vsx, double vsy, double vsz, const Vec3* periodicBoxVectors, bool usePeriodic) :
+    VoxelHash(double vsx, double vsy, double vsz, const Vec3* periodicBoxVectors, const AtomLocationList& atomLocations, bool usePeriodic) :
             voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxVectors(periodicBoxVectors), usePeriodic(usePeriodic) {
         if (usePeriodic) {
             nx = (int) floor(periodicBoxVectors[0][0]/voxelSizeX+0.5);
@@ -94,6 +94,19 @@ public:
             voxelSizeY = periodicBoxVectors[1][1]/ny;
             voxelSizeZ = periodicBoxVectors[2][2]/nz;
         }
+        else {
+            Vec3 minPos = atomLocations[0];
+            Vec3 maxPos = atomLocations[0];
+            for (int i = 1; i < atomLocations.size(); i++) {
+                for (int j = 0; j < 3; j++) {
+                    minPos[j] = min(minPos[j], atomLocations[i][j]);
+                    maxPos[j] = max(maxPos[j], atomLocations[i][j]);
+                }
+            }
+            nx = max(1, min(500, (int) round((maxPos[0]-minPos[0])/voxelSizeY)));
+            ny = max(1, min(500, (int) round((maxPos[1]-minPos[1])/voxelSizeY)));
+            nz = max(1, min(500, (int) round((maxPos[2]-minPos[2])/voxelSizeY)));
+        }
     }
 
     void insert(const AtomIndex& item, const Vec3& location)
@@ -240,7 +253,7 @@ void OPENMM_EXPORT computeNeighborListVoxelHash(
         edgeSizeY = 0.5*periodicBoxVectors[1][1]/floor(periodicBoxVectors[1][1]/maxDistance);
         edgeSizeZ = 0.5*periodicBoxVectors[2][2]/floor(periodicBoxVectors[2][2]/maxDistance);
     }
-    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxVectors, usePeriodic);
+    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxVectors, atomLocations, usePeriodic);
     for (AtomIndex atomJ = 0; atomJ < (AtomIndex) nAtoms; ++atomJ) // use "j", because j > i for pairs
     {
         // 1) Find other atoms that are close to this one

platforms/reference/src/SimTKReference/SimTKOpenMMUtilities.cpp

@@ -311,6 +311,7 @@ void SimTKOpenMMUtilities::setRandomNumberSeed(uint32_t seed) {
     else
         _randomNumberSeed = seed;
    _randomInitialized = false;
+   nextGaussianIsValid = false;
 }
 
 void SimTKOpenMMUtilities::createCheckpoint(std::ostream& stream) {

platforms/reference/tests/TestReferenceDPDIntegrator.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2025 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "ReferenceTests.h"
+#include "TestDPDIntegrator.h"
+
+void runPlatformTests() {
+}

serialization/include/openmm/serialization/DPDIntegratorProxy.h

+#ifndef OPENMM_DPD_INTEGRATOR_PROXY_H_
+#define OPENMM_DPD_INTEGRATOR_PROXY_H_
+
+#include "openmm/serialization/XmlSerializer.h"
+
+namespace OpenMM {
+
+class DPDIntegratorProxy : public SerializationProxy {
+public:
+    DPDIntegratorProxy();
+    void serialize(const void* object, SerializationNode& node) const;
+    void* deserialize(const SerializationNode& node) const;
+};
+
+}
+
+#endif /*OPENMM_DPD_INTEGRATOR_PROXY_H_*/
\ No newline at end of file