Skip to content

GitLab

Commit ddf1216d authored by peastman's avatar peastman
Browse files

Parallelized LangevinIntegrator for CPU platform

parent da998c51
Hide whitespace changes
Inline Side-by-side
Showing with 779 additions and 3 deletions
platforms/cpu/include/CpuKernels.h
View file @ ddf1216d
......@@ -33,6 +33,7 @@
* -------------------------------------------------------------------------- */
#include "CpuGBSAOBCForce.h"
#include "CpuLangevinDynamics.h"
#include "CpuNeighborList.h"
#include "CpuNonbondedForce.h"
#include "CpuPlatform.h"
......@@ -171,6 +172,43 @@ private:
CpuGBSAOBCForce obc;
};
/**
* This kernel is invoked by LangevinIntegrator to take one time step.
*/
class CpuIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public:
CpuIntegrateLangevinStepKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform),
data(data), dynamics(NULL) {
}
~CpuIntegrateLangevinStepKernel();
/**
* Initialize the kernel, setting up the particle masses.
*
* @param system the System this kernel will be applied to
* @param integrator the LangevinIntegrator this kernel will be used for
*/
void initialize(const System& system, const LangevinIntegrator& integrator);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the LangevinIntegrator this kernel is being used for
*/
void execute(ContextImpl& context, const LangevinIntegrator& integrator);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the LangevinIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator);
private:
CpuPlatform::PlatformData& data;
CpuLangevinDynamics* dynamics;
std::vector<RealOpenMM> masses;
double prevTemp, prevFriction, prevStepSize;
};
} // namespace OpenMM
#endif /*OPENMM_CPUKERNELS_H_*/
......
platforms/cpu/include/CpuLangevinDynamics.h 0 → 100644
View file @ ddf1216d
/* Portions copyright (c) 2013 Stanford University and Simbios.
* 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.
*/
#ifndef __CPU_LANGEVIN_DYNAMICS_H__
#define __CPU_LANGEVIN_DYNAMICS_H__
#include "ReferenceStochasticDynamics.h"
#include "CpuRandom.h"
#include "openmm/internal/ThreadPool.h"
#include "sfmt/SFMT.h"
// ---------------------------------------------------------------------------------------
class CpuLangevinDynamics : public ReferenceStochasticDynamics {
public:
class Update1Task;
class Update2Task;
/**
* Constructor.
*
* @param numberOfAtoms number of atoms
* @param deltaT delta t for dynamics
* @param tau viscosity
* @param temperature temperature
* @param threads thread pool for parallelizing computation
* @param random random number generator
*/
CpuLangevinDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM tau, RealOpenMM temperature, OpenMM::ThreadPool& threads, OpenMM::CpuRandom& random);
/**
* Destructor.
*/
~CpuLangevinDynamics();
/**
* First update step.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates atom coordinates
* @param velocities velocities
* @param forces forces
* @param inverseMasses inverse atom masses
* @param xPrime xPrime
*/
void updatePart1(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& inverseMasses, std::vector<OpenMM::RealVec>& xPrime);
/**
* Second update step.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates atom coordinates
* @param velocities velocities
* @param forces forces
* @param inverseMasses inverse atom masses
* @param xPrime xPrime
*/
void updatePart2(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& inverseMasses, std::vector<OpenMM::RealVec>& xPrime);
private:
void threadUpdate1(int threadIndex);
void threadUpdate2(int threadIndex);
OpenMM::ThreadPool& threads;
OpenMM::CpuRandom& random;
std::vector<OpenMM_SFMT::SFMT> threadRandom;
// The following variables are used to make information accessible to the individual threads.
int numberOfAtoms;
OpenMM::RealVec* atomCoordinates;
OpenMM::RealVec* velocities;
OpenMM::RealVec* forces;
RealOpenMM* inverseMasses;
OpenMM::RealVec* xPrime;
};
// ---------------------------------------------------------------------------------------
#endif // __CPU_LANGEVIN_DYNAMICS_H__
platforms/cpu/include/CpuPlatform.h
View file @ ddf1216d
......@@ -33,6 +33,7 @@
* -------------------------------------------------------------------------- */
#include "AlignedArray.h"
#include "CpuRandom.h"
#include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/ThreadPool.h"
......@@ -74,6 +75,7 @@ public:
std::vector<AlignedArray<float> > threadForce;
ThreadPool threads;
bool isPeriodic;
CpuRandom random;
};
} // namespace OpenMM
......
platforms/cpu/include/CpuRandom.h 0 → 100644
View file @ ddf1216d
#ifndef OPENMM_CPURANDOM_H_
#define OPENMM_CPURANDOM_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) 2013 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 "sfmt/SFMT.h"
#include <vector>
namespace OpenMM {
/**
* This class provides a multithreaded random number generator.
*/
class OPENMM_EXPORT CpuRandom {
public:
CpuRandom();
~CpuRandom();
void initialize(int seed, int numThreads);
float getGaussianRandom(int threadIndex);
float getUniformRandom(int threadIndex);
private:
bool hasInitialized;
int randomSeed;
std::vector<OpenMM_SFMT::SFMT*> threadRandom;
std::vector<float> nextGaussian;
std::vector<int> nextGaussianIsValid;
};
} // namespace OpenMM
#endif /*OPENMM_CPURANDOM_H_*/
platforms/cpu/src/CpuKernelFactory.cpp
View file @ ddf1216d
......@@ -45,5 +45,7 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
return new CpuCalcNonbondedForceKernel(name, platform, data);
if (name == CalcGBSAOBCForceKernel::Name())
return new CpuCalcGBSAOBCForceKernel(name, platform, data);
if (name == IntegrateLangevinStepKernel::Name())
return new CpuIntegrateLangevinStepKernel(name, platform, data);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());
}
platforms/cpu/src/CpuKernels.cpp
View file @ ddf1216d
......@@ -31,6 +31,7 @@
#include "CpuKernels.h"
#include "ReferenceBondForce.h"
#include "ReferenceConstraints.h"
#include "ReferenceKernelFactory.h"
#include "ReferenceKernels.h"
#include "ReferenceLJCoulomb14.h"
......@@ -64,6 +65,47 @@ static RealVec& extractBoxSize(ContextImpl& context) {
return *(RealVec*) data->periodicBoxSize;
}
static ReferenceConstraints& extractConstraints(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *(ReferenceConstraints*) data->constraints;
}
/**
* Compute the kinetic energy of the system, possibly shifting the velocities in time to account
* for a leapfrog integrator.
*/
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& masses, double timeShift) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context);
int numParticles = context.getSystem().getNumParticles();
// Compute the shifted velocities.
vector<RealVec> shiftedVel(numParticles);
for (int i = 0; i < numParticles; ++i) {
if (masses[i] > 0)
shiftedVel[i] = velData[i]+forceData[i]*(timeShift/masses[i]);
else
shiftedVel[i] = velData[i];
}
// Apply constraints to them.
vector<double> inverseMasses(numParticles);
for (int i = 0; i < numParticles; i++)
inverseMasses[i] = (masses[i] == 0 ? 0 : 1/masses[i]);
extractConstraints(context).applyToVelocities(posData, shiftedVel, inverseMasses, 1e-4);
// Compute the kinetic energy.
double energy = 0.0;
for (int i = 0; i < numParticles; ++i)
if (masses[i] > 0)
energy += masses[i]*(shiftedVel[i].dot(shiftedVel[i]));
return 0.5*energy;
}
CpuCalcForcesAndEnergyKernel::CpuCalcForcesAndEnergyKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data, ContextImpl& context) :
CalcForcesAndEnergyKernel(name, platform), data(data) {
// Create a Reference platform version of this kernel.
......@@ -460,3 +502,45 @@ void CpuCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, co
}
obc.setParticleParameters(particleParams);
}
CpuIntegrateLangevinStepKernel::~CpuIntegrateLangevinStepKernel() {
if (dynamics)
delete dynamics;
}
void CpuIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
int numParticles = system.getNumParticles();
masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
data.random.initialize(integrator.getRandomNumberSeed(), data.threads.getNumThreads());
}
void CpuIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
double temperature = integrator.getTemperature();
double friction = integrator.getFriction();
double stepSize = integrator.getStepSize();
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context);
if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Recreate the computation objects with the new parameters.
if (dynamics)
delete dynamics;
RealOpenMM tau = (friction == 0.0 ? 0.0 : 1.0/friction);
dynamics = new CpuLangevinDynamics(context.getSystem().getNumParticles(), stepSize, tau, temperature, data.threads, data.random);
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature;
prevFriction = friction;
prevStepSize = stepSize;
}
dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
ReferencePlatform::PlatformData* refData = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
refData->time += stepSize;
refData->stepCount++;
}
double CpuIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
}
platforms/cpu/src/CpuLangevinDynamics.cpp 0 → 100644
View file @ ddf1216d
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* 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 "SimTKOpenMMCommon.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
#include "CpuLangevinDynamics.h"
using namespace OpenMM;
using namespace std;
class CpuLangevinDynamics::Update1Task : public ThreadPool::Task {
public:
Update1Task(CpuLangevinDynamics& owner) : owner(owner) {
}
void execute(ThreadPool& threads, int threadIndex) {
owner.threadUpdate1(threadIndex);
}
CpuLangevinDynamics& owner;
};
class CpuLangevinDynamics::Update2Task : public ThreadPool::Task {
public:
Update2Task(CpuLangevinDynamics& owner) : owner(owner) {
}
void execute(ThreadPool& threads, int threadIndex) {
owner.threadUpdate2(threadIndex);
}
CpuLangevinDynamics& owner;
};
CpuLangevinDynamics::CpuLangevinDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM tau, RealOpenMM temperature, ThreadPool& threads, CpuRandom& random) :
ReferenceStochasticDynamics(numberOfAtoms, deltaT, tau, temperature), threads(threads), random(random) {
}
CpuLangevinDynamics::~CpuLangevinDynamics() {
}
void CpuLangevinDynamics::updatePart1(int numberOfAtoms, vector<RealVec>& atomCoordinates, vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& inverseMasses, vector<RealVec>& xPrime) {
// Record the parameters for the threads.
this->numberOfAtoms = numberOfAtoms;
this->atomCoordinates = &atomCoordinates[0];
this->velocities = &velocities[0];
this->forces = &forces[0];
this->inverseMasses = &inverseMasses[0];
this->xPrime = &xPrime[0];
// Signal the threads to start running and wait for them to finish.
Update1Task task(*this);
threads.execute(task);
threads.waitForThreads();
}
void CpuLangevinDynamics::updatePart2(int numberOfAtoms, vector<RealVec>& atomCoordinates, vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& inverseMasses, vector<RealVec>& xPrime) {
// Record the parameters for the threads.
this->numberOfAtoms = numberOfAtoms;
this->atomCoordinates = &atomCoordinates[0];
this->velocities = &velocities[0];
this->forces = &forces[0];
this->inverseMasses = &inverseMasses[0];
this->xPrime = &xPrime[0];
// Signal the threads to start running and wait for them to finish.
Update2Task task(*this);
threads.execute(task);
threads.waitForThreads();
}
void CpuLangevinDynamics::threadUpdate1(int threadIndex) {
const RealOpenMM tau = getTau();
const RealOpenMM vscale = EXP(-getDeltaT()/tau);
const RealOpenMM fscale = (1-vscale)*tau;
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
int start = threadIndex*numberOfAtoms/threads.getNumThreads();
int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
for (int i = start; i < end; i++) {
if (inverseMasses[i] != 0.0) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[i]);
RealVec noise(random.getGaussianRandom(threadIndex), random.getGaussianRandom(threadIndex), random.getGaussianRandom(threadIndex));
velocities[i] = velocities[i]*vscale + forces[i]*(fscale*inverseMasses[i]) + noise*(noisescale*sqrtInvMass);
}
}
}
void CpuLangevinDynamics::threadUpdate2(int threadIndex) {
const RealOpenMM dt = getDeltaT();
int start = threadIndex*numberOfAtoms/threads.getNumThreads();
int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
for (int i = start; i < end; i++) {
if (inverseMasses[i] != 0.0) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[i]);
xPrime[i] = atomCoordinates[i]+velocities[i]*dt;
}
}
}
platforms/cpu/src/CpuPlatform.cpp
View file @ ddf1216d
......@@ -53,6 +53,7 @@ CpuPlatform::CpuPlatform() {
registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
}
double CpuPlatform::getSpeed() const {
......
platforms/cpu/src/CpuRandom.cpp 0 → 100644
View file @ ddf1216d
/* Portions copyright (c) 2013 Stanford University and Simbios.
* 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 "CpuRandom.h"
#include "openmm/OpenMMException.h"
#include <cmath>
using namespace std;
using namespace OpenMM;
CpuRandom::CpuRandom() : hasInitialized(false) {
}
CpuRandom::~CpuRandom() {
for (int i = 0; i < threadRandom.size(); i++)
delete threadRandom[i];
}
void CpuRandom::initialize(int seed, int numThreads) {
if (hasInitialized) {
if (seed == randomSeed)
return; // Already initialized with the same seed.
throw OpenMMException("Random number generator initialized twice with different seeds");
}
randomSeed = seed;
hasInitialized = true;
threadRandom.resize(numThreads);
nextGaussian.resize(numThreads);
nextGaussianIsValid.resize(numThreads, false);
// Use a quick and dirty RNG to pick seeds for the real random number generator.
unsigned int r = (unsigned int) seed;
for (int i = 0; i < numThreads; i++) {
r = (1664525*r + 1013904223) & 0xFFFFFFFF;
threadRandom[i] = new OpenMM_SFMT::SFMT();
init_gen_rand(r, *threadRandom[i]);
}
}
float CpuRandom::getGaussianRandom(int threadIndex) {
if (nextGaussianIsValid[threadIndex]) {
nextGaussianIsValid[threadIndex] = false;
return nextGaussian[threadIndex];
}
// Use the polar form of the Box-Muller transformation to generate two Gaussian random numbers.
float x, y, r2;
do {
x = 2.0f*(float) genrand_real2(*threadRandom[threadIndex])-1.0f;
y = 2.0f*(float) genrand_real2(*threadRandom[threadIndex])-1.0f;
r2 = x*x + y*y;
} while (r2 >= 1.0f || r2 == 0.0f);
float multiplier = sqrtf((-2.0f*logf(r2))/r2);
nextGaussian[threadIndex] = y*multiplier;
nextGaussianIsValid[threadIndex] = true;
return x*multiplier;
}
float CpuRandom::getUniformRandom(int threadIndex) {
return genrand_real2(*threadRandom[threadIndex]);
}
platforms/cpu/tests/TestCpuLangevinIntegrator.cpp 0 → 100644
View file @ ddf1216d
/* -------------------------------------------------------------------------- *
* 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-2013 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 "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond() {
CpuPlatform platform;
System system;
system.addParticle(2.0);
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
Context 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);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 10000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= 10000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 6/std::sqrt(10000.0));
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-5);
}
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
CpuPlatform platform;
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
LangevinIntegrator integrator(300.0, 2.0, 0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities | State::Positions);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
integrator.setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
integrator.setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT_EQUAL_TOL(state1.getPositions()[i][j], state2.getPositions()[i][j], 1e-5);
ASSERT_EQUAL_TOL(state3.getPositions()[i][j], state4.getPositions()[i][j], 1e-5);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
int main() {
try {
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/include/ReferenceStochasticDynamics.h
View file @ ddf1216d
......@@ -31,7 +31,7 @@
class ReferenceStochasticDynamics : public ReferenceDynamics {
private:
protected:
std::vector<OpenMM::RealVec> xPrime;
std::vector<RealOpenMM> inverseMasses;
......@@ -99,7 +99,7 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void updatePart1( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
virtual void updatePart1( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& inverseMasses, std::vector<OpenMM::RealVec>& xPrime );
/**---------------------------------------------------------------------------------------
......@@ -114,7 +114,7 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void updatePart2( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
virtual void updatePart2( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& inverseMasses, std::vector<OpenMM::RealVec>& xPrime );
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment