Skip to content

GitLab

Commit 37a87d8a authored by Peter Eastman's avatar Peter Eastman
Browse files

Langevin integrator supports setting friction to 0

parent 9f717609
Show whitespace changes
Inline Side-by-side
Showing with 75 additions and 109 deletions
platforms/cpu/include/CpuLangevinDynamics.h
View file @ 37a87d8a
/* Portions copyright (c) 2013-2015 Stanford University and Simbios.
/* Portions copyright (c) 2013-2016 Stanford University and Simbios.
* Authors: Peter Eastman
* Contributors:
*
......@@ -43,12 +43,12 @@ public:
*
* @param numberOfAtoms number of atoms
* @param deltaT delta t for dynamics
* @param tau viscosity
* @param friction friction coefficient
* @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);
CpuLangevinDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM friction, RealOpenMM temperature, OpenMM::ThreadPool& threads, OpenMM::CpuRandom& random);
/**
* Destructor.
......
platforms/cpu/src/CpuKernels.cpp
View file @ 37a87d8a
......@@ -1285,8 +1285,7 @@ void CpuIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langevi
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 = new CpuLangevinDynamics(context.getSystem().getNumParticles(), stepSize, friction, temperature, data.threads, data.random);
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature;
prevFriction = friction;
......
platforms/cpu/src/CpuLangevinDynamics.cpp
View file @ 37a87d8a
/* Portions copyright (c) 2006-2015 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Authors: Peter Eastman
* Contributors:
*
......@@ -59,8 +59,8 @@ public:
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(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM friction, RealOpenMM temperature, ThreadPool& threads, CpuRandom& random) :
ReferenceStochasticDynamics(numberOfAtoms, deltaT, friction, temperature), threads(threads), random(random) {
}
CpuLangevinDynamics::~CpuLangevinDynamics() {
......@@ -120,11 +120,12 @@ void CpuLangevinDynamics::updatePart3(int numberOfAtoms, vector<RealVec>& atomCo
}
void CpuLangevinDynamics::threadUpdate1(int threadIndex) {
const RealOpenMM tau = getTau();
const RealOpenMM vscale = EXP(-getDeltaT()/tau);
const RealOpenMM fscale = (1-vscale)*tau;
RealOpenMM dt = getDeltaT();
RealOpenMM friction = getFriction();
const RealOpenMM vscale = EXP(-dt*friction);
const RealOpenMM fscale = (friction == 0 ? dt : (1-vscale)/friction);
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
const RealOpenMM noisescale = SQRT(kT*(1-vscale*vscale));
int start = threadIndex*numberOfAtoms/threads.getNumThreads();
int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
......
platforms/cuda/src/CudaKernels.cpp
View file @ 37a87d8a
......@@ -6497,11 +6497,10 @@ void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langev
if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Calculate the integration parameters.
double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
double kT = BOLTZ*temperature;
double vscale = exp(-stepSize/tau);
double fscale = (1-vscale)*tau;
double noisescale = sqrt(2*kT/tau)*sqrt(0.5*(1-vscale*vscale)*tau);
double vscale = exp(-stepSize*friction);
double fscale = (friction == 0 ? stepSize : (1-vscale)/friction);
double noisescale = sqrt(kT*(1-vscale*vscale));
if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
vector<double> p(params->getSize());
p[0] = vscale;
......@@ -6719,14 +6718,14 @@ double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
float maxStepSizeFloat = (float) maxStepSize;
double tol = integrator.getErrorTolerance();
float tolFloat = (float) tol;
double tau = integrator.getFriction() == 0.0 ? 0.0 : 1.0/integrator.getFriction();
float tauFloat = (float) tau;
double friction = integrator.getFriction();
float frictionFloat = (float) friction;
double kT = BOLTZ*integrator.getTemperature();
float kTFloat = (float) kT;
bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
void* argsSelect[] = {&numAtoms, &paddedNumAtoms, useDouble ? (void*) &maxStepSize : (void*) &maxStepSizeFloat,
useDouble ? (void*) &tol : (void*) &tolFloat,
useDouble ? (void*) &tau : (void*) &tauFloat,
useDouble ? (void*) &friction : (void*) &frictionFloat,
useDouble ? (void*) &kT : (void*) &kTFloat,
&cu.getIntegrationUtilities().getStepSize().getDevicePointer(),
&cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &params->getDevicePointer()};
......
platforms/cuda/src/kernels/langevin.cu
View file @ 37a87d8a
......@@ -78,7 +78,7 @@ extern "C" __global__ void integrateLangevinPart2(int numAtoms, real4* __restric
* Select the step size to use for the next step.
*/
extern "C" __global__ void selectLangevinStepSize(int numAtoms, int paddedNumAtoms, mixed maxStepSize, mixed errorTol, mixed tau, mixed kT, mixed2* __restrict__ dt,
extern "C" __global__ void selectLangevinStepSize(int numAtoms, int paddedNumAtoms, mixed maxStepSize, mixed errorTol, mixed friction, mixed kT, mixed2* __restrict__ dt,
const mixed4* __restrict__ velm, const long long* __restrict__ force, mixed* __restrict__ paramBuffer) {
// Calculate the error.
......@@ -119,9 +119,9 @@ extern "C" __global__ void selectLangevinStepSize(int numAtoms, int paddedNumAto
// Recalculate the integration parameters.
mixed vscale = EXP(-newStepSize/tau);
mixed fscale = (1-vscale)*tau;
mixed noisescale = SQRT(2*kT/tau)*SQRT(0.5f*(1-vscale*vscale)*tau);
mixed vscale = EXP(-newStepSize*friction);
mixed fscale = (friction == 0 ? newStepSize : (1-vscale)/friction);
mixed noisescale = SQRT(kT*(1-vscale*vscale));
params[VelScale] = vscale;
params[ForceScale] = fscale;
params[NoiseScale] = noisescale;
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 37a87d8a
......@@ -6733,11 +6733,10 @@ void OpenCLIntegrateLangevinStepKernel::execute(ContextImpl& context, const Lang
if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Calculate the integration parameters.
double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
double kT = BOLTZ*temperature;
double vscale = exp(-stepSize/tau);
double fscale = (1-vscale)*tau;
double noisescale = sqrt(2*kT/tau)*sqrt(0.5*(1-vscale*vscale)*tau);
double vscale = exp(-stepSize*friction);
double fscale = (friction == 0 ? stepSize : (1-vscale)/friction);
double noisescale = sqrt(kT*(1-vscale*vscale));
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
vector<cl_double> p(params->getSize());
p[0] = vscale;
......@@ -7015,13 +7014,13 @@ double OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context,
if (useDouble) {
selectSizeKernel.setArg<cl_double>(0, maxStepSize);
selectSizeKernel.setArg<cl_double>(1, integrator.getErrorTolerance());
selectSizeKernel.setArg<cl_double>(2, integrator.getFriction() == 0.0 ? 0.0 : 1.0/integrator.getFriction());
selectSizeKernel.setArg<cl_double>(2, integrator.getFriction());
selectSizeKernel.setArg<cl_double>(3, BOLTZ*integrator.getTemperature());
}
else {
selectSizeKernel.setArg<cl_float>(0, maxStepSizeFloat);
selectSizeKernel.setArg<cl_float>(1, (cl_float) integrator.getErrorTolerance());
selectSizeKernel.setArg<cl_float>(2, (cl_float) (integrator.getFriction() == 0.0 ? 0.0 : 1.0/integrator.getFriction()));
selectSizeKernel.setArg<cl_float>(2, (cl_float) integrator.getFriction());
selectSizeKernel.setArg<cl_float>(3, (cl_float) (BOLTZ*integrator.getTemperature()));
}
cl.executeKernel(selectSizeKernel, blockSize, blockSize);
......
platforms/opencl/src/kernels/langevin.cl
View file @ 37a87d8a
......@@ -72,7 +72,7 @@ __kernel void integrateLangevinPart2(__global real4* restrict posq, __global rea
* Select the step size to use for the next step.
*/
__kernel void selectLangevinStepSize(mixed maxStepSize, mixed errorTol, mixed tau, mixed kT, __global mixed2* restrict dt,
__kernel void selectLangevinStepSize(mixed maxStepSize, mixed errorTol, mixed friction, mixed kT, __global mixed2* restrict dt,
__global const mixed4* restrict velm, __global const real4* restrict force, __global mixed* restrict paramBuffer, __local mixed* restrict params, __local mixed* restrict error) {
// Calculate the error.
......@@ -110,9 +110,9 @@ __kernel void selectLangevinStepSize(mixed maxStepSize, mixed errorTol, mixed ta
// Recalculate the integration parameters.
mixed vscale = exp(-newStepSize/tau);
mixed fscale = (1-vscale)*tau;
mixed noisescale = sqrt(2*kT/tau)*sqrt(0.5f*(1-vscale*vscale)*tau);
mixed vscale = EXP(-newStepSize*friction);
mixed fscale = (friction == 0 ? newStepSize : (1-vscale)/friction);
mixed noisescale = SQRT(kT*(1-vscale*vscale));
params[VelScale] = vscale;
params[ForceScale] = fscale;
params[NoiseScale] = noisescale;
......
platforms/reference/include/ReferenceStochasticDynamics.h
View file @ 37a87d8a
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -36,7 +36,7 @@ class OPENMM_EXPORT ReferenceStochasticDynamics : public ReferenceDynamics {
std::vector<OpenMM::RealVec> xPrime;
std::vector<RealOpenMM> inverseMasses;
RealOpenMM _tau;
RealOpenMM friction;
public:
......@@ -46,12 +46,12 @@ class OPENMM_EXPORT ReferenceStochasticDynamics : public ReferenceDynamics {
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param tau viscosity
@param friction friction coefficient
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceStochasticDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM tau, RealOpenMM temperature);
ReferenceStochasticDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM friction, RealOpenMM temperature);
/**---------------------------------------------------------------------------------------
......@@ -63,13 +63,11 @@ class OPENMM_EXPORT ReferenceStochasticDynamics : public ReferenceDynamics {
/**---------------------------------------------------------------------------------------
Get tau
@return tau
Get friction coefficient
--------------------------------------------------------------------------------------- */
RealOpenMM getTau() const;
RealOpenMM getFriction() const;
/**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceVariableStochasticDynamics.h
View file @ 37a87d8a
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -35,7 +35,7 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
std::vector<OpenMM::RealVec> xPrime;
std::vector<RealOpenMM> inverseMasses;
RealOpenMM _tau, _accuracy;
RealOpenMM friction, _accuracy;
public:
......@@ -44,13 +44,13 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
Constructor
@param numberOfAtoms number of atoms
@param tau viscosity
@param friction friction coefficient
@param temperature temperature
@param accuracy required accuracy
--------------------------------------------------------------------------------------- */
ReferenceVariableStochasticDynamics(int numberOfAtoms, RealOpenMM tau, RealOpenMM temperature, RealOpenMM accuracy);
ReferenceVariableStochasticDynamics(int numberOfAtoms, RealOpenMM friction, RealOpenMM temperature, RealOpenMM accuracy);
/**---------------------------------------------------------------------------------------
......@@ -62,13 +62,11 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
/**---------------------------------------------------------------------------------------
Get tau
@return tau
Get friction coefficient
--------------------------------------------------------------------------------------- */
RealOpenMM getTau() const;
RealOpenMM getFriction() const;
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 37a87d8a
......@@ -2053,11 +2053,10 @@ void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const L
if (dynamics)
delete dynamics;
RealOpenMM tau = static_cast<RealOpenMM>(friction == 0.0 ? 0.0 : 1.0/friction);
dynamics = new ReferenceStochasticDynamics(
context.getSystem().getNumParticles(),
static_cast<RealOpenMM>(stepSize),
static_cast<RealOpenMM>(tau),
static_cast<RealOpenMM>(friction),
static_cast<RealOpenMM>(temperature));
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature;
......@@ -2142,8 +2141,7 @@ double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& contex
if (dynamics)
delete dynamics;
RealOpenMM tau = static_cast<RealOpenMM>(friction == 0.0 ? 0.0 : 1.0/friction);
dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), (RealOpenMM) tau, (RealOpenMM) temperature, (RealOpenMM) errorTol);
dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), (RealOpenMM) friction, (RealOpenMM) temperature, (RealOpenMM) errorTol);
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature;
prevFriction = friction;
......
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp
View file @ 37a87d8a
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -41,20 +41,15 @@ using namespace OpenMM;
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param tau viscosity(?)
@param friction friction coefficient
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceStochasticDynamics::ReferenceStochasticDynamics(int numberOfAtoms,
RealOpenMM deltaT, RealOpenMM tau,
RealOpenMM deltaT, RealOpenMM friction,
RealOpenMM temperature) :
ReferenceDynamics(numberOfAtoms, deltaT, temperature), _tau(tau) {
if (tau <= 0) {
std::stringstream message;
message << "illegal tau value: " << tau;
throw OpenMMException(message.str());
}
ReferenceDynamics(numberOfAtoms, deltaT, temperature), friction(friction) {
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
}
......@@ -77,21 +72,12 @@ ReferenceStochasticDynamics::~ReferenceStochasticDynamics() {
/**---------------------------------------------------------------------------------------
Get tau
@return tau
Get friction coefficient
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceStochasticDynamics::getTau() const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceStochasticDynamics::getTau";
// ---------------------------------------------------------------------------------------
return _tau;
RealOpenMM ReferenceStochasticDynamics::getFriction() const {
return friction;
}
/**---------------------------------------------------------------------------------------
......@@ -120,11 +106,12 @@ void ReferenceStochasticDynamics::updatePart1(int numberOfAtoms, vector<RealVec>
// perform first update
RealOpenMM tau = getTau();
const RealOpenMM vscale = EXP(-getDeltaT()/tau);
const RealOpenMM fscale = (1-vscale)*tau;
RealOpenMM dt = getDeltaT();
RealOpenMM friction = getFriction();
const RealOpenMM vscale = EXP(-dt*friction);
const RealOpenMM fscale = (friction == 0 ? dt : (1-vscale)/friction);
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
const RealOpenMM noisescale = SQRT(kT*(1-vscale*vscale));
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (inverseMasses[ii] != 0.0) {
......
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp
View file @ 37a87d8a
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -42,21 +42,16 @@ using namespace OpenMM;
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param tau viscosity(?)
@param friction friction coefficient
@param temperature temperature
@param accuracy required accuracy
--------------------------------------------------------------------------------------- */
ReferenceVariableStochasticDynamics::ReferenceVariableStochasticDynamics(int numberOfAtoms,
RealOpenMM tau, RealOpenMM temperature,
RealOpenMM friction, RealOpenMM temperature,
RealOpenMM accuracy) :
ReferenceDynamics(numberOfAtoms, 0.0f, temperature), _tau(tau), _accuracy(accuracy) {
if (tau <= 0) {
std::stringstream message;
message << "illegal tau value: " << tau;
throw OpenMMException(message.str());
}
ReferenceDynamics(numberOfAtoms, 0.0f, temperature), friction(friction), _accuracy(accuracy) {
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
}
......@@ -101,21 +96,12 @@ void ReferenceVariableStochasticDynamics::setAccuracy(RealOpenMM accuracy) {
/**---------------------------------------------------------------------------------------
Get tau
@return tau
Get friction coefficient
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceVariableStochasticDynamics::getTau() const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceVariableStochasticDynamics::getTau";
// ---------------------------------------------------------------------------------------
return _tau;
RealOpenMM ReferenceVariableStochasticDynamics::getFriction() const {
return friction;
}
/**---------------------------------------------------------------------------------------
......@@ -178,11 +164,12 @@ void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<
// perform first update
RealOpenMM tau = getTau();
const RealOpenMM vscale = EXP(-getDeltaT()/tau);
const RealOpenMM fscale = (1-vscale)*tau;
RealOpenMM dt = getDeltaT();
RealOpenMM friction = getFriction();
const RealOpenMM vscale = EXP(-dt*friction);
const RealOpenMM fscale = (friction == 0 ? dt : (1-vscale)/friction);
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
const RealOpenMM noisescale = SQRT(kT*(1-vscale*vscale));
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] != 0) {
......
tests/TestLangevinIntegrator.h
View file @ 37a87d8a
......@@ -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) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -74,9 +74,9 @@ void testSingleBond() {
integrator.step(1);
}
// Not set the friction to a tiny value and see if it conserves energy.
// Now set the friction to 0 and see if it conserves energy.
integrator.setFriction(5e-5);
integrator.setFriction(0.0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
......
tests/TestVariableLangevinIntegrator.h
View file @ 37a87d8a
......@@ -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) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -74,9 +74,9 @@ void testSingleBond() {
integrator.step(1);
}
// Now set the friction to a tiny value and see if it conserves energy.
// Now set the friction to 0 and see if it conserves energy.
integrator.setFriction(5e-5);
integrator.setFriction(0.0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
......
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