Radius of gyration force (#5031)

* Reference implementation of RGForce

* GPU implementation of RGForce

* Serialization

* Documentation

* Fix compilation error

* Fixed error building API docs

docs-source/api-c++/forces.rst

@@ -70,6 +70,7 @@ be combined in arbitrary ways.
     generated/MonteCarloBarostat
     generated/MonteCarloFlexibleBarostat
     generated/MonteCarloMembraneBarostat
+    generated/RGForce
     generated/RMSDForce
     generated/RPMDMonteCarloBarostat
 

docs-source/usersguide/theory/02_standard_forces.rst

@@ -790,3 +790,21 @@ This force is normally used with a CustomCVForce (see Section :numref:`customcvf
 One rarely wants a force whose energy exactly equals the RMSD, but there are many
 situations where it is useful to have a restraining or biasing force that depends
 on the RMSD in some way.
+
+RGForce
+*********
+
+RGForce computes the radius of gyration (Rg) of a set of particle positions:
+
+.. math::
+   \text{Rg} = \sqrt{\frac{\sum_{i} \| \mathbf{x}_i - \mathbf{x}_c \|^2}{N}}
+
+where :math:`\mathbf{x}_c` is the center position,
+
+.. math::
+   \mathbf{x}_c = \frac{\sum_{i} \mathbf{x}_i}{N}
+
+This force is normally used with a CustomCVForce (see Section :numref:`customcvforce`).
+One rarely wants a force whose energy exactly equals Rg, but there are many
+situations where it is useful to have a restraining or biasing force that depends
+on Rg in some way.

olla/include/openmm/kernels.h

@@ -59,6 +59,7 @@
 #include "openmm/PeriodicTorsionForce.h"
 #include "openmm/QTBIntegrator.h"
 #include "openmm/RBTorsionForce.h"
+#include "openmm/RGForce.h"
 #include "openmm/RMSDForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
@@ -1071,6 +1072,34 @@ public:
     virtual void copyParametersToContext(ContextImpl& context, const RMSDForce& force) = 0;
 };
 
+/**
+ * This kernel is invoked by RGForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcRGForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcRGForce";
+    }
+    CalcRGForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RGForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const RGForce& force) = 0;
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/OpenMM.h

@@ -70,6 +70,7 @@
 #include "openmm/PeriodicTorsionForce.h"
 #include "openmm/QTBIntegrator.h"
 #include "openmm/RBTorsionForce.h"
+#include "openmm/RGForce.h"
 #include "openmm/RMSDForce.h"
 #include "openmm/State.h"
 #include "openmm/System.h"

openmmapi/include/openmm/RGForce.h

+#ifndef OPENMM_RGFORCE_H_
+#define OPENMM_RGFORCE_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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 "Force.h"
+#include "Vec3.h"
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This is a force whose energy equals the radius of gyration (Rg) of a set of
+ * particles.  It is intended for use with CustomCVForce.  You will not normally
+ * want a force that exactly equals the radius of gyration, but there are many
+ * situations where it is useful to have a restraining or biasing force that
+ * depends on Rg in some way.
+ */
+
+class OPENMM_EXPORT RGForce : public Force {
+public:
+    /**
+     * Create an RGForce.
+     *
+     * @param particles    the indices of the particles to use when computing
+     *                     Rg.  If this is empty (the default), all particles in
+     *                     the system will be used.
+     */
+    explicit RGForce(const std::vector<int>& particles=std::vector<int>());
+    /**
+     * Get the indices of the particles to use when computing Rg.  If this
+     * is empty, all particles in the system will be used.
+     */
+    const std::vector<int>& getParticles() const {
+        return particles;
+    }
+    /**
+     * Set the indices of the particles to use when computing Rg.  If this
+     * is empty, all particles in the system will be used.
+     */
+    void setParticles(const std::vector<int>& particles);
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return false;
+    }
+protected:
+    ForceImpl* createImpl() const;
+private:
+    std::vector<int> particles;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_RGFORCE_H_*/

openmmapi/include/openmm/internal/RGForceImpl.h

+#ifndef OPENMM_RGFORCEIMPL_H_
+#define OPENMM_RGFORCEIMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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 "ForceImpl.h"
+#include "openmm/RGForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of RGForce.
+ */
+
+class RGForceImpl : public ForceImpl {
+public:
+    RGForceImpl(const RGForce& owner);
+    ~RGForceImpl();
+    void initialize(ContextImpl& context);
+    const RGForce& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context, bool& forcesInvalid) {
+        // This force field doesn't update the state directly.
+    }
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
+    std::map<std::string, double> getDefaultParameters() {
+        return std::map<std::string, double>(); // This force doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+private:
+    const RGForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_RGFORCEIMPL_H_*/

openmmapi/src/RGForce.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 "openmm/RGForce.h"
+#include "openmm/internal/RGForceImpl.h"
+
+using namespace OpenMM;
+using namespace std;
+
+RGForce::RGForce(const vector<int>& particles) : particles(particles) {
+}
+
+void RGForce::setParticles(const std::vector<int>& particles) {
+    this->particles = particles;
+}
+
+ForceImpl* RGForce::createImpl() const {
+    return new RGForceImpl(*this);
+}

openmmapi/src/RGForceImpl.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 "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/RGForceImpl.h"
+#include "openmm/kernels.h"
+#include <set>
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+RGForceImpl::RGForceImpl(const RGForce& owner) : owner(owner) {
+    forceGroup = owner.getForceGroup();
+}
+
+RGForceImpl::~RGForceImpl() {
+}
+
+void RGForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcRGForceKernel::Name(), context);
+
+    // Check for errors in the specification of particles.
+    const System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
+    set<int> particles;
+    for (int i : owner.getParticles()) {
+        if (i < 0 || i >= numParticles) {
+            stringstream msg;
+            msg << "RGForce: Illegal particle index for Rg: ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+        if (particles.find(i) != particles.end()) {
+            stringstream msg;
+            msg << "RGForce: Duplicated particle index for Rg: ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+        particles.insert(i);
+    }
+    kernel.getAs<CalcRGForceKernel>().initialize(context.getSystem(), owner);
+}
+
+double RGForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<forceGroup)) != 0)
+        return kernel.getAs<CalcRGForceKernel>().execute(context, includeForces, includeEnergy);
+    return 0.0;
+}
+
+vector<string> RGForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcRGForceKernel::Name());
+    return names;
+}

platforms/common/include/openmm/common/CommonKernels.h

@@ -1158,6 +1158,37 @@ private:
     ComputeKernel kernel1, kernel2;
 };
 
+/**
+ * This kernel is invoked by RGForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CommonCalcRGForceKernel : public CalcRGForceKernel {
+public:
+    CommonCalcRGForceKernel(std::string name, const Platform& platform, ComputeContext& cc) : CalcRGForceKernel(name, platform), cc(cc) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RGForce this kernel will be used for
+     */
+    void initialize(const System& system, const RGForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+private:
+    class ReorderListener;
+    ComputeContext& cc;
+    int blockSize;
+    ComputeArray particles, centerBuffer, rgBuffer;
+    ComputeKernel centerKernel, rgKernel, forceKernel;
+};
+
 /**
  * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
  */

platforms/common/src/CommonKernels.cpp

@@ -3592,6 +3592,90 @@ void CommonCalcRMSDForceKernel::copyParametersToContext(ContextImpl& context, co
     cc.invalidateMolecules(info);
 }
 
+class CommonCalcRGForceKernel::ReorderListener : public ComputeContext::ReorderListener {
+public:
+    ReorderListener(ComputeContext& cc, const vector<int>& particleIndices, ArrayInterface& particles) : cc(cc),
+            particleIndices(particleIndices), particles(particles) {
+    }
+    void execute() {
+        vector<int> particleVec(particles.getSize());
+        const vector<int>& order = cc.getAtomIndex();
+        vector<int> invOrder(cc.getPaddedNumAtoms());
+        for (int i = 0; i < order.size(); i++)
+            invOrder[order[i]] = i;
+        for (int i = 0; i < particleIndices.size(); i++)
+            particleVec[i] = invOrder[particleIndices[i]];
+        particles.upload(particleVec);
+    }
+private:
+    ComputeContext& cc;
+    const vector<int>& particleIndices;
+    ArrayInterface& particles;
+};
+
+void CommonCalcRGForceKernel::initialize(const System& system, const RGForce& force) {
+    // Create data structures.
+
+    ContextSelector selector(cc);
+    bool useDouble = cc.getUseDoublePrecision();
+    int elementSize = (useDouble ? sizeof(double) : sizeof(float));
+    int numParticles = force.getParticles().size();
+    if (numParticles == 0)
+        numParticles = system.getNumParticles();
+    particles.initialize<int>(cc, numParticles, "particles");
+    centerBuffer.initialize(cc, 3*(cc.getNumThreadBlocks()+1), elementSize, "centerBuffer");
+    rgBuffer.initialize(cc, cc.getNumThreadBlocks(), elementSize, "rgBuffer");
+
+    // Create the kernels.
+
+    blockSize = min(256, cc.getMaxThreadBlockSize());
+    map<string, string> defines;
+    defines["THREAD_BLOCK_SIZE"] = cc.intToString(blockSize);
+    defines["PADDED_NUM_ATOMS"] = cc.intToString(cc.getPaddedNumAtoms());
+    ComputeProgram program = cc.compileProgram(CommonKernelSources::rg, defines);
+    centerKernel = program->createKernel("computeCenterPosition");
+    rgKernel = program->createKernel("computeRg");
+    forceKernel = program->createKernel("computeForces");
+    centerKernel->addArg(numParticles);
+    centerKernel->addArg(cc.getPosq());
+    centerKernel->addArg(particles);
+    centerKernel->addArg(centerBuffer);
+    rgKernel->addArg(numParticles);
+    rgKernel->addArg(cc.getPosq());
+    rgKernel->addArg(particles);
+    rgKernel->addArg(centerBuffer);
+    rgKernel->addArg(rgBuffer);
+    forceKernel->addArg(numParticles);
+    forceKernel->addArg(cc.getPosq());
+    forceKernel->addArg(particles);
+    forceKernel->addArg(centerBuffer);
+    forceKernel->addArg(rgBuffer);
+    forceKernel->addArg(cc.getLongForceBuffer());
+    forceKernel->addArg(cc.getEnergyBuffer());
+
+    // Create the listener for updating the list of particles.
+
+    if (force.getParticles().size() == 0) {
+        vector<int> particleVec(numParticles);
+        for (int i = 0; i < numParticles; i++)
+            particleVec[i] = i;
+        particles.upload(particleVec);
+    }
+    else {
+        ReorderListener* listener = new ReorderListener(cc, force.getParticles(), particles);
+        cc.addReorderListener(listener);
+        listener->execute();
+    }
+}
+
+double CommonCalcRGForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    ContextSelector selector(cc);
+    centerKernel->execute(particles.getSize(), blockSize);
+    rgKernel->execute(particles.getSize(), blockSize);
+    forceKernel->execute(particles.getSize(), blockSize);
+    return 0.0;
+}
+
 void CommonApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
     ContextSelector selector(cc);
     randomSeed = thermostat.getRandomNumberSeed();

platforms/common/src/kernels/rg.cc

+/**
+ * Sum a value over all threads.
+ */
+DEVICE real reduceValue(real value, LOCAL_ARG volatile real* temp) {
+    const int thread = LOCAL_ID;
+    SYNC_THREADS;
+    temp[thread] = value;
+    SYNC_THREADS;
+    for (int step = 1; step < 32; step *= 2) {
+        if (thread+step < LOCAL_SIZE && thread%(2*step) == 0)
+            temp[thread] = temp[thread] + temp[thread+step];
+        SYNC_WARPS;
+    }
+    for (int step = 32; step < LOCAL_SIZE; step *= 2) {
+        if (thread+step < LOCAL_SIZE && thread%(2*step) == 0)
+            temp[thread] = temp[thread] + temp[thread+step];
+        SYNC_THREADS;
+    }
+    return temp[0];
+}
+
+/**
+ * First step of computing the radius of gyration: each thread block computes a contribution
+ * to the center position.
+ */
+KERNEL void computeCenterPosition(int numParticles, GLOBAL const real4* RESTRICT posq, GLOBAL const int* RESTRICT particles,
+        GLOBAL real* centerBuffer) {
+    LOCAL volatile real temp[THREAD_BLOCK_SIZE];
+    real3 center = make_real3(0);
+    for (int i = GLOBAL_ID; i < numParticles; i += GLOBAL_SIZE)
+        center += trimTo3(posq[particles[i]]);
+    center.x = reduceValue(center.x, temp);
+    center.y = reduceValue(center.y, temp);
+    center.z = reduceValue(center.z, temp);
+    if (LOCAL_ID == 0) {
+        centerBuffer[3*GROUP_ID] = center.x;
+        centerBuffer[3*GROUP_ID+1] = center.y;
+        centerBuffer[3*GROUP_ID+2] = center.z;
+    }
+}
+
+/**
+ * Second step of computing the radius of gyration: each thread block computes a contribution to Rg^2.
+ */
+KERNEL void computeRg(int numParticles, GLOBAL const real4* RESTRICT posq, GLOBAL const int* RESTRICT particles,
+        GLOBAL real* centerBuffer, GLOBAL real* rgBuffer) {
+    LOCAL volatile real temp[THREAD_BLOCK_SIZE];
+
+    // Sum the contributions computed in the previous kernel to find the center position.
+
+    real3 center = make_real3(0);
+    for (int i = LOCAL_ID; i < NUM_GROUPS; i += LOCAL_SIZE) {
+        center.x += centerBuffer[3*i];
+        center.y += centerBuffer[3*i+1];
+        center.z += centerBuffer[3*i+2];
+    }
+    center.x = reduceValue(center.x, temp)/numParticles;
+    center.y = reduceValue(center.y, temp)/numParticles;
+    center.z = reduceValue(center.z, temp)/numParticles;
+    if (GLOBAL_ID == 0) {
+        // Save the result so we don't need to compute it again in the next kernel.
+
+        centerBuffer[3*NUM_GROUPS] = center.x;
+        centerBuffer[3*NUM_GROUPS+1] = center.y;
+        centerBuffer[3*NUM_GROUPS+2] = center.z;
+    }
+
+    // Compute the contributions to Rg^2.
+
+    real rg2 = 0;
+    for (int i = GLOBAL_ID; i < numParticles; i += GLOBAL_SIZE) {
+        real3 delta = trimTo3(posq[particles[i]]) - center;
+        rg2 += dot(delta, delta);
+    }
+    rg2 = reduceValue(rg2, temp);
+    if (LOCAL_ID == 0)
+        rgBuffer[GROUP_ID] = rg2;
+}
+
+/**
+ * Third step of computing the radius of gyration: compute the forces.
+ */
+KERNEL void computeForces(int numParticles, GLOBAL const real4* RESTRICT posq, GLOBAL const int* RESTRICT particles,
+        GLOBAL const real* centerBuffer, GLOBAL real* rgBuffer, GLOBAL mm_ulong* RESTRICT forceBuffers, GLOBAL mixed* RESTRICT energyBuffer) {
+    LOCAL volatile real temp[THREAD_BLOCK_SIZE];
+    real3 center = make_real3(centerBuffer[3*NUM_GROUPS], centerBuffer[3*NUM_GROUPS+1], centerBuffer[3*NUM_GROUPS+2]);
+
+    // Sum the contributions computed in the previous kernel to find Rg.
+
+    real rg2 = 0;
+    for (int i = LOCAL_ID; i < NUM_GROUPS; i += LOCAL_SIZE)
+        rg2 += rgBuffer[i];
+    real rg = SQRT(reduceValue(rg2, temp)/numParticles);
+    if (GLOBAL_ID == 0)
+        energyBuffer[0] += rg;
+
+    // Compute the forces.
+
+    real scale = 1/(rg*numParticles);
+    for (int i = GLOBAL_ID; i < numParticles; i += GLOBAL_SIZE) {
+        int particle = particles[i];
+        real3 delta = trimTo3(posq[particle]) - center;
+        ATOMIC_ADD(&forceBuffers[particle], (mm_ulong) realToFixedPoint(-scale*delta.x));
+        ATOMIC_ADD(&forceBuffers[particle+PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(-scale*delta.y));
+        ATOMIC_ADD(&forceBuffers[particle+2*PADDED_NUM_ATOMS], (mm_ulong) realToFixedPoint(-scale*delta.z));
+    }
+}

platforms/cuda/src/CudaKernelFactory.cpp

@@ -123,6 +123,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcATMForceKernel(name, platform, cu);
     if (name == CalcCustomCPPForceKernel::Name())
         return new CommonCalcCustomCPPForceKernel(name, platform, context, cu);
+    if (name == CalcRGForceKernel::Name())
+        return new CommonCalcRGForceKernel(name, platform, cu);
     if (name == CalcRMSDForceKernel::Name())
         return new CommonCalcRMSDForceKernel(name, platform, cu);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/cuda/src/CudaPlatform.cpp

@@ -94,6 +94,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
     registerKernelFactory(CalcATMForceKernel::Name(), factory);
+    registerKernelFactory(CalcRGForceKernel::Name(), factory);
     registerKernelFactory(CalcRMSDForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);

platforms/cuda/tests/TestCudaRGForce.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 "CudaTests.h"
+#include "TestRGForce.h"
+
+void runPlatformTests() {
+}

platforms/hip/src/HipKernelFactory.cpp

@@ -122,6 +122,8 @@ KernelImpl* HipKernelFactory::createKernelImpl(std::string name, const Platform&
         return new HipCalcCustomCVForceKernel(name, platform, cu);
     if (name == CalcCustomCPPForceKernel::Name())
         return new CommonCalcCustomCPPForceKernel(name, platform, context, cu);
+    if (name == CalcRGForceKernel::Name())
+        return new CommonCalcRGForceKernel(name, platform, cu);
     if (name == CalcRMSDForceKernel::Name())
         return new CommonCalcRMSDForceKernel(name, platform, cu);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/hip/src/HipPlatform.cpp

@@ -94,6 +94,7 @@ HipPlatform::HipPlatform() {
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
+    registerKernelFactory(CalcRGForceKernel::Name(), factory);
     registerKernelFactory(CalcRMSDForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);

platforms/hip/tests/TestHipRGForce.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 "TestRGForce.h"
+
+void runPlatformTests() {
+}

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -121,6 +121,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcATMForceKernel(name, platform, cl);
     if (name == CalcCustomCPPForceKernel::Name())
         return new CommonCalcCustomCPPForceKernel(name, platform, context, cl);
+    if (name == CalcRGForceKernel::Name())
+        return new CommonCalcRGForceKernel(name, platform, cl);
     if (name == CalcRMSDForceKernel::Name())
         return new CommonCalcRMSDForceKernel(name, platform, cl);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/opencl/src/OpenCLPlatform.cpp

@@ -85,6 +85,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
     registerKernelFactory(CalcATMForceKernel::Name(), factory);
+    registerKernelFactory(CalcRGForceKernel::Name(), factory);
     registerKernelFactory(CalcRMSDForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);

platforms/opencl/tests/TestOpenCLRGForce.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 "TestRGForce.h"
+
+void runPlatformTests() {
+}