Skip to content

GitLab

"vscode:/vscode.git/clone" did not exist on "2c4bf882689fda969c84dfc86c4158032de10e1b"
Commit 178aa003 authored by Yutong Zhao's avatar Yutong Zhao
Browse files

Merge pull request #8 from peastman/DrudeSCFIntegrator 

Completed implementation of DrudeSCFIntegrator.  Also cleaned up other details of Drude plugin.
parents 08f1039f a5dc836e
Show whitespace changes
Inline Side-by-side
Showing with 1113 additions and 10 deletions
docs/DoxyfileC++.in
View file @ 178aa003
...@@ -78,7 +78,12 @@ WARN_LOGFILE = ...@@ -78,7 +78,12 @@ WARN_LOGFILE =
#--------------------------------------------------------------------------- #---------------------------------------------------------------------------
# configuration options related to the input files # configuration options related to the input files
#--------------------------------------------------------------------------- #---------------------------------------------------------------------------
INPUT = "@CMAKE_SOURCE_DIR@/openmmapi" "@CMAKE_SOURCE_DIR@/olla" "@CMAKE_SOURCE_DIR@/serialization/include/openmm/serialization/XmlSerializer.h" INPUT = "@CMAKE_SOURCE_DIR@/openmmapi" \
"@CMAKE_SOURCE_DIR@/olla" \
"@CMAKE_SOURCE_DIR@/serialization/include/openmm/serialization/XmlSerializer.h" \
"@CMAKE_SOURCE_DIR@/plugins/drude/openmmapi/include" \
"@CMAKE_SOURCE_DIR@/plugins/rpmd/openmmapi/include" \
"@CMAKE_SOURCE_DIR@/plugins/amoeba/openmmapi/include"
INPUT_ENCODING = UTF-8 INPUT_ENCODING = UTF-8
FILE_PATTERNS = FILE_PATTERNS =
RECURSIVE = YES RECURSIVE = YES
...@@ -87,7 +92,10 @@ EXCLUDE_SYMLINKS = NO ...@@ -87,7 +92,10 @@ EXCLUDE_SYMLINKS = NO
EXCLUDE_PATTERNS = */tests/* \ EXCLUDE_PATTERNS = */tests/* \
*/openmmapi/src/* \ */openmmapi/src/* \
*/.svn/* \ */.svn/* \
*/olla/include/openmm/kernels.h */olla/include/openmm/kernels.h \
*/DrudeKernels.h \
*/RpmdKernels.h \
*/amoebaKernels.h \
EXCLUDE_SYMBOLS = EXCLUDE_SYMBOLS =
EXAMPLE_PATH = EXAMPLE_PATH =
EXAMPLE_PATTERNS = EXAMPLE_PATTERNS =
......
plugins/drude/openmmapi/include/OpenMMDrude.h 0 → 100644
View file @ 178aa003
#ifndef OPENMM_DRUDE_H_
#define OPENMM_DRUDE_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/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/DrudeSCFIntegrator.h"
#endif /*OPENMM_DRUDE_H_*/
plugins/drude/openmmapi/include/openmm/DrudeForce.h
View file @ 178aa003
...@@ -125,7 +125,7 @@ public: ...@@ -125,7 +125,7 @@ public:
*/ */
void setParticleParameters(int index, int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34); void setParticleParameters(int index, int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34);
/** /**
* Add an interaction to the list of screenedPairs. * Add an interaction to the list of screened pairs.
* *
* @param particle1 the index within this Force of the first particle involved in the interaction * @param particle1 the index within this Force of the first particle involved in the interaction
* @param particle2 the index within this Force of the second particle involved in the interaction * @param particle2 the index within this Force of the second particle involved in the interaction
...@@ -152,10 +152,10 @@ public: ...@@ -152,10 +152,10 @@ public:
*/ */
void setScreenedPairParameters(int index, int particle1, int particle2, double thole); void setScreenedPairParameters(int index, int particle1, int particle2, double thole);
/** /**
* Update the particle and screenedPair parameters in a Context to match those stored in this Force object. This method * Update the particle and screened pair parameters in a Context to match those stored in this Force object. This method
* provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it. * provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
* Simply call setParticleParameters() and setScreenedPairParameters() to modify this object's parameters, then call * Simply call setParticleParameters() and setScreenedPairParameters() to modify this object's parameters, then call
* updateParametersInState() to copy them over to the Context. * updateParametersInContext() to copy them over to the Context.
* *
* This method has several limitations. It can be used to modify the numeric parameters associated with a particle or * This method has several limitations. It can be used to modify the numeric parameters associated with a particle or
* screened pair (polarizability, thole, etc.), but not the identities of the particles they involve. It also cannot * screened pair (polarizability, thole, etc.), but not the identities of the particles they involve. It also cannot
......
plugins/drude/platforms/cuda/src/CudaDrudeKernelFactory.cpp
View file @ 178aa003
...@@ -43,6 +43,7 @@ extern "C" void registerKernelFactories() { ...@@ -43,6 +43,7 @@ extern "C" void registerKernelFactories() {
CudaDrudeKernelFactory* factory = new CudaDrudeKernelFactory(); CudaDrudeKernelFactory* factory = new CudaDrudeKernelFactory();
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);
} }
catch (std::exception ex) { catch (std::exception ex) {
// Ignore // Ignore
...@@ -65,5 +66,7 @@ KernelImpl* CudaDrudeKernelFactory::createKernelImpl(std::string name, const Pla ...@@ -65,5 +66,7 @@ KernelImpl* CudaDrudeKernelFactory::createKernelImpl(std::string name, const Pla
return new CudaCalcDrudeForceKernel(name, platform, cu); return new CudaCalcDrudeForceKernel(name, platform, cu);
if (name == IntegrateDrudeLangevinStepKernel::Name()) if (name == IntegrateDrudeLangevinStepKernel::Name())
return new CudaIntegrateDrudeLangevinStepKernel(name, platform, cu); return new CudaIntegrateDrudeLangevinStepKernel(name, platform, cu);
if (name == IntegrateDrudeSCFStepKernel::Name())
return new CudaIntegrateDrudeSCFStepKernel(name, platform, cu);
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/cuda/src/CudaDrudeKernels.cpp
View file @ 178aa003
...@@ -179,7 +179,59 @@ double CudaCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForce ...@@ -179,7 +179,59 @@ double CudaCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForce
} }
void CudaCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) { void CudaCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) {
int numContexts = cu.getPlatformData().contexts.size();
// Set the particle parameters.
int startParticleIndex = cu.getContextIndex()*force.getNumParticles()/numContexts;
int endParticleIndex = (cu.getContextIndex()+1)*force.getNumParticles()/numContexts;
int numParticles = endParticleIndex-startParticleIndex;
if (numParticles > 0) {
if (particleParams == NULL || numParticles != particleParams->getSize())
throw OpenMMException("updateParametersInContext: The number of Drude particles has changed");
vector<float4> paramVector(numParticles);
for (int i = 0; i < numParticles; i++) {
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(startParticleIndex+i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
double a1 = (p2 == -1 ? 1 : aniso12);
double a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34);
double a3 = 3-a1-a2;
double k3 = charge*charge/(polarizability*a3);
double k1 = charge*charge/(polarizability*a1) - k3;
double k2 = charge*charge/(polarizability*a2) - k3;
if (p2 == -1)
k1 = 0;
if (p3 == -1 || p4 == -1)
k2 = 0;
paramVector[i] = make_float4((float) k1, (float) k2, (float) k3, 0.0f);
}
particleParams->upload(paramVector);
}
// Set the pair parameters.
int startPairIndex = cu.getContextIndex()*force.getNumScreenedPairs()/numContexts;
int endPairIndex = (cu.getContextIndex()+1)*force.getNumScreenedPairs()/numContexts;
int numPairs = endPairIndex-startPairIndex;
if (numPairs > 0) {
if (pairParams == NULL || numPairs != pairParams->getSize())
throw OpenMMException("updateParametersInContext: The number of screened pairs has changed");
vector<float2> paramVector(numPairs);
for (int i = 0; i < numPairs; i++) {
int drude1, drude2;
double thole;
force.getScreenedPairParameters(startPairIndex+i, drude1, drude2, thole);
int p, p1, p2, p3, p4;
double charge1, charge2, polarizability1, polarizability2, aniso12, aniso34;
force.getParticleParameters(drude1, p, p1, p2, p3, p4, charge1, polarizability1, aniso12, aniso34);
force.getParticleParameters(drude2, p, p1, p2, p3, p4, charge2, polarizability2, aniso12, aniso34);
double screeningScale = thole/pow(polarizability1*polarizability2, 1.0/6.0);
double energyScale = ONE_4PI_EPS0*charge1*charge2;
paramVector[i] = make_float2((float) screeningScale, (float) energyScale);
}
pairParams->upload(paramVector);
}
} }
CudaIntegrateDrudeLangevinStepKernel::~CudaIntegrateDrudeLangevinStepKernel() { CudaIntegrateDrudeLangevinStepKernel::~CudaIntegrateDrudeLangevinStepKernel() {
...@@ -311,3 +363,214 @@ void CudaIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const D ...@@ -311,3 +363,214 @@ void CudaIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const D
double CudaIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) { double CudaIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
return cu.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize()); return cu.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
} }
class CudaIntegrateDrudeSCFStepKernel::ReorderListener : public CudaContext::ReorderListener {
public:
ReorderListener(CudaContext& cu, const vector<int>& drudeParticles, vector<int>& reorderedDrudeParticles) :
cu(cu), drudeParticles(drudeParticles), reorderedDrudeParticles(reorderedDrudeParticles) {
}
void execute() {
const vector<int>& order = cu.getAtomIndex();
int numParticles = order.size();
vector<int> inverseOrder(numParticles);
for (int i = 0; i < numParticles; i++)
inverseOrder[order[i]] = i;
int numDrudeParticles = drudeParticles.size();
for (int i = 0; i < numDrudeParticles; i++)
reorderedDrudeParticles[i] = inverseOrder[drudeParticles[i]];
}
private:
CudaContext& cu;
const vector<int>& drudeParticles;
vector<int>& reorderedDrudeParticles;
};
CudaIntegrateDrudeSCFStepKernel::~CudaIntegrateDrudeSCFStepKernel() {
if (minimizerPos != NULL)
lbfgs_free(minimizerPos);
}
void CudaIntegrateDrudeSCFStepKernel::initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) {
cu.getPlatformData().initializeContexts(system);
cu.setAsCurrent();
// Identify Drude particles.
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);
reorderedDrudeParticles.push_back(p);
}
cu.addReorderListener(new ReorderListener(cu, drudeParticles, reorderedDrudeParticles));
// Initialize the energy minimizer.
minimizerPos = lbfgs_malloc(drudeParticles.size()*3);
if (minimizerPos == NULL)
throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory");
lbfgs_parameter_init(&minimizerParams);
minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
// Create the kernels.
map<string, string> defines;
defines["NUM_ATOMS"] = cu.intToString(cu.getNumAtoms());
defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
CUmodule module = cu.createModule(CudaKernelSources::verlet, defines, "");
kernel1 = cu.getKernel(module, "integrateVerletPart1");
kernel2 = cu.getKernel(module, "integrateVerletPart2");
prevStepSize = -1.0;
}
void CudaIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
cu.setAsCurrent();
CudaIntegrationUtilities& integration = cu.getIntegrationUtilities();
int numAtoms = cu.getNumAtoms();
double dt = integrator.getStepSize();
if (dt != prevStepSize) {
if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
vector<double2> stepSizeVec(1);
stepSizeVec[0] = make_double2(dt, dt);
cu.getIntegrationUtilities().getStepSize().upload(stepSizeVec);
}
else {
vector<float2> stepSizeVec(1);
stepSizeVec[0] = make_float2((float) dt, (float) dt);
cu.getIntegrationUtilities().getStepSize().upload(stepSizeVec);
}
prevStepSize = dt;
}
// Call the first integration kernel.
CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
void* args1[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
&cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
cu.executeKernel(kernel1, args1, numAtoms);
// Apply constraints.
integration.applyConstraints(integrator.getConstraintTolerance());
// Call the second integration kernel.
void* args2[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
&cu.getVelm().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
cu.executeKernel(kernel2, args2, numAtoms);
// Update the positions of virtual sites and Drude particles.
integration.computeVirtualSites();
minimize(context, integrator.getMinimizationErrorTolerance());
// Update the time and step count.
cu.setTime(cu.getTime()+dt);
cu.setStepCount(cu.getStepCount()+1);
}
double CudaIntegrateDrudeSCFStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
return cu.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
}
struct MinimizerData {
ContextImpl& context;
CudaContext& cu;
vector<int>& reorderedDrudeParticles;
MinimizerData(ContextImpl& context, CudaContext& cu, vector<int>& reorderedDrudeParticles) : context(context), cu(cu), reorderedDrudeParticles(reorderedDrudeParticles) {}
};
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;
CudaContext& cu = data->cu;
vector<int>& reorderedDrudeParticles = data->reorderedDrudeParticles;
int numDrudeParticles = reorderedDrudeParticles.size();
// Set the particle positions.
cu.getPosq().download(cu.getPinnedBuffer());
double4 periodicBoxSize = cu.getPeriodicBoxSize();
if (cu.getUseDoublePrecision()) {
double4* posq = (double4*) cu.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
double4& p = posq[reorderedDrudeParticles[i]];
int4 offset = cu.getPosCellOffsets()[reorderedDrudeParticles[i]];
p.x = x[3*i]+offset.x*periodicBoxSize.x;
p.y = x[3*i+1]+offset.y*periodicBoxSize.y;
p.z = x[3*i+2]+offset.z*periodicBoxSize.z;
}
}
else {
float4* posq = (float4*) cu.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
float4& p = posq[reorderedDrudeParticles[i]];
int4 offset = cu.getPosCellOffsets()[reorderedDrudeParticles[i]];
p.x = x[3*i]+offset.x*periodicBoxSize.x;
p.y = x[3*i+1]+offset.y*periodicBoxSize.y;
p.z = x[3*i+2]+offset.z*periodicBoxSize.z;
}
}
cu.getPosq().upload(cu.getPinnedBuffer());
// Compute the forces and energy for this configuration.
double energy = context.calcForcesAndEnergy(true, true);
long long* force = (long long*) cu.getPinnedBuffer();
cu.getForce().download(force);
double forceScale = -1.0/0x100000000;
int paddedNumAtoms = cu.getPaddedNumAtoms();
for (int i = 0; i < numDrudeParticles; ++i) {
int index = reorderedDrudeParticles[i];
g[3*i] = forceScale*force[index];
g[3*i+1] = forceScale*force[index+paddedNumAtoms];
g[3*i+2] = forceScale*force[index+paddedNumAtoms*2];
}
return energy;
}
void CudaIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double tolerance) {
// Record the initial positions.
int numDrudeParticles = reorderedDrudeParticles.size();
cu.getPosq().download(cu.getPinnedBuffer());
double4 periodicBoxSize = cu.getPeriodicBoxSize();
if (cu.getUseDoublePrecision()) {
double4* posq = (double4*) cu.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
double4 p = posq[reorderedDrudeParticles[i]];
int4 offset = cu.getPosCellOffsets()[reorderedDrudeParticles[i]];
minimizerPos[3*i] = p.x-offset.x*periodicBoxSize.x;
minimizerPos[3*i+1] = p.y-offset.y*periodicBoxSize.y;
minimizerPos[3*i+2] = p.z-offset.z*periodicBoxSize.z;
}
}
else {
float4* posq = (float4*) cu.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
float4 p = posq[reorderedDrudeParticles[i]];
int4 offset = cu.getPosCellOffsets()[reorderedDrudeParticles[i]];
minimizerPos[3*i] = p.x-offset.x*periodicBoxSize.x;
minimizerPos[3*i+1] = p.y-offset.y*periodicBoxSize.y;
minimizerPos[3*i+2] = p.z-offset.z*periodicBoxSize.z;
}
minimizerParams.xtol = 1e-7;
}
// Determine a normalization constant for scaling the tolerance.
double norm = 0.0;
for (int i = 0; i < 3*numDrudeParticles; i++)
norm += minimizerPos[i]*minimizerPos[i];
norm /= numDrudeParticles;
norm = (norm < 1 ? 1 : sqrt(norm));
minimizerParams.epsilon = tolerance/norm;
// Perform the minimization.
lbfgsfloatval_t fx;
MinimizerData data(context, cu, reorderedDrudeParticles);
lbfgs(numDrudeParticles*3, minimizerPos, &fx, evaluate, NULL, &data, &minimizerParams);
}
\ No newline at end of file
plugins/drude/platforms/cuda/src/CudaDrudeKernels.h
View file @ 178aa003
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include "openmm/DrudeKernels.h" #include "openmm/DrudeKernels.h"
#include "CudaContext.h" #include "CudaContext.h"
#include "CudaArray.h" #include "CudaArray.h"
#include "lbfgs.h"
namespace OpenMM { namespace OpenMM {
...@@ -115,6 +116,49 @@ private: ...@@ -115,6 +116,49 @@ private:
CUfunction kernel1, kernel2; CUfunction kernel1, kernel2;
}; };
/**
* This kernel is invoked by DrudeSCFIntegrator to take one time step
*/
class CudaIntegrateDrudeSCFStepKernel : public IntegrateDrudeSCFStepKernel {
public:
CudaIntegrateDrudeSCFStepKernel(std::string name, const Platform& platform, CudaContext& cu) :
IntegrateDrudeSCFStepKernel(name, platform), cu(cu), minimizerPos(NULL) {
}
~CudaIntegrateDrudeSCFStepKernel();
/**
* 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:
class ReorderListener;
void minimize(ContextImpl& context, double tolerance);
CudaContext& cu;
double prevStepSize;
std::vector<int> drudeParticles;
std::vector<int> reorderedDrudeParticles;
lbfgsfloatval_t *minimizerPos;
lbfgs_parameter_t minimizerParams;
CUfunction kernel1, kernel2;
};
} // namespace OpenMM } // namespace OpenMM
#endif /*CUDA_DRUDE_KERNELS_H_*/ #endif /*CUDA_DRUDE_KERNELS_H_*/
plugins/drude/platforms/cuda/tests/TestCudaDrudeForce.cpp
View file @ 178aa003
...@@ -156,6 +156,43 @@ void testThole() { ...@@ -156,6 +156,43 @@ void testThole() {
validateForce(system, positions, energySpring1+energySpring2+energyDipole); validateForce(system, positions, energySpring1+energySpring2+energyDipole);
} }
void testChangingParameters() {
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
Platform& platform = Platform::getPlatformByName("CUDA");
// Create the system.
System system;
system.addParticle(1.0);
system.addParticle(1.0);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
system.addForce(drude);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(2, 0, 0);
// Check the energy.
VerletIntegrator integ(1.0);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*k*3*3, state.getPotentialEnergy(), 1e-5);
// Modify the parameters.
const double k2 = 2.2;
const double charge2 = 0.3;
const double alpha2 = charge2*charge2/k2;
drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
drude->updateParametersInContext(context);
state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*k2*3*3, state.getPotentialEnergy(), 1e-5);
}
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
registerDrudeCudaKernelFactories(); registerDrudeCudaKernelFactories();
...@@ -164,6 +201,7 @@ int main(int argc, char* argv[]) { ...@@ -164,6 +201,7 @@ int main(int argc, char* argv[]) {
testSingleParticle(); testSingleParticle();
testAnisotropicParticle(); testAnisotropicParticle();
testThole(); testThole();
testChangingParameters();
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/drude/platforms/cuda/tests/TestCudaDrudeSCFIntegrator.cpp 0 → 100644
View file @ 178aa003
/* -------------------------------------------------------------------------- *
* 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 CUDA 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;
extern "C" OPENMM_EXPORT void registerDrudeCudaKernelFactories();
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("CUDA");
Context context(system, integ, platform);
context.setPositions(positions);
context.applyConstraints(1e-5);
context.setVelocitiesToTemperature(300.0);
State state = context.getState(State::Energy);
double initialEnergy;
int numSteps = 1000;
double maxNorm = (platform.getPropertyValue(context, "CudaPrecision") == "double" ? 1.0 : 5.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
state = context.getState(State::Energy | State::Forces);
if (i == 0)
initialEnergy = state.getPotentialEnergy()+state.getKineticEnergy();
else
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 < maxNorm);
}
}
int main(int argc, char* argv[]) {
try {
registerDrudeCudaKernelFactories();
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", string(argv[1]));
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/platforms/opencl/src/OpenCLDrudeKernelFactory.cpp
View file @ 178aa003
...@@ -43,6 +43,7 @@ extern "C" void registerKernelFactories() { ...@@ -43,6 +43,7 @@ extern "C" void registerKernelFactories() {
OpenCLDrudeKernelFactory* factory = new OpenCLDrudeKernelFactory(); OpenCLDrudeKernelFactory* factory = new OpenCLDrudeKernelFactory();
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);
} }
catch (std::exception ex) { catch (std::exception ex) {
// Ignore // Ignore
...@@ -65,5 +66,7 @@ KernelImpl* OpenCLDrudeKernelFactory::createKernelImpl(std::string name, const P ...@@ -65,5 +66,7 @@ KernelImpl* OpenCLDrudeKernelFactory::createKernelImpl(std::string name, const P
return new OpenCLCalcDrudeForceKernel(name, platform, cl); return new OpenCLCalcDrudeForceKernel(name, platform, cl);
if (name == IntegrateDrudeLangevinStepKernel::Name()) if (name == IntegrateDrudeLangevinStepKernel::Name())
return new OpenCLIntegrateDrudeLangevinStepKernel(name, platform, cl); return new OpenCLIntegrateDrudeLangevinStepKernel(name, platform, cl);
if (name == IntegrateDrudeSCFStepKernel::Name())
return new OpenCLIntegrateDrudeSCFStepKernel(name, platform, cl);
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/opencl/src/OpenCLDrudeKernels.cpp
View file @ 178aa003
...@@ -42,6 +42,13 @@ ...@@ -42,6 +42,13 @@
using namespace OpenMM; using namespace OpenMM;
using namespace std; using namespace std;
static void setPosqCorrectionArg(OpenCLContext& cl, cl::Kernel& kernel, int index) {
if (cl.getUseMixedPrecision())
kernel.setArg<cl::Buffer>(index, cl.getPosqCorrection().getDeviceBuffer());
else
kernel.setArg<void*>(index, NULL);
}
class OpenCLDrudeForceInfo : public OpenCLForceInfo { class OpenCLDrudeForceInfo : public OpenCLForceInfo {
public: public:
OpenCLDrudeForceInfo(const DrudeForce& force) : OpenCLForceInfo(0), force(force) { OpenCLDrudeForceInfo(const DrudeForce& force) : OpenCLForceInfo(0), force(force) {
...@@ -177,7 +184,59 @@ double OpenCLCalcDrudeForceKernel::execute(ContextImpl& context, bool includeFor ...@@ -177,7 +184,59 @@ double OpenCLCalcDrudeForceKernel::execute(ContextImpl& context, bool includeFor
} }
void OpenCLCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) { void OpenCLCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) {
int numContexts = cl.getPlatformData().contexts.size();
// Set the particle parameters.
int startParticleIndex = cl.getContextIndex()*force.getNumParticles()/numContexts;
int endParticleIndex = (cl.getContextIndex()+1)*force.getNumParticles()/numContexts;
int numParticles = endParticleIndex-startParticleIndex;
if (numParticles > 0) {
if (particleParams == NULL || numParticles != particleParams->getSize())
throw OpenMMException("updateParametersInContext: The number of Drude particles has changed");
vector<mm_float4> paramVector(numParticles);
for (int i = 0; i < numParticles; i++) {
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(startParticleIndex+i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
double a1 = (p2 == -1 ? 1 : aniso12);
double a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34);
double a3 = 3-a1-a2;
double k3 = charge*charge/(polarizability*a3);
double k1 = charge*charge/(polarizability*a1) - k3;
double k2 = charge*charge/(polarizability*a2) - k3;
if (p2 == -1)
k1 = 0;
if (p3 == -1 || p4 == -1)
k2 = 0;
paramVector[i] = mm_float4((float) k1, (float) k2, (float) k3, 0.0f);
}
particleParams->upload(paramVector);
}
// Set the pair parameters.
int startPairIndex = cl.getContextIndex()*force.getNumScreenedPairs()/numContexts;
int endPairIndex = (cl.getContextIndex()+1)*force.getNumScreenedPairs()/numContexts;
int numPairs = endPairIndex-startPairIndex;
if (numPairs > 0) {
if (pairParams == NULL || numPairs != pairParams->getSize())
throw OpenMMException("updateParametersInContext: The number of screened pairs has changed");
vector<mm_float2> paramVector(numPairs);
for (int i = 0; i < numPairs; i++) {
int drude1, drude2;
double thole;
force.getScreenedPairParameters(startPairIndex+i, drude1, drude2, thole);
int p, p1, p2, p3, p4;
double charge1, charge2, polarizability1, polarizability2, aniso12, aniso34;
force.getParticleParameters(drude1, p, p1, p2, p3, p4, charge1, polarizability1, aniso12, aniso34);
force.getParticleParameters(drude2, p, p1, p2, p3, p4, charge2, polarizability2, aniso12, aniso34);
double screeningScale = thole/pow(polarizability1*polarizability2, 1.0/6.0);
double energyScale = ONE_4PI_EPS0*charge1*charge2;
paramVector[i] = mm_float2((float) screeningScale, (float) energyScale);
}
pairParams->upload(paramVector);
}
} }
OpenCLIntegrateDrudeLangevinStepKernel::~OpenCLIntegrateDrudeLangevinStepKernel() { OpenCLIntegrateDrudeLangevinStepKernel::~OpenCLIntegrateDrudeLangevinStepKernel() {
...@@ -310,3 +369,235 @@ void OpenCLIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const ...@@ -310,3 +369,235 @@ void OpenCLIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const
double OpenCLIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) { double OpenCLIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
return cl.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize()); return cl.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
} }
class OpenCLIntegrateDrudeSCFStepKernel::ReorderListener : public OpenCLContext::ReorderListener {
public:
ReorderListener(OpenCLContext& cl, const vector<int>& drudeParticles, vector<int>& reorderedDrudeParticles) :
cl(cl), drudeParticles(drudeParticles), reorderedDrudeParticles(reorderedDrudeParticles) {
}
void execute() {
const vector<int>& order = cl.getAtomIndex();
int numParticles = order.size();
vector<int> inverseOrder(numParticles);
for (int i = 0; i < numParticles; i++)
inverseOrder[order[i]] = i;
int numDrudeParticles = drudeParticles.size();
for (int i = 0; i < numDrudeParticles; i++)
reorderedDrudeParticles[i] = inverseOrder[drudeParticles[i]];
}
private:
OpenCLContext& cl;
const vector<int>& drudeParticles;
vector<int>& reorderedDrudeParticles;
};
OpenCLIntegrateDrudeSCFStepKernel::~OpenCLIntegrateDrudeSCFStepKernel() {
if (minimizerPos != NULL)
lbfgs_free(minimizerPos);
}
void OpenCLIntegrateDrudeSCFStepKernel::initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) {
cl.getPlatformData().initializeContexts(system);
// Identify Drude particles.
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);
reorderedDrudeParticles.push_back(p);
}
cl.addReorderListener(new ReorderListener(cl, drudeParticles, reorderedDrudeParticles));
// Initialize the energy minimizer.
minimizerPos = lbfgs_malloc(drudeParticles.size()*3);
if (minimizerPos == NULL)
throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory");
lbfgs_parameter_init(&minimizerParams);
minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
// Create the kernels.
cl::Program program = cl.createProgram(OpenCLKernelSources::verlet, "");
kernel1 = cl::Kernel(program, "integrateVerletPart1");
kernel2 = cl::Kernel(program, "integrateVerletPart2");
prevStepSize = -1.0;
}
void OpenCLIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
int numAtoms = cl.getNumAtoms();
double dt = integrator.getStepSize();
if (!hasInitializedKernels) {
hasInitializedKernels = true;
kernel1.setArg<cl_int>(0, numAtoms);
kernel1.setArg<cl::Buffer>(1, cl.getIntegrationUtilities().getStepSize().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
setPosqCorrectionArg(cl, kernel1, 3);
kernel1.setArg<cl::Buffer>(4, cl.getVelm().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(5, cl.getForce().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(6, integration.getPosDelta().getDeviceBuffer());
kernel2.setArg<cl_int>(0, numAtoms);
kernel2.setArg<cl::Buffer>(1, cl.getIntegrationUtilities().getStepSize().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
setPosqCorrectionArg(cl, kernel2, 3);
kernel2.setArg<cl::Buffer>(4, cl.getVelm().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(5, integration.getPosDelta().getDeviceBuffer());
}
if (dt != prevStepSize) {
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
vector<mm_double2> stepSizeVec(1);
stepSizeVec[0] = mm_double2(dt, dt);
cl.getIntegrationUtilities().getStepSize().upload(stepSizeVec);
}
else {
vector<mm_float2> stepSizeVec(1);
stepSizeVec[0] = mm_float2((cl_float) dt, (cl_float) dt);
cl.getIntegrationUtilities().getStepSize().upload(stepSizeVec);
}
prevStepSize = dt;
}
// Call the first integration kernel.
cl.executeKernel(kernel1, numAtoms);
// Apply constraints.
integration.applyConstraints(integrator.getConstraintTolerance());
// Call the second integration kernel.
cl.executeKernel(kernel2, numAtoms);
// Update the positions of virtual sites and Drude particles.
integration.computeVirtualSites();
minimize(context, integrator.getMinimizationErrorTolerance());
// Update the time and step count.
cl.setTime(cl.getTime()+dt);
cl.setStepCount(cl.getStepCount()+1);
// Reduce UI lag.
#ifdef WIN32
cl.getQueue().flush();
#endif
}
double OpenCLIntegrateDrudeSCFStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
return cl.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
}
struct MinimizerData {
ContextImpl& context;
OpenCLContext& cl;
vector<int>& reorderedDrudeParticles;
MinimizerData(ContextImpl& context, OpenCLContext& cl, vector<int>& reorderedDrudeParticles) : context(context), cl(cl), reorderedDrudeParticles(reorderedDrudeParticles) {}
};
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;
OpenCLContext& cl = data->cl;
vector<int>& reorderedDrudeParticles = data->reorderedDrudeParticles;
int numDrudeParticles = reorderedDrudeParticles.size();
// Set the particle positions.
cl.getPosq().download(cl.getPinnedBuffer());
mm_double4 periodicBoxSize = cl.getPeriodicBoxSizeDouble();
if (cl.getUseDoublePrecision()) {
mm_double4* posq = (mm_double4*) cl.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
mm_double4& p = posq[reorderedDrudeParticles[i]];
mm_int4 offset = cl.getPosCellOffsets()[reorderedDrudeParticles[i]];
p.x = x[3*i]+offset.x*periodicBoxSize.x;
p.y = x[3*i+1]+offset.y*periodicBoxSize.y;
p.z = x[3*i+2]+offset.z*periodicBoxSize.z;
}
}
else {
mm_float4* posq = (mm_float4*) cl.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
mm_float4& p = posq[reorderedDrudeParticles[i]];
mm_int4 offset = cl.getPosCellOffsets()[reorderedDrudeParticles[i]];
p.x = x[3*i]+offset.x*periodicBoxSize.x;
p.y = x[3*i+1]+offset.y*periodicBoxSize.y;
p.z = x[3*i+2]+offset.z*periodicBoxSize.z;
}
}
cl.getPosq().upload(cl.getPinnedBuffer());
// Compute the forces and energy for this configuration.
double energy = context.calcForcesAndEnergy(true, true);
cl.getForce().download(cl.getPinnedBuffer());
if (cl.getUseDoublePrecision()) {
mm_double4* force = (mm_double4*) cl.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
int index = reorderedDrudeParticles[i];
g[3*i] = -force[index].x;
g[3*i+1] = -force[index].y;
g[3*i+2] = -force[index].z;
}
}
else {
mm_float4* force = (mm_float4*) cl.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
int index = reorderedDrudeParticles[i];
g[3*i] = -force[index].x;
g[3*i+1] = -force[index].y;
g[3*i+2] = -force[index].z;
}
}
return energy;
}
void OpenCLIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double tolerance) {
// Record the initial positions.
int numDrudeParticles = reorderedDrudeParticles.size();
cl.getPosq().download(cl.getPinnedBuffer());
mm_double4 periodicBoxSize = cl.getPeriodicBoxSizeDouble();
if (cl.getUseDoublePrecision()) {
mm_double4* posq = (mm_double4*) cl.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
mm_double4 p = posq[reorderedDrudeParticles[i]];
mm_int4 offset = cl.getPosCellOffsets()[reorderedDrudeParticles[i]];
minimizerPos[3*i] = p.x-offset.x*periodicBoxSize.x;
minimizerPos[3*i+1] = p.y-offset.y*periodicBoxSize.y;
minimizerPos[3*i+2] = p.z-offset.z*periodicBoxSize.z;
}
}
else {
mm_float4* posq = (mm_float4*) cl.getPinnedBuffer();
for (int i = 0; i < numDrudeParticles; ++i) {
mm_float4 p = posq[reorderedDrudeParticles[i]];
mm_int4 offset = cl.getPosCellOffsets()[reorderedDrudeParticles[i]];
minimizerPos[3*i] = p.x-offset.x*periodicBoxSize.x;
minimizerPos[3*i+1] = p.y-offset.y*periodicBoxSize.y;
minimizerPos[3*i+2] = p.z-offset.z*periodicBoxSize.z;
}
minimizerParams.xtol = 1e-7;
}
// Determine a normalization constant for scaling the tolerance.
double norm = 0.0;
for (int i = 0; i < 3*numDrudeParticles; i++)
norm += minimizerPos[i]*minimizerPos[i];
norm /= numDrudeParticles;
norm = (norm < 1 ? 1 : sqrt(norm));
minimizerParams.epsilon = tolerance/norm;
// Perform the minimization.
lbfgsfloatval_t fx;
MinimizerData data(context, cl, reorderedDrudeParticles);
lbfgs(numDrudeParticles*3, minimizerPos, &fx, evaluate, NULL, &data, &minimizerParams);
}
\ No newline at end of file
plugins/drude/platforms/opencl/src/OpenCLDrudeKernels.h
View file @ 178aa003
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include "openmm/DrudeKernels.h" #include "openmm/DrudeKernels.h"
#include "OpenCLContext.h" #include "OpenCLContext.h"
#include "OpenCLArray.h" #include "OpenCLArray.h"
#include "lbfgs.h"
namespace OpenMM { namespace OpenMM {
...@@ -116,6 +117,50 @@ private: ...@@ -116,6 +117,50 @@ private:
cl::Kernel kernel1, kernel2; cl::Kernel kernel1, kernel2;
}; };
/**
* This kernel is invoked by DrudeSCFIntegrator to take one time step
*/
class OpenCLIntegrateDrudeSCFStepKernel : public IntegrateDrudeSCFStepKernel {
public:
OpenCLIntegrateDrudeSCFStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) :
IntegrateDrudeSCFStepKernel(name, platform), cl(cl), hasInitializedKernels(false), minimizerPos(NULL) {
}
~OpenCLIntegrateDrudeSCFStepKernel();
/**
* 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:
class ReorderListener;
void minimize(ContextImpl& context, double tolerance);
OpenCLContext& cl;
bool hasInitializedKernels;
double prevStepSize;
std::vector<int> drudeParticles;
std::vector<int> reorderedDrudeParticles;
lbfgsfloatval_t *minimizerPos;
lbfgs_parameter_t minimizerParams;
cl::Kernel kernel1, kernel2;
};
} // namespace OpenMM } // namespace OpenMM
#endif /*OPENCL_DRUDE_KERNELS_H_*/ #endif /*OPENCL_DRUDE_KERNELS_H_*/
plugins/drude/platforms/opencl/tests/TestOpenCLDrudeForce.cpp
View file @ 178aa003
...@@ -156,6 +156,43 @@ void testThole() { ...@@ -156,6 +156,43 @@ void testThole() {
validateForce(system, positions, energySpring1+energySpring2+energyDipole); validateForce(system, positions, energySpring1+energySpring2+energyDipole);
} }
void testChangingParameters() {
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
Platform& platform = Platform::getPlatformByName("OpenCL");
// Create the system.
System system;
system.addParticle(1.0);
system.addParticle(1.0);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
system.addForce(drude);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(2, 0, 0);
// Check the energy.
VerletIntegrator integ(1.0);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*k*3*3, state.getPotentialEnergy(), 1e-5);
// Modify the parameters.
const double k2 = 2.2;
const double charge2 = 0.3;
const double alpha2 = charge2*charge2/k2;
drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
drude->updateParametersInContext(context);
state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*k2*3*3, state.getPotentialEnergy(), 1e-5);
}
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
registerDrudeOpenCLKernelFactories(); registerDrudeOpenCLKernelFactories();
...@@ -164,6 +201,7 @@ int main(int argc, char* argv[]) { ...@@ -164,6 +201,7 @@ int main(int argc, char* argv[]) {
testSingleParticle(); testSingleParticle();
testAnisotropicParticle(); testAnisotropicParticle();
testThole(); testThole();
testChangingParameters();
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/drude/platforms/opencl/tests/TestOpenCLDrudeSCFIntegrator.cpp 0 → 100644
View file @ 178aa003
/* -------------------------------------------------------------------------- *
* 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 OpenCL 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;
extern "C" OPENMM_EXPORT void registerDrudeOpenCLKernelFactories();
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("OpenCL");
Context context(system, integ, platform);
context.setPositions(positions);
context.applyConstraints(1e-5);
context.setVelocitiesToTemperature(300.0);
State state = context.getState(State::Energy);
double initialEnergy;
int numSteps = 1000;
double maxNorm = (platform.getPropertyValue(context, "OpenCLPrecision") == "double" ? 1.0 : 5.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
state = context.getState(State::Energy | State::Forces);
if (i == 0)
initialEnergy = state.getPotentialEnergy()+state.getKineticEnergy();
else
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 < maxNorm);
}
}
int main(int argc, char* argv[]) {
try {
registerDrudeOpenCLKernelFactories();
if (argc > 1)
Platform::getPlatformByName("OpenCL").setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
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/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @ 178aa003
...@@ -418,8 +418,9 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons ...@@ -418,8 +418,9 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons
if (minimizerPos == NULL) if (minimizerPos == NULL)
throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory"); throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory");
lbfgs_parameter_init(&minimizerParams); lbfgs_parameter_init(&minimizerParams);
minimizerParams.max_iterations = 0;
minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE; minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
if (sizeof(RealOpenMM) < 8)
minimizerParams.xtol = 1e-7;
} }
void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) { void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
......
plugins/drude/platforms/reference/tests/TestReferenceDrudeForce.cpp
View file @ 178aa003
...@@ -154,11 +154,49 @@ void testThole() { ...@@ -154,11 +154,49 @@ void testThole() {
validateForce(system, positions, energySpring1+energySpring2+energyDipole); validateForce(system, positions, energySpring1+energySpring2+energyDipole);
} }
void testChangingParameters() {
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
Platform& platform = Platform::getPlatformByName("Reference");
// Create the system.
System system;
system.addParticle(1.0);
system.addParticle(1.0);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
system.addForce(drude);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(2, 0, 0);
// Check the energy.
VerletIntegrator integ(1.0);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*k*3*3, state.getPotentialEnergy(), 1e-5);
// Modify the parameters.
const double k2 = 2.2;
const double charge2 = 0.3;
const double alpha2 = charge2*charge2/k2;
drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
drude->updateParametersInContext(context);
state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*k2*3*3, state.getPotentialEnergy(), 1e-5);
}
int main() { int main() {
try { try {
testSingleParticle(); testSingleParticle();
testAnisotropicParticle(); testAnisotropicParticle();
testThole(); testThole();
testChangingParameters();
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/drude/platforms/reference/tests/TestReferenceDrudeSCFIntegrator.cpp
View file @ 178aa003
...@@ -106,11 +106,14 @@ void testWater() { ...@@ -106,11 +106,14 @@ void testWater() {
context.applyConstraints(1e-5); context.applyConstraints(1e-5);
context.setVelocitiesToTemperature(300.0); context.setVelocitiesToTemperature(300.0);
State state = context.getState(State::Energy); State state = context.getState(State::Energy);
double initialEnergy = state.getPotentialEnergy()+state.getKineticEnergy(); double initialEnergy;
int numSteps = 1000; int numSteps = 1000;
for (int i = 0; i < numSteps; i++) { for (int i = 0; i < numSteps; i++) {
integ.step(1); integ.step(1);
state = context.getState(State::Energy | State::Forces); state = context.getState(State::Energy | State::Forces);
if (i == 0)
initialEnergy = state.getPotentialEnergy()+state.getKineticEnergy();
else
ASSERT_EQUAL_TOL(initialEnergy, state.getPotentialEnergy()+state.getKineticEnergy(), 0.01); ASSERT_EQUAL_TOL(initialEnergy, state.getPotentialEnergy()+state.getKineticEnergy(), 0.01);
const vector<Vec3>& force = state.getForces(); const vector<Vec3>& force = state.getForces();
double norm = 0.0; double norm = 0.0;
......
wrappers/python/src/swig_doxygen/OpenMM.i
View file @ 178aa003
...@@ -45,8 +45,7 @@ namespace std { ...@@ -45,8 +45,7 @@ namespace std {
#include "OpenMM.h" #include "OpenMM.h"
#include "OpenMMAmoeba.h" #include "OpenMMAmoeba.h"
#include "openmm/RPMDIntegrator.h" #include "openmm/RPMDIntegrator.h"
#include "openmm/DrudeForce.h" #include "OpenMMDrude.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/serialization/SerializationNode.h" #include "openmm/serialization/SerializationNode.h"
#include "openmm/serialization/SerializationProxy.h" #include "openmm/serialization/SerializationProxy.h"
#include "openmm/serialization/XmlSerializer.h" #include "openmm/serialization/XmlSerializer.h"
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 178aa003
...@@ -141,6 +141,7 @@ SKIP_METHODS = [('State',), ...@@ -141,6 +141,7 @@ SKIP_METHODS = [('State',),
('DrudeForceImpl',), ('DrudeForceImpl',),
('CalcDrudeForceKernel',), ('CalcDrudeForceKernel',),
('IntegrateDrudeLangevinStepKernel',), ('IntegrateDrudeLangevinStepKernel',),
('IntegrateDrudeSCFStepKernel',),
('XmlSerializer', 'serialize'), ('XmlSerializer', 'serialize'),
('XmlSerializer', 'deserialize'), ('XmlSerializer', 'deserialize'),
] ]
...@@ -444,5 +445,6 @@ UNITS = { ...@@ -444,5 +445,6 @@ UNITS = {
("System", "getVirtualSite") : (None, ()), ("System", "getVirtualSite") : (None, ()),
("DrudeLangevinIntegrator", "getDrudeTemperature") : ("unit.kelvin", ()), ("DrudeLangevinIntegrator", "getDrudeTemperature") : ("unit.kelvin", ()),
("DrudeLangevinIntegrator", "getDrudeFriction") : ("1/unit.picosecond", ()), ("DrudeLangevinIntegrator", "getDrudeFriction") : ("1/unit.picosecond", ()),
("DrudeSCFIntegrator", "getMinimizationErrorTolerance") : ("unit.kilojoules_per_mole/unit.nanometer", ()),
} }
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