Continuing to implement ReferencePlatform (implemented kinetic energy...

Continuing to implement ReferencePlatform (implemented kinetic energy calculation and LangevinIntegrator)

olla/include/kernels.h

@@ -273,7 +273,7 @@ public:
      * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
      * @return the kinetic energy of the system
      */
-    virtual double execute(const Stream& positions) = 0;
+    virtual double execute(const Stream& velocities) = 0;
 };
 
 } // namespace OpenMM

openmmapi/src/OpenMMContextImpl.cpp

@@ -60,7 +60,13 @@ OpenMMContextImpl::OpenMMContextImpl(OpenMMContext& owner, System& system, Integ
     positions = platform->createStream("atomPositions", system.getNumAtoms(), Stream::Double3);
     velocities = platform->createStream("atomVelocities", system.getNumAtoms(), Stream::Double3);
     forces = platform->createStream("atomForces", system.getNumAtoms(), Stream::Double3);
+    double zero[] = {0.0, 0.0, 0.0};
+    velocities.fillWithValue(&zero);
     kineticEnergyKernel = platform->createKernel(CalcKineticEnergyKernel::Name());
+    vector<double> masses(system.getNumAtoms());
+    for (int i = 0; i < masses.size(); ++i)
+        masses[i] = system.getAtomMass(i);
+    dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).initialize(masses);
     integrator.initialize(*this);
 }
 

openmmapi/src/State.cpp

@@ -39,32 +39,32 @@ double State::getTime() const {
     return time;
 }
 const vector<Vec3>& State::getPositions() const {
-    if (types&Positions == 0)
+    if ((types&Positions) == 0)
         throw OpenMMException("Invoked getPositions() on a State which does not contain positions.");
     return positions;
 }
 const vector<Vec3>& State::getVelocities() const {
-    if (types&Velocities == 0)
+    if ((types&Velocities) == 0)
         throw OpenMMException("Invoked getVelocities() on a State which does not contain velocities.");
     return velocities;
 }
 const vector<Vec3>& State::getForces() const {
-    if (types&Forces== 0)
+    if ((types&Forces) == 0)
         throw OpenMMException("Invoked getForces() on a State which does not contain forces.");
     return forces;
 }
 double State::getKineticEnergy() const {
-    if (types&Energy== 0)
+    if ((types&Energy) == 0)
         throw OpenMMException("Invoked getKineticEnergy() on a State which does not contain energies.");
     return ke;
 }
 double State::getPotentialEnergy() const {
-    if (types&Energy== 0)
+    if ((types&Energy) == 0)
         throw OpenMMException("Invoked getPotentialEnergy() on a State which does not contain energies.");
     return pe;
 }
 const map<string, double>& State::getParameters() const {
-    if (types&Parameters== 0)
+    if ((types&Parameters) == 0)
         throw OpenMMException("Invoked getParameters() on a State which does not contain parameters.");
     return parameters;
 }

platforms/reference/src/ReferenceKernels.cpp

@@ -38,7 +38,10 @@
 #include "SimTKReference/ReferenceLJCoulombIxn.h"
 #include "SimTKReference/ReferenceProperDihedralBond.h"
 #include "SimTKReference/ReferenceRbDihedralBond.h"
+#include "SimTKReference/ReferenceStochasticDynamics.h"
+#include "SimTKReference/ReferenceShakeAlgorithm.h"
 #include <cmath>
+#include <limits>
 
 using namespace OpenMM;
 using namespace std;
@@ -155,7 +158,7 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const vector<vector<int
 }
 
 void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
-    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData());
+    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData()); // Reference code needs to be made const correct
     RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) forces.getImpl()).getData();
     ReferenceBondForce refBondForce;
     ReferenceHarmonicBondIxn harmonicBond;
@@ -173,7 +176,7 @@ void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positi
 }
 
 double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions) {
-    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData());
+    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData()); // Reference code needs to be made const correct
     RealOpenMM** forceData = allocateRealArray(numAtoms, 3);
     int arraySize = max(max(max(max(numAtoms, numBonds), numAngles), numPeriodicTorsions), numRBTorsions);
     RealOpenMM* energyArray = new RealOpenMM[arraySize];
@@ -227,14 +230,56 @@ void ReferenceIntegrateVerletStepKernel::initialize(const vector<double>& masses
 void ReferenceIntegrateVerletStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double stepSize) {
     
 }
+#include <iostream>
+ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
+    if (dynamics)
+        delete dynamics;
+    if (shake)
+        delete shake;
+    if (masses)
+        delete[] masses;
+    if (constraintIndices)
+        disposeIntArray(constraintIndices, numConstraints);
+    if (shakeParameters)
+        disposeRealArray(shakeParameters, numConstraints);
+}
 
 void ReferenceIntegrateLangevinStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
         const vector<double>& constraintLengths) {
-    
+    this->masses = new RealOpenMM[masses.size()];
+    for (int i = 0; i < masses.size(); ++i)
+        this->masses[i] = masses[i];
+    numConstraints = constraintIndices.size();
+    this->constraintIndices = allocateIntArray(numConstraints, 2);
+    for (int i = 0; i < numConstraints; ++i) {
+        this->constraintIndices[i][0] = constraintIndices[i][0];
+        this->constraintIndices[i][1] = constraintIndices[i][1];
+    }
+    shakeParameters = allocateRealArray(constraintLengths.size(), 1);
+    for (int i = 0; i < constraintLengths.size(); ++i)
+        shakeParameters[i][0] = constraintLengths[i];
 }
 
 void ReferenceIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
-    
+    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) positions.getImpl()).getData();
+    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
+    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) forces.getImpl()).getData()); // Reference code needs to be made const correct
+    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
+        // Recreate the computation objects with the new parameters.
+        
+        if (dynamics) {
+            delete dynamics;
+            delete shake;
+        }
+        RealOpenMM tau = (friction == 0.0 ? 0.0 : 1.0/friction);
+        dynamics = new ReferenceStochasticDynamics(positions.getSize(), stepSize, tau, temperature);
+        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
+        dynamics->setReferenceShakeAlgorithm(shake);
+        prevTemp = temperature;
+        prevFriction = friction;
+        prevStepSize = stepSize;
+    }
+    dynamics->update(positions.getSize(), posData, velData, forceData, masses);
 }
 
 void ReferenceIntegrateBrownianStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
@@ -255,9 +300,13 @@ void ReferenceApplyAndersenThermostatKernel::execute(Stream& velocities, double
 }
 
 void ReferenceCalcKineticEnergyKernel::initialize(const vector<double>& masses) {
-    
+    this->masses = masses;
 }
 
-double ReferenceCalcKineticEnergyKernel::execute(const Stream& positions) {
-    return 0.0; // TODO implement correctly
+double ReferenceCalcKineticEnergyKernel::execute(const Stream& velocities) {
+    RealOpenMM** velData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) velocities.getImpl()).getData()); // Reference code needs to be made const correct
+    double energy = 0.0;
+    for (int i = 0; i < masses.size(); ++i)
+        energy += masses[i]*(velData[i][0]*velData[i][0]+velData[i][1]*velData[i][1]+velData[i][2]*velData[i][2]);
+    return 0.5*energy;
 }

platforms/reference/src/ReferenceKernels.h

@@ -35,6 +35,9 @@
 #include "kernels.h"
 #include "SimTKUtilities/SimTKOpenMMRealType.h"
 
+class ReferenceStochasticDynamics;
+class ReferenceShakeAlgorithm;
+
 namespace OpenMM {
 
 /**
@@ -158,8 +161,10 @@ public:
  */
 class ReferenceIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
 public:
-    ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform) : IntegrateLangevinStepKernel(name, platform) {
+    ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform) : IntegrateLangevinStepKernel(name, platform),
+        dynamics(0), shake(0), masses(0), shakeParameters(0), constraintIndices(0) {
     }
+    ~ReferenceIntegrateLangevinStepKernel();
     /**
      * Initialize the kernel, setting up all parameters related to integrator.
      * 
@@ -180,6 +185,14 @@ public:
      * @param stepSize           the integration step size
      */
     void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
+private:
+    ReferenceStochasticDynamics* dynamics;
+    ReferenceShakeAlgorithm* shake;
+    RealOpenMM* masses;
+    RealOpenMM** shakeParameters;
+    int** constraintIndices;
+    int numConstraints;
+    double prevTemp, prevFriction, prevStepSize;
 };
 
 /**
@@ -254,7 +267,9 @@ public:
      * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
      * @return the kinetic energy of the system
      */
-    double execute(const Stream& positions);
+    double execute(const Stream& velocities);
+private:
+    std::vector<double> masses;
 };
 
 } // namespace OpenMM

platforms/reference/src/SimTKReference/ReferenceDynamics.h

@@ -26,6 +26,7 @@
 #define __ReferenceDynamics_H__
 
 #include "ReferenceShakeAlgorithm.h"
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
 
 // ---------------------------------------------------------------------------------------
 

platforms/reference/src/SimTKReference/ReferenceShakeAlgorithm.cpp

@@ -68,13 +68,14 @@ ReferenceShakeAlgorithm::ReferenceShakeAlgorithm( int numberOfConstraints,
 
    // work arrays
 
+   if (_numberOfConstraints > 0) {
        _r_ij                       = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfConstraints, threeI, NULL, 
                                                                                         1, zero, "r_ij" );
     
        _d_ij2                      = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "dij_2" );
        _distanceTolerance          = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "distanceTolerance" );
        _reducedMasses              = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "reducedMasses" );
-
+   }
 }
 
 /**---------------------------------------------------------------------------------------
@@ -91,11 +92,13 @@ ReferenceShakeAlgorithm::~ReferenceShakeAlgorithm( ){
 
    // ---------------------------------------------------------------------------------------
 
+    if (_numberOfConstraints > 0) {
        SimTKOpenMMUtilities::freeTwoDRealOpenMMArray( _r_ij,  "r_ij" );
     
        SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _d_ij2, "d_ij2" );
        SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _distanceTolerance, "distanceTolerance" );
        SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _reducedMasses, "reducedMasses" );
+    }
 }
 
 /**---------------------------------------------------------------------------------------

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

@@ -567,6 +567,7 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
 
       referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime,
                                            inverseMasses );
+   }
 
    // copy xPrime -> atomCoordinates
 
@@ -575,7 +576,6 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
       atomCoordinates[ii][1] = xPrime[ii][1];
       atomCoordinates[ii][2] = xPrime[ii][2];
    }
-   }
 
    incrementTimeStep();
 

platforms/reference/src/SimTKUtilities/SimTKOpenMMUtilities.cpp

@@ -832,7 +832,7 @@ void SimTKOpenMMUtilities::Xfree( const char* name, char* fileName, int line, vo
 #ifdef UseGromacsMalloc
    return save_free( name, fileName, line, ptr );
 #else
-   delete (char*) ptr;
+   delete[] (char*) ptr;
    return;
 #endif
 }

platforms/reference/tests/TestReferenceKineticEnergy.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) 2008 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 the kernel to calculate kinetic energy.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "System.h"
+#include "VerletIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testCalcKE() {
+    ReferencePlatform platform;
+    System system(4, 0);
+    for (int i = 0; i < 4; ++i)
+        system.setAtomMass(i, i+1);
+    VerletIntegrator integrator(0.01);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> velocities(4);
+    velocities[0] = Vec3(1, 0, 0);
+    velocities[1] = Vec3(0, 1, 0);
+    velocities[2] = Vec3(0, 0, 2);
+    velocities[3] = Vec3(std::sqrt(2.0), 0, std::sqrt(2.0));
+    context.setVelocities(velocities);
+    State state = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(0.5*(1+2+4*3+4*4), state.getKineticEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testCalcKE();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/reference/tests/TestReferenceLangevinIntegrator.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) 2008 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 LangevinIntegrator.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "StandardMMForceField.h"
+#include "System.h"
+#include "LangevinIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testSingleBond() {
+    ReferencePlatform platform;
+    System system(2, 0);
+    system.setAtomMass(0, 2.0);
+    system.setAtomMass(1, 2.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
+    forceField->setBondParameters(0, 0, 1, 1.5, 1);
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    
+    // This is simply a damped harmonic oscillator, so compare it to the analytical solution.
+    
+    double freq = std::sqrt(1-0.05*0.05);
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Positions | State::Velocities);
+        double time = state.getTime();
+        double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
+        ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
+        ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
+        double expectedSpeed = -0.5*std::exp(-0.05*time)*(0.05*std::cos(freq*time)+freq*std::sin(freq*time));
+        ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
+        ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
+        integrator.step(1);
+    }
+    
+    // Not set the friction to a tiny value and see if it conserves energy.
+    
+    integrator.setFriction(5e-5);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
+    for (int i = 0; i < 1000; ++i) {
+        state = context.getState(State::Energy);
+        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
+//        std::cout << state.getKineticEnergy()+state.getPotentialEnergy() << std::endl;
+//        std::cout << state.getPositions()[1]<<" "<<state.getKineticEnergy()<<" "<<state.getPotentialEnergy()<<" "<<(state.getKineticEnergy()+state.getPotentialEnergy()) << std::endl;
+        integrator.step(1);
+    }
+}
+
+int main() {
+    try {
+        testSingleBond();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}