Skip to content

GitLab

Commit 3976ad91 authored by Peter Eastman's avatar Peter Eastman
Browse files

Began implementing nonbonded forces

parent 4c489434
Hide whitespace changes
Inline Side-by-side
Showing with 1392 additions and 201 deletions
platforms/opencl/src/OpenCLContext.cpp
View file @ 3976ad91
...@@ -28,6 +28,7 @@ ...@@ -28,6 +28,7 @@
#include "OpenCLArray.h" #include "OpenCLArray.h"
#include "OpenCLForceInfo.h" #include "OpenCLForceInfo.h"
#include "OpenCLIntegrationUtilities.h" #include "OpenCLIntegrationUtilities.h"
#include "OpenCLNonbondedUtilities.h"
#include "openmm/Platform.h" #include "openmm/Platform.h"
#include "openmm/System.h" #include "openmm/System.h"
#include <fstream> #include <fstream>
...@@ -61,8 +62,10 @@ OpenCLContext::OpenCLContext(int numParticles, int deviceIndex) : time(0.0), ste ...@@ -61,8 +62,10 @@ OpenCLContext::OpenCLContext(int numParticles, int deviceIndex) : time(0.0), ste
numAtoms = numParticles; numAtoms = numParticles;
paddedNumAtoms = TileSize*((numParticles+TileSize-1)/TileSize); paddedNumAtoms = TileSize*((numParticles+TileSize-1)/TileSize);
numAtomBlocks = (paddedNumAtoms+(TileSize-1))/TileSize; numAtomBlocks = (paddedNumAtoms+(TileSize-1))/TileSize;
numTiles = numAtomBlocks*(numAtomBlocks+1)/2;
numThreadBlocks = device.getInfo<CL_DEVICE_MAX_WORK_ITEM_SIZES>()[0]/ThreadBlockSize; numThreadBlocks = device.getInfo<CL_DEVICE_MAX_WORK_ITEM_SIZES>()[0]/ThreadBlockSize;
nonbonded = new OpenCLNonbondedUtilities(*this);
posq = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "posq", true);
velm = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "velm", true);
// Create utility kernels that are used in multiple places. // Create utility kernels that are used in multiple places.
...@@ -74,26 +77,25 @@ OpenCLContext::OpenCLContext(int numParticles, int deviceIndex) : time(0.0), ste ...@@ -74,26 +77,25 @@ OpenCLContext::OpenCLContext(int numParticles, int deviceIndex) : time(0.0), ste
OpenCLContext::~OpenCLContext() { OpenCLContext::~OpenCLContext() {
for (int i = 0; i < (int) forces.size(); i++) for (int i = 0; i < (int) forces.size(); i++)
delete forces[i]; delete forces[i];
delete posq; if (posq != NULL)
delete velm; delete posq;
delete force; if (velm != NULL)
delete forceBuffers; delete velm;
delete energyBuffer; if (force != NULL)
delete atomIndex; delete force;
delete integration; if (forceBuffers != NULL)
delete forceBuffers;
if (energyBuffer != NULL)
delete energyBuffer;
if (atomIndex != NULL)
delete atomIndex;
if (integration != NULL)
delete integration;
if (nonbonded != NULL)
delete nonbonded;
} }
void OpenCLContext::initialize(const System& system) { void OpenCLContext::initialize(const System& system) {
// forceBufferPerWarp = true;
// numForceBuffers = numThreadBlocks*ThreadBlockSize/TileSize;
// if (numForceBuffers >= numAtomBlocks) {
// // For small systems, it is more efficient to have one force buffer per block of 32 atoms instead of one per warp.
//
// forceBufferPerWarp = false;
// numForceBuffers = numAtomBlocks;
// }
posq = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "posq", true);
velm = new OpenCLArray<mm_float4>(*this, paddedNumAtoms, "velm", true);
for (int i = 0; i < numAtoms; i++) for (int i = 0; i < numAtoms; i++)
(*velm)[i].w = (float) (1.0/system.getParticleMass(i)); (*velm)[i].w = (float) (1.0/system.getParticleMass(i));
velm->upload(); velm->upload();
...@@ -108,6 +110,7 @@ void OpenCLContext::initialize(const System& system) { ...@@ -108,6 +110,7 @@ void OpenCLContext::initialize(const System& system) {
(*atomIndex)[i] = i; (*atomIndex)[i] = i;
atomIndex->upload(); atomIndex->upload();
integration = new OpenCLIntegrationUtilities(*this, system); integration = new OpenCLIntegrationUtilities(*this, system);
nonbonded->initialize(system);
} }
void OpenCLContext::addForce(OpenCLForceInfo* force) { void OpenCLContext::addForce(OpenCLForceInfo* force) {
......
platforms/opencl/src/OpenCLContext.h
View file @ 3976ad91
...@@ -36,6 +36,7 @@ template <class T> ...@@ -36,6 +36,7 @@ template <class T>
class OpenCLArray; class OpenCLArray;
class OpenCLForceInfo; class OpenCLForceInfo;
class OpenCLIntegrationUtilities; class OpenCLIntegrationUtilities;
class OpenCLNonbondedUtilities;
class System; class System;
/** /**
...@@ -210,12 +211,6 @@ public: ...@@ -210,12 +211,6 @@ public:
int getNumThreadBlocks() const { int getNumThreadBlocks() const {
return numThreadBlocks; return numThreadBlocks;
} }
/**
* Get the total number of tiles used for nonbonded computation.
*/
int getNumTiles() const {
return numTiles;
}
/** /**
* Get the number of force buffers. * Get the number of force buffers.
*/ */
...@@ -228,13 +223,18 @@ public: ...@@ -228,13 +223,18 @@ public:
OpenCLIntegrationUtilities& getIntegrationUtilties() { OpenCLIntegrationUtilities& getIntegrationUtilties() {
return *integration; return *integration;
} }
/**
* Get the OpenCLNonbondedUtilities for this context.
*/
OpenCLNonbondedUtilities& getNonbondedUtilties() {
return *nonbonded;
}
private: private:
double time; double time;
int stepCount; int stepCount;
int numAtoms; int numAtoms;
int paddedNumAtoms; int paddedNumAtoms;
int numAtomBlocks; int numAtomBlocks;
int numTiles;
int numThreadBlocks; int numThreadBlocks;
int numForceBuffers; int numForceBuffers;
cl::Context context; cl::Context context;
...@@ -251,6 +251,7 @@ private: ...@@ -251,6 +251,7 @@ private:
OpenCLArray<cl_float>* energyBuffer; OpenCLArray<cl_float>* energyBuffer;
OpenCLArray<cl_int>* atomIndex; OpenCLArray<cl_int>* atomIndex;
OpenCLIntegrationUtilities* integration; OpenCLIntegrationUtilities* integration;
OpenCLNonbondedUtilities* nonbonded;
}; };
} // namespace OpenMM } // namespace OpenMM
......
platforms/opencl/src/OpenCLKernelFactory.cpp
View file @ 3976ad91
...@@ -46,8 +46,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo ...@@ -46,8 +46,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, cl, context.getSystem()); return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, cl, context.getSystem());
if (name == CalcRBTorsionForceKernel::Name()) if (name == CalcRBTorsionForceKernel::Name())
return new OpenCLCalcRBTorsionForceKernel(name, platform, cl, context.getSystem()); return new OpenCLCalcRBTorsionForceKernel(name, platform, cl, context.getSystem());
// if (name == CalcNonbondedForceKernel::Name()) if (name == CalcNonbondedForceKernel::Name())
// return new OpenCLCalcNonbondedForceKernel(name, platform, cl, context.getSystem()); return new OpenCLCalcNonbondedForceKernel(name, platform, cl, context.getSystem());
// if (name == CalcCustomNonbondedForceKernel::Name()) // if (name == CalcCustomNonbondedForceKernel::Name())
// return new OpenCLCalcCustomNonbondedForceKernel(name, platform, cl, context.getSystem()); // return new OpenCLCalcCustomNonbondedForceKernel(name, platform, cl, context.getSystem());
// if (name == CalcGBSAOBCForceKernel::Name()) // if (name == CalcGBSAOBCForceKernel::Name())
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 3976ad91
...@@ -30,6 +30,7 @@ ...@@ -30,6 +30,7 @@
#include "openmm/Context.h" #include "openmm/Context.h"
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
#include "OpenCLIntegrationUtilities.h" #include "OpenCLIntegrationUtilities.h"
#include "OpenCLNonbondedUtilities.h"
#include <cmath> #include <cmath>
using namespace OpenMM; using namespace OpenMM;
...@@ -46,10 +47,12 @@ void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) { ...@@ -46,10 +47,12 @@ void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
void OpenCLCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) { void OpenCLCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
cl.clearBuffer(cl.getForceBuffers()); cl.clearBuffer(cl.getForceBuffers());
cl.getNonbondedUtilties().prepareInteractions();
} }
void OpenCLCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) { void OpenCLCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
cl.reduceBuffer(cl.getForceBuffers(), cl.getNumForceBuffers()); cl.reduceBuffer(cl.getForceBuffers(), cl.getNumForceBuffers());
cl.getNonbondedUtilties().prepareInteractions();
} }
void OpenCLCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) { void OpenCLCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
...@@ -452,150 +455,127 @@ double OpenCLCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) { ...@@ -452,150 +455,127 @@ double OpenCLCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
return 0.0; return 0.0;
} }
//OpenCLCalcNonbondedForceKernel::~OpenCLCalcNonbondedForceKernel() { OpenCLCalcNonbondedForceKernel::~OpenCLCalcNonbondedForceKernel() {
//} if (sigmaEpsilon != NULL)
// delete sigmaEpsilon;
//void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) { }
// if (data.primaryKernel == NULL)
// data.primaryKernel = this; void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
// data.hasNonbonded = true;
// numParticles = force.getNumParticles(); // Identify which exceptions are 1-4 interactions.
// _gpuContext* gpu = data.gpu;
// vector<pair<int, int> > exclusions;
// // Identify which exceptions are 1-4 interactions. vector<int> exceptions;
// for (int i = 0; i < force.getNumExceptions(); i++) {
// vector<pair<int, int> > exclusions; int particle1, particle2;
// vector<int> exceptions; double chargeProd, sigma, epsilon;
// for (int i = 0; i < force.getNumExceptions(); i++) { force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
// int particle1, particle2; exclusions.push_back(pair<int, int>(particle1, particle2));
// double chargeProd, sigma, epsilon; if (chargeProd != 0.0 || epsilon != 0.0)
// force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon); exceptions.push_back(i);
// exclusions.push_back(pair<int, int>(particle1, particle2)); }
// if (chargeProd != 0.0 || epsilon != 0.0)
// exceptions.push_back(i); // Initialize nonbonded interactions.
int numParticles = force.getNumParticles();
sigmaEpsilon = new OpenCLArray<mm_float2>(cl, numParticles, "sigmaEpsilon");
OpenCLArray<mm_float4>& posq = cl.getPosq();
vector<mm_float2> sigmaEpsilonVector(numParticles);
vector<vector<int> > exclusionList(numParticles);
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);
}
for (int i = 0; i < (int) exclusions.size(); i++) {
exclusionList[exclusions[i].first].push_back(exclusions[i].second);
exclusionList[exclusions[i].second].push_back(exclusions[i].first);
}
posq.upload();
sigmaEpsilon->upload(sigmaEpsilonVector);
bool useCutoff = (force.getNonbondedMethod() != NonbondedForce::NoCutoff);
bool usePeriodic = (force.getNonbondedMethod() != NonbondedForce::NoCutoff && force.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic);
// if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
// method = CUTOFF;
// } // }
// // if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
// // Initialize nonbonded interactions. // method = PERIODIC;
// // }
// { // if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
// vector<int> particle(numParticles); // double ewaldErrorTol = force.getEwaldErrorTolerance();
// vector<float> c6(numParticles); // double alpha = (1.0/force.getCutoffDistance())*std::sqrt(-std::log(ewaldErrorTol));
// vector<float> c12(numParticles); // double mx = boxVectors[0][0]/force.getCutoffDistance();
// vector<float> q(numParticles); // double my = boxVectors[1][1]/force.getCutoffDistance();
// vector<char> symbol; // double mz = boxVectors[2][2]/force.getCutoffDistance();
// vector<vector<int> > exclusionList(numParticles); // double pi = 3.1415926535897932385;
// for (int i = 0; i < numParticles; i++) { // int kmaxx = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
// double charge, radius, depth; // int kmaxy = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
// force.getParticleParameters(i, charge, radius, depth); // int kmaxz = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
// particle[i] = i; // if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
// q[i] = (float) charge; // if (kmaxx%2 == 0)
// c6[i] = (float) (4*depth*pow(radius, 6.0)); // kmaxx++;
// c12[i] = (float) (4*depth*pow(radius, 12.0)); // if (kmaxy%2 == 0)
// exclusionList[i].push_back(i); // kmaxy++;
// } // if (kmaxz%2 == 0)
// for (int i = 0; i < (int)exclusions.size(); i++) { // kmaxz++;
// exclusionList[exclusions[i].first].push_back(exclusions[i].second); // gpuSetEwaldParameters(gpu, (float) alpha, kmaxx, kmaxy, kmaxz);
// exclusionList[exclusions[i].second].push_back(exclusions[i].first); // method = EWALD;
// }
// Vec3 boxVectors[3];
// system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
// gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
// OpenCLNonbondedMethod method = NO_CUTOFF;
// if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
// gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), force.getReactionFieldDielectric());
// 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);
//// printf("%d %d\n", gridSizeX, (int) (kmaxx*std::sqrt(-std::log(ewaldErrorTol))));
//// gridSizeX = 0.02*mx*std::pow(-std::log(1.5*ewaldErrorTol), 3);
//// gridSizeY = 0.02*my*std::pow(-std::log(1.5*ewaldErrorTol), 3);
//// gridSizeZ = 0.02*mz*std::pow(-std::log(1.5*ewaldErrorTol), 3);
//// double scale = 0.698*std::pow(ewaldErrorTol, -0.312);
//// double scale = 0.713*std::pow(ewaldErrorTol, -0.261);
//// printf("%f\n", scale);
// gridSizeX = mx*NonbondedForce::PMEscale;
// gridSizeY = my*NonbondedForce::PMEscale;
// gridSizeZ = mz*NonbondedForce::PMEscale;
//// printf("%d %d %d\n", gridSizeX, gridSizeY, gridSizeZ);
//// gridSizeX = mx*scale;
//// gridSizeY = my*scale;
//// gridSizeZ = mz*scale;
// gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
// method = PARTICLE_MESH_EWALD;
// }
// } // }
// data.nonbondedMethod = method; // else {
// gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method); // int gridSizeX = -0.5*kmaxx*std::log(ewaldErrorTol);
// // int gridSizeY = -0.5*kmaxy*std::log(ewaldErrorTol);
// // Compute the Ewald self energy. // int gridSizeZ = -0.5*kmaxz*std::log(ewaldErrorTol);
// // gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
// data.ewaldSelfEnergy = 0.0; // method = PARTICLE_MESH_EWALD;
// 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++)
// data.ewaldSelfEnergy -= selfEnergyScale*q[i]*q[i];
// } // }
// } // }
// // data.nonbondedMethod = method;
// // Initialize 1-4 nonbonded interactions. // gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
// cl.getNonbondedUtilties().addInteraction(useCutoff, usePeriodic, force.getCutoffDistance(), exclusionList);
// { cl.getNonbondedUtilties().addParameter("sigmaEpsilon", "float2", 8, sigmaEpsilon->getDeviceBuffer());
// int numExceptions = exceptions.size();
// vector<int> particle1(numExceptions); // Compute the Ewald self energy.
// vector<int> particle2(numExceptions);
// vector<float> c6(numExceptions); ewaldSelfEnergy = 0.0;
// vector<float> c12(numExceptions); if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
// vector<float> q1(numExceptions); // double selfEnergyScale = gpu->sim.epsfac*gpu->sim.alphaEwald/std::sqrt(PI);
// vector<float> q2(numExceptions); // for (int i = 0; i < numParticles; i++)
// for (int i = 0; i < numExceptions; i++) { // ewaldSelfEnergy -= selfEnergyScale*q[i]*q[i];
// double charge, sig, eps; }
// force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
// c6[i] = (float) (4*eps*pow(sig, 6.0)); // Initialize 1-4 nonbonded interactions.
// c12[i] = (float) (4*eps*pow(sig, 12.0));
// q1[i] = (float) charge; {
// q2[i] = 1.0f; int numExceptions = exceptions.size();
// } vector<int> particle1(numExceptions);
vector<int> particle2(numExceptions);
vector<float> c6(numExceptions);
vector<float> c12(numExceptions);
vector<float> q1(numExceptions);
vector<float> q2(numExceptions);
for (int i = 0; i < numExceptions; i++) {
double charge, sig, eps;
force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
c6[i] = (float) (4*eps*pow(sig, 6.0));
c12[i] = (float) (4*eps*pow(sig, 12.0));
q1[i] = (float) charge;
q2[i] = 1.0f;
}
// gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2); // gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
// } }
//} }
//
//void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) { void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
// if (data.primaryKernel == this) cl.getNonbondedUtilties().computeInteractions();
// calcForces(context, data); }
//}
// double OpenCLCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
//double OpenCLCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) { executeForces(context);
// if (data.primaryKernel == this) return ewaldSelfEnergy;
// return calcEnergy(context, data, system); }
// return 0.0;
//}
//
//OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() { //OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
//} //}
// //
......
platforms/opencl/src/OpenCLKernels.h
View file @ 3976ad91
...@@ -294,40 +294,40 @@ private: ...@@ -294,40 +294,40 @@ private:
cl::Kernel kernel; cl::Kernel kernel;
}; };
///** /**
// * This kernel is invoked by NonbondedForce to calculate the forces acting on the system. * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
// */ */
//class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel { class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
//public: public:
// OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform), cl(cl), system(system) { OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform), cl(cl) {
// } }
// ~OpenCLCalcNonbondedForceKernel(); ~OpenCLCalcNonbondedForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the NonbondedForce this kernel will be used for * @param force the NonbondedForce this kernel will be used for
// */ */
// void initialize(const System& system, const NonbondedForce& force); void initialize(const System& system, const NonbondedForce& force);
// /** /**
// * Execute the kernel to calculate the forces. * Execute the kernel to calculate the forces.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// */ */
// void executeForces(ContextImpl& context); void executeForces(ContextImpl& context);
// /** /**
// * Execute the kernel to calculate the energy. * Execute the kernel to calculate the energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @return the potential energy due to the NonbondedForce * @return the potential energy due to the NonbondedForce
// */ */
// double executeEnergy(ContextImpl& context); double executeEnergy(ContextImpl& context);
//private: private:
// OpenCLContext& cl; OpenCLContext& cl;
// int numParticles; OpenCLArray<mm_float2>* sigmaEpsilon;
// System& system; double ewaldSelfEnergy;
//}; };
//
///** ///**
// * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system. // * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
// */ // */
......
platforms/opencl/src/OpenCLNonbondedUtilities.cpp 0 → 100644
View file @ 3976ad91
/* -------------------------------------------------------------------------- *
* 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) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "OpenCLNonbondedUtilities.h"
#include "OpenCLArray.h"
#include <map>
using namespace OpenMM;
using namespace std;
OpenCLNonbondedUtilities::OpenCLNonbondedUtilities(OpenCLContext& context) : context(context), cutoff(-1.0), useCutoff(false),
numForceBuffers(0), tiles(NULL), exclusionIndex(NULL), exclusions(NULL) {
}
OpenCLNonbondedUtilities::~OpenCLNonbondedUtilities() {
if (tiles != NULL)
delete tiles;
if (exclusionIndex != NULL)
delete exclusionIndex;
if (exclusions != NULL)
delete exclusions;
}
void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList) {
if (cutoff != -1.0) {
if (usesCutoff != useCutoff)
throw OpenMMException("All Forces must agree on whether to use a cutoff");
if (usesPeriodic != usePeriodic)
throw OpenMMException("All Forces must agree on whether to use periodic boundary conditions");
if (cutoffDistance != cutoff)
throw OpenMMException("All Forces must use the same cutoff distance");
bool sameExclusions = (exclusionList.size() == atomExclusions.size());
for (int i = 0; i < exclusionList.size() && sameExclusions; i++) {
if (exclusionList[i].size() != atomExclusions[i].size())
sameExclusions = false;
for (int j = 0; j < exclusionList[i].size(); j++)
if (exclusionList[i][j] != atomExclusions[i][j])
sameExclusions = false;
}
if (!sameExclusions)
throw OpenMMException("All Forces must have identical exceptions");
}
else {
useCutoff = usesCutoff;
usePeriodic = usesPeriodic;
cutoff = cutoffDistance;
atomExclusions = exclusionList;
}
}
void OpenCLNonbondedUtilities::addParameter(const string& name, const string& type, int size, cl::Buffer& buffer) {
parameters.push_back(ParameterInfo(name, type, size, buffer));
}
void OpenCLNonbondedUtilities::initialize(const System& system) {
if (cutoff == -1.0)
return; // There are no nonbonded interactions in the System.
// Create the list of tiles.
int numAtomBlocks = context.getNumAtomBlocks();
int numTiles = numAtomBlocks*(numAtomBlocks+1)/2;
tiles = new OpenCLArray<cl_uint>(context, numTiles, "tiles");
vector<cl_uint> tileVec(tiles->getSize());
unsigned int count = 0;
for (unsigned int y = 0; y < numAtomBlocks; y++)
for (unsigned int x = y; x < numAtomBlocks; x++)
tileVec[count++] = (x << 17) | (y << 2);
// Decide how many force buffers to use.
bool forceBufferPerAtomBlock = false;
numForceBuffers = context.getNumThreadBlocks()*OpenCLContext::ThreadBlockSize/OpenCLContext::TileSize;
if (numForceBuffers >= numAtomBlocks) {
// For small systems, it is more efficient to have one force buffer per block of 32 atoms instead of one per warp.
forceBufferPerAtomBlock = true;
numForceBuffers = numAtomBlocks;
}
// Create kernels.
cl::Program forceProgram = context.createProgram(context.loadSourceFromFile("nonbonded.cl"));
forceKernel = cl::Kernel(forceProgram, "computeNonbonded");
// Mark which tiles have exclusions.
for (int atom1 = 0; atom1 < (int) atomExclusions.size(); ++atom1) {
int x = atom1/OpenCLContext::TileSize;
for (int j = 0; j < (int) atomExclusions[atom1].size(); ++j) {
int atom2 = atomExclusions[atom1][j];
int y = atom2/OpenCLContext::TileSize;
int index = (x > y ? x+y*numAtomBlocks-y*(y+1)/2 : y+x*numAtomBlocks-x*(x+1)/2);
tileVec[index] |= 1;
}
}
if (context.getPaddedNumAtoms() > context.getNumAtoms()) {
int lastTile = context.getNumAtoms()/OpenCLContext::TileSize;
for (int i = 0; i < numTiles; ++i) {
int x = tileVec[i]>>17;
int y = (tileVec[i]>>2)&0x7FFF;
if (x == lastTile || y == lastTile)
tileVec[i] |= 1;
}
}
// Build a list of indices for the tiles with exclusions.
exclusionIndex = new OpenCLArray<cl_uint>(context, numTiles, "exclusionIndex");
vector<cl_uint> exclusionIndexVec(exclusionIndex->getSize());
int numWithExclusions = 0;
for (int i = 0; i < numTiles; ++i)
if ((tileVec[i]&1) == 1)
exclusionIndexVec[i] = (numWithExclusions++)*OpenCLContext::TileSize;
// Record the exclusion data.
exclusions = new OpenCLArray<cl_uint>(context, numWithExclusions*OpenCLContext::TileSize, "exclusions");
vector<cl_uint> exclusionVec(exclusions->getSize());
for (int i = 0; i < exclusions->getSize(); ++i)
exclusionVec[i] = 0xFFFFFFFF;
for (int atom1 = 0; atom1 < (int) atomExclusions.size(); ++atom1) {
int x = atom1/OpenCLContext::TileSize;
int offset1 = atom1-x*OpenCLContext::TileSize;
for (int j = 0; j < (int) atomExclusions[atom1].size(); ++j) {
int atom2 = atomExclusions[atom1][j];
int y = atom2/OpenCLContext::TileSize;
int offset2 = atom2-y*OpenCLContext::TileSize;
if (x > y) {
int tile = x+y*numAtomBlocks-y*(y+1)/2;
exclusionVec[exclusionIndexVec[tile]+offset1] &= 0xFFFFFFFF-(1<<offset2);
}
else {
int tile = y+x*numAtomBlocks-x*(x+1)/2;
exclusionVec[exclusionIndexVec[tile]+offset2] &= 0xFFFFFFFF-(1<<offset1);
}
}
}
// Mark all interactions that involve a padding atom as being excluded.
for (int atom1 = context.getNumAtoms(); atom1 < context.getPaddedNumAtoms(); ++atom1) {
int x = atom1/OpenCLContext::TileSize;
int offset1 = atom1-x*OpenCLContext::TileSize;
for (int atom2 = 0; atom2 < context.getPaddedNumAtoms(); ++atom2) {
int y = atom2/OpenCLContext::TileSize;
int offset2 = atom2-y*OpenCLContext::TileSize;
if (x >= y) {
int tile = x+y*numAtomBlocks-y*(y+1)/2;
exclusionVec[exclusionIndexVec[tile]+offset1] &= 0xFFFFFFFF-(1<<offset2);
}
if (y >= x) {
int tile = y+x*numAtomBlocks-x*(x+1)/2;
exclusionVec[exclusionIndexVec[tile]+offset2] &= 0xFFFFFFFF-(1<<offset1);
}
}
}
atomExclusions.clear(); // We won't use this again, so free the memory it used
tiles->upload(tileVec);
exclusions->upload(exclusionVec);
exclusionIndex->upload(exclusionIndexVec);
Vec3 boxVectors[3];
system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
periodicBoxSize = (mm_float4) {(float) boxVectors[0][0], (float) boxVectors[1][1], (float) boxVectors[2][2], 0.0f};
}
void OpenCLNonbondedUtilities::prepareInteractions() {
hasComputedInteractions = false;
if (!useCutoff)
return;
// TODO compute the neighbor list
}
void OpenCLNonbondedUtilities::computeInteractions() {
if (hasComputedInteractions)
return;
hasComputedInteractions = true;
forceKernel.setArg<cl_int>(0, tiles->getSize());
forceKernel.setArg<cl_int>(1, context.getPaddedNumAtoms());
forceKernel.setArg<cl_float>(2, cutoff*cutoff);
forceKernel.setArg<mm_float4>(3, periodicBoxSize);
forceKernel.setArg<cl::Buffer>(4, context.getForceBuffers().getDeviceBuffer());
forceKernel.setArg<cl::Buffer>(5, context.getEnergyBuffer().getDeviceBuffer());
forceKernel.setArg<cl::Buffer>(6, context.getPosq().getDeviceBuffer());
forceKernel.setArg<cl::Buffer>(7, tiles->getDeviceBuffer());
forceKernel.setArg<cl::Buffer>(8, exclusions->getDeviceBuffer());
forceKernel.setArg<cl::Buffer>(9, exclusionIndex->getDeviceBuffer());
forceKernel.setArg(10, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
forceKernel.setArg(11, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
for (int i = 0; i < (int) parameters.size(); i++) {
forceKernel.setArg<cl::Buffer>(i*2+12, *parameters[i].buffer);
forceKernel.setArg(i*2+13, OpenCLContext::ThreadBlockSize*parameters[i].size, NULL);
}
context.executeKernel(forceKernel, tiles->getSize()*OpenCLContext::TileSize);
}
platforms/opencl/src/OpenCLNonbondedUtilities.h 0 → 100644
View file @ 3976ad91
#ifndef OPENMM_OPENCLNONBONDEDUTILITIES_H_
#define OPENMM_OPENCLNONBONDEDUTILITIES_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) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "OpenCLContext.h"
#include "openmm/System.h"
#include <string>
#include <vector>
namespace OpenMM {
/**
* This class implements features that are used by several different force. It provides
* a generic interface for calculating nonbonded interactions.
*/
class OpenCLNonbondedUtilities {
public:
OpenCLNonbondedUtilities(OpenCLContext& context);
~OpenCLNonbondedUtilities();
/**
* Add a nonbonded interaction.
*
* @param usesCutoff specifies whether a cutoff should be applied to this interaction
* @param usesPeriodic specifies whether periodic boundary conditions should be applied to this interaction
* @param cutoffDistance the cutoff distance for this interaction (ignored if usesCutoff is false)
* @param exclusionList for each atom, specifies the list of other atoms whose interactions should be excluded
*/
void addInteraction(bool usesCutoff, bool usesPeriodic, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList);
/**
* Add a per-atom parameter that interactions may depend on.
*
* @param name the name of the parameter
* @param type the data type of the parameter
* @param size the size of the parameter in bytes
* @param buffer the buffer containing the parameter values
*/
void addParameter(const std::string& name, const std::string& type, int size, cl::Buffer& buffer);
/**
* Initialize this object in preparation for a simulation.
*/
void initialize(const System& system);
/**
* Get the number of force buffers required for nonbonded forces.
*/
int getNumForceBuffers() {
return numForceBuffers;
}
/**
* Prepare to compute interactions. This updates the neighbor list.
*/
void prepareInteractions();
/**
* Compute the nonbonded interactions. This will only be executed once after each call to
* prepareInteractions(). Additional calls return immediately without doing anything.
*/
void computeInteractions();
private:
class ParameterInfo;
OpenCLContext& context;
cl::Kernel forceKernel;
OpenCLArray<cl_uint>* tiles;
OpenCLArray<cl_uint>* exclusionIndex;
OpenCLArray<cl_uint>* exclusions;
std::vector<std::vector<int> > atomExclusions;
std::vector<ParameterInfo> parameters;
double cutoff;
bool useCutoff, usePeriodic, hasComputedInteractions;
int numForceBuffers;
mm_float4 periodicBoxSize;
};
class OpenCLNonbondedUtilities::ParameterInfo {
public:
ParameterInfo(const std::string& name, const std::string& type, int size, cl::Buffer& buffer) :
name(name), type(type), size(size), buffer(&buffer) {
}
std::string name;
std::string type;
int size;
cl::Buffer* buffer;
};
} // namespace OpenMM
#endif /*OPENMM_OPENCLNONBONDEDUTILITIES_H_*/
platforms/opencl/src/OpenCLPlatform.cpp
View file @ 3976ad91
...@@ -51,7 +51,7 @@ OpenCLPlatform::OpenCLPlatform() { ...@@ -51,7 +51,7 @@ OpenCLPlatform::OpenCLPlatform() {
registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory); registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory); registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory); registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
// registerKernelFactory(CalcNonbondedForceKernel::Name(), factory); registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
// registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory); // registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
// registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory); // registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory); registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
......
platforms/opencl/src/kernels/nonbonded.cl 0 → 100644
View file @ 3976ad91
const unsigned int TileSize = 32;
const float EpsilonFactor = 138.935485f;
/**
* Compute nonbonded interactions.
*/
__kernel void computeNonbonded(int numTiles, int paddedNumAtoms, float cutoffSquared, float4 periodicBoxSize,
__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global unsigned int* tiles,
__global unsigned int* exclusions, __global unsigned int* exclusionIndices, __local float4* local_posq, __local float4* local_force,
__global float2* sigmaEpsilon, __local float2* local_sigmaEpsilon) {
unsigned int totalWarps = get_global_size(0)/TileSize;
unsigned int warp = get_global_id(0)/TileSize;
unsigned int pos = warp*numTiles/totalWarps;
unsigned int end = (warp+1)*numTiles/totalWarps;
float energy = 0.0f;
#ifdef USE_CUTOFF
float3* tempBuffer = (float3*) &sA[cSim.nonbond_threads_per_block];
#endif
unsigned int lasty = 0xFFFFFFFF;
while (pos < end) {
// Extract tile coordinates from appropriate work unit
unsigned int x = tiles[pos];
unsigned int y = ((x >> 2) & 0x7fff)*TileSize;
bool hasExclusions = (x & 0x1);
x = (x>>17)*TileSize;
float4 apos; // Local atom x, y, z, q
float4 af = 0.0f; // Local atom fx, fy, fz
unsigned int tgx = get_local_id(0) & (TileSize-1);
unsigned int tbx = get_local_id(0) - tgx;
unsigned int tj = tgx;
unsigned int i = x + tgx;
apos = posq[i];
float2 a = sigmaEpsilon[i];
if (x == y) {
// Handle diagonals uniquely at 50% efficiency
// Read fixed atom data into registers and GRF
local_posq[get_local_id(0)] = apos;
local_sigmaEpsilon[get_local_id(0)] = a;
barrier(CLK_LOCAL_MEM_FENCE);
apos.w *= EpsilonFactor;
unsigned int xi = x/TileSize;
unsigned int tile = xi+xi*paddedNumAtoms/TileSize-xi*(xi+1)/2;
unsigned int excl = exclusions[exclusionIndices[tile]+tgx];
for (unsigned int j = 0; j < TileSize; j++) {
bool isExcluded = !(excl & 0x1);
float4 delta = (float4) (local_posq[tbx+j].xyz - apos.xyz, 0.0f);
#ifdef USE_PERIODIC
delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
#endif
float r2 = delta.x * delta.x + delta.y * delta.y + delta.z * delta.z;
float invR = 1.0f / sqrt(r2);
float sig = a.x + local_sigmaEpsilon[tbx+j].x;
float sig2 = invR * sig;
sig2 *= sig2;
float sig6 = sig2 * sig2 * sig2;
float eps = a.y * local_sigmaEpsilon[tbx+j].y;
float dEdR = eps * (12.0f * sig6 - 6.0f) * sig6;
float tempEnergy = eps * (sig6 - 1.0f) * sig6;
#ifdef USE_CUTOFF
dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * cSim.reactionFieldK * r2);
tempEnergy += apos.w * local_posq[tbx+j].w * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
#else
dEdR += apos.w * local_posq[tbx+j].w * invR;
tempEnergy += apos.w * local_posq[tbx+j].w * invR;
#endif
dEdR *= invR * invR;
#ifdef USE_CUTOFF
if (isExcluded || r2 > cutoffSquared) {
#else
if (isExcluded) {
#endif
dEdR = 0.0f;
tempEnergy = 0.0f;
}
energy += 0.5f*tempEnergy;
delta.xyz *= dEdR;
af.xyz -= delta.xyz;
excl >>= 1;
}
// Write results
float4 of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
unsigned int offset = x + tgx + warp*paddedNumAtoms;
of = forceBuffers[offset];
of.xyz += af.xyz;
forceBuffers[offset] = of;
#else
of.xyz = af.xyz;
of.w = 0.0f;
unsigned int offset = x + tgx + (x/TileSize) * paddedNumAtoms;
forceBuffers[offset] = of;
#endif
}
else {
// 100% utilization
// Read fixed atom data into registers and GRF
if (lasty != y) {
unsigned int j = y + tgx;
float2 temp1 = sigmaEpsilon[j];
local_posq[get_local_id(0)] = posq[j];
local_sigmaEpsilon[get_local_id(0)] = sigmaEpsilon[j];
}
local_force[get_local_id(0)] = 0.0f;
barrier(CLK_LOCAL_MEM_FENCE);
apos.w *= EpsilonFactor;
#ifdef USE_CUTOFF
unsigned int flags = cSim.pInteractionFlag[pos];
if (!hasExclusions && flags != 0xFFFFFFFF) {
if (flags == 0) {
// No interactions in this tile.
}
else {
// Compute only a subset of the interactions in this tile.
for (unsigned int j = 0; j < TileSize; j++) {
if ((flags&(1<<j)) != 0) {
bool isExcluded = false;
float4 delta = (float4) (local_posq[tbx+j].xyz - apos.xyz, 0.0f);
#ifdef USE_PERIODIC
delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
#endif
float r2 = delta.x * delta.x + delta.y * delta.y + delta.z * delta.z;
float invR = 1.0f / sqrt(r2);
float sig = a.x + local_sigmaEpsilon[tbx+j].x;
float sig2 = invR * sig;
sig2 *= sig2;
float sig6 = sig2 * sig2 * sig2;
float eps = a.y * local_sigmaEpsilon[tbx+j].Y;
float dEdR = eps * (12.0f * sig6 - 6.0f) * sig6;
float tempEnergy = eps * (sig6 - 1.0f) * sig6;
#ifdef USE_CUTOFF
dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * cSim.reactionFieldK * r2);
tempEnergy += apos.w * local_posq[tbx+j].w * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
#else
dEdR += apos.w * local_posq[tbx+j].w * invR;
tempEnergy += apos.w * local_posq[tbx+j].w * invR;
#endif
dEdR *= invR * invR;
#ifdef USE_CUTOFF
if (r2 > cutoffSquared) {
dEdR = 0.0f;
tempEnergy = 0.0f;
}
#endif
energy += tempEnergy;
delta.xyz *= dEdR;
af.xyz -= delta.xyz;
tempBuffer[get_local_id(0)] = delta;
// Sum the forces on atom j.
if (tgx % 2 == 0)
tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+1].xyz;
barrier(CLK_LOCAL_MEM_FENCE);
if (tgx % 4 == 0)
tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+2].xyz;
barrier(CLK_LOCAL_MEM_FENCE);
if (tgx % 8 == 0)
tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+4].xyz;
barrier(CLK_LOCAL_MEM_FENCE);
if (tgx % 16 == 0)
tempBuffer[get_local_id(0)].xyz += tempBuffer[get_local_id(0)+8].xyz;
barrier(CLK_LOCAL_MEM_FENCE);
if (tgx == 0)
local_force[tbx+j].xyz += tempBuffer[get_local_id(0)].xyz + tempBuffer[get_local_id(0)+16].xyz;
barrier(CLK_LOCAL_MEM_FENCE);
}
}
}
}
else // bExclusion
#endif
{
// Read fixed atom data into registers and GRF
unsigned int xi = x/TileSize;
unsigned int yi = y/TileSize;
unsigned int tile = xi+yi*paddedNumAtoms/TileSize-yi*(yi+1)/2;
unsigned int excl = exclusions[exclusionIndices[tile]+tgx];
excl = (excl >> tgx) | (excl << (TileSize - tgx));
for (unsigned int j = 0; j < TileSize; j++) {
bool isExcluded = !(excl & 0x1);
float4 delta = (float4) (local_posq[tbx+tj].xyz - apos.xyz, 0.0f);
#ifdef USE_PERIODIC
delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
#endif
float r2 = delta.x * delta.x + delta.y * delta.y + delta.z * delta.z;
float invR = 1.0f / sqrt(r2);
float sig = a.x + local_sigmaEpsilon[tbx+j].x;
float sig2 = invR * sig;
sig2 *= sig2;
float sig6 = sig2 * sig2 * sig2;
float eps = a.y * local_sigmaEpsilon[tbx+j].y;
float dEdR = eps * (12.0f * sig6 - 6.0f) * sig6;
float tempEnergy = eps * (sig6 - 1.0f) * sig6;
#ifdef USE_CUTOFF
dEdR += apos.w * local_posq[tbx+j].w * (invR - 2.0f * cSim.reactionFieldK * r2);
tempEnergy += apos.w * local_posq[tbx+j].w * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
#else
dEdR += apos.w * local_posq[tbx+j].w * invR;
tempEnergy += apos.w * local_posq[tbx+j].w * invR;
#endif
dEdR *= invR * invR;
#ifdef USE_CUTOFF
if (isExcluded || r2 > cutoffSquared) {
#else
if (isExcluded) {
#endif
dEdR = 0.0f;
tempEnergy = 0.0f;
}
energy += tempEnergy;
delta.xyz *= dEdR;
af.xyz -= delta.xyz;
local_force[tbx+tj].xyz += delta.xyz;
barrier(CLK_LOCAL_MEM_FENCE);
excl >>= 1;
tj = (tj + 1) & (TileSize - 1);
}
}
// Write results
float4 of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
unsigned int offset = x + tgx + warp*paddedNumAtoms;
of = forceBuffers[offset];
of.xyz += af.xyz;
forceBuffers[offset] = of;
offset = y + tgx + warp*paddedNumAtoms;
of = forceBuffers[offset];
of.xyz += local_foce[get_local_id(0)].xyz;
forceBuffers[offset] = of;
#else
of.xyz = af.xyz;
of.w = 0.0f;
unsigned int offset = x + tgx + (y/TileSize) * paddedNumAtoms;
forceBuffers[offset] = of;
of = local_force[get_local_id(0)];
offset = y + tgx + (x/TileSize) * paddedNumAtoms;
forceBuffers[offset] = of;
#endif
lasty = y;
}
pos++;
}
energyBuffer[get_global_id(0)] += energy;
}
platforms/opencl/tests/TestOpenCLNonbondedForce.cpp 0 → 100644
View file @ 3976ad91
/* -------------------------------------------------------------------------- *
* 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 all the different force terms in the reference implementation of NonbondedForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
//#include "kernels/gputypes.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 testCoulomb() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0.5, 1, 0);
forceField->addParticle(-1.5, 1, 0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double force = 138.935485*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
}
void testLJ() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0, 1.2, 1);
forceField->addParticle(0, 1.4, 2);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.3/2.0;
double eps = SQRT_TWO;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
}
void testExclusionsAnd14() {
OpenCLPlatform platform;
System system;
LangevinIntegrator integrator(0.0, 0.1, 0.01);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < 5; ++i) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
nonbonded->setParticleParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
nonbonded->setParticleParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.5/r;
double eps = 1.0;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
nonbonded->setParticleParameters(0, 2, 1.5, 0);
nonbonded->setParticleParameters(i, 2, 1.5, 0);
nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
Context context2(system, integrator, platform);
context2.setPositions(positions);
state = context2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r);
energy = 138.935485*4/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testCutoff() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(1.0, 1, 0);
forceField->addParticle(1.0, 1, 0);
forceField->addParticle(1.0, 1, 0);
forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
const double eps = 50.0;
forceField->setReactionFieldDielectric(eps);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 2, 0);
positions[2] = Vec3(0, 3, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0);
const double force2 = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
const double energy1 = 138.935485*(1.0)*(0.5+krf*4.0-crf);
const double energy2 = 138.935485*(1.0)*(1.0+krf*1.0-crf);
ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
}
void testCutoff14() {
OpenCLPlatform platform;
System system;
LangevinIntegrator integrator(0.0, 0.1, 0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
for (int i = 0; i < 5; ++i) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
const double cutoff = 3.5;
nonbonded->setCutoffDistance(cutoff);
const double eps = 30.0;
nonbonded->setReactionFieldDielectric(eps);
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(5);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(2, 0, 0);
positions[3] = Vec3(3, 0, 0);
positions[4] = Vec3(4, 0, 0);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
nonbonded->setParticleParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
nonbonded->setParticleParameters(j, 0, 1.5, 0);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double r = positions[i][0];
double x = 1.5/r;
double e = 1.0;
double force = 4.0*e*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*e*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3 || r > cutoff) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
const double q = 0.7;
nonbonded->setParticleParameters(0, q, 1.5, 0);
nonbonded->setParticleParameters(i, q, 1.5, 0);
nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
force = 138.935485*q*q*(1.0/(r*r)-2.0*krf*r);
energy = 138.935485*q*q*(1.0/r+krf*r*r-crf);
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3 || r > cutoff) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testPeriodic() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addException(0, 1, 0.0, 1.0, 0.0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
positions[2] = Vec3(3, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testLargeSystem() {
const int numMolecules = 600;
const int numParticles = numMolecules*2;
const double cutoff = 2.0;
const double boxSize = 20.0;
const double tol = 1e-3;
OpenCLPlatform cl;
ReferencePlatform reference;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
HarmonicBondForce* bonds = new HarmonicBondForce();
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
nonbonded->addParticle(-1.0, 0.2, 0.1);
nonbonded->addParticle(1.0, 0.1, 0.1);
}
else {
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(1.0, 0.1, 0.2);
}
positions[2*i] = Vec3(boxSize*genrand_real2(), boxSize*genrand_real2(), boxSize*genrand_real2());
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
bonds->addBond(2*i, 2*i+1, 1.0, 0.1);
nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
}
// Try with cutoffs but not periodic boundary conditions, and make sure the cl and Reference
// platforms agree.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
system.addForce(bonds);
Context clContext(system, integrator, cl);
Context referenceContext(system, integrator, reference);
clContext.setPositions(positions);
clContext.setVelocities(velocities);
referenceContext.setPositions(positions);
referenceContext.setVelocities(velocities);
State clState = clContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(clState.getPositions()[i], referenceState.getPositions()[i], tol);
ASSERT_EQUAL_VEC(clState.getVelocities()[i], referenceState.getVelocities()[i], tol);
ASSERT_EQUAL_VEC(clState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(clState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// Now do the same thing with periodic boundary conditions.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
clContext.reinitialize();
referenceContext.reinitialize();
clContext.setPositions(positions);
clContext.setVelocities(velocities);
referenceContext.setPositions(positions);
referenceContext.setVelocities(velocities);
clState = clContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_TOL(fmod(clState.getPositions()[i][0]-referenceState.getPositions()[i][0], boxSize), 0, tol);
ASSERT_EQUAL_TOL(fmod(clState.getPositions()[i][1]-referenceState.getPositions()[i][1], boxSize), 0, tol);
ASSERT_EQUAL_TOL(fmod(clState.getPositions()[i][2]-referenceState.getPositions()[i][2], boxSize), 0, tol);
ASSERT_EQUAL_VEC(clState.getVelocities()[i], referenceState.getVelocities()[i], tol);
ASSERT_EQUAL_VEC(clState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(clState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
}
//void testBlockInteractions(bool periodic) {
// const int blockSize = 32;
// const int numBlocks = 100;
// const int numParticles = blockSize*numBlocks;
// const double cutoff = 1.0;
// const double boxSize = (periodic ? 5.1 : 1.1);
// OpenCLPlatform cl;
// System system;
// VerletIntegrator integrator(0.01);
// NonbondedForce* nonbonded = new NonbondedForce();
// vector<Vec3> positions(numParticles);
// init_gen_rand(0);
// for (int i = 0; i < numParticles; i++) {
// system.addParticle(1.0);
// nonbonded->addParticle(1.0, 0.2, 0.2);
// positions[i] = Vec3(boxSize*(3*genrand_real2()-1), boxSize*(3*genrand_real2()-1), boxSize*(3*genrand_real2()-1));
// }
// nonbonded->setNonbondedMethod(periodic ? NonbondedForce::CutoffPeriodic : NonbondedForce::CutoffNonPeriodic);
// nonbonded->setCutoffDistance(cutoff);
// system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
// system.addForce(nonbonded);
// Context context(system, integrator, cl);
// context.setPositions(positions);
// State state = context.getState(State::Positions | State::Velocities | State::Forces);
// ContextImpl* contextImpl = *reinterpret_cast<ContextImpl**>(&context);
// OpenCLPlatform::PlatformData& data = *static_cast<OpenCLPlatform::PlatformData*>(contextImpl->getPlatformData());
//
// // Verify that the bounds of each block were calculated correctly.
//
// data.gpu->psPosq4->Download();
// data.gpu->psGridBoundingBox->Download();
// data.gpu->psGridCenter->Download();
// for (int i = 0; i < numBlocks; i++) {
// float4 gridSize = (*data.gpu->psGridBoundingBox)[i];
// float4 center = (*data.gpu->psGridCenter)[i];
// if (periodic) {
// ASSERT(gridSize.x < 0.5*boxSize);
// ASSERT(gridSize.y < 0.5*boxSize);
// ASSERT(gridSize.z < 0.5*boxSize);
// }
// float minx = 0.0, maxx = 0.0, miny = 0.0, maxy = 0.0, minz = 0.0, maxz = 0.0, radius = 0.0;
// for (int j = 0; j < blockSize; j++) {
// float4 pos = (*data.gpu->psPosq4)[i*blockSize+j];
// float dx = pos.x-center.x;
// float dy = pos.y-center.y;
// float dz = pos.z-center.z;
// if (periodic) {
// dx -= (float)(floor(0.5+dx/boxSize)*boxSize);
// dy -= (float)(floor(0.5+dy/boxSize)*boxSize);
// dz -= (float)(floor(0.5+dz/boxSize)*boxSize);
// }
// ASSERT(abs(dx) < gridSize.x+TOL);
// ASSERT(abs(dy) < gridSize.y+TOL);
// ASSERT(abs(dz) < gridSize.z+TOL);
// minx = min(minx, dx);
// maxx = max(maxx, dx);
// miny = min(miny, dy);
// maxy = max(maxy, dy);
// minz = min(minz, dz);
// maxz = max(maxz, dz);
// }
// ASSERT_EQUAL_TOL(-minx, gridSize.x, TOL);
// ASSERT_EQUAL_TOL(maxx, gridSize.x, TOL);
// ASSERT_EQUAL_TOL(-miny, gridSize.y, TOL);
// ASSERT_EQUAL_TOL(maxy, gridSize.y, TOL);
// ASSERT_EQUAL_TOL(-minz, gridSize.z, TOL);
// ASSERT_EQUAL_TOL(maxz, gridSize.z, TOL);
// }
//
// // Verify that interactions were identified correctly.
//
// data.gpu->psInteractionCount->Download();
// int numWithInteractions = (*data.gpu->psInteractionCount)[0];
// vector<bool> hasInteractions(data.gpu->sim.workUnits, false);
// data.gpu->psInteractingWorkUnit->Download();
// data.gpu->psInteractionFlag->Download();
// const unsigned int atoms = data.gpu->sim.paddedNumberOfAtoms;
// const unsigned int grid = data.gpu->grid;
// const unsigned int dim = (atoms+(grid-1))/grid;
// for (int i = 0; i < numWithInteractions; i++) {
// unsigned int workUnit = (*data.gpu->psInteractingWorkUnit)[i];
// unsigned int x = (workUnit >> 17);
// unsigned int y = ((workUnit >> 2) & 0x7fff);
// int tile = (x > y ? x+y*dim-y*(y+1)/2 : y+x*dim-x*(x+1)/2);
// hasInteractions[tile] = true;
//
// // Make sure this tile really should have been flagged based on bounding volumes.
//
// float4 gridSize1 = (*data.gpu->psGridBoundingBox)[x];
// float4 gridSize2 = (*data.gpu->psGridBoundingBox)[y];
// float4 center1 = (*data.gpu->psGridCenter)[x];
// float4 center2 = (*data.gpu->psGridCenter)[y];
// float dx = center1.x-center2.x;
// float dy = center1.y-center2.y;
// float dz = center1.z-center2.z;
// if (periodic) {
// dx -= (float)(floor(0.5+dx/boxSize)*boxSize);
// dy -= (float)(floor(0.5+dy/boxSize)*boxSize);
// dz -= (float)(floor(0.5+dz/boxSize)*boxSize);
// }
// dx = max(0.0f, abs(dx)-gridSize1.x-gridSize2.x);
// dy = max(0.0f, abs(dy)-gridSize1.y-gridSize2.y);
// dz = max(0.0f, abs(dz)-gridSize1.z-gridSize2.z);
// ASSERT(sqrt(dx*dx+dy*dy+dz*dz) < cutoff+TOL);
//
// // Check the interaction flags.
//
// unsigned int flags = (*data.gpu->psInteractionFlag)[i];
// for (int atom2 = 0; atom2 < 32; atom2++) {
// if ((flags & 1) == 0) {
// float4 pos2 = (*data.gpu->psPosq4)[y*blockSize+atom2];
// for (int atom1 = 0; atom1 < blockSize; ++atom1) {
// float4 pos1 = (*data.gpu->psPosq4)[x*blockSize+atom1];
// float dx = pos2.x-pos1.x;
// float dy = pos2.y-pos1.y;
// float dz = pos2.z-pos1.z;
// if (periodic) {
// dx -= (float)(floor(0.5+dx/boxSize)*boxSize);
// dy -= (float)(floor(0.5+dy/boxSize)*boxSize);
// dz -= (float)(floor(0.5+dz/boxSize)*boxSize);
// }
// if (dx*dx+dy*dy+dz*dz < cutoff*cutoff) {
// dx = pos2.x-center1.x;
// dy = pos2.y-center1.y;
// dz = pos2.z-center1.z;
// dx = max(0.0f, abs(dx)-gridSize1.x);
// dy = max(0.0f, abs(dy)-gridSize1.y);
// dz = max(0.0f, abs(dz)-gridSize1.z);
// }
// ASSERT(dx*dx+dy*dy+dz*dz > cutoff*cutoff);
// }
// }
// flags >>= 1;
// }
// }
//
// // Check the tiles that did not have interactions to make sure all atoms are beyond the cutoff.
//
// data.gpu->psWorkUnit->Download();
// for (int i = 0; i < (int)hasInteractions.size(); i++)
// if (!hasInteractions[i]) {
// unsigned int workUnit = (*data.gpu->psWorkUnit)[i];
// unsigned int x = (workUnit >> 17);
// unsigned int y = ((workUnit >> 2) & 0x7fff);
// for (int atom1 = 0; atom1 < blockSize; ++atom1) {
// float4 pos1 = (*data.gpu->psPosq4)[x*blockSize+atom1];
// for (int atom2 = 0; atom2 < blockSize; ++atom2) {
// float4 pos2 = (*data.gpu->psPosq4)[y*blockSize+atom2];
// float dx = pos1.x-pos2.x;
// float dy = pos1.y-pos2.y;
// float dz = pos1.z-pos2.z;
// if (periodic) {
// dx -= (float)(floor(0.5+dx/boxSize)*boxSize);
// dy -= (float)(floor(0.5+dy/boxSize)*boxSize);
// dz -= (float)(floor(0.5+dz/boxSize)*boxSize);
// }
// ASSERT(dx*dx+dy*dy+dz*dz > cutoff*cutoff);
// }
// }
// }
//}
int main() {
try {
testCoulomb();
testLJ();
// testExclusionsAnd14();
// testCutoff();
// testCutoff14();
// testPeriodic();
// testLargeSystem();
// testBlockInteractions(false);
// testBlockInteractions(true);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
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