Skip to content

GitLab

Commit 0634e3cc authored by Peter Eastman's avatar Peter Eastman
Browse files

Created OpenCL implementations of HarmonicBondForce, HarmonicAngleForce, and PeriodicTorsionForce.

parent 9532c446
Hide whitespace changes
Inline Side-by-side
Showing with 550 additions and 155 deletions
platforms/opencl/src/OpenCLContext.h
View file @ 0634e3cc
......@@ -55,6 +55,9 @@ typedef struct {
typedef struct {
cl_int x, y, z, w;
} mm_int4;
typedef struct {
cl_int s0, s1, s2, s3, s4, s5, s6, s7;
} mm_int8;
/**
* This class contains the information associated with a Context by the OpenCL Platform.
......
platforms/opencl/src/OpenCLKernelFactory.cpp
View file @ 0634e3cc
......@@ -39,10 +39,10 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
return new OpenCLUpdateStateDataKernel(name, platform, data);
if (name == CalcHarmonicBondForceKernel::Name())
return new OpenCLCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
// if (name == CalcHarmonicAngleForceKernel::Name())
// return new OpenCLCalcHarmonicAngleForceKernel(name, platform, data, context.getSystem());
// if (name == CalcPeriodicTorsionForceKernel::Name())
// return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, data, context.getSystem());
if (name == CalcHarmonicAngleForceKernel::Name())
return new OpenCLCalcHarmonicAngleForceKernel(name, platform, data, context.getSystem());
if (name == CalcPeriodicTorsionForceKernel::Name())
return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, data, context.getSystem());
// if (name == CalcRBTorsionForceKernel::Name())
// return new OpenCLCalcRBTorsionForceKernel(name, platform, data, context.getSystem());
// if (name == CalcNonbondedForceKernel::Name())
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 0634e3cc
......@@ -137,13 +137,6 @@ void OpenCLUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Ve
}
}
OpenCLCalcHarmonicBondForceKernel::~OpenCLCalcHarmonicBondForceKernel() {
if (params != NULL)
delete params;
if (indices != NULL)
delete indices;
}
class OpenCLBondForceInfo : public OpenCLForceInfo {
public:
OpenCLBondForceInfo(int requiredBuffers, const HarmonicBondForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
......@@ -170,10 +163,17 @@ private:
const HarmonicBondForce& force;
};
OpenCLCalcHarmonicBondForceKernel::~OpenCLCalcHarmonicBondForceKernel() {
if (params != NULL)
delete params;
if (indices != NULL)
delete indices;
}
void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
numBonds = force.getNumBonds();
params = new OpenCLArray<mm_float2>(*data.context, numBonds, "BondParams");
indices = new OpenCLArray<mm_int4>(*data.context, numBonds, "BondIndices");
params = new OpenCLArray<mm_float2>(*data.context, numBonds, "bondParams");
indices = new OpenCLArray<mm_int4>(*data.context, numBonds, "bondIndices");
vector<int> forceBufferCounter(system.getNumParticles(), 0);
vector<mm_float2> paramVector(numBonds);
vector<mm_int4> indicesVector(numBonds);
......@@ -211,77 +211,162 @@ double OpenCLCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
executeForces(context);
return 0.0;
}
//
//OpenCLCalcHarmonicAngleForceKernel::~OpenCLCalcHarmonicAngleForceKernel() {
//}
//
//void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
// if (data.primaryKernel == NULL)
// data.primaryKernel = this;
// data.hasAngles = true;
// numAngles = force.getNumAngles();
// const float RadiansToDegrees = (float) (180.0/3.14159265);
// vector<int> particle1(numAngles);
// vector<int> particle2(numAngles);
// vector<int> particle3(numAngles);
// vector<float> angle(numAngles);
// vector<float> k(numAngles);
// for (int i = 0; i < numAngles; i++) {
// double angleValue, kValue;
// force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
// angle[i] = (float) (angleValue*RadiansToDegrees);
// k[i] = (float) kValue;
// }
// gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
//}
//
//void OpenCLCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
// if (data.primaryKernel == this)
// calcForces(context, data);
//}
//
//double OpenCLCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
// if (data.primaryKernel == this)
// return calcEnergy(context, data, system);
// return 0.0;
//}
//
//OpenCLCalcPeriodicTorsionForceKernel::~OpenCLCalcPeriodicTorsionForceKernel() {
//}
//
//void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
// if (data.primaryKernel == NULL)
// data.primaryKernel = this;
// data.hasPeriodicTorsions = true;
// numTorsions = force.getNumTorsions();
// const float RadiansToDegrees = (float)(180.0/3.14159265);
// vector<int> particle1(numTorsions);
// vector<int> particle2(numTorsions);
// vector<int> particle3(numTorsions);
// vector<int> particle4(numTorsions);
// vector<float> k(numTorsions);
// vector<float> phase(numTorsions);
// vector<int> periodicity(numTorsions);
// for (int i = 0; i < numTorsions; i++) {
// double kValue, phaseValue;
// force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
// k[i] = (float) kValue;
// phase[i] = (float) (phaseValue*RadiansToDegrees);
// }
// gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
//}
//
//void OpenCLCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
// if (data.primaryKernel == this)
// calcForces(context, data);
//}
//
//double OpenCLCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
// if (data.primaryKernel == this)
// return calcEnergy(context, data, system);
// return 0.0;
//}
//
class OpenCLAngleForceInfo : public OpenCLForceInfo {
public:
OpenCLAngleForceInfo(int requiredBuffers, const HarmonicAngleForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
}
int getNumParticleGroups() {
return force.getNumAngles();
}
void getParticlesInGroup(int index, std::vector<int>& particles) {
int particle1, particle2, particle3;
double angle, k;
force.getAngleParameters(index, particle1, particle2, particle3, angle, k);
particles.resize(3);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3;
double angle1, angle2, k1, k2;
force.getAngleParameters(group1, particle1, particle2, particle3, angle1, k1);
force.getAngleParameters(group2, particle1, particle2, particle3, angle2, k2);
return (angle1 == angle2 && k1 == k2);
}
private:
const HarmonicAngleForce& force;
};
OpenCLCalcHarmonicAngleForceKernel::~OpenCLCalcHarmonicAngleForceKernel() {
if (params != NULL)
delete params;
if (indices != NULL)
delete indices;
}
void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
numAngles = force.getNumAngles();
params = new OpenCLArray<mm_float2>(*data.context, numAngles, "angleParams");
indices = new OpenCLArray<mm_int8>(*data.context, numAngles, "angleIndices");
vector<int> forceBufferCounter(system.getNumParticles(), 0);
vector<mm_float2> paramVector(numAngles);
vector<mm_int8> indicesVector(numAngles);
for (int i = 0; i < numAngles; i++) {
int particle1, particle2, particle3;
double angle, k;
force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
paramVector[i] = (mm_float2) {angle, k};
indicesVector[i] = (mm_int8) {particle1, particle2, particle3,
forceBufferCounter[particle1]++, forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, 0, 0};
}
params->upload(paramVector);
indices->upload(indicesVector);
int maxBuffers = 1;
for (int i = 0; i < forceBufferCounter.size(); i++) {
maxBuffers = max(maxBuffers, forceBufferCounter[i]);
}
data.context->addForce(new OpenCLAngleForceInfo(maxBuffers, force));
cl::Program program = data.context->createProgram(data.context->loadSourceFromFile("harmonicAngleForce.cl"));
kernel = cl::Kernel(program, "calcHarmonicAngleForce");
}
void OpenCLCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
kernel.setArg<cl_int>(0, data.context->getPaddedNumAtoms());
kernel.setArg<cl_int>(1, numAngles);
kernel.setArg<cl::Buffer>(2, data.context->getForceBuffers().getDeviceBuffer());
kernel.setArg<cl::Buffer>(3, data.context->getEnergyBuffer().getDeviceBuffer());
kernel.setArg<cl::Buffer>(4, data.context->getPosq().getDeviceBuffer());
kernel.setArg<cl::Buffer>(5, params->getDeviceBuffer());
kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
data.context->executeKernel(kernel, numAngles);
}
double OpenCLCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
executeForces(context);
return 0.0;
}
class OpenCLPeriodicTorsionForceInfo : public OpenCLForceInfo {
public:
OpenCLPeriodicTorsionForceInfo(int requiredBuffers, const PeriodicTorsionForce& force) : OpenCLForceInfo(requiredBuffers, false, 0.0), force(force) {
}
int getNumParticleGroups() {
return force.getNumTorsions();
}
void getParticlesInGroup(int index, std::vector<int>& particles) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(index, particle1, particle2, particle3, particle4, periodicity, phase, k);
particles.resize(4);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
particles[3] = particle4;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3, particle4, periodicity1, periodicity2;
double phase1, phase2, k1, k2;
force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, periodicity1, phase1, k1);
force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, periodicity2, phase2, k2);
return (periodicity1 == periodicity2 && phase1 == phase2 && k1 == k2);
}
private:
const PeriodicTorsionForce& force;
};
OpenCLCalcPeriodicTorsionForceKernel::~OpenCLCalcPeriodicTorsionForceKernel() {
if (params != NULL)
delete params;
if (indices != NULL)
delete indices;
}
void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
numTorsions = force.getNumTorsions();
params = new OpenCLArray<mm_float4>(*data.context, numTorsions, "periodicTorsionParams");
indices = new OpenCLArray<mm_int8>(*data.context, numTorsions, "periodicTorsionIndices");
vector<int> forceBufferCounter(system.getNumParticles(), 0);
vector<mm_float4> paramVector(numTorsions);
vector<mm_int8> indicesVector(numTorsions);
for (int i = 0; i < numTorsions; i++) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
paramVector[i] = (mm_float4) {k, phase, (float) periodicity};
indicesVector[i] = (mm_int8) {particle1, particle2, particle3, particle4,
forceBufferCounter[particle1]++, forceBufferCounter[particle2]++, forceBufferCounter[particle3]++, forceBufferCounter[particle4]++};
}
params->upload(paramVector);
indices->upload(indicesVector);
int maxBuffers = 1;
for (int i = 0; i < forceBufferCounter.size(); i++) {
maxBuffers = max(maxBuffers, forceBufferCounter[i]);
}
data.context->addForce(new OpenCLPeriodicTorsionForceInfo(maxBuffers, force));
cl::Program program = data.context->createProgram(data.context->loadSourceFromFile("periodicTorsionForce.cl"));
kernel = cl::Kernel(program, "calcPeriodicTorsionForce");
}
void OpenCLCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
kernel.setArg<cl_int>(0, data.context->getPaddedNumAtoms());
kernel.setArg<cl_int>(1, numTorsions);
kernel.setArg<cl::Buffer>(2, data.context->getForceBuffers().getDeviceBuffer());
kernel.setArg<cl::Buffer>(3, data.context->getEnergyBuffer().getDeviceBuffer());
kernel.setArg<cl::Buffer>(4, data.context->getPosq().getDeviceBuffer());
kernel.setArg<cl::Buffer>(5, params->getDeviceBuffer());
kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
data.context->executeKernel(kernel, numTorsions);
}
double OpenCLCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
executeForces(context);
return 0.0;
}
//OpenCLCalcRBTorsionForceKernel::~OpenCLCalcRBTorsionForceKernel() {
//}
//
......
platforms/opencl/src/OpenCLKernels.h
View file @ 0634e3cc
......@@ -183,74 +183,80 @@ private:
cl::Kernel kernel;
};
///**
// * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
// */
//class OpenCLCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
//public:
// OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcHarmonicAngleForceKernel(name, platform), data(data), system(system) {
// }
// ~OpenCLCalcHarmonicAngleForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the HarmonicAngleForce this kernel will be used for
// */
// void initialize(const System& system, const HarmonicAngleForce& force);
// /**
// * Execute the kernel to calculate the forces.
// *
// * @param context the context in which to execute this kernel
// */
// void executeForces(ContextImpl& context);
// /**
// * Execute the kernel to calculate the energy.
// *
// * @param context the context in which to execute this kernel
// * @return the potential energy due to the HarmonicAngleForce
// */
// double executeEnergy(ContextImpl& context);
//private:
// int numAngles;
// OpenCLPlatform::PlatformData& data;
// System& system;
//};
//
///**
// * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
// */
//class OpenCLCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
//public:
// OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcPeriodicTorsionForceKernel(name, platform), data(data), system(system) {
// }
// ~OpenCLCalcPeriodicTorsionForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the PeriodicTorsionForce this kernel will be used for
// */
// void initialize(const System& system, const PeriodicTorsionForce& force);
// /**
// * Execute the kernel to calculate the forces.
// *
// * @param context the context in which to execute this kernel
// */
// void executeForces(ContextImpl& context);
// /**
// * Execute the kernel to calculate the energy.
// *
// * @param context the context in which to execute this kernel
// * @return the potential energy due to the PeriodicTorsionForce
// */
// double executeEnergy(ContextImpl& context);
//private:
// int numTorsions;
// OpenCLPlatform::PlatformData& data;
// System& system;
//};
//
/**
* This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
*/
class OpenCLCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
public:
OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcHarmonicAngleForceKernel(name, platform), data(data), system(system) {
}
~OpenCLCalcHarmonicAngleForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the HarmonicAngleForce this kernel will be used for
*/
void initialize(const System& system, const HarmonicAngleForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(ContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the HarmonicAngleForce
*/
double executeEnergy(ContextImpl& context);
private:
int numAngles;
OpenCLPlatform::PlatformData& data;
System& system;
OpenCLArray<mm_float2>* params;
OpenCLArray<mm_int8>* indices;
cl::Kernel kernel;
};
/**
* This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class OpenCLCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
public:
OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, System& system) : CalcPeriodicTorsionForceKernel(name, platform), data(data), system(system) {
}
~OpenCLCalcPeriodicTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the PeriodicTorsionForce this kernel will be used for
*/
void initialize(const System& system, const PeriodicTorsionForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(ContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the PeriodicTorsionForce
*/
double executeEnergy(ContextImpl& context);
private:
int numTorsions;
OpenCLPlatform::PlatformData& data;
System& system;
OpenCLArray<mm_float4>* params;
OpenCLArray<mm_int8>* indices;
cl::Kernel kernel;
};
///**
// * This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
// */
......
platforms/opencl/src/OpenCLPlatform.cpp
View file @ 0634e3cc
......@@ -48,8 +48,8 @@ OpenCLPlatform::OpenCLPlatform() {
registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
registerKernelFactory(UpdateStateDataKernel::Name(), factory);
registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
// registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
// registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
// registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
// registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
// registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
......
platforms/opencl/src/kernels/harmonicAngleForce.cl 0 → 100644
View file @ 0634e3cc
/**
* Evaluate the forces due to harmonic angles.
*/
__kernel void calcHarmonicAngleForce(int numAtoms, int numAngles, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global float2* params, __global int8* indices) {
int index = get_global_id(0);
float energy = 0.0f;
while (index < numAngles) {
// Look up the data for this angle.
int8 atoms = indices[index];
float2 angleParams = params[index];
float4 a1 = posq[atoms.x];
float4 a2 = posq[atoms.y];
float4 a3 = posq[atoms.z];
// Compute the force.
float4 v0 = a2-a1;
float4 v1 = a2-a3;
float4 cp = cross(v0, v1);
float rp = dot(cp.xyz, cp.xyz);
rp = max(sqrt(rp), 1.0e-06f);
float r21 = dot(v0.xyz, v0.xyz);
float r23 = dot(v1.xyz, v1.xyz);
float dot = dot(v0.xyz, v1.xyz);
float cosine = dot/sqrt(r21*r23);
float deltaIdeal = acos(cosine)-angleParams.x;
energy += 0.5f*angleParams.y*deltaIdeal*deltaIdeal;
float dEdR = angleParams.y*deltaIdeal;
float4 c21 = cross(v0, cp)*(dEdR/(r21*rp));
float4 c23 = cross(cp, v1)*(dEdR/(r23*rp));
// Record the force on each of the three atoms.
unsigned int offsetA = atoms.s0+atoms.s3*numAtoms;
unsigned int offsetB = atoms.s1+atoms.s4*numAtoms;
unsigned int offsetC = atoms.s2+atoms.s5*numAtoms;
float4 forceA = forceBuffers[offsetA];
float4 forceB = forceBuffers[offsetB];
float4 forceC = forceBuffers[offsetC];
forceA.xyz += c21.xyz;
forceB.xyz -= c21.xyz+c23.xyz;
forceC.xyz += c23.xyz;
forceBuffers[offsetA] = forceA;
forceBuffers[offsetB] = forceB;
forceBuffers[offsetC] = forceC;
index += get_global_size(0);
}
energyBuffer[get_global_id(0)] += energy;
}
platforms/opencl/src/kernels/harmonicBondForce.cl
View file @ 0634e3cc
......@@ -6,15 +6,23 @@ __kernel void calcHarmonicBondForce(int numAtoms, int numBonds, __global float4*
int index = get_global_id(0);
float energy = 0.0f;
while (index < numBonds) {
// Look up the data for this bonds.
int4 atoms = indices[index];
float4 delta = posq[atoms.y]-posq[atoms.x];
float2 bondParams = params[index];
float r = sqrt(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
// Compute the force.
float r = sqrt(dot(delta.xyz, delta.xyz));
float deltaIdeal = r-bondParams.x;
energy += 0.5f * bondParams.y*deltaIdeal*deltaIdeal;
float dEdR = bondParams.y * deltaIdeal;
dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f;
delta.xyz *= dEdR;
// Record the force on each of the two atoms.
unsigned int offsetA = atoms.x+atoms.z*numAtoms;
unsigned int offsetB = atoms.y+atoms.w*numAtoms;
float4 forceA = forceBuffers[offsetA];
......
platforms/opencl/src/kernels/periodicTorsionForce.cl 0 → 100644
View file @ 0634e3cc
/**
* Evaluate the forces due to periodic torsions.
*/
__kernel void calcPeriodicTorsionForce(int numAtoms, int numTorsions, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global float4* params, __global int8* indices) {
int index = get_global_id(0);
float energy = 0.0f;
while (index < numTorsions) {
// Look up the data for this torsion.
int8 atoms = indices[index];
float4 torsionParams = params[index];
float4 a1 = posq[atoms.x];
float4 a2 = posq[atoms.y];
float4 a3 = posq[atoms.z];
float4 a4 = posq[atoms.w];
// Compute the force.
float4 v0 = (float4) (a1.xyz-a2.xyz, 0.0f);
float4 v1 = (float4) (a3.xyz-a2.xyz, 0.0f);
float4 v2 = (float4) (a3.xyz-a4.xyz, 0.0f);
float4 cp0 = cross(v0, v1);
float4 cp1 = cross(v1, v2);
float cosangle = dot(normalize(cp0), normalize(cp1));
cosangle = clamp(cosangle, -1.0f, 1.0f);
float dihedralAngle = acos(cosangle);
float deltaAngle = torsionParams.z*dihedralAngle-torsionParams.y;
energy += torsionParams.x*(1.0f+cos(deltaAngle));
float sinDeltaAngle = sin(deltaAngle);
float dEdAngle = -torsionParams.x*torsionParams.z*sinDeltaAngle;
float normCross1 = dot(cp0, cp0);
float normBC = sqrt(dot(v1, v1));
float normCross2 = dot(cp1, cp1);
float dp = 1.0f/normBC;
float4 ff = (float4) ((-dEdAngle*normBC)/normCross1, dot(v0, v1)*dp, dot(v2, v1)*dp, (dEdAngle*normBC)/normCross2);
float4 internalF0 = ff.x*cp0;
float4 internalF3 = ff.w*cp1;
float4 s = ff.y*internalF0 - ff.z*internalF3;
// Record the force on each of the four atoms.
unsigned int offsetA = atoms.s0+atoms.s4*numAtoms;
unsigned int offsetB = atoms.s1+atoms.s5*numAtoms;
unsigned int offsetC = atoms.s2+atoms.s6*numAtoms;
unsigned int offsetD = atoms.s3+atoms.s7*numAtoms;
float4 forceA = forceBuffers[offsetA];
float4 forceB = forceBuffers[offsetB];
float4 forceC = forceBuffers[offsetC];
float4 forceD = forceBuffers[offsetD];
forceA.xyz += internalF0.xyz;
forceB.xyz += s.xyz-internalF0.xyz;
forceC.xyz += -s.xyz-internalF3.xyz;
forceD.xyz += internalF3.xyz;
forceBuffers[offsetA] = forceA;
forceBuffers[offsetB] = forceB;
forceBuffers[offsetC] = forceC;
forceBuffers[offsetD] = forceD;
index += get_global_size(0);
}
energyBuffer[get_global_id(0)] += energy;
}
platforms/opencl/tests/TestOpenCLHarmonicAngleForce.cpp 0 → 100644
View file @ 0634e3cc
/* -------------------------------------------------------------------------- *
* 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 OpenCL implementation of HarmonicAngleForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testAngles() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
HarmonicAngleForce* forceField = new HarmonicAngleForce();
forceField->addAngle(0, 1, 2, PI_M/3, 1.1);
forceField->addAngle(1, 2, 3, PI_M/2, 1.2);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(2, 1, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double torque1 = 1.1*PI_M/6;
double torque2 = 1.2*PI_M/4;
ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
}
int main() {
try {
testAngles();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/opencl/tests/TestOpenCLPeriodicTorsionForce.cpp 0 → 100644
View file @ 0634e3cc
/* -------------------------------------------------------------------------- *
* 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 OpenCL implementation of PeriodicTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testPeriodicTorsions() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
}
int main() {
try {
testPeriodicTorsions();
}
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