Continuing reference implementation of RPMD

plugins/rpmd/CMakeLists.txt

@@ -148,7 +148,7 @@ IF(OPENMM_BUILD_STATIC_LIB)
     TARGET_LINK_LIBRARIES( ${STATIC_RPMD_TARGET} ${STATIC_TARGET} )
 ENDIF(OPENMM_BUILD_STATIC_LIB)
 
-#ADD_SUBDIRECTORY(platforms/reference/tests)
+ADD_SUBDIRECTORY(platforms/reference/tests)
 
 # Which hardware platforms to build
 IF(CUDA_FOUND)

plugins/rpmd/openmmapi/include/openmm/RPMDIntegrator.h

@@ -41,7 +41,17 @@
 namespace OpenMM {
 
 /**
- * This is an Integrator which simulates a System using RPMD dynamics.
+ * This is an Integrator which simulates a System using ring polymer molecular dynamics (RPMD).
+ * It simulates many copies of the System, with successive copies connected by harmonic
+ * springs to form a ring.  This allows certain quantum mechanical effects to be efficiently
+ * simulated.
+ * 
+ * Because this Integrator simulates many copies of the System at once, it must be used
+ * differently from other Integrators.  Instead of setting positions and velocities by
+ * calling methods of the Context, you should use the corresponding methods of the Integrator
+ * to set them for specific copies of the System.  Similarly, you should retrieve state information
+ * for particular copies by calling getState() on the Integrator.  Do not query the Context for
+ * state information.
  */
 
 class OPENMM_EXPORT RPMDIntegrator : public Integrator {

plugins/rpmd/platforms/reference/include/ReferenceRpmdKernelFactory.h

+#ifndef OPENMM_REFERENCERPMDKERNELFACTORY_H_
+#define OPENMM_REFERENCERPMDKERNELFACTORY_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) 2011 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/KernelFactory.h"
+
+namespace OpenMM {
+
+/**
+ * This KernelFactory creates kernels for the reference implementation of RPMDIntegrator.
+ */
+
+class ReferenceRpmdKernelFactory : public KernelFactory {
+public:
+    KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_REFERENCERPMDKERNELFACTORY_H_*/

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernelFactory.cpp

+/* -------------------------------------------------------------------------- *
+ *                              OpenMMAmoeba                                  *
+ * -------------------------------------------------------------------------- *
+ * 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) 2011 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "ReferenceRpmdKernelFactory.h"
+#include "ReferenceRpmdKernels.h"
+#include "ReferencePlatform.h"
+#include "openmm/internal/windowsExport.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/OpenMMException.h"
+
+using namespace OpenMM;
+
+extern "C" void registerPlatforms() {
+}
+
+#if defined(WIN32)
+    #include <windows.h>
+    extern "C" void registerKernelFactories();
+    BOOL WINAPI DllMain(HANDLE hModule, DWORD  ul_reason_for_call, LPVOID lpReserved) {
+        if (ul_reason_for_call == DLL_PROCESS_ATTACH)
+            registerKernelFactories();
+        return TRUE;
+    }
+#else
+    extern "C" void __attribute__((constructor)) registerKernelFactories();
+#endif
+
+extern "C" void registerKernelFactories() {
+    Platform& platform = Platform::getPlatformByName("Reference");
+    ReferenceRpmdKernelFactory* factory = new ReferenceRpmdKernelFactory();
+    platform.registerKernelFactory(IntegrateRPMDStepKernel::Name(), factory);
+}
+
+KernelImpl* ReferenceRpmdKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
+    if (name == IntegrateRPMDStepKernel::Name())
+        return new ReferenceIntegrateRPMDStepKernel(name, platform);
+    throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
+}

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp

@@ -66,7 +66,7 @@ void ReferenceIntegrateRPMDStepKernel::initialize(const System& system, const RP
         velocities[i].resize(numParticles);
         forces[i].resize(numParticles);
     }
-    fftpack_init_1d(&fft, numParticles);
+    fftpack_init_1d(&fft, numCopies);
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
 
@@ -115,7 +115,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     vector<t_complex> p(numCopies);
     vector<t_complex> q(numCopies);
     const RealOpenMM scale = 1.0/sqrt(numCopies);
-    const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/1000;
+    const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
     const RealOpenMM nkT = numCopies*BOLTZ*integrator.getTemperature();
     const RealOpenMM twown = 2.0*nkT/hbar;
     for (int particle = 0; particle < numParticles; particle++) {
@@ -132,7 +132,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
                 const RealOpenMM wt = wk*dt;
                 const RealOpenMM sinwt2 = sin(wt/2);
                 const RealOpenMM wm = wk*system.getParticleMass(particle);
-                const t_complex pprime = p[k] - q[k]*2.0*wm*sinwt2; // Advance momentum from t-dt/2 to t+dt/2
+                const t_complex pprime = p[k] - q[k]*(2.0*wm*sinwt2); // Advance momentum from t-dt/2 to t+dt/2
                 q[k] = p[k]*(2.0*sinwt2/wm) + q[k]*(1.0-4.0*sinwt2*sinwt2); // Advance position from t to t+dt
                 p[k] = pprime;
             }
@@ -150,7 +150,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     const RealOpenMM c1_0 = exp(-dt*integrator.getFriction());
     const RealOpenMM c2_0 = sqrt(1.0-c1_0*c1_0);
     for (int particle = 0; particle < numParticles; particle++) {
-        const RealOpenMM localc3 = c2_0*sqrt(system.getParticleMass(particle));
+        const RealOpenMM c3_0 = c2_0*sqrt(system.getParticleMass(particle)*nkT);
         for (int component = 0; component < 3; component++) {
             for (int k = 0; k < numCopies; k++)
                 p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
@@ -158,11 +158,11 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
             
             // Apply a local Langevin thermostat to the centroid mode.
 
-            p[0].re = p[0].re*c1_0 + localc3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+            p[0].re = p[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
 
             // Use critical damping white noise for the remaining modes.
             
-            for (int k = 1; k < numCopies/2; k++) {
+            for (int k = 1; k <= numCopies/2; k++) {
                 const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
                 const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
                 const RealOpenMM c1 = exp(-2.0*wk*dt);

plugins/rpmd/platforms/reference/tests/CMakeLists.txt

+#
+# Testing
+#
+ENABLE_TESTING()
+#INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/include)
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/openmmapi/include/openmm)
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src)
+#INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src/kernels)
+
+Set( SHARED_OPENMM_RPMD_TARGET OpenMMRPMD)
+Set( SHARED_CUDA_TARGET OpenMMRPMDCuda )
+
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d)
+    SET(SHARED_OPENMM_RPMD_TARGET ${SHARED_OPENMM_RPMD_TARGET}_d)
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+
+#LINK_DIRECTORIES
+
+# Automatically create tests using files named "Test*.cpp"
+FILE(GLOB TEST_PROGS "*Test*.cpp")
+FOREACH(TEST_PROG ${TEST_PROGS})
+    GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
+
+    # Link with shared library
+
+    ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
+    TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_OPENMM_RPMD_TARGET} )
+    ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
+ENDFOREACH(TEST_PROG ${TEST_PROGS})

plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.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) 2011 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Reference implementation of RPMDIntegrator.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/Platform.h"
+#include "openmm/System.h"
+#include "openmm/RPMDIntegrator.h"
+#include "SimTKUtilities/SimTKOpenMMUtilities.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testIntegration() {
+    const int numParticles = 5;
+    const int numCopies = 50;
+    const double temperature = 300.0;
+    System system;
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(i+1);
+        positions[i] = Vec3(i, 0, 0);
+    }
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    for (int i = 0; i < numParticles-1; i++)
+        bonds->addBond(i, i+1, 1.0, 100.0);
+    RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
+    Context context(system, integ);
+    context.setPositions(positions);
+    const int numSteps = 5000;
+    integ.step(1000);
+    vector<double> ke(numCopies, 0.0);
+    vector<double> rg(numParticles, 0.0);
+    const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
+    const double wn = numCopies*BOLTZ*temperature/hbar;
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++)
+            state[j] = integ.getState(j, State::Positions | State::Velocities | State::Energy);
+        for (int j = 0; j < numCopies; j++)
+            ke[j] += state[j].getKineticEnergy();
+        double totalEnergy = 0.0;
+        for (int j = 0; j < numCopies; j++) {
+            totalEnergy += state[j].getKineticEnergy()+state[j].getPotentialEnergy();
+            for (int k = 0; k < numParticles; k++) {
+                Vec3 delta = state[j].getPositions()[k]-state[j].getPositions()[j == 0 ? numCopies-1 : j-1];
+                totalEnergy += 0.5*system.getParticleMass(k)*wn*wn*delta.dot(delta);
+            }
+        }
+        for (int j = 0; j < numParticles; j++) {
+            double rg2 = 0.0;
+            for (int k = 0; k < numCopies; k++)
+                for (int m = 0; m < numCopies; m++) {
+                    Vec3 delta = state[k].getPositions()[j]-state[m].getPositions()[j];
+                    rg2 += delta.dot(delta);
+                }
+            rg[j] += sqrt(rg2/(2*numCopies*numCopies));
+        }
+    }
+    for (int i = 0; i < numCopies; i++) {
+        double value = ke[i]/numSteps;
+        double expected = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
+        printf("%d: %g %g %g\n", i, value, expected, value/expected);
+    }
+    for (int i = 0; i < numParticles; i++) {
+        double value = rg[i]/numSteps;
+        double expected = hbar/(2*sqrt(system.getParticleMass(i)*BOLTZ*temperature));
+        printf("%d: %g %g %g\n", i, value, expected, value/expected);
+    }
+}
+
+int main() {
+    try {
+        Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
+        testIntegration();
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "Done" << std::endl;
+    return 0;
+}