Skip to content

GitLab

Commit aa96846e authored by Peter Eastman's avatar Peter Eastman
Browse files

Created OpenCL implementation of Ewald

parent 014f3406
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1358 additions and 57 deletions
platforms/opencl/src/OpenCLKernels.cpp
View file @ aa96846e
......@@ -29,6 +29,7 @@
#include "openmm/LangevinIntegrator.h"
#include "openmm/Context.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/NonbondedForceImpl.h"
#include "OpenCLIntegrationUtilities.h"
#include "OpenCLNonbondedUtilities.h"
#include <cmath>
......@@ -42,6 +43,19 @@ static const double AVOGADRO = 6.0221367e23; // ()
static const double RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
static const double BOLTZ = (RGAS/KILO); // (kJ/(mol K))
static string doubleToString(double value) {
stringstream s;
s.precision(8);
s << scientific << value << "f";
return s.str();
}
static string intToString(int value) {
stringstream s;
s << value;
return s.str();
}
void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
}
......@@ -522,12 +536,14 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
OpenCLArray<mm_float4>& posq = cl.getPosq();
vector<mm_float2> sigmaEpsilonVector(numParticles);
vector<vector<int> > exclusionList(numParticles);
double sumSquaredCharges = 0.0;
for (int i = 0; i < numParticles; i++) {
double charge, sigma, epsilon;
force.getParticleParameters(i, charge, sigma, epsilon);
posq[i].w = (float) charge;
sigmaEpsilonVector[i] = (mm_float2) {(float) (0.5*sigma), (float) (2.0*sqrt(epsilon))};
exclusionList[i].push_back(i);
sumSquaredCharges += charge*charge;
}
for (int i = 0; i < (int) exclusions.size(); i++) {
exclusionList[exclusions[i].first].push_back(exclusions[i].second);
......@@ -537,68 +553,56 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
sigmaEpsilon->upload(sigmaEpsilonVector);
bool useCutoff = (force.getNonbondedMethod() != NonbondedForce::NoCutoff);
bool usePeriodic = (force.getNonbondedMethod() != NonbondedForce::NoCutoff && force.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic);
Vec3 boxVectors[3];
system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
map<string, string> defines;
if (useCutoff) {
// Compute the reaction field constants.
double reactionFieldK = pow(force.getCutoffDistance(), -3.0)*(force.getReactionFieldDielectric()-1.0)/(2.0*force.getReactionFieldDielectric()+1.0);
double reactionFieldC = (1.0 / force.getCutoffDistance())*(3.0*force.getReactionFieldDielectric())/(2.0*force.getReactionFieldDielectric()+1.0);
stringstream k, c;
k.precision(8);
c.precision(8);
k << scientific << reactionFieldK << "f";
c << scientific << reactionFieldC << "f";
defines["REACTION_FIELD_K"] = k.str();
defines["REACTION_FIELD_C"] = c.str();
defines["REACTION_FIELD_K"] = doubleToString(reactionFieldK);
defines["REACTION_FIELD_C"] = doubleToString(reactionFieldC);
}
// if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
// method = CUTOFF;
// }
// if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
// method = PERIODIC;
// }
// if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
// double ewaldErrorTol = force.getEwaldErrorTolerance();
// double alpha = (1.0/force.getCutoffDistance())*std::sqrt(-std::log(ewaldErrorTol));
// double mx = boxVectors[0][0]/force.getCutoffDistance();
// double my = boxVectors[1][1]/force.getCutoffDistance();
// double mz = boxVectors[2][2]/force.getCutoffDistance();
// double pi = 3.1415926535897932385;
// int kmaxx = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
// int kmaxy = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
// int kmaxz = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
// if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
// if (kmaxx%2 == 0)
// kmaxx++;
// if (kmaxy%2 == 0)
// kmaxy++;
// if (kmaxz%2 == 0)
// kmaxz++;
// gpuSetEwaldParameters(gpu, (float) alpha, kmaxx, kmaxy, kmaxz);
// method = EWALD;
// }
// else {
// int gridSizeX = -0.5*kmaxx*std::log(ewaldErrorTol);
// int gridSizeY = -0.5*kmaxy*std::log(ewaldErrorTol);
// int gridSizeZ = -0.5*kmaxz*std::log(ewaldErrorTol);
// gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
// method = PARTICLE_MESH_EWALD;
// }
// }
// data.nonbondedMethod = method;
// gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
// Compute the Ewald parameters.
double alpha;
int kmaxx, kmaxy, kmaxz;
NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmaxx, kmaxy, kmaxz);
defines["EWALD_ALPHA"] = doubleToString(alpha);
defines["TWO_OVER_SQRT_PI"] = doubleToString(2.0/sqrt(M_PI));
defines["USE_EWALD"] = "1";
double selfEnergyScale = 138.935485*alpha/std::sqrt(M_PI);
ewaldSelfEnergy = - 138.935485*alpha*sumSquaredCharges/std::sqrt(M_PI);
// Create the reciprocal space kernels.
map<string, string> replacements;
replacements["NUM_ATOMS"] = intToString(numParticles);
replacements["KMAX_X"] = intToString(kmaxx);
replacements["KMAX_Y"] = intToString(kmaxy);
replacements["KMAX_Z"] = intToString(kmaxz);
replacements["RECIPROCAL_BOX_SIZE_X"] = doubleToString(2.0*M_PI/boxVectors[0][0]);
replacements["RECIPROCAL_BOX_SIZE_Y"] = doubleToString(2.0*M_PI/boxVectors[1][1]);
replacements["RECIPROCAL_BOX_SIZE_Z"] = doubleToString(2.0*M_PI/boxVectors[2][2]);
replacements["RECIPROCAL_COEFFICIENT"] = doubleToString(138.935485*4*M_PI/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]));
replacements["EXP_COEFFICIENT"] = doubleToString(-1.0/(4.0*alpha*alpha));
cl::Program program = cl.createProgram(cl.loadSourceFromFile("ewald.cl"), replacements);
ewaldSumsKernel = cl::Kernel(program, "calculateEwaldCosSinSums");
ewaldForcesKernel = cl::Kernel(program, "calculateEwaldForces");
cosSinSums = new OpenCLArray<mm_float2>(cl, (2*kmaxx-1)*(2*kmaxy-1)*(2*kmaxz-1), "cosSinSums");
}
else
ewaldSelfEnergy = 0.0;
// Add the interaction to the default nonbonded kernel.
string source = cl.loadSourceFromFile("coulombLennardJones.cl", defines);
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float2", sizeof(cl_float2), sigmaEpsilon->getDeviceBuffer()));
cutoffSquared = force.getCutoffDistance()*force.getCutoffDistance();
// Compute the Ewald self energy.
ewaldSelfEnergy = 0.0;
if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
// double selfEnergyScale = gpu->sim.epsfac*gpu->sim.alphaEwald/std::sqrt(PI);
// for (int i = 0; i < numParticles; i++)
// ewaldSelfEnergy -= selfEnergyScale*q[i]*q[i];
}
// Initialize the exceptions.
int numExceptions = exceptions.size();
......@@ -645,6 +649,16 @@ void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
exceptionsKernel.setArg<cl::Buffer>(8, exceptionIndices->getDeviceBuffer());
cl.executeKernel(exceptionsKernel, numExceptions);
}
if (cosSinSums != NULL) {
ewaldSumsKernel.setArg<cl::Buffer>(0, cl.getEnergyBuffer().getDeviceBuffer());
ewaldSumsKernel.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
ewaldSumsKernel.setArg<cl::Buffer>(2, cosSinSums->getDeviceBuffer());
cl.executeKernel(ewaldSumsKernel, cosSinSums->getSize());
ewaldForcesKernel.setArg<cl::Buffer>(0, cl.getForceBuffers().getDeviceBuffer());
ewaldForcesKernel.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
ewaldForcesKernel.setArg<cl::Buffer>(2, cosSinSums->getDeviceBuffer());
cl.executeKernel(ewaldForcesKernel, cl.getNumAtoms());
}
}
double OpenCLCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
......
platforms/opencl/src/OpenCLKernels.h
View file @ aa96846e
......@@ -328,7 +328,10 @@ private:
OpenCLArray<mm_float2>* sigmaEpsilon;
OpenCLArray<mm_float4>* exceptionParams;
OpenCLArray<mm_int4>* exceptionIndices;
OpenCLArray<mm_float2>* cosSinSums;
cl::Kernel exceptionsKernel;
cl::Kernel ewaldSumsKernel;
cl::Kernel ewaldForcesKernel;
double cutoffSquared, ewaldSelfEnergy;
};
......
platforms/opencl/src/kernels/coulombLennardJones.cl
View file @ aa96846e
#if USE_EWALD
if (r2 < CUTOFF_SQUARED) {
bool needCorrection = isExcluded && atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS;
if (!isExcluded || needCorrection) {
const float prefactor = 138.935485f*posq1.w*posq2.w*invR;
float alphaR = EWALD_ALPHA*r;
float erfcAlphaR = erfc(alphaR);
float tempForce;
if (needCorrection) {
// Subtract off the part of this interaction that was included in the reciprocal space contribution.
tempForce = -prefactor*((1.0f-erfcAlphaR)-alphaR*exp(-alphaR*alphaR)*TWO_OVER_SQRT_PI);
tempEnergy += -prefactor*(1.0f-erfcAlphaR);
}
else {
float sig = sigmaEpsilon1.x + sigmaEpsilon2.x;
float sig2 = invR*sig;
sig2 *= sig2;
float sig6 = sig2*sig2*sig2;
float eps = sigmaEpsilon1.y*sigmaEpsilon2.y;
tempForce = eps*(12.0f*sig6 - 6.0f)*sig6 + prefactor*(erfcAlphaR+alphaR*exp(-alphaR*alphaR)*TWO_OVER_SQRT_PI);
tempEnergy += eps*(sig6 - 1.0f)*sig6 + prefactor*erfcAlphaR;
}
dEdR += tempForce*invR*invR;
}
}
#else
#ifdef USE_CUTOFF
if (!isExcluded && r2 < CUTOFF_SQUARED) {
#else
......@@ -10,13 +37,15 @@ if (!isExcluded) {
float eps = sigmaEpsilon1.y*sigmaEpsilon2.y;
float tempForce = eps*(12.0f*sig6 - 6.0f)*sig6;
tempEnergy += eps*(sig6 - 1.0f)*sig6;
const float EpsilonFactor = 138.935485f;
#ifdef USE_CUTOFF
tempForce += EpsilonFactor*posq1.w*posq2.w*(invR - 2.0f*REACTION_FIELD_K*r2);
tempEnergy += EpsilonFactor*posq1.w*posq2.w*(invR + REACTION_FIELD_K*r2 - REACTION_FIELD_C);
const float prefactor = 138.935485f*posq1.w*posq2.w;
tempForce += prefactor*(invR - 2.0f*REACTION_FIELD_K*r2);
tempEnergy += prefactor*(invR + REACTION_FIELD_K*r2 - REACTION_FIELD_C);
#else
tempForce += EpsilonFactor*posq1.w*posq2.w*invR;
tempEnergy += EpsilonFactor*posq1.w*posq2.w*invR;
const float prefactor = 138.935485f*posq1.w*posq2.w*invR;
tempForce += prefactor;
tempEnergy += prefactor;
#endif
dEdR += tempForce*invR*invR;
}
#endif
\ No newline at end of file
platforms/opencl/src/kernels/ewald.cl 0 → 100644
View file @ aa96846e
float2 multofFloat2(float2 a, float2 b) {
return (float2) (a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}
/**
* Precompute the cosine and sine sums which appear in each force term.
*/
__kernel void calculateEwaldCosSinSums(__global float* energyBuffer, __global float4* posq, __global float2* cosSinSum) {
const unsigned int ksizex = 2*KMAX_X-1;
const unsigned int ksizey = 2*KMAX_Y-1;
const unsigned int ksizez = 2*KMAX_Z-1;
const unsigned int totalK = ksizex*ksizey*ksizez;
unsigned int index = get_global_id(0);
float energy = 0.0f;
while (index < (KMAX_Y-1)*ksizez+KMAX_Z)
index += get_global_size(0);
while (index < totalK) {
// Find the wave vector (kx, ky, kz) this index corresponds to.
int rx = index/(ksizey*ksizez);
int remainder = index - rx*ksizey*ksizez;
int ry = remainder/ksizez;
int rz = remainder - ry*ksizez - KMAX_Z + 1;
ry += -KMAX_Y + 1;
float kx = rx*RECIPROCAL_BOX_SIZE_X;
float ky = ry*RECIPROCAL_BOX_SIZE_Y;
float kz = rz*RECIPROCAL_BOX_SIZE_Z;
// Compute the sum for this wave vector.
float2 sum = 0.0f;
for (int atom = 0; atom < NUM_ATOMS; atom++) {
float4 apos = posq[atom];
float phase = apos.x*kx;
float2 structureFactor = (float2) (cos(phase), sin(phase));
phase = apos.y*ky;
structureFactor = multofFloat2(structureFactor, (float2) (cos(phase), sin(phase)));
phase = apos.z*kz;
structureFactor = multofFloat2(structureFactor, (float2) (cos(phase), sin(phase)));
sum += apos.w*structureFactor;
}
cosSinSum[index] = sum;
// Compute the contribution to the energy.
float k2 = kx*kx + ky*ky + kz*kz;
float ak = exp(k2*EXP_COEFFICIENT) / k2;
energy += RECIPROCAL_COEFFICIENT*ak*(sum.x*sum.x + sum.y*sum.y);
index += get_global_size(0);
}
energyBuffer[get_global_id(0)] += energy;
}
/**
* Compute the reciprocal space part of the Ewald force, using the precomputed sums from the
* previous routine.
*/
__kernel void calculateEwaldForces(__global float4* forceBuffers, __global float4* posq, __global float2* cosSinSum) {
unsigned int atom = get_global_id(0);
while (atom < NUM_ATOMS) {
float4 force = forceBuffers[atom];
float4 apos = posq[atom];
// Loop over all wave vectors.
int lowry = 0;
int lowrz = 1;
for (int rx = 0; rx < KMAX_X; rx++) {
float kx = rx*RECIPROCAL_BOX_SIZE_X;
for (int ry = lowry; ry < KMAX_Y; ry++) {
float ky = ry*RECIPROCAL_BOX_SIZE_Y;
float phase = apos.x*kx;
float2 tab_xy = (float2) (cos(phase), sin(phase));
phase = apos.y*ky;
tab_xy = multofFloat2(tab_xy, (float2) (cos(phase), sin(phase)));
for (int rz = lowrz; rz < KMAX_Z; rz++) {
float kz = rz*RECIPROCAL_BOX_SIZE_Z;
// Compute the force contribution of this wave vector.
int index = rx*(KMAX_Y*2-1)*(KMAX_Z*2-1) + (ry+KMAX_Y-1)*(KMAX_Z*2-1) + (rz+KMAX_Z-1);
float k2 = kx*kx + ky*ky + kz*kz;
float ak = exp(k2*EXP_COEFFICIENT)/k2;
phase = apos.z*kz;
float2 structureFactor = multofFloat2(tab_xy, (float2) (cos(phase), sin(phase)));
float2 sum = cosSinSum[index];
float dEdR = 2*RECIPROCAL_COEFFICIENT*ak*apos.w*(sum.x*structureFactor.y - sum.y*structureFactor.x);
force.x += dEdR*kx;
force.y += dEdR*ky;
force.z += dEdR*kz;
lowrz = 1 - KMAX_Z;
}
lowry = 1 - KMAX_Y;
}
}
// Record the force on the atom.
forceBuffers[atom] = force;
atom += get_global_size(0);
}
}
platforms/opencl/tests/TestOpenCLEwald.cpp 0 → 100644
View file @ aa96846e
/* -------------------------------------------------------------------------- *
* 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-2009 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 Ewald summation method OpenCL implementation of NonbondedForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testEwaldPME() {
// Use amorphous NaCl system for the tests
const int numParticles = 894;
const double cutoff = 1.2;
const double boxSize = 3.00646;
double tol = 1e-5;
OpenCLPlatform cl;
ReferencePlatform reference;
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(tol);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
// (1) Check whether the Reference and OpenCL platforms agree when using Ewald Method
Context clContext(system, integrator, cl);
Context referenceContext(system, integrator, reference);
clContext.setPositions(positions);
referenceContext.setPositions(positions);
State clState = clContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], clState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), clState.getPotentialEnergy(), tol);
// (2) Check whether Ewald method in OpenCL is self-consistent
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = clState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = clState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
Context clContext2(system, integrator, cl);
clContext2.setPositions(positions);
tol = 1e-3;
State clState2 = clContext2.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (clState2.getPotentialEnergy()-clState.getPotentialEnergy())/delta, tol)
// (3) Check whether the Reference and OpenCL platforms agree when using PME
/*
nonbonded->setNonbondedMethod(NonbondedForce::PME);
clContext.reinitialize();
referenceContext.reinitialize();
clContext.setPositions(positions);
referenceContext.setPositions(positions);
clState = clContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], clState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), clState.getPotentialEnergy(), tol);
// (4) Check whether PME method in OpenCL is self-consistent
norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = clState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = clState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
Context clContext3(system, integrator, cl);
clContext3.setPositions(positions);
tol = 1e-3;
State clState3 = clContext3.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (clState3.getPotentialEnergy()-clState.getPotentialEnergy())/delta, tol)*/
}
void testEwald2Ions() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
}
void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 5.0;
System system;
system.setPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
init_gen_rand(0);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(), boxWidth*genrand_real2(), boxWidth*genrand_real2());
}
force->setNonbondedMethod(method);
OpenCLPlatform platform;
VerletIntegrator integrator(0.01);
// For various values of the cutoff and error tolerance, see if the actual error is reasonable.
for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff += 0.2) {
force->setCutoffDistance(cutoff);
vector<Vec3> refForces;
double norm = 0.0;
for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
force->setEwaldErrorTolerance(tol);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces);
if (refForces.size() == 0) {
refForces = state.getForces();
for (int i = 0; i < numParticles; i++)
norm += refForces[i].dot(refForces[i]);
norm = sqrt(norm);
}
else {
double diff = 0.0;
for (int i = 0; i < numParticles; i++) {
Vec3 delta = refForces[i]-state.getForces()[i];
diff += delta.dot(delta);
}
diff = sqrt(diff)/norm;
ASSERT(diff < 5*tol);
}
}
}
}
int main() {
try {
testEwaldPME();
// testEwald2Ions();
testErrorTolerance(NonbondedForce::Ewald);
// testErrorTolerance(NonbondedForce::PME);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/opencl/tests/nacl_amorph.dat 0 → 100644
View file @ aa96846e
This diff is collapsed.
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