Skip to content

GitLab

Commit f2276667 authored by Yutong Zhao's avatar Yutong Zhao
Browse files

Merge branch 'master' of https://github.com/SimTk/openmm

parents ec39f6ff 7521da4a
Hide whitespace changes
Inline Side-by-side
Showing with 815 additions and 71 deletions
README.md 0 → 100644
View file @ f2276667
## OpenMM: A High Performance Molecular Dynamics Library
Introduction
------------
[OpenMM](https://simtk.org/home/openmm) is a toolkit for molecular simulation. It can be used either as a stand-alone application for running simulations, or as a library you call from your own code. It
provides a combination of extreme flexibility (through custom forces and integrators), openness, and high performance (especially on recent GPUs) that make it truly unique among simulation codes.
Getting Help
------------
Need Help? Check out the [documentation](https://simtk.org/docman/?group_id=161) and [discussion forums](https://simtk.org/forums/viewforum.php?f=161).
[C++ API Reference](https://simtk.org/api_docs/openmm/api5_0/c++/)
[Python API Reference](https://simtk.org/api_docs/openmm/api5_0/python/)
libraries/lbfgs/include/lbfgs.h
View file @ f2276667
...@@ -33,6 +33,16 @@ ...@@ -33,6 +33,16 @@
extern "C" { extern "C" {
#endif/*__cplusplus*/ #endif/*__cplusplus*/
#ifdef _MSC_VER
#if defined(OPENMM_BUILDING_SHARED_LIBRARY)
#define WINDOWS_EXPORT __declspec(dllexport)
#else
#define WINDOWS_EXPORT
#endif
#else
#define WINDOWS_EXPORT
#endif
/* /*
* The default precision of floating point values is 64bit (double). * The default precision of floating point values is 64bit (double).
*/ */
...@@ -474,7 +484,7 @@ In this formula, ||.|| denotes the Euclidean norm. ...@@ -474,7 +484,7 @@ In this formula, ||.|| denotes the Euclidean norm.
* minimization process terminates without an error. A * minimization process terminates without an error. A
* non-zero value indicates an error. * non-zero value indicates an error.
*/ */
int lbfgs( int WINDOWS_EXPORT lbfgs(
int n, int n,
lbfgsfloatval_t *x, lbfgsfloatval_t *x,
lbfgsfloatval_t *ptr_fx, lbfgsfloatval_t *ptr_fx,
...@@ -492,7 +502,7 @@ int lbfgs( ...@@ -492,7 +502,7 @@ int lbfgs(
* *
* @param param The pointer to the parameter structure. * @param param The pointer to the parameter structure.
*/ */
void lbfgs_parameter_init(lbfgs_parameter_t *param); void WINDOWS_EXPORT lbfgs_parameter_init(lbfgs_parameter_t *param);
/** /**
* Allocate an array for variables. * Allocate an array for variables.
...@@ -505,7 +515,7 @@ void lbfgs_parameter_init(lbfgs_parameter_t *param); ...@@ -505,7 +515,7 @@ void lbfgs_parameter_init(lbfgs_parameter_t *param);
* *
* @param n The number of variables. * @param n The number of variables.
*/ */
lbfgsfloatval_t* lbfgs_malloc(int n); lbfgsfloatval_t WINDOWS_EXPORT * lbfgs_malloc(int n);
/** /**
* Free an array of variables. * Free an array of variables.
...@@ -513,7 +523,7 @@ lbfgsfloatval_t* lbfgs_malloc(int n); ...@@ -513,7 +523,7 @@ lbfgsfloatval_t* lbfgs_malloc(int n);
* @param x The array of variables allocated by ::lbfgs_malloc * @param x The array of variables allocated by ::lbfgs_malloc
* function. * function.
*/ */
void lbfgs_free(lbfgsfloatval_t *x); void WINDOWS_EXPORT lbfgs_free(lbfgsfloatval_t *x);
/** @} */ /** @} */
......
platforms/cuda/src/CudaKernels.cpp
View file @ f2276667
...@@ -4554,7 +4554,7 @@ void CudaIntegrateCustomStepKernel::initialize(const System& system, const Custo ...@@ -4554,7 +4554,7 @@ void CudaIntegrateCustomStepKernel::initialize(const System& system, const Custo
numGlobalVariables = integrator.getNumGlobalVariables(); numGlobalVariables = integrator.getNumGlobalVariables();
int elementSize = (cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float)); int elementSize = (cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
globalValues = new CudaArray(cu, max(1, numGlobalVariables), elementSize, "globalVariables"); globalValues = new CudaArray(cu, max(1, numGlobalVariables), elementSize, "globalVariables");
sumBuffer = new CudaArray(cu, 3*system.getNumParticles(), elementSize, "sumBuffer"); sumBuffer = new CudaArray(cu, ((3*system.getNumParticles()+3)/4)*4, elementSize, "sumBuffer");
potentialEnergy = new CudaArray(cu, 1, cu.getEnergyBuffer().getElementSize(), "potentialEnergy"); potentialEnergy = new CudaArray(cu, 1, cu.getEnergyBuffer().getElementSize(), "potentialEnergy");
kineticEnergy = new CudaArray(cu, 1, elementSize, "kineticEnergy"); kineticEnergy = new CudaArray(cu, 1, elementSize, "kineticEnergy");
perDofValues = new CudaParameterSet(cu, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cu.getUseDoublePrecision() || cu.getUseMixedPrecision()); perDofValues = new CudaParameterSet(cu, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cu.getUseDoublePrecision() || cu.getUseMixedPrecision());
...@@ -5103,6 +5103,7 @@ void CudaIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrat ...@@ -5103,6 +5103,7 @@ void CudaIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrat
kernelArgs[i][0][10] = &randomIndex; kernelArgs[i][0][10] = &randomIndex;
if (requiredUniform[i] > 0) if (requiredUniform[i] > 0)
cu.executeKernel(randomKernel, &randomArgs[0], numAtoms); cu.executeKernel(randomKernel, &randomArgs[0], numAtoms);
cu.clearBuffer(*sumBuffer);
cu.executeKernel(kernels[i][0], &kernelArgs[i][0][0], numAtoms); cu.executeKernel(kernels[i][0], &kernelArgs[i][0][0], numAtoms);
cu.executeKernel(kernels[i][1], &kernelArgs[i][1][0], CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize); cu.executeKernel(kernels[i][1], &kernelArgs[i][1][0], CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);
} }
...@@ -5150,6 +5151,7 @@ double CudaIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& context, ...@@ -5150,6 +5151,7 @@ double CudaIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& context,
int randomIndex = 0; int randomIndex = 0;
kineticEnergyArgs[1] = &posCorrection; kineticEnergyArgs[1] = &posCorrection;
kineticEnergyArgs[10] = &randomIndex; kineticEnergyArgs[10] = &randomIndex;
cu.clearBuffer(*sumBuffer);
cu.executeKernel(kineticEnergyKernel, &kineticEnergyArgs[0], cu.getNumAtoms()); cu.executeKernel(kineticEnergyKernel, &kineticEnergyArgs[0], cu.getNumAtoms());
void* args[] = {&sumBuffer->getDevicePointer(), &kineticEnergy->getDevicePointer()}; void* args[] = {&sumBuffer->getDevicePointer(), &kineticEnergy->getDevicePointer()};
cu.executeKernel(sumKineticEnergyKernel, args, CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize); cu.executeKernel(sumKineticEnergyKernel, args, CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);
......
platforms/cuda/tests/TestCudaCustomIntegrator.cpp
View file @ f2276667
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. * * Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -345,20 +345,20 @@ void testSum() { ...@@ -345,20 +345,20 @@ void testSum() {
OpenMM_SFMT::SFMT sfmt; OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt); init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
system.addParticle(1.5); system.addParticle(i%10 == 0 ? 0.0 : 1.5);
nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1); nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
bool close = true; bool close = true;
while (close) { while (close) {
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt)); positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
close = false; close = false;
for (int j = 0; j < i; ++j) { for (int j = 0; j < i; ++j) {
Vec3 delta = positions[i]-positions[j]; Vec3 delta = positions[i]-positions[j];
if (delta.dot(delta) < 0.1) if (delta.dot(delta) < 1)
close = true; close = true;
} }
} }
} }
CustomIntegrator integrator(0.01); CustomIntegrator integrator(0.005);
integrator.addGlobalVariable("ke", 0); integrator.addGlobalVariable("ke", 0);
integrator.addComputePerDof("v", "v+dt*f/m"); integrator.addComputePerDof("v", "v+dt*f/m");
integrator.addComputePerDof("x", "x+dt*v"); integrator.addComputePerDof("x", "x+dt*v");
...@@ -368,10 +368,8 @@ void testSum() { ...@@ -368,10 +368,8 @@ void testSum() {
// See if the sum is being computed correctly. // See if the sum is being computed correctly.
State state = context.getState(State::Energy);
const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 100; ++i) { for (int i = 0; i < 100; ++i) {
state = context.getState(State::Energy); State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5); ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
integrator.step(1); integrator.step(1);
} }
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ f2276667
...@@ -4786,7 +4786,7 @@ void OpenCLIntegrateCustomStepKernel::initialize(const System& system, const Cus ...@@ -4786,7 +4786,7 @@ void OpenCLIntegrateCustomStepKernel::initialize(const System& system, const Cus
numGlobalVariables = integrator.getNumGlobalVariables(); numGlobalVariables = integrator.getNumGlobalVariables();
int elementSize = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision() ? sizeof(double) : sizeof(float)); int elementSize = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
globalValues = new OpenCLArray(cl, max(1, numGlobalVariables), elementSize, "globalVariables"); globalValues = new OpenCLArray(cl, max(1, numGlobalVariables), elementSize, "globalVariables");
sumBuffer = new OpenCLArray(cl, 3*system.getNumParticles(), elementSize, "sumBuffer"); sumBuffer = new OpenCLArray(cl, ((3*system.getNumParticles()+3)/4)*4, elementSize, "sumBuffer");
potentialEnergy = new OpenCLArray(cl, 1, cl.getEnergyBuffer().getElementSize(), "potentialEnergy"); potentialEnergy = new OpenCLArray(cl, 1, cl.getEnergyBuffer().getElementSize(), "potentialEnergy");
kineticEnergy = new OpenCLArray(cl, 1, elementSize, "kineticEnergy"); kineticEnergy = new OpenCLArray(cl, 1, elementSize, "kineticEnergy");
perDofValues = new OpenCLParameterSet(cl, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cl.getUseDoublePrecision() || cl.getUseMixedPrecision()); perDofValues = new OpenCLParameterSet(cl, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cl.getUseDoublePrecision() || cl.getUseMixedPrecision());
...@@ -5329,6 +5329,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr ...@@ -5329,6 +5329,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
kernels[i][0].setArg<cl_uint>(10, integration.prepareRandomNumbers(requiredGaussian[i])); kernels[i][0].setArg<cl_uint>(10, integration.prepareRandomNumbers(requiredGaussian[i]));
if (requiredUniform[i] > 0) if (requiredUniform[i] > 0)
cl.executeKernel(randomKernel, numAtoms); cl.executeKernel(randomKernel, numAtoms);
cl.clearBuffer(*sumBuffer);
cl.executeKernel(kernels[i][0], numAtoms); cl.executeKernel(kernels[i][0], numAtoms);
cl.executeKernel(kernels[i][1], OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize); cl.executeKernel(kernels[i][1], OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
} }
...@@ -5375,6 +5376,7 @@ double OpenCLIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& contex ...@@ -5375,6 +5376,7 @@ double OpenCLIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& contex
cl.executeKernel(sumPotentialEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize); cl.executeKernel(sumPotentialEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
forcesAreValid = true; forcesAreValid = true;
} }
cl.clearBuffer(*sumBuffer);
cl.executeKernel(kineticEnergyKernel, cl.getNumAtoms()); cl.executeKernel(kineticEnergyKernel, cl.getNumAtoms());
cl.executeKernel(sumKineticEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize); cl.executeKernel(sumKineticEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) { if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
......
platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp
View file @ f2276667
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2011 Stanford University and the Authors. * * Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -345,20 +345,20 @@ void testSum() { ...@@ -345,20 +345,20 @@ void testSum() {
OpenMM_SFMT::SFMT sfmt; OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt); init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
system.addParticle(1.5); system.addParticle(i%10 == 0 ? 0.0 : 1.5);
nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1); nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
bool close = true; bool close = true;
while (close) { while (close) {
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt)); positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
close = false; close = false;
for (int j = 0; j < i; ++j) { for (int j = 0; j < i; ++j) {
Vec3 delta = positions[i]-positions[j]; Vec3 delta = positions[i]-positions[j];
if (delta.dot(delta) < 0.1) if (delta.dot(delta) < 1)
close = true; close = true;
} }
} }
} }
CustomIntegrator integrator(0.01); CustomIntegrator integrator(0.005);
integrator.addGlobalVariable("ke", 0); integrator.addGlobalVariable("ke", 0);
integrator.addComputePerDof("v", "v+dt*f/m"); integrator.addComputePerDof("v", "v+dt*f/m");
integrator.addComputePerDof("x", "x+dt*v"); integrator.addComputePerDof("x", "x+dt*v");
...@@ -368,10 +368,8 @@ void testSum() { ...@@ -368,10 +368,8 @@ void testSum() {
// See if the sum is being computed correctly. // See if the sum is being computed correctly.
State state = context.getState(State::Energy);
const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 100; ++i) { for (int i = 0; i < 100; ++i) {
state = context.getState(State::Energy); State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5); ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
integrator.step(1); integrator.step(1);
} }
......
platforms/reference/tests/TestReferenceCustomIntegrator.cpp
View file @ f2276667
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2011 Stanford University and the Authors. * * Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -336,20 +336,20 @@ void testSum() { ...@@ -336,20 +336,20 @@ void testSum() {
OpenMM_SFMT::SFMT sfmt; OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt); init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
system.addParticle(1.5); system.addParticle(i%10 == 0 ? 0.0 : 1.5);
nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1); nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
bool close = true; bool close = true;
while (close) { while (close) {
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt)); positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
close = false; close = false;
for (int j = 0; j < i; ++j) { for (int j = 0; j < i; ++j) {
Vec3 delta = positions[i]-positions[j]; Vec3 delta = positions[i]-positions[j];
if (delta.dot(delta) < 0.1) if (delta.dot(delta) < 1)
close = true; close = true;
} }
} }
} }
CustomIntegrator integrator(0.01); CustomIntegrator integrator(0.005);
integrator.addGlobalVariable("ke", 0); integrator.addGlobalVariable("ke", 0);
integrator.addComputePerDof("v", "v+dt*f/m"); integrator.addComputePerDof("v", "v+dt*f/m");
integrator.addComputePerDof("x", "x+dt*v"); integrator.addComputePerDof("x", "x+dt*v");
...@@ -359,10 +359,8 @@ void testSum() { ...@@ -359,10 +359,8 @@ void testSum() {
// See if the sum is being computed correctly. // See if the sum is being computed correctly.
State state = context.getState(State::Energy);
const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 100; ++i) { for (int i = 0; i < 100; ++i) {
state = context.getState(State::Energy); State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5); ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
integrator.step(1); integrator.step(1);
} }
......
plugins/drude/CMakeLists.txt
View file @ f2276667
...@@ -205,9 +205,9 @@ ELSE (EXECUTABLE_OUTPUT_PATH) ...@@ -205,9 +205,9 @@ ELSE (EXECUTABLE_OUTPUT_PATH)
SET (TEST_PATH .) SET (TEST_PATH .)
ENDIF (EXECUTABLE_OUTPUT_PATH) ENDIF (EXECUTABLE_OUTPUT_PATH)
#IF (OPENMM_BUILD_SERIALIZATION_SUPPORT) IF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
# ADD_SUBDIRECTORY(serialization) ADD_SUBDIRECTORY(serialization)
#ENDIF (OPENMM_BUILD_SERIALIZATION_SUPPORT) ENDIF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
#INCLUDE(ApiDoxygen.cmake) #INCLUDE(ApiDoxygen.cmake)
......
plugins/drude/openmmapi/include/openmm/DrudeKernels.h
View file @ f2276667
...@@ -34,6 +34,7 @@ ...@@ -34,6 +34,7 @@
#include "openmm/DrudeForce.h" #include "openmm/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h" #include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/DrudeSCFIntegrator.h"
#include "openmm/Platform.h" #include "openmm/Platform.h"
#include "openmm/System.h" #include "openmm/System.h"
#include "openmm/Vec3.h" #include "openmm/Vec3.h"
...@@ -108,6 +109,37 @@ public: ...@@ -108,6 +109,37 @@ public:
virtual double computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) = 0; virtual double computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) = 0;
}; };
/**
* This kernel is invoked by DrudeSCFIntegrator to take one time step.
*/
class IntegrateDrudeSCFStepKernel : public KernelImpl {
public:
static std::string Name() {
return "IntegrateDrudeSCFStep";
}
IntegrateDrudeSCFStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the DrudeSCFIntegrator this kernel will be used for
* @param force the DrudeForce to get particle parameters from
*/
virtual void initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) = 0;
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the DrudeSCFIntegrator this kernel is being used for
*/
virtual void execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) = 0;
/**
* Compute the kinetic energy.
*/
virtual double computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) = 0;
};
} // namespace OpenMM } // namespace OpenMM
#endif /*DRUDE_KERNELS_H_*/ #endif /*DRUDE_KERNELS_H_*/
plugins/drude/openmmapi/include/openmm/DrudeSCFIntegrator.h 0 → 100644
View file @ f2276667
#ifndef OPENMM_DRUDESCFINTEGRATOR_H_
#define OPENMM_DRUDESCFINTEGRATOR_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) 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. *
* -------------------------------------------------------------------------- */
#include "openmm/Integrator.h"
#include "openmm/Kernel.h"
#include "openmm/internal/windowsExportDrude.h"
namespace OpenMM {
/**
* This is a leap-frog Verlet Integrator that simulates systems with Drude particles. It uses the
* self-consistent field (SCF) method: at every time step, the positions of Drude particles are
* adjusted to minimize the potential energy.
*
* This Integrator requires the System to include a DrudeForce, which it uses to identify the Drude
* particles.
*/
class OPENMM_EXPORT_DRUDE DrudeSCFIntegrator : public Integrator {
public:
/**
* Create a DrudeSCFIntegrator.
*
* @param stepSize the step size with which to integrator the system (in picoseconds)
*/
DrudeSCFIntegrator(double stepSize);
/**
* Get the error tolerance to use when minimizing the potential energy. This roughly corresponds
* to the maximum allowed force magnitude on the Drude particles after minimization.
*
* @return the error tolerance to use, measured in kJ/mol/nm
*/
double getMinimizationErrorTolerance() const {
return tolerance;
}
/**
* Set the error tolerance to use when minimizing the potential energy. This roughly corresponds
* to the maximum allowed force magnitude on the Drude particles after minimization.
*
* @param tol the error tolerance to use, measured in kJ/mol/nm
*/
void setMinimizationErrorTolerance(double tol) {
tolerance = tol;
}
/**
* Advance a simulation through time by taking a series of time steps.
*
* @param steps the number of time steps to take
*/
void step(int steps);
protected:
/**
* This will be called by the Context when it is created. It informs the Integrator
* of what context it will be integrating, and gives it a chance to do any necessary initialization.
* It will also get called again if the application calls reinitialize() on the Context.
*/
void initialize(ContextImpl& context);
/**
* This will be called by the Context when it is destroyed to let the Integrator do any necessary
* cleanup. It will also get called again if the application calls reinitialize() on the Context.
*/
void cleanup();
/**
* When the user modifies the state, we need to mark that the forces need to be recalculated.
*/
void stateChanged(State::DataType changed);
/**
* Get the names of all Kernels used by this Integrator.
*/
std::vector<std::string> getKernelNames();
/**
* Compute the kinetic energy of the system at the current time.
*/
double computeKineticEnergy();
private:
double tolerance;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_DRUDESCFINTEGRATOR_H_*/
plugins/drude/openmmapi/src/DrudeSCFIntegrator.cpp 0 → 100644
View file @ f2276667
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
#include "openmm/DrudeSCFIntegrator.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/DrudeKernels.h"
#include <cmath>
#include <ctime>
#include <string>
using namespace OpenMM;
using std::string;
using std::vector;
DrudeSCFIntegrator::DrudeSCFIntegrator(double stepSize) {
setStepSize(stepSize);
setMinimizationErrorTolerance(0.1);
setConstraintTolerance(1e-5);
}
void DrudeSCFIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
const DrudeForce* force = NULL;
const System& system = contextRef.getSystem();
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
if (force == NULL)
force = dynamic_cast<const DrudeForce*>(&system.getForce(i));
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (force == NULL)
throw OpenMMException("The System does not contain a DrudeForce");
context = &contextRef;
owner = &contextRef.getOwner();
kernel = context->getPlatform().createKernel(IntegrateDrudeSCFStepKernel::Name(), contextRef);
kernel.getAs<IntegrateDrudeSCFStepKernel>().initialize(contextRef.getSystem(), *this, *force);
}
void DrudeSCFIntegrator::cleanup() {
kernel = Kernel();
}
void DrudeSCFIntegrator::stateChanged(State::DataType changed) {
}
vector<string> DrudeSCFIntegrator::getKernelNames() {
std::vector<std::string> names;
names.push_back(IntegrateDrudeSCFStepKernel::Name());
return names;
}
double DrudeSCFIntegrator::computeKineticEnergy() {
return kernel.getAs<IntegrateDrudeSCFStepKernel>().computeKineticEnergy(*context, *this);
}
void DrudeSCFIntegrator::step(int steps) {
for (int i = 0; i < steps; ++i) {
context->updateContextState();
context->calcForcesAndEnergy(true, false);
kernel.getAs<IntegrateDrudeSCFStepKernel>().execute(*context, *this);
}
}
plugins/drude/platforms/cuda/src/CudaDrudeKernelFactory.cpp
View file @ f2276667
/* -------------------------------------------------------------------------- * /* -------------------------------------------------------------------------- *
* OpenMMAmoeba * * OpenMMDrude *
* -------------------------------------------------------------------------- * * -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from * * This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of * * Simbios, the NIH National Center for Physics-Based Simulation of *
......
plugins/drude/platforms/opencl/src/OpenCLDrudeKernelFactory.cpp
View file @ f2276667
/* -------------------------------------------------------------------------- * /* -------------------------------------------------------------------------- *
* OpenMMAmoeba * * OpenMMDrude *
* -------------------------------------------------------------------------- * * -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from * * This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of * * Simbios, the NIH National Center for Physics-Based Simulation of *
......
plugins/drude/platforms/opencl/src/kernels/drudeLangevin.cl 100644 → 100755
View file @ f2276667
...@@ -13,10 +13,12 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob ...@@ -13,10 +13,12 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob
int index = normalParticles[i]; int index = normalParticles[i];
mixed4 velocity = velm[index]; mixed4 velocity = velm[index];
if (velocity.w != 0) { if (velocity.w != 0) {
mixed sqrtInvMass = SQRT(velocity.w); mixed sqrtInvMass = sqrt(velocity.w);
float4 rand = random[randomIndex+index]; float4 rand = random[randomIndex+index];
real4 f = force[index]; real4 f = force[index];
velocity.xyz = vscale*velocity.xyz + fscale*velocity.w*f.xyz + noisescale*sqrtInvMass*rand.xyz; velocity.x = vscale*velocity.x + fscale*velocity.w*f.x + noisescale*sqrtInvMass*rand.x;
velocity.y = vscale*velocity.y + fscale*velocity.w*f.y + noisescale*sqrtInvMass*rand.y;
velocity.z = vscale*velocity.z + fscale*velocity.w*f.z + noisescale*sqrtInvMass*rand.z;
velm[index] = velocity; velm[index] = velocity;
posDelta[index] = (mixed4) (stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0); posDelta[index] = (mixed4) (stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0);
} }
...@@ -29,24 +31,28 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob ...@@ -29,24 +31,28 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob
int2 particles = pairParticles[i]; int2 particles = pairParticles[i];
mixed4 velocity1 = velm[particles.x]; mixed4 velocity1 = velm[particles.x];
mixed4 velocity2 = velm[particles.y]; mixed4 velocity2 = velm[particles.y];
mixed mass1 = RECIP(velocity1.w); mixed mass1 = 1/velocity1.w;
mixed mass2 = RECIP(velocity2.w); mixed mass2 = 1/velocity2.w;
mixed invTotalMass = RECIP(mass1+mass2); mixed invTotalMass = 1/(mass1+mass2);
mixed invReducedMass = (mass1+mass2)*velocity1.w*velocity2.w; mixed invReducedMass = (mass1+mass2)*velocity1.w*velocity2.w;
mixed mass1fract = invTotalMass*mass1; mixed mass1fract = invTotalMass*mass1;
mixed mass2fract = invTotalMass*mass2; mixed mass2fract = invTotalMass*mass2;
mixed sqrtInvTotalMass = SQRT(invTotalMass); mixed sqrtInvTotalMass = sqrt(invTotalMass);
mixed sqrtInvReducedMass = SQRT(invReducedMass); mixed sqrtInvReducedMass = sqrt(invReducedMass);
mixed4 cmVel = velocity1*mass1fract+velocity2*mass2fract; mixed4 cmVel = velocity1*mass1fract+velocity2*mass2fract;
mixed4 relVel = velocity2-velocity1; mixed4 relVel = velocity2-velocity1;
mixed4 force1 = force[particles.x]; mixed4 force1 = convert_mixed4(force[particles.x]);
mixed4 force2 = force[particles.y]; mixed4 force2 = convert_mixed4(force[particles.y]);
mixed4 cmForce = force1+force2; mixed4 cmForce = force1+force2;
mixed4 relForce = force2*mass1fract - force1*mass2fract; mixed4 relForce = force2*mass1fract - force1*mass2fract;
float4 rand1 = random[randomIndex+2*i]; float4 rand1 = random[randomIndex+2*i];
float4 rand2 = random[randomIndex+2*i+1]; float4 rand2 = random[randomIndex+2*i+1];
cmVel.xyz = vscale*cmVel.xyz + fscale*invTotalMass*cmForce.xyz + noisescale*sqrtInvTotalMass*rand1.xyz; cmVel.x = vscale*cmVel.x + fscale*invTotalMass*cmForce.x + noisescale*sqrtInvTotalMass*rand1.x;
relVel.xyz = vscaleDrude*relVel.xyz + fscaleDrude*invReducedMass*relForce.xyz + noisescaleDrude*sqrtInvReducedMass*rand2.xyz; cmVel.y = vscale*cmVel.y + fscale*invTotalMass*cmForce.y + noisescale*sqrtInvTotalMass*rand1.y;
cmVel.z = vscale*cmVel.z + fscale*invTotalMass*cmForce.z + noisescale*sqrtInvTotalMass*rand1.z;
relVel.x = vscaleDrude*relVel.x + fscaleDrude*invReducedMass*relForce.x + noisescaleDrude*sqrtInvReducedMass*rand2.x;
relVel.y = vscaleDrude*relVel.y + fscaleDrude*invReducedMass*relForce.y + noisescaleDrude*sqrtInvReducedMass*rand2.y;
relVel.z = vscaleDrude*relVel.z + fscaleDrude*invReducedMass*relForce.z + noisescaleDrude*sqrtInvReducedMass*rand2.z;
velocity1.xyz = cmVel.xyz-relVel.xyz*mass2fract; velocity1.xyz = cmVel.xyz-relVel.xyz*mass2fract;
velocity2.xyz = cmVel.xyz+relVel.xyz*mass1fract; velocity2.xyz = cmVel.xyz+relVel.xyz*mass1fract;
velm[particles.x] = velocity1; velm[particles.x] = velocity1;
......
plugins/drude/platforms/reference/src/ReferenceDrudeKernelFactory.cpp
View file @ f2276667
/* -------------------------------------------------------------------------- * /* -------------------------------------------------------------------------- *
* OpenMMAmoeba * * OpenMMDrude *
* -------------------------------------------------------------------------- * * -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from * * This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of * * Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2011 Stanford University and the Authors. * * Portions copyright (c) 2011-2013 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -49,6 +49,7 @@ extern "C" void initDrudeReferenceKernels() { ...@@ -49,6 +49,7 @@ extern "C" void initDrudeReferenceKernels() {
ReferenceDrudeKernelFactory* factory = new ReferenceDrudeKernelFactory(); ReferenceDrudeKernelFactory* factory = new ReferenceDrudeKernelFactory();
platform.registerKernelFactory(CalcDrudeForceKernel::Name(), factory); platform.registerKernelFactory(CalcDrudeForceKernel::Name(), factory);
platform.registerKernelFactory(IntegrateDrudeLangevinStepKernel::Name(), factory); platform.registerKernelFactory(IntegrateDrudeLangevinStepKernel::Name(), factory);
platform.registerKernelFactory(IntegrateDrudeSCFStepKernel::Name(), factory);
} }
KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const { KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
...@@ -57,5 +58,7 @@ KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, cons ...@@ -57,5 +58,7 @@ KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, cons
return new ReferenceCalcDrudeForceKernel(name, platform); return new ReferenceCalcDrudeForceKernel(name, platform);
if (name == IntegrateDrudeLangevinStepKernel::Name()) if (name == IntegrateDrudeLangevinStepKernel::Name())
return new ReferenceIntegrateDrudeLangevinStepKernel(name, platform, data); return new ReferenceIntegrateDrudeLangevinStepKernel(name, platform, data);
if (name == IntegrateDrudeSCFStepKernel::Name())
return new ReferenceIntegrateDrudeSCFStepKernel(name, platform, data);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str()); throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
} }
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @ f2276667
...@@ -70,6 +70,40 @@ static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorith ...@@ -70,6 +70,40 @@ static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorith
} }
} }
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift, ReferenceConstraintAlgorithm* constraints) {
const System& system = context.getSystem();
int numParticles = system.getNumParticles();
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context);
// Compute the shifted velocities.
vector<RealVec> shiftedVel(numParticles);
for (int i = 0; i < numParticles; ++i) {
if (inverseMasses[i] > 0)
shiftedVel[i] = velData[i]+forceData[i]*(timeShift*inverseMasses[i]);
else
shiftedVel[i] = velData[i];
}
// Apply constraints to them.
if (constraints != NULL) {
constraints->setTolerance(1e-4);
constraints->applyToVelocities(numParticles, posData, shiftedVel, inverseMasses);
}
// Compute the kinetic energy.
double energy = 0.0;
for (int i = 0; i < numParticles; ++i)
if (inverseMasses[i] > 0)
energy += (shiftedVel[i].dot(shiftedVel[i]))/inverseMasses[i];
return 0.5*energy;
}
void ReferenceCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) { void ReferenceCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) {
// Initialize particle parameters. // Initialize particle parameters.
...@@ -330,35 +364,156 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co ...@@ -330,35 +364,156 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
} }
double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) { double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
const System& system = context.getSystem(); return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize(), constraints);
int numParticles = system.getNumParticles(); }
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context); ReferenceIntegrateDrudeSCFStepKernel::~ReferenceIntegrateDrudeSCFStepKernel() {
vector<RealVec>& forceData = extractForces(context); if (constraints != NULL)
delete constraints;
if (minimizerPos != NULL)
lbfgs_free(minimizerPos);
}
void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) {
// Identify Drude particles.
// Compute the shifted velocities. for (int i = 0; i < force.getNumParticles(); i++) {
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
drudeParticles.push_back(p);
}
// Record particle masses.
vector<RealOpenMM> particleMass;
for (int i = 0; i < system.getNumParticles(); i++) {
double mass = system.getParticleMass(i);
particleMass.push_back(mass);
particleInvMass.push_back(mass == 0.0 ? 0.0 : 1.0/mass);
}
vector<RealVec> shiftedVel(numParticles); // Prepare constraints.
double timeShift = 0.5*integrator.getStepSize();
for (int i = 0; i < numParticles; ++i) { int numConstraints = system.getNumConstraints();
if (particleInvMass[i] > 0) if (numConstraints > 0) {
shiftedVel[i] = velData[i]+forceData[i]*(timeShift*particleInvMass[i]); vector<pair<int, int> > constraintIndices(numConstraints);
else vector<RealOpenMM> constraintDistances(numConstraints);
shiftedVel[i] = velData[i]; for (int i = 0; i < numConstraints; ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
constraintIndices[i].first = particle1;
constraintIndices[i].second = particle2;
constraintDistances[i] = static_cast<RealOpenMM>(distance);
}
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, particleMass, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
// Apply constraints to them. // Initialize the energy minimizer.
if (constraints != NULL) { minimizerPos = lbfgs_malloc(drudeParticles.size()*3);
constraints->setTolerance(1e-4); if (minimizerPos == NULL)
constraints->applyToVelocities(numParticles, posData, shiftedVel, particleInvMass); throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory");
lbfgs_parameter_init(&minimizerParams);
minimizerParams.max_iterations = 0;
minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
}
void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
vector<RealVec>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context);
vector<RealVec>& force = extractForces(context);
// Update the positions and velocities.
int numParticles = particleInvMass.size();
vector<RealVec> xPrime(numParticles);
RealOpenMM dt = integrator.getStepSize();
for (int i = 0; i < numParticles; i++) {
if (particleInvMass[i] != 0.0) {
vel[i] += force[i]*particleInvMass[i]*dt;
xPrime[i] = pos[i]+vel[i]*dt;
}
} }
// Apply constraints.
// Compute the kinetic energy. if (constraints != NULL)
constraints->apply(numParticles, pos, xPrime, particleInvMass);
double energy = 0.0; // Record the constrained positions and velocities.
for (int i = 0; i < numParticles; ++i)
if (particleInvMass[i] > 0) RealOpenMM dtInv = 1.0/dt;
energy += (shiftedVel[i].dot(shiftedVel[i]))/particleInvMass[i]; for (int i = 0; i < numParticles; i++) {
return 0.5*energy; if (particleInvMass[i] != 0.0) {
vel[i] = (xPrime[i]-pos[i])*dtInv;
pos[i] = xPrime[i];
}
}
// Update the positions of virtual sites and Drude particles.
ReferenceVirtualSites::computePositions(context.getSystem(), pos);
minimize(context, integrator.getMinimizationErrorTolerance());
data.time += integrator.getStepSize();
data.stepCount++;
}
double ReferenceIntegrateDrudeSCFStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize(), constraints);
}
struct MinimizerData {
ContextImpl& context;
vector<int>& drudeParticles;
MinimizerData(ContextImpl& context, vector<int>& drudeParticles) : context(context), drudeParticles(drudeParticles) {}
};
static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsfloatval_t *g, const int n, const lbfgsfloatval_t step) {
MinimizerData* data = reinterpret_cast<MinimizerData*>(instance);
ContextImpl& context = data->context;
vector<int>& drudeParticles = data->drudeParticles;
int numDrudeParticles = drudeParticles.size();
// Compute the force and energy for this configuration.
vector<RealVec>& pos = extractPositions(context);
vector<RealVec>& force = extractForces(context);
for (int i = 0; i < numDrudeParticles; i++)
pos[drudeParticles[i]] = RealVec(x[3*i], x[3*i+1], x[3*i+2]);
double energy = context.calcForcesAndEnergy(true, true);
for (int i = 0; i < numDrudeParticles; i++) {
RealVec f = force[drudeParticles[i]];
g[3*i] = -f[0];
g[3*i+1] = -f[1];
g[3*i+2] = -f[2];
}
return energy;
} }
void ReferenceIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double tolerance) {
// Record the initial positions and determine a normalization constant for scaling the tolerance.
vector<RealVec>& pos = extractPositions(context);
int numDrudeParticles = drudeParticles.size();
double norm = 0.0;
for (int i = 0; i < numDrudeParticles; i++) {
RealVec p = pos[drudeParticles[i]];
minimizerPos[3*i] = p[0];
minimizerPos[3*i+1] = p[1];
minimizerPos[3*i+2] = p[2];
norm += p.dot(p);
}
norm /= numDrudeParticles;
norm = (norm < 1 ? 1 : sqrt(norm));
minimizerParams.epsilon = tolerance/norm;
// Perform the minimization.
lbfgsfloatval_t fx;
MinimizerData data(context, drudeParticles);
lbfgs(numDrudeParticles*3, minimizerPos, &fx, evaluate, NULL, &data, &minimizerParams);
}
\ No newline at end of file
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h
View file @ f2276667
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include "ReferencePlatform.h" #include "ReferencePlatform.h"
#include "openmm/DrudeKernels.h" #include "openmm/DrudeKernels.h"
#include "SimTKUtilities/RealVec.h" #include "SimTKUtilities/RealVec.h"
#include "lbfgs.h"
#include <utility> #include <utility>
#include <vector> #include <vector>
...@@ -120,6 +121,47 @@ private: ...@@ -120,6 +121,47 @@ private:
ReferenceConstraintAlgorithm* constraints; ReferenceConstraintAlgorithm* constraints;
}; };
/**
* This kernel is invoked by DrudeSCFIntegrator to take one time step
*/
class ReferenceIntegrateDrudeSCFStepKernel : public IntegrateDrudeSCFStepKernel {
public:
ReferenceIntegrateDrudeSCFStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
IntegrateDrudeSCFStepKernel(name, platform), data(data), constraints(NULL), minimizerPos(NULL) {
}
~ReferenceIntegrateDrudeSCFStepKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the DrudeSCFIntegrator this kernel will be used for
* @param force the DrudeForce to get particle parameters from
*/
void initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the DrudeSCFIntegrator this kernel is being used for
*/
void execute(ContextImpl& context, const DrudeSCFIntegrator& integrator);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the DrudeSCFIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator);
private:
void minimize(ContextImpl& context, double tolerance);
ReferencePlatform::PlatformData& data;
std::vector<int> drudeParticles;
std::vector<double> particleInvMass;
ReferenceConstraintAlgorithm* constraints;
lbfgsfloatval_t *minimizerPos;
lbfgs_parameter_t minimizerParams;
};
} // namespace OpenMM } // namespace OpenMM
#endif /*REFERENCE_DRUDE_KERNELS_H_*/ #endif /*REFERENCE_DRUDE_KERNELS_H_*/
plugins/drude/platforms/reference/tests/TestReferenceDrudeSCFIntegrator.cpp 0 → 100644
View file @ f2276667
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of DrudeSCFIntegrator.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/VirtualSite.h"
#include "openmm/DrudeForce.h"
#include "openmm/DrudeSCFIntegrator.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testWater() {
// Create a box of SWM4-NDP water molecules. This involves constraints, virtual sites,
// and Drude particles.
const int gridSize = 4;
const int numMolecules = gridSize*gridSize*gridSize;
const double spacing = 0.6;
const double boxSize = spacing*(gridSize+1);
System system;
NonbondedForce* nonbonded = new NonbondedForce();
DrudeForce* drude = new DrudeForce();
system.addForce(nonbonded);
system.addForce(drude);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(1.0);
for (int i = 0; i < numMolecules; i++) {
int startIndex = system.getNumParticles();
system.addParticle(15.6); // O
system.addParticle(0.4); // D
system.addParticle(1.0); // H1
system.addParticle(1.0); // H2
system.addParticle(0.0); // M
nonbonded->addParticle(1.71636, 0.318395, 0.21094*4.184);
nonbonded->addParticle(-1.71636, 1, 0);
nonbonded->addParticle(0.55733, 1, 0);
nonbonded->addParticle(0.55733, 1, 0);
nonbonded->addParticle(-1.11466, 1, 0);
for (int j = 0; j < 5; j++)
for (int k = 0; k < j; k++)
nonbonded->addException(startIndex+j, startIndex+k, 0, 1, 0);
system.addConstraint(startIndex, startIndex+2, 0.09572);
system.addConstraint(startIndex, startIndex+3, 0.09572);
system.addConstraint(startIndex+2, startIndex+3, 0.15139);
system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
}
vector<Vec3> positions;
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos(i*spacing, j*spacing, k*spacing);
positions.push_back(pos);
positions.push_back(pos);
positions.push_back(pos+Vec3(0.09572, 0, 0));
positions.push_back(pos+Vec3(-0.023999, 0.092663, 0));
positions.push_back(pos);
}
// Simulate it and check energy conservation and the total force on the Drude particles.
DrudeSCFIntegrator integ(0.0005);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
context.setPositions(positions);
context.applyConstraints(1e-5);
context.setVelocitiesToTemperature(300.0);
State state = context.getState(State::Energy);
double initialEnergy = state.getPotentialEnergy()+state.getKineticEnergy();
int numSteps = 1000;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
state = context.getState(State::Energy | State::Forces);
ASSERT_EQUAL_TOL(initialEnergy, state.getPotentialEnergy()+state.getKineticEnergy(), 0.01);
const vector<Vec3>& force = state.getForces();
double norm = 0.0;
for (int j = 1; j < (int) force.size(); j += 5)
norm += sqrt(force[j].dot(force[j]));
norm = (norm/numMolecules);
ASSERT(norm < 1.0);
}
}
int main() {
try {
testWater();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
}
plugins/drude/serialization/CMakeLists.txt 0 → 100644
View file @ f2276667
#---------------------------------------------------
# OpenMMDrude Serialization Library
#
# Creates OpenMMDrude serializatin library, base name=OpenMMDrudeSerialization.
# Default libraries are shared & optimized. Variants
# are created for static (_static) and debug (_d).
#
# Windows:
# OpenMMDrudeSerialization[_d].dll
# OpenMMDrudeSerialization[_d].lib
# Unix:
# libOpenMMDrudeSerialization[_d].so
#----------------------------------------------------
# The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
SET(OPENMM_SOURCE_SUBDIRS . ../../../serialization)
SET(OPENMM_DRUDE_SOURCE_SUBDIRS . )
# Collect up information about the version of the OpenMM library we're building
# and make it available to the code so it can be built into the binaries.
SET(OPENMM_DRUDE_LIBRARY_NAME OpenMMDrude)
SET(OPENMM_SERIALIZATION_LIBRARY_NAME OpenMMSerialization)
SET(OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME OpenMMDrudeSerialization)
# Ensure that debug libraries have "_d" appended to their names.
# CMake gets this right on Windows automatically with this definition.
IF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
SET(CMAKE_DEBUG_POSTFIX "_d" CACHE INTERNAL "" FORCE)
ENDIF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
# But on Unix or Cygwin we have to add the suffix manually
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(OPENMM_DRUDE_LIBRARY_NAME ${OPENMM_DRUDE_LIBRARY_NAME}_d)
SET(OPENMM_SERIALIZATION_LIBRARY_NAME ${OPENMM_SERIALIZATION_LIBRARY_NAME}_d)
SET(OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME ${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME}_d)
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
# These are all the places to search for header files which are
# to be part of the API.
SET(API_INCLUDE_DIRS) # start empty
FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
# append
SET(API_INCLUDE_DIRS ${API_INCLUDE_DIRS}
${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include
${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include/internal)
ENDFOREACH(subdir)
# We'll need both *relative* path names, starting with their API_INCLUDE_DIRS,
# and absolute pathnames.
SET(API_REL_INCLUDE_FILES) # start these out empty
SET(API_ABS_INCLUDE_FILES)
FOREACH(dir ${API_INCLUDE_DIRS})
FILE(GLOB fullpaths ${dir}/*.h) # returns full pathnames
SET(API_ABS_INCLUDE_FILES ${API_ABS_INCLUDE_FILES} ${fullpaths})
FOREACH(pathname ${fullpaths})
GET_FILENAME_COMPONENT(filename ${pathname} NAME)
SET(API_REL_INCLUDE_FILES ${API_REL_INCLUDE_FILES} ${dir}/${filename})
ENDFOREACH(pathname)
ENDFOREACH(dir)
# collect up source files
SET(SOURCE_FILES) # empty
SET(SOURCE_INCLUDE_FILES)
FOREACH(subdir ${OPENMM_DRUDE_SOURCE_SUBDIRS})
FILE(GLOB_RECURSE src_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.c)
FILE(GLOB incl_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.h ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.hpp)
SET(SOURCE_FILES ${SOURCE_FILES} ${src_files}) #append
SET(SOURCE_INCLUDE_FILES ${SOURCE_INCLUDE_FILES} ${incl_files})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/../../../serialization/${subdir}/include)
ENDFOREACH(subdir)
#Message( "API_REL_INCLUDE_FILES=${API_REL_INCLUDE_FILES}\n" )
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src )
# Create the library
ADD_LIBRARY(${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
TARGET_LINK_LIBRARIES(${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} ${OPENMM_DRUDE_LIBRARY_NAME} ${OPENMM_SERIALIZATION_LIBRARY_NAME} ${SHARED_TARGET})
SET_TARGET_PROPERTIES(${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME} PROPERTIES COMPILE_FLAGS "-DOPENMM_DRUDE_SERIALIZATION_BUILDING_SHARED_LIBRARY -DTIXML_USE_STL -DIEEE_8087")
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${OPENMM_DRUDE_SERIALIZATION_LIBRARY_NAME})
ADD_SUBDIRECTORY(tests)
plugins/drude/serialization/include/openmm/serialization/DrudeForceProxy.h 0 → 100644
View file @ f2276667
#ifndef OPENMM_DRUDE_FORCE_PROXY_H_
#define OPENMM_DRUDE_FORCE_PROXY_H_
/* -------------------------------------------------------------------------- *
* OpenMMDrude *
* -------------------------------------------------------------------------- *
* 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 "openmm/serialization/internal/windowsExportDrudeSerialization.h"
#include "openmm/serialization/SerializationProxy.h"
namespace OpenMM {
/**
* This is a proxy for serializing DrudeForce objects.
*/
class OPENMM_EXPORT_DRUDE_SERIALIZATION DrudeForceProxy : public SerializationProxy {
public:
DrudeForceProxy();
void serialize(const void* object, SerializationNode& node) const;
void* deserialize(const SerializationNode& node) const;
};
} // namespace OpenMM
#endif /*OPENMM_DRUDE_FORCE_PROXY_H_*/
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