Commit 8ee2c644 authored by peastman's avatar peastman
Browse files

CPU implementation of BAOABLangevinIntegrator

parent 11a479b0
/* Portions copyright (c) 2013-2019 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_BAOAB_DYNAMICS_H__
#define __CPU_BAOAB_DYNAMICS_H__
#include "ReferenceBAOABDynamics.h"
#include "CpuRandom.h"
#include "openmm/internal/ThreadPool.h"
#include "sfmt/SFMT.h"
namespace OpenMM {
class CpuBAOABDynamics : public ReferenceBAOABDynamics {
public:
/**
* Constructor.
*
* @param numberOfAtoms number of atoms
* @param deltaT delta t for dynamics
* @param friction friction coefficient
* @param temperature temperature
* @param threads thread pool for parallelizing computation
* @param random random number generator
*/
CpuBAOABDynamics(int numberOfAtoms, double deltaT, double friction, double temperature, OpenMM::ThreadPool& threads, OpenMM::CpuRandom& random);
/**
* Destructor.
*/
~CpuBAOABDynamics();
/**
* 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::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
/**
* Second update step.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates atom coordinates
* @param velocities velocities
* @param inverseMasses inverse atom masses
* @param xPrime xPrime
*/
void updatePart2(int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
/**
* Third update
*
* @param context the context this integrator is updating
* @param numberOfAtoms number of atoms
* @param atomCoordinates atom coordinates
* @param velocities velocities
* @param inverseMasses inverse atom masses
* @param xPrime xPrime
*/
void updatePart3(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
private:
void threadUpdate1(int threadIndex);
void threadUpdate2(int threadIndex);
void threadUpdate3(int threadIndex);
void threadUpdate4(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::Vec3* atomCoordinates;
OpenMM::Vec3* velocities;
OpenMM::Vec3* forces;
double* inverseMasses;
OpenMM::Vec3* xPrime;
};
} // namespace OpenMM
#endif // __CPU_BAOAB_DYNAMICS_H__
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2018 Stanford University and the Authors. *
* Portions copyright (c) 2013-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "CpuBAOABDynamics.h"
#include "CpuBondForce.h"
#include "CpuCustomGBForce.h"
#include "CpuCustomManyParticleForce.h"
......@@ -536,6 +537,47 @@ private:
double prevTemp, prevFriction, prevStepSize;
};
/**
* This kernel is invoked by BAOABLangevinIntegrator to take one time step.
*/
class CpuIntegrateBAOABStepKernel : public IntegrateBAOABStepKernel {
public:
CpuIntegrateBAOABStepKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : IntegrateBAOABStepKernel(name, platform),
data(data), dynamics(0) {
}
~CpuIntegrateBAOABStepKernel();
/**
* Initialize the kernel, setting up the particle masses.
*
* @param system the System this kernel will be applied to
* @param integrator the BAOABLangevinIntegrator this kernel will be used for
*/
void initialize(const System& system, const BAOABLangevinIntegrator& integrator);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the BAOABLangevinIntegrator this kernel is being used for
* @param forcesAreValid if the context has been modified since the last time step, this will be
* false to show that cached forces are invalid and must be recalculated.
* On exit, this should specify whether the cached forces are valid at the
* end of the step.
*/
void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the BAOABLangevinIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator);
private:
CpuPlatform::PlatformData& data;
CpuBAOABDynamics* dynamics;
std::vector<double> masses;
double prevTemp, prevFriction, prevStepSize;
};
} // namespace OpenMM
#endif /*OPENMM_CPUKERNELS_H_*/
......
/* Portions copyright (c) 2006-2019 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 "SimTKOpenMMUtilities.h"
#include "CpuBAOABDynamics.h"
using namespace OpenMM;
using namespace std;
CpuBAOABDynamics::CpuBAOABDynamics(int numberOfAtoms, double deltaT, double friction, double temperature, ThreadPool& threads, CpuRandom& random) :
ReferenceBAOABDynamics(numberOfAtoms, deltaT, friction, temperature), threads(threads), random(random) {
}
CpuBAOABDynamics::~CpuBAOABDynamics() {
}
void CpuBAOABDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates, vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& inverseMasses, vector<Vec3>& 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.
threads.execute([&] (ThreadPool& threads, int threadIndex) { threadUpdate1(threadIndex); });
threads.waitForThreads();
}
void CpuBAOABDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCoordinates, vector<Vec3>& velocities,
vector<double>& inverseMasses, vector<Vec3>& xPrime) {
// Record the parameters for the threads.
this->numberOfAtoms = numberOfAtoms;
this->atomCoordinates = &atomCoordinates[0];
this->velocities = &velocities[0];
this->inverseMasses = &inverseMasses[0];
this->xPrime = &xPrime[0];
// Signal the threads to start running and wait for them to finish.
threads.execute([&] (ThreadPool& threads, int threadIndex) { threadUpdate2(threadIndex); });
threads.waitForThreads();
}
void CpuBAOABDynamics::updatePart3(ContextImpl& context, int numberOfAtoms, vector<Vec3>& atomCoordinates, vector<Vec3>& velocities,
vector<double>& inverseMasses, vector<Vec3>& xPrime) {
// Record the parameters for the threads.
this->numberOfAtoms = numberOfAtoms;
this->atomCoordinates = &atomCoordinates[0];
this->velocities = &velocities[0];
this->inverseMasses = &inverseMasses[0];
this->xPrime = &xPrime[0];
// Signal the threads to start running and wait for them to finish.
threads.execute([&] (ThreadPool& threads, int threadIndex) { threadUpdate3(threadIndex); });
threads.waitForThreads();
context.calcForcesAndEnergy(true, false);
threads.execute([&] (ThreadPool& threads, int threadIndex) { threadUpdate4(threadIndex); });
threads.waitForThreads();
}
void CpuBAOABDynamics::threadUpdate1(int threadIndex) {
const double halfdt = 0.5*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) {
velocities[i] += (halfdt*inverseMasses[i])*forces[i];
xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
oldx[i] = xPrime[i];
}
}
}
void CpuBAOABDynamics::threadUpdate2(int threadIndex) {
const double halfdt = 0.5*getDeltaT();
const double kT = BOLTZ*getTemperature();
const double friction = getFriction();
const double vscale = exp(-getDeltaT()*friction);
const double noisescale = sqrt(1-vscale*vscale);
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) {
velocities[i] += (xPrime[i]-oldx[i])/halfdt;
Vec3 noise(random.getGaussianRandom(threadIndex), random.getGaussianRandom(threadIndex), random.getGaussianRandom(threadIndex));
velocities[i] = vscale*velocities[i] + noisescale*sqrt(kT*inverseMasses[i])*noise;
atomCoordinates[i] = xPrime[i];
xPrime[i] = atomCoordinates[i] + velocities[i]*halfdt;
oldx[i] = xPrime[i];
}
}
}
void CpuBAOABDynamics::threadUpdate3(int threadIndex) {
const double halfdt = 0.5*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) {
velocities[i] += (xPrime[i]-oldx[i])/halfdt;
atomCoordinates[i] = xPrime[i];
}
}
void CpuBAOABDynamics::threadUpdate4(int threadIndex) {
const double halfdt = 0.5*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)
velocities[i] += (halfdt*inverseMasses[i])*forces[i];
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2016 Stanford University and the Authors. *
* Portions copyright (c) 2013-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -61,5 +61,7 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
return new CpuCalcGayBerneForceKernel(name, platform, data);
if (name == IntegrateLangevinStepKernel::Name())
return new CpuIntegrateLangevinStepKernel(name, platform, data);
if (name == IntegrateBAOABStepKernel::Name())
return new CpuIntegrateBAOABStepKernel(name, platform, data);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2018 Stanford University and the Authors. *
* Portions copyright (c) 2013-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -1352,3 +1352,43 @@ void CpuIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langevi
double CpuIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
}
CpuIntegrateBAOABStepKernel::~CpuIntegrateBAOABStepKernel() {
if (dynamics)
delete dynamics;
}
void CpuIntegrateBAOABStepKernel::initialize(const System& system, const BAOABLangevinIntegrator& integrator) {
int numParticles = system.getNumParticles();
masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
data.random.initialize(integrator.getRandomNumberSeed(), data.threads.getNumThreads());
}
void CpuIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) {
double temperature = integrator.getTemperature();
double friction = integrator.getFriction();
double stepSize = integrator.getStepSize();
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& velData = extractVelocities(context);
if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Recreate the computation objects with the new parameters.
if (dynamics)
delete dynamics;
dynamics = new CpuBAOABDynamics(context.getSystem().getNumParticles(), stepSize, friction, temperature, data.threads, data.random);
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature;
prevFriction = friction;
prevStepSize = stepSize;
}
dynamics->update(context, posData, velData, masses, forcesAreValid, integrator.getConstraintTolerance());
ReferencePlatform::PlatformData* refData = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
refData->time += stepSize;
refData->stepCount++;
}
double CpuIntegrateBAOABStepKernel::computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.0);
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2018 Stanford University and the Authors. *
* Portions copyright (c) 2013-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -74,6 +74,7 @@ CpuPlatform::CpuPlatform() {
registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateBAOABStepKernel::Name(), factory);
platformProperties.push_back(CpuThreads());
platformProperties.push_back(CpuDeterministicForces());
int threads = getNumProcessors();
......
/* -------------------------------------------------------------------------- *
* 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) 2019 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 "CpuTests.h"
#include "TestBAOABLangevinIntegrator.h"
void runPlatformTests() {
}
......@@ -95,13 +95,13 @@ class OPENMM_EXPORT ReferenceBAOABDynamics : public ReferenceDynamics {
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param inverseMasses inverse atom masses
@param xPrime xPrime
--------------------------------------------------------------------------------------- */
virtual void updatePart1(int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, std::vector<OpenMM::Vec3>& xPrime);
std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
/**---------------------------------------------------------------------------------------
......@@ -110,12 +110,13 @@ class OPENMM_EXPORT ReferenceBAOABDynamics : public ReferenceDynamics {
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param masses atom masses
@param inverseMasses inverse atom masses
@param xPrime xPrime
--------------------------------------------------------------------------------------- */
virtual void updatePart2(int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
std::vector<double>& masses, std::vector<OpenMM::Vec3>& xPrime);
std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
/**---------------------------------------------------------------------------------------
......@@ -126,12 +127,13 @@ class OPENMM_EXPORT ReferenceBAOABDynamics : public ReferenceDynamics {
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param inverseMasses inverse atom masses
@param xPrime xPrime
--------------------------------------------------------------------------------------- */
virtual void updatePart3(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& velocities,
std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, std::vector<OpenMM::Vec3>& xPrime);
std::vector<OpenMM::Vec3>& forces, std::vector<double>& inverseMasses, std::vector<OpenMM::Vec3>& xPrime);
};
} // namespace OpenMM
......
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