Skip to content

GitLab

Commit 5c4a033f authored by Peter Eastman's avatar Peter Eastman
Browse files

Incorporating Cuda code into OpenMM

parent 170e493a
Show whitespace changes
Inline Side-by-side
Showing with 890 additions and 235 deletions
olla/include/Stream.h
View file @ 5c4a033f
......@@ -107,7 +107,7 @@ public:
* the values should be packed into a single array: all the values for the first element, followed by all the values
* for the next element, etc.
*/
void saveToArray(void* array);
void saveToArray(void* array) const;
/**
* Set every element of this stream to the same value.
*
......
olla/src/Stream.cpp
View file @ 5c4a033f
......@@ -75,7 +75,7 @@ void Stream::loadFromArray(const void* array) {
impl->loadFromArray(array);
}
void Stream::saveToArray(void* array) {
void Stream::saveToArray(void* array) const {
impl->saveToArray(array);
}
......
openmmapi/src/OpenMMContextImpl.cpp
View file @ 5c4a033f
......@@ -66,11 +66,6 @@ OpenMMContextImpl::OpenMMContextImpl(OpenMMContext& owner, System& system, Integ
else if (!platform->supportsKernels(kernelNames))
throw OpenMMException("Specified a Platform for an OpenMMContext which does not support all required kernels");
platform->contextCreated(*this);
positions = platform->createStream("atomPositions", system.getNumAtoms(), Stream::Double3, *this);
velocities = platform->createStream("atomVelocities", system.getNumAtoms(), Stream::Double3, *this);
forces = platform->createStream("atomForces", system.getNumAtoms(), Stream::Double3, *this);
double zero[] = {0.0, 0.0, 0.0};
velocities.fillWithValue(&zero);
kineticEnergyKernel = platform->createKernel(CalcKineticEnergyKernel::Name(), *this);
vector<double> masses(system.getNumAtoms());
for (size_t i = 0; i < masses.size(); ++i)
......@@ -79,6 +74,11 @@ OpenMMContextImpl::OpenMMContextImpl(OpenMMContext& owner, System& system, Integ
for (size_t i = 0; i < forceImpls.size(); ++i)
forceImpls[i]->initialize(*this);
integrator.initialize(*this);
positions = platform->createStream("atomPositions", system.getNumAtoms(), Stream::Double3, *this);
velocities = platform->createStream("atomVelocities", system.getNumAtoms(), Stream::Double3, *this);
forces = platform->createStream("atomForces", system.getNumAtoms(), Stream::Double3, *this);
double zero[] = {0.0, 0.0, 0.0};
velocities.fillWithValue(&zero);
}
OpenMMContextImpl::~OpenMMContextImpl() {
......
platforms/cuda/src/CudaKernelFactory.cpp
View file @ 5c4a033f
......@@ -38,19 +38,19 @@ using namespace OpenMM;
KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform& platform, OpenMMContextImpl& context) const {
_gpuContext* gpu = static_cast<_gpuContext*>(context.getPlatformData());
if (name == CalcStandardMMForceFieldKernel::Name())
return new CudaCalcStandardMMForceFieldKernel(name, platform, gpu);
return new CudaCalcStandardMMForceFieldKernel(name, platform, gpu, context.getSystem());
// if (name == CalcGBSAOBCForceFieldKernel::Name())
// return new CudaCalcGBSAOBCForceFieldKernel(name, platform);
// if (name == IntegrateVerletStepKernel::Name())
// return new CudaIntegrateVerletStepKernel(name, platform);
// if (name == IntegrateLangevinStepKernel::Name())
// return new CudaIntegrateLangevinStepKernel(name, platform);
if (name == IntegrateLangevinStepKernel::Name())
return new CudaIntegrateLangevinStepKernel(name, platform, gpu);
// if (name == IntegrateBrownianStepKernel::Name())
// return new CudaIntegrateBrownianStepKernel(name, platform);
// if (name == ApplyAndersenThermostatKernel::Name())
// return new CudaApplyAndersenThermostatKernel(name, platform);
// if (name == CalcKineticEnergyKernel::Name())
// return new CudaCalcKineticEnergyKernel(name, platform);
if (name == CalcKineticEnergyKernel::Name())
return new CudaCalcKineticEnergyKernel(name, platform);
// if (name == RemoveCMMotionKernel::Name())
// return new CudaRemoveCMMotionKernel(name, platform);
}
platforms/cuda/src/CudaKernels.cpp
View file @ 5c4a033f
......@@ -29,9 +29,14 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "OpenMMContext.h"
#include "CudaKernels.h"
#include "CudaStreamImpl.h"
#include "LangevinIntegrator.h"
#include "ReferencePlatform.h"
#include "kernels/gputypes.h"
#include "kernels/cudaKernels.h"
#include <cmath>
extern "C" int gpuSetConstants( gpuContext gpu );
......@@ -55,169 +60,166 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
numPeriodicTorsions = periodicTorsionIndices.size();
numRBTorsions = rbTorsionIndices.size();
num14 = bonded14Indices.size();
const float RadiansToDegrees = 180.0/3.14159265;
// Initialize bonds.
gpu->sim.bonds = numBonds;
CUDAStream<int4>* psBondID = new CUDAStream<int4>(numBonds, 1);
gpu->psBondID = psBondID;
gpu->sim.pBondID = psBondID->_pDevStream[0];
CUDAStream<float2>* psBondParameter = new CUDAStream<float2>(numBonds, 1);
gpu->psBondParameter = psBondParameter;
gpu->sim.pBondParameter = psBondParameter->_pDevStream[0];
for (int i = 0; i < numBonds; i++ ) {
psBondID->_pSysStream[0][i].x = bondIndices[i][0];
psBondID->_pSysStream[0][i].y = bondIndices[i][1];
psBondID->_pSysStream[0][i].z = gpu->pOutputBufferCounter[psBondID->_pSysStream[0][i].x]++;
psBondID->_pSysStream[0][i].w = gpu->pOutputBufferCounter[psBondID->_pSysStream[0][i].y]++;
psBondParameter->_pSysStream[0][i].x = (float) bondParameters[i][0];
psBondParameter->_pSysStream[0][i].y = (float) bondParameters[i][1];
{
vector<int> atom1(numBonds);
vector<int> atom2(numBonds);
vector<float> length(numBonds);
vector<float> k(numBonds);
for (int i = 0; i < numBonds; i++) {
atom1[i] = bondIndices[i][0];
atom2[i] = bondIndices[i][1];
length[i] = (float) bondParameters[i][0];
k[i] = (float) bondParameters[i][1];
}
gpuSetBondParameters(gpu, atom1, atom2, length, k);
}
psBondID->Upload();
psBondParameter->Upload();
// Initialize angles.
gpu->sim.bond_angles = numAngles;
CUDAStream<int4>* psBondAngleID1 = new CUDAStream<int4>(numAngles, 1);
gpu->psBondAngleID1 = psBondAngleID1;
gpu->sim.pBondAngleID1 = psBondAngleID1->_pDevStream[0];
CUDAStream<int2>* psBondAngleID2 = new CUDAStream<int2>(numAngles, 1);
gpu->psBondAngleID2 = psBondAngleID2;
gpu->sim.pBondAngleID2 = psBondAngleID2->_pDevStream[0];
CUDAStream<float2>* psBondAngleParameter = new CUDAStream<float2>(numAngles, 1);
gpu->psBondAngleParameter = psBondAngleParameter;
gpu->sim.pBondAngleParameter = psBondAngleParameter->_pDevStream[0];
{
vector<int> atom1(numAngles);
vector<int> atom2(numAngles);
vector<int> atom3(numAngles);
vector<float> angle(numAngles);
vector<float> k(numAngles);
for (int i = 0; i < numAngles; i++) {
psBondAngleID1->_pSysStream[0][i].x = angleIndices[i][0];
psBondAngleID1->_pSysStream[0][i].y = angleIndices[i][1];
psBondAngleID1->_pSysStream[0][i].z = angleIndices[i][2];
psBondAngleID1->_pSysStream[0][i].w = gpu->pOutputBufferCounter[psBondAngleID1->_pSysStream[0][i].x]++;
psBondAngleID2->_pSysStream[0][i].x = gpu->pOutputBufferCounter[psBondAngleID1->_pSysStream[0][i].y]++;
psBondAngleID2->_pSysStream[0][i].y = gpu->pOutputBufferCounter[psBondAngleID1->_pSysStream[0][i].z]++;
psBondAngleParameter->_pSysStream[0][i].x = (float) (angleParameters[i][0]*180.0/3.14159265);
psBondAngleParameter->_pSysStream[0][i].y = (float) angleParameters[i][1];
atom1[i] = angleIndices[i][0];
atom2[i] = angleIndices[i][1];
atom3[i] = angleIndices[i][2];
angle[i] = (float) (angleParameters[i][0]*RadiansToDegrees);
k[i] = (float) angleParameters[i][1];
}
gpuSetBondAngleParameters(gpu, atom1, atom2, atom3, angle, k);
}
psBondAngleID1->Upload();
psBondAngleID2->Upload();
psBondAngleParameter->Upload();
// Initialize periodic torsions.
gpu->sim.dihedrals = numPeriodicTorsions;
CUDAStream<int4>* psDihedralID1 = new CUDAStream<int4>(numPeriodicTorsions, 1);
gpu->psDihedralID1 = psDihedralID1;
gpu->sim.pDihedralID1 = psDihedralID1->_pDevStream[0];
CUDAStream<int4>* psDihedralID2 = new CUDAStream<int4>(numPeriodicTorsions, 1);
gpu->psDihedralID2 = psDihedralID2;
gpu->sim.pDihedralID2 = psDihedralID2->_pDevStream[0];
CUDAStream<float4>* psDihedralParameter = new CUDAStream<float4>(numPeriodicTorsions, 1);
gpu->psDihedralParameter = psDihedralParameter;
gpu->sim.pDihedralParameter = psDihedralParameter->_pDevStream[0];
{
vector<int> atom1(numPeriodicTorsions);
vector<int> atom2(numPeriodicTorsions);
vector<int> atom3(numPeriodicTorsions);
vector<int> atom4(numPeriodicTorsions);
vector<float> k(numPeriodicTorsions);
vector<float> phase(numPeriodicTorsions);
vector<int> periodicity(numPeriodicTorsions);
for (int i = 0; i < numPeriodicTorsions; i++) {
psDihedralID1->_pSysStream[0][i].x = periodicTorsionIndices[i][0];
psDihedralID1->_pSysStream[0][i].y = periodicTorsionIndices[i][1];
psDihedralID1->_pSysStream[0][i].z = periodicTorsionIndices[i][2];
psDihedralID1->_pSysStream[0][i].w = periodicTorsionIndices[i][3];
psDihedralID2->_pSysStream[0][i].x = gpu->pOutputBufferCounter[psDihedralID1->_pSysStream[0][i].x]++;
psDihedralID2->_pSysStream[0][i].y = gpu->pOutputBufferCounter[psDihedralID1->_pSysStream[0][i].y]++;
psDihedralID2->_pSysStream[0][i].z = gpu->pOutputBufferCounter[psDihedralID1->_pSysStream[0][i].z]++;
psDihedralID2->_pSysStream[0][i].w = gpu->pOutputBufferCounter[psDihedralID1->_pSysStream[0][i].w]++;
psDihedralParameter->_pSysStream[0][i].x = (float) periodicTorsionParameters[i][0];
psDihedralParameter->_pSysStream[0][i].y = (float) periodicTorsionParameters[i][1];
psDihedralParameter->_pSysStream[0][i].z = (float) periodicTorsionParameters[i][2];
psDihedralParameter->_pSysStream[0][i].w = 0.0f;
atom1[i] = periodicTorsionIndices[i][0];
atom2[i] = periodicTorsionIndices[i][1];
atom3[i] = periodicTorsionIndices[i][2];
atom4[i] = periodicTorsionIndices[i][3];
k[i] = (float) periodicTorsionParameters[i][0];
phase[i] = (float) (periodicTorsionParameters[i][1]*RadiansToDegrees);
periodicity[i] = (int) periodicTorsionParameters[i][2];
}
gpuSetDihedralParameters(gpu, atom1, atom2, atom3, atom4, k, phase, periodicity);
}
psDihedralID1->Upload();
psDihedralID2->Upload();
psDihedralParameter->Upload();
// Initialize Ryckaert-Bellemans torsions.
gpu->sim.rb_dihedrals = numRBTorsions;
CUDAStream<int4>* psRbDihedralID1 = new CUDAStream<int4>(numRBTorsions, 1);
gpu->psRbDihedralID1 = psRbDihedralID1;
gpu->sim.pRbDihedralID1 = psRbDihedralID1->_pDevStream[0];
CUDAStream<int4>* psRbDihedralID2 = new CUDAStream<int4>(numRBTorsions, 1);
gpu->psRbDihedralID2 = psRbDihedralID2;
gpu->sim.pRbDihedralID2 = psRbDihedralID2->_pDevStream[0];
CUDAStream<float4>* psRbDihedralParameter1 = new CUDAStream<float4>(numRBTorsions, 1);
gpu->psRbDihedralParameter1 = psRbDihedralParameter1;
gpu->sim.pRbDihedralParameter1 = psRbDihedralParameter1->_pDevStream[0];
CUDAStream<float2>* psRbDihedralParameter2 = new CUDAStream<float2>(numRBTorsions, 1);
gpu->psRbDihedralParameter2 = psRbDihedralParameter2;
gpu->sim.pRbDihedralParameter2 = psRbDihedralParameter2->_pDevStream[0];
{
vector<int> atom1(numRBTorsions);
vector<int> atom2(numRBTorsions);
vector<int> atom3(numRBTorsions);
vector<int> atom4(numRBTorsions);
vector<float> c0(numRBTorsions);
vector<float> c1(numRBTorsions);
vector<float> c2(numRBTorsions);
vector<float> c3(numRBTorsions);
vector<float> c4(numRBTorsions);
vector<float> c5(numRBTorsions);
for (int i = 0; i < numRBTorsions; i++) {
psRbDihedralID1->_pSysStream[0][i].x = rbTorsionIndices[i][0];
psRbDihedralID1->_pSysStream[0][i].y = rbTorsionIndices[i][1];
psRbDihedralID1->_pSysStream[0][i].z = rbTorsionIndices[i][2];
psRbDihedralID1->_pSysStream[0][i].w = rbTorsionIndices[i][3];
psRbDihedralID2->_pSysStream[0][i].x = gpu->pOutputBufferCounter[psRbDihedralID1->_pSysStream[0][i].x]++;
psRbDihedralID2->_pSysStream[0][i].y = gpu->pOutputBufferCounter[psRbDihedralID1->_pSysStream[0][i].y]++;
psRbDihedralID2->_pSysStream[0][i].z = gpu->pOutputBufferCounter[psRbDihedralID1->_pSysStream[0][i].z]++;
psRbDihedralID2->_pSysStream[0][i].w = gpu->pOutputBufferCounter[psRbDihedralID1->_pSysStream[0][i].w]++;
psRbDihedralParameter1->_pSysStream[0][i].x = (float) rbTorsionParameters[i][0];
psRbDihedralParameter1->_pSysStream[0][i].y = (float) rbTorsionParameters[i][1];
psRbDihedralParameter1->_pSysStream[0][i].z = (float) rbTorsionParameters[i][2];
psRbDihedralParameter1->_pSysStream[0][i].w = (float) rbTorsionParameters[i][3];
psRbDihedralParameter2->_pSysStream[0][i].x = (float) rbTorsionParameters[i][4];
psRbDihedralParameter2->_pSysStream[0][i].y = (float) rbTorsionParameters[i][5];
atom1[i] = rbTorsionIndices[i][0];
atom2[i] = rbTorsionIndices[i][1];
atom3[i] = rbTorsionIndices[i][2];
atom4[i] = rbTorsionIndices[i][3];
c0[i] = (float) rbTorsionParameters[i][0];
c1[i] = (float) rbTorsionParameters[i][1];
c2[i] = (float) rbTorsionParameters[i][2];
c3[i] = (float) rbTorsionParameters[i][3];
c4[i] = (float) rbTorsionParameters[i][4];
c5[i] = (float) rbTorsionParameters[i][5];
}
gpuSetRbDihedralParameters(gpu, atom1, atom2, atom3, atom4, c0, c1, c2, c3, c4, c5);
}
psRbDihedralID1->Upload();
psRbDihedralID2->Upload();
psRbDihedralParameter1->Upload();
psRbDihedralParameter2->Upload();
// Initialize nonbonded interactions.
{
vector<int> atom(numAtoms);
vector<float> c6(numAtoms);
vector<float> c12(numAtoms);
vector<float> q(numAtoms);
vector<char> symbol;
vector<vector<int> > exclusionList(numAtoms);
for (int i = 0; i < numAtoms; i++) {
gpu->psPosq4->_pSysStream[0][i].w = (float) nonbondedParameters[i][0];
gpu->psSigEps2->_pSysStream[0][i].x = (float) nonbondedParameters[i][1];
gpu->psSigEps2->_pSysStream[0][i].y = (float) nonbondedParameters[i][2];
atom[i] = i;
q[i] = (float) nonbondedParameters[i][0];
c6[i] = (float) (4*nonbondedParameters[i][2]*pow(nonbondedParameters[i][1], 6.0));
c12[i] = (float) (4*nonbondedParameters[i][2]*pow(nonbondedParameters[i][1], 12.0));
exclusionList[i] = vector<int>(exclusions[i].begin(), exclusions[i].end());
exclusionList[i].push_back(i);
}
gpuSetCoulombParameters(gpu, 138.935485f, atom, c6, c12, q, symbol, exclusionList);
}
gpu->psPosq4->Upload();
gpu->psSigEps2->Upload();
// Initialize 1-4 nonbonded interactions.
gpu->sim.LJ14s = num14;
CUDAStream<int4>* psLJ14ID = new CUDAStream<int4>(num14, 1);
gpu->psLJ14ID = psLJ14ID;
gpu->sim.pLJ14ID = psLJ14ID->_pDevStream[0];
CUDAStream<float4>* psLJ14Parameter = new CUDAStream<float4>(num14, 1);
gpu->psLJ14Parameter = psLJ14Parameter;
gpu->sim.pLJ14Parameter = psLJ14Parameter->_pDevStream[0];
double sqrtEps = std::sqrt(138.935485);
{
vector<int> atom1(num14);
vector<int> atom2(num14);
vector<float> c6(num14);
vector<float> c12(num14);
vector<float> q1(num14);
vector<float> q2(num14);
for (int i = 0; i < num14; i++) {
int atom1 = bonded14Indices[i][0];
int atom2 = bonded14Indices[i][1];
double atom1params[] = {0.5*nonbondedParameters[atom1][1], 2.0*sqrt(nonbondedParameters[atom1][2]), nonbondedParameters[atom1][0]*sqrtEps};
double atom2params[] = {0.5*nonbondedParameters[atom2][1], 2.0*sqrt(nonbondedParameters[atom2][2]), nonbondedParameters[atom2][0]*sqrtEps};
psLJ14ID->_pSysStream[0][i].x = atom1;
psLJ14ID->_pSysStream[0][i].y = atom2;
psLJ14ID->_pSysStream[0][i].z = gpu->pOutputBufferCounter[psLJ14ID->_pSysStream[0][i].x]++;
psLJ14ID->_pSysStream[0][i].w = gpu->pOutputBufferCounter[psLJ14ID->_pSysStream[0][i].y]++;
psLJ14Parameter->_pSysStream[0][i].x = (float) (atom1params[0]+atom2params[0]);
psLJ14Parameter->_pSysStream[0][i].y = (float) (lj14Scale*(atom1params[1]*atom2params[1]));
psLJ14Parameter->_pSysStream[0][i].z = (float) (coulomb14Scale*(atom1params[2]*atom2params[2]));
atom1[i] = bonded14Indices[i][0];
atom2[i] = bonded14Indices[i][1];
double sig = 0.5*(nonbondedParameters[atom1[i]][1]+nonbondedParameters[atom2[i]][1]);
double eps = sqrt(nonbondedParameters[atom1[i]][2]*nonbondedParameters[atom2[i]][2]);
c6[i] = (float) (4*eps*pow(sig, 6.0)*lj14Scale);
c12[i] = (float) (4*eps*pow(sig, 12.0)*lj14Scale);
q1[i] = (float) nonbondedParameters[atom1[i]][0];
q2[i] = (float) nonbondedParameters[atom2[i]][0];
}
gpuSetLJ14Parameters(gpu, 138.935485f, (float) coulomb14Scale, atom1, atom2, c6, c12, q1, q2);
}
psLJ14ID->Upload();
psLJ14Parameter->Upload();
// Initialize exclusions.
// TODO
// Finish initialization.
gpuBuildThreadBlockWorkList(gpu);
gpuBuildExclusionList(gpu);
gpuBuildOutputBuffers(gpu);
gpuSetConstants(gpu);
kClearBornForces(gpu);
kClearForces(gpu);
}
void CudaCalcStandardMMForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
kClearForces(gpu);
kCalculateCDLJForces(gpu);
kCalculateLocalForces(gpu);
kReduceBornSumAndForces(gpu);
}
double CudaCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions) {
return 0.0;
// We don't currently have GPU kernels to calculate energy, so instead we have the reference
// platform do it. This is VERY slow.
LangevinIntegrator integrator(0.0, 1.0, 0.0);
ReferencePlatform platform;
OpenMMContext context(system, integrator, platform);
double* posData = new double[positions.getSize()*3];
positions.saveToArray(posData);
vector<Vec3> pos(positions.getSize());
for (int i = 0; i < pos.size(); i++)
pos[i] = Vec3(posData[3*i], posData[3*i+1], posData[3*i+2]);
delete[] posData;
context.setPositions(pos);
return context.getState(State::Energy).getPotentialEnergy();
}
//CudaCalcGBSAOBCForceFieldKernel::~CudaCalcGBSAOBCForceFieldKernel() {
......@@ -240,16 +242,57 @@ double CudaCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions
//
//void CudaIntegrateVerletStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double stepSize) {
//}
//
//CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
//}
//
//void CudaIntegrateLangevinStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
// const vector<double>& constraintLengths) {
//}
//
//void CudaIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
//}
CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}
void CudaIntegrateLangevinStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
const vector<double>& constraintLengths) {
// Set masses.
vector<float> mass(masses.size());
for (int i = 0; i < (int) mass.size(); i++)
mass[i] = (float) masses[i];
gpuSetMass(gpu, mass);
// Set constraints.
int numConstraints = constraintLengths.size();
vector<int> atom1(numConstraints);
vector<int> atom2(numConstraints);
vector<float> distance(numConstraints);
vector<float> invMass1(numConstraints);
vector<float> invMass2(numConstraints);
for (int i = 0; i < numConstraints; i++) {
atom1[i] = constraintIndices[i][0];
atom2[i] = constraintIndices[i][1];
distance[i] = (float) constraintLengths[i];
invMass1[i] = 1.0f/mass[atom1[i]];
invMass2[i] = 1.0f/mass[atom2[i]];
}
gpuSetShakeParameters(gpu, atom1, atom2, distance, invMass1, invMass2);
gpuSetConstants(gpu);
prevStepSize = -1.0;
}
void CudaIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Initialize the GPU parameters.
double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
gpuSetIntegrationParameters(gpu, tau, stepSize, temperature);
gpuSetConstants(gpu);
kGenerateRandoms(gpu);
prevTemp = temperature;
prevFriction = friction;
prevStepSize = stepSize;
}
kUpdatePart1(gpu);
kApplyFirstShake(gpu);
kUpdatePart2(gpu);
kApplySecondShake(gpu);
}
//
//CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
//}
......@@ -269,12 +312,23 @@ double CudaCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions
//
//void CudaApplyAndersenThermostatKernel::execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize) {
//}
//
//void CudaCalcKineticEnergyKernel::initialize(const vector<double>& masses) {
//}
//
//double CudaCalcKineticEnergyKernel::execute(const Stream& velocities) {
//}
void CudaCalcKineticEnergyKernel::initialize(const vector<double>& masses) {
this->masses = masses;
}
double CudaCalcKineticEnergyKernel::execute(const Stream& velocities) {
// We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
// on the CPU.
double* v = new double[velocities.getSize()*3];
velocities.saveToArray(v);
double energy = 0.0;
for (size_t i = 0; i < masses.size(); ++i)
energy += masses[i]*(v[i*3]*v[i*3]+v[i*3+1]*v[i*3+1]+v[i*3+2]*v[i*3+2]);
delete v;
return 0.5*energy;
}
//
//void CudaRemoveCMMotionKernel::initialize(const vector<double>& masses) {
//}
......
platforms/cuda/src/CudaKernels.h
View file @ 5c4a033f
......@@ -34,6 +34,7 @@
#include "kernels.h"
#include "kernels/gpuTypes.h"
#include "System.h"
class CudaAndersenThermostat;
class CudaBrownianDynamics;
......@@ -48,7 +49,7 @@ namespace OpenMM {
*/
class CudaCalcStandardMMForceFieldKernel : public CalcStandardMMForceFieldKernel {
public:
CudaCalcStandardMMForceFieldKernel(std::string name, const Platform& platform, _gpuContext* gpu) : CalcStandardMMForceFieldKernel(name, platform), gpu(gpu) {
CudaCalcStandardMMForceFieldKernel(std::string name, const Platform& platform, _gpuContext* gpu, System& system) : CalcStandardMMForceFieldKernel(name, platform), gpu(gpu), system(system) {
}
~CudaCalcStandardMMForceFieldKernel();
/**
......@@ -96,8 +97,7 @@ public:
private:
_gpuContext* gpu;
int numAtoms, numBonds, numAngles, numPeriodicTorsions, numRBTorsions, num14;
// int **bondIndexArray, **angleIndexArray, **periodicTorsionIndexArray, **rbTorsionIndexArray, **exclusionArray, **bonded14IndexArray;
// RealOpenMM **bondParamArray, **angleParamArray, **periodicTorsionParamArray, **rbTorsionParamArray, **atomParamArray, **bonded14ParamArray;
System& system;
};
///**
......@@ -172,45 +172,39 @@ private:
// int numConstraints;
// double prevStepSize;
//};
//
///**
// * This kernel is invoked by LangevinIntegrator to take one time step.
// */
//class CudaIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
//public:
// CudaIntegrateLangevinStepKernel(std::string name, const Platform& platform) : IntegrateLangevinStepKernel(name, platform),
// dynamics(0), shake(0), masses(0), shakeParameters(0), constraintIndices(0) {
// }
// ~CudaIntegrateLangevinStepKernel();
// /**
// * Initialize the kernel, setting up all parameters related to integrator.
// *
// * @param masses the mass of each atom
// * @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
// * @param constraintLengths the required distance between each pair of constrained atoms
// */
// void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
// const std::vector<double>& constraintLengths);
// /**
// * Execute the kernel.
// *
// * @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
// * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
// * @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
// * @param temperature the temperature of the heat bath
// * @param friction the friction coefficient coupling the system to the heat bath
// * @param stepSize the integration step size
// */
// void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
//private:
// CudaStochasticDynamics* dynamics;
// CudaShakeAlgorithm* shake;
// RealOpenMM* masses;
// RealOpenMM** shakeParameters;
// int** constraintIndices;
// int numConstraints;
// double prevTemp, prevFriction, prevStepSize;
//};
/**
* This kernel is invoked by LangevinIntegrator to take one time step.
*/
class CudaIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public:
CudaIntegrateLangevinStepKernel(std::string name, const Platform& platform, _gpuContext* gpu) : IntegrateLangevinStepKernel(name, platform), gpu(gpu) {
}
~CudaIntegrateLangevinStepKernel();
/**
* Initialize the kernel, setting up all parameters related to integrator.
*
* @param masses the mass of each atom
* @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
* @param constraintLengths the required distance between each pair of constrained atoms
*/
void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
const std::vector<double>& constraintLengths);
/**
* Execute the kernel.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
* @param temperature the temperature of the heat bath
* @param friction the friction coefficient coupling the system to the heat bath
* @param stepSize the integration step size
*/
void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
private:
_gpuContext* gpu;
double prevTemp, prevFriction, prevStepSize;
};
//
///**
// * This kernel is invoked by BrownianIntegrator to take one time step.
......@@ -278,30 +272,30 @@ private:
// CudaAndersenThermostat* thermostat;
// RealOpenMM* masses;
//};
//
///**
// * This kernel is invoked to calculate the kinetic energy of the system.
// */
//class CudaCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
//public:
// CudaCalcKineticEnergyKernel(std::string name, const Platform& platform) : CalcKineticEnergyKernel(name, platform) {
// }
// /**
// * Initialize the kernel, setting up the atomic masses.
// *
// * @param masses the mass of each atom
// */
// void initialize(const std::vector<double>& masses);
// /**
// * Execute the kernel.
// *
// * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
// * @return the kinetic energy of the system
// */
// double execute(const Stream& velocities);
//private:
// std::vector<double> masses;
//};
/**
* This kernel is invoked to calculate the kinetic energy of the system.
*/
class CudaCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
public:
CudaCalcKineticEnergyKernel(std::string name, const Platform& platform) : CalcKineticEnergyKernel(name, platform) {
}
/**
* Initialize the kernel, setting up the atomic masses.
*
* @param masses the mass of each atom
*/
void initialize(const std::vector<double>& masses);
/**
* Execute the kernel.
*
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @return the kinetic energy of the system
*/
double execute(const Stream& velocities);
private:
std::vector<double> masses;
};
//
///**
// * This kernel is invoked to remove center of mass motion from the system.
......
platforms/cuda/src/CudaPlatform.cpp
View file @ 5c4a033f
......@@ -43,10 +43,10 @@ CudaPlatform::CudaPlatform() {
registerKernelFactory(CalcStandardMMForceFieldKernel::Name(), factory);
// registerKernelFactory(CalcGBSAOBCForceFieldKernel::Name(), factory);
// registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
// registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
// registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
// registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
// registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
// registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
}
......
platforms/cuda/src/CudaStreamImpl.h
View file @ 5c4a033f
......@@ -157,6 +157,7 @@ void CudaStreamImpl<T>::loadFromArray(const void* array) {
template <class T>
void CudaStreamImpl<T>::saveToArray(void* array) {
stream->Download();
float* data = reinterpret_cast<float*>(stream->_pSysData);
if (baseType == Stream::Float) {
float* arrayData = (float*) array;
......@@ -176,7 +177,6 @@ void CudaStreamImpl<T>::saveToArray(void* array) {
for (int j = 0; j < width; ++j)
arrayData[i*width+j] = (int) data[i*rowOffset+j];
}
stream->Download();
}
template <class T>
......
platforms/cuda/tests/TestCudaLangevinIntegrator.cpp 0 → 100644
View file @ 5c4a033f
/* -------------------------------------------------------------------------- *
* 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 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of LangevinIntegrator.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "CudaPlatform.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "LangevinIntegrator.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 testSingleBond() {
CudaPlatform platform;
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// This is simply a damped harmonic oscillator, so compare it to the analytical solution.
double freq = std::sqrt(1-0.05*0.05);
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities);
double time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
double expectedSpeed = -0.5*std::exp(-0.05*time)*(0.05*std::cos(freq*time)+freq*std::sin(freq*time));
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
integrator.step(1);
}
// Not set the friction to a tiny value and see if it conserves energy.
integrator.setFriction(5e-5);
context.setPositions(positions);
State state = context.getState(State::Energy);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Energy);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testTemperature() {
const int numAtoms = 8;
const double temp = 100.0;
CudaPlatform platform;
System system(numAtoms, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
for (int i = 0; i < numAtoms; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= 1000;
double expected = 0.5*numAtoms*3*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
}
void testConstraints() {
const int numAtoms = 8;
const int numConstraints = numAtoms-1;
const double temp = 100.0;
CudaPlatform platform;
System system(numAtoms, numConstraints);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numConstraints; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numConstraints; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-4);
}
integrator.step(1);
}
}
int main() {
try {
testSingleBond();
testTemperature();
testConstraints();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda/tests/TestCudaStandardMMForceField.cpp 0 → 100644
View file @ 5c4a033f
/* -------------------------------------------------------------------------- *
* 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 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 StandardMMForceField.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "CudaPlatform.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBonds() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 2, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 0.8);
forceField->setBondParameters(1, 1, 2, 1.2, 0.7);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
void testAngles() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 2, 0, 0);
forceField->setAngleParameters(0, 0, 1, 2, PI_M/3, 1.1);
forceField->setAngleParameters(1, 1, 2, 3, PI_M/2, 1.2);
system.addForce(forceField);
OpenMMContext 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);
}
void testPeriodicTorsions() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 1, 0);
forceField->setPeriodicTorsionParameters(0, 0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
OpenMMContext 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);
}
void testRBTorsions() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 0, 1);
forceField->setRBTorsionParameters(0, 0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
OpenMMContext 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, 1, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.1*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), 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);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.1*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
void testCoulomb() {
CudaPlatform platform;
System system(2, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0);
forceField->setAtomParameters(0, 0.5, 1, 0);
forceField->setAtomParameters(1, -1.5, 1, 0);
system.addForce(forceField);
OpenMMContext 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() {
CudaPlatform platform;
System system(2, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0);
forceField->setAtomParameters(0, 0, 1.2, 1);
forceField->setAtomParameters(1, 0, 1.4, 2);
system.addForce(forceField);
OpenMMContext 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() {
CudaPlatform platform;
System system(5, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1, 0);
forceField->setBondParameters(1, 1, 2, 1, 0);
forceField->setBondParameters(2, 2, 3, 1, 0);
forceField->setBondParameters(3, 3, 4, 1, 0);
system.addForce(forceField);
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) {
forceField->setAtomParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
forceField->setAtomParameters(0, 0, 1.5, 1);
forceField->setAtomParameters(i, 0, 1.5, 1);
positions[i] = Vec3(r, 0, 0);
OpenMMContext 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
forceField->setAtomParameters(0, 2, 1.5, 0);
forceField->setAtomParameters(i, 2, 1.5, 0);
OpenMMContext 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() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 0, 0, 0, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
system.addForce(forceField);
OpenMMContext 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 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 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() {
CudaPlatform platform;
System system(5, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1, 0);
forceField->setBondParameters(1, 1, 2, 1, 0);
forceField->setBondParameters(2, 2, 3, 1, 0);
forceField->setBondParameters(3, 3, 4, 1, 0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
const double cutoff = 3.5;
forceField->setCutoffDistance(cutoff);
system.addForce(forceField);
OpenMMContext 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
forceField->setAtomParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
forceField->setAtomParameters(j, 0, 1.5, 0);
forceField->setAtomParameters(i, 0, 1.5, 1);
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;
forceField->setAtomParameters(0, q, 1.5, 0);
forceField->setAtomParameters(i, q, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = 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);
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() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 1, 0, 0, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setBondParameters(0, 0, 1, 1.0, 0.0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffPeriodic);
const double cutoff = 2.0;
forceField->setCutoffDistance(cutoff);
forceField->setPeriodicBoxSize(4.0, 4.0, 4.0);
system.addForce(forceField);
OpenMMContext 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);
}
int main() {
try {
testBonds();
testAngles();
testPeriodicTorsions();
testRBTorsions();
testCoulomb();
testLJ();
testExclusionsAnd14();
// testCutoff();
// testCutoff14();
// testPeriodic();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceStandardMMForceField.cpp
View file @ 5c4a033f
......@@ -206,21 +206,19 @@ void testExclusionsAnd14() {
forceField->setBondParameters(3, 3, 4, 1, 0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(5);
const double r = 1.0;
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(r, 0, 0);
positions[2] = Vec3(r, 0, 0);
positions[3] = Vec3(r, 0, 0);
positions[4] = Vec3(r, 0, 0);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
forceField->setAtomParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
forceField->setAtomParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
forceField->setAtomParameters(0, 0, 1.5, 1);
forceField->setAtomParameters(i, 0, 1.5, 1);
positions[i] = Vec3(r, 0, 0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......
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