Skip to content

GitLab

Commit c88213f8 authored by peastman's avatar peastman
Browse files

Created CUDA and OpenCL implementations of LocalCoordinatesSite

parent 4be5daf1
Hide whitespace changes
Inline Side-by-side
Showing with 548 additions and 19 deletions
openmmapi/include/openmm/VirtualSite.h
View file @ c88213f8
......@@ -187,7 +187,7 @@ private:
* the x axis points from particle 1 toward the midpoint between particles 2 and 3.
*
* The z direction is computed as zdir = xdir x ydir. To ensure the axes are all orthogonal,
* ydir is then recomputed as ydir = zdir x xdir. All three axis vectors are then normalize, and
* ydir is then recomputed as ydir = zdir x xdir. All three axis vectors are then normalized, and
* the virtual site location is set to
*
* origin + x*xdir + y*ydir + z*zdir
......
platforms/cuda/include/CudaIntegrationUtilities.h
View file @ c88213f8
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009-2013 Stanford University and the Authors. *
* Portions copyright (c) 2009-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -150,6 +150,8 @@ private:
CudaArray* vsite3AvgWeights;
CudaArray* vsiteOutOfPlaneAtoms;
CudaArray* vsiteOutOfPlaneWeights;
CudaArray* vsiteLocalCoordsAtoms;
CudaArray* vsiteLocalCoordsParams;
int randomPos;
int lastSeed, numVsites;
struct ShakeCluster;
......
platforms/cuda/src/CudaIntegrationUtilities.cpp
View file @ c88213f8
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009-2013 Stanford University and the Authors. *
* Portions copyright (c) 2009-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -101,7 +101,7 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
ccmaReducedMass(NULL), ccmaAtomConstraints(NULL), ccmaNumAtomConstraints(NULL), ccmaConstraintMatrixColumn(NULL),
ccmaConstraintMatrixValue(NULL), ccmaDelta1(NULL), ccmaDelta2(NULL), ccmaConverged(NULL), ccmaConvergedMemory(NULL),
vsite2AvgAtoms(NULL), vsite2AvgWeights(NULL), vsite3AvgAtoms(NULL), vsite3AvgWeights(NULL),
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL) {
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsAtoms(NULL), vsiteLocalCoordsParams(NULL) {
// Create workspace arrays.
if (context.getUseDoublePrecision() || context.getUseMixedPrecision()) {
......@@ -553,6 +553,8 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
vector<double4> vsite3AvgWeightVec;
vector<int4> vsiteOutOfPlaneAtomVec;
vector<double4> vsiteOutOfPlaneWeightVec;
vector<int4> vsiteLocalCoordsAtomVec;
vector<double> vsiteLocalCoordsParamVec;
for (int i = 0; i < numAtoms; i++) {
if (system.isVirtualSite(i)) {
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -576,35 +578,65 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
vsiteOutOfPlaneAtomVec.push_back(make_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
vsiteOutOfPlaneWeightVec.push_back(make_double4(site.getWeight12(), site.getWeight13(), site.getWeightCross(), 0.0));
}
else if (dynamic_cast<const LocalCoordinatesSite*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
vsiteLocalCoordsAtomVec.push_back(make_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
Vec3 origin = site.getOriginWeights();
Vec3 x = site.getXWeights();
Vec3 y = site.getYWeights();
Vec3 pos = site.getLocalPosition();
vsiteLocalCoordsParamVec.push_back(origin[0]);
vsiteLocalCoordsParamVec.push_back(origin[1]);
vsiteLocalCoordsParamVec.push_back(origin[2]);
vsiteLocalCoordsParamVec.push_back(x[0]);
vsiteLocalCoordsParamVec.push_back(x[1]);
vsiteLocalCoordsParamVec.push_back(x[2]);
vsiteLocalCoordsParamVec.push_back(y[0]);
vsiteLocalCoordsParamVec.push_back(y[1]);
vsiteLocalCoordsParamVec.push_back(y[2]);
vsiteLocalCoordsParamVec.push_back(pos[0]);
vsiteLocalCoordsParamVec.push_back(pos[1]);
vsiteLocalCoordsParamVec.push_back(pos[2]);
}
}
}
int num2Avg = vsite2AvgAtomVec.size();
int num3Avg = vsite3AvgAtomVec.size();
int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
int numLocalCoords = vsiteLocalCoordsAtomVec.size();
vsite2AvgAtoms = CudaArray::create<int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
vsite3AvgAtoms = CudaArray::create<int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
vsiteOutOfPlaneAtoms = CudaArray::create<int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
vsiteLocalCoordsAtoms = CudaArray::create<int4>(context, max(1, numLocalCoords), "vsiteLocalCoordinatesAtoms");
if (num2Avg > 0)
vsite2AvgAtoms->upload(vsite2AvgAtomVec);
if (num3Avg > 0)
vsite3AvgAtoms->upload(vsite3AvgAtomVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneAtoms->upload(vsiteOutOfPlaneAtomVec);
if (numLocalCoords > 0)
vsiteLocalCoordsAtoms->upload(vsiteLocalCoordsAtomVec);
if (context.getUseDoublePrecision()) {
vsite2AvgWeights = CudaArray::create<double2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = CudaArray::create<double4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = CudaArray::create<double4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = CudaArray::create<double>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
if (num2Avg > 0)
vsite2AvgWeights->upload(vsite2AvgWeightVec);
if (num3Avg > 0)
vsite3AvgWeights->upload(vsite3AvgWeightVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
if (numLocalCoords > 0)
vsiteLocalCoordsParams->upload(vsiteLocalCoordsParamVec);
}
else {
vsite2AvgWeights = CudaArray::create<float2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = CudaArray::create<float4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = CudaArray::create<float4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = CudaArray::create<float>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
if (num2Avg > 0) {
vector<float2> floatWeights(num2Avg);
for (int i = 0; i < num2Avg; i++)
......@@ -623,6 +655,12 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
floatWeights[i] = make_float4((float) vsiteOutOfPlaneWeightVec[i].x, (float) vsiteOutOfPlaneWeightVec[i].y, (float) vsiteOutOfPlaneWeightVec[i].z, 0.0f);
vsiteOutOfPlaneWeights->upload(floatWeights);
}
if (numLocalCoords > 0) {
vector<float> floatParams(vsiteLocalCoordsParamVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsParamVec.size(); i++)
floatParams[i] = (float) vsiteLocalCoordsParamVec[i];
vsiteLocalCoordsParams->upload(floatParams);
}
}
// Create the kernels used by this class.
......@@ -633,6 +671,7 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
defines["NUM_2_AVERAGE"] = context.intToString(num2Avg);
defines["NUM_3_AVERAGE"] = context.intToString(num3Avg);
defines["NUM_OUT_OF_PLANE"] = context.intToString(numOutOfPlane);
defines["NUM_LOCAL_COORDS"] = context.intToString(numLocalCoords);
defines["PADDED_NUM_ATOMS"] = context.intToString(context.getPaddedNumAtoms());
CUmodule module = context.createModule(CudaKernelSources::vectorOps+CudaKernelSources::integrationUtilities, defines);
settlePosKernel = context.getKernel(module, "applySettleToPositions");
......@@ -647,7 +686,7 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
CHECK_RESULT2(cuEventCreate(&ccmaEvent, CU_EVENT_DISABLE_TIMING), "Error creating event for CCMA");
vsitePositionKernel = context.getKernel(module, "computeVirtualSites");
vsiteForceKernel = context.getKernel(module, "distributeVirtualSiteForces");
numVsites = num2Avg+num3Avg+numOutOfPlane;
numVsites = num2Avg+num3Avg+numOutOfPlane+numLocalCoords;
randomKernel = context.getKernel(module, "generateRandomNumbers");
timeShiftKernel = context.getKernel(module, "timeShiftVelocities");
}
......@@ -704,6 +743,10 @@ CudaIntegrationUtilities::~CudaIntegrationUtilities() {
delete vsiteOutOfPlaneAtoms;
if (vsiteOutOfPlaneWeights != NULL)
delete vsiteOutOfPlaneWeights;
if (vsiteLocalCoordsAtoms != NULL)
delete vsiteLocalCoordsAtoms;
if (vsiteLocalCoordsParams != NULL)
delete vsiteLocalCoordsParams;
}
void CudaIntegrationUtilities::applyConstraints(double tol) {
......@@ -779,7 +822,8 @@ void CudaIntegrationUtilities::computeVirtualSites() {
CUdeviceptr posCorrection = (context.getUseMixedPrecision() ? context.getPosqCorrection().getDevicePointer() : 0);
void* args[] = {&context.getPosq().getDevicePointer(), &posCorrection, &vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
&vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer()};
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer(),
&vsiteLocalCoordsAtoms->getDevicePointer(), &vsiteLocalCoordsParams->getDevicePointer()};
context.executeKernel(vsitePositionKernel, args, numVsites);
}
}
......@@ -790,7 +834,8 @@ void CudaIntegrationUtilities::distributeForcesFromVirtualSites() {
void* args[] = {&context.getPosq().getDevicePointer(), &posCorrection, &context.getForce().getDevicePointer(),
&vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
&vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer()};
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer(),
&vsiteLocalCoordsAtoms->getDevicePointer(), &vsiteLocalCoordsParams->getDevicePointer()};
context.executeKernel(vsiteForceKernel, args, numVsites);
}
}
......
platforms/cuda/src/kernels/integrationUtilities.cu
View file @ c88213f8
......@@ -685,7 +685,8 @@ extern "C" __global__ void updateCCMAAtomPositions(const int* __restrict__ numAt
*/
extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, real4* __restrict__ posqCorrection, const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights) {
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights,
const int4* __restrict__ localCoordsAtoms, const real* __restrict__ localCoordsParams) {
// Two particle average sites.
......@@ -733,6 +734,35 @@ extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, real4*
pos.z = pos1.z + v12.z*weights.x + v13.z*weights.y + cr.z*weights.z;
storePos(posq, posqCorrection, atoms.x, pos);
}
// Local coordinates sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_LOCAL_COORDS; index += blockDim.x*gridDim.x) {
int4 atoms = localCoordsAtoms[index];
const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed3 pos1 = make_mixed3(pos1_4.x, pos1_4.y, pos1_4.z);
mixed3 pos2 = make_mixed3(pos2_4.x, pos2_4.y, pos2_4.z);
mixed3 pos3 = make_mixed3(pos3_4.x, pos3_4.y, pos3_4.z);
mixed3 originWeights = make_mixed3(params[0], params[1], params[2]);
mixed3 xWeights = make_mixed3(params[3], params[4], params[5]);
mixed3 yWeights = make_mixed3(params[6], params[7], params[8]);
mixed3 localPosition = make_mixed3(params[9], params[10], params[11]);
mixed3 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed3 xdir = pos1*xWeights.x + pos2*xWeights.y + pos3*xWeights.z;
mixed3 ydir = pos1*yWeights.x + pos2*yWeights.y + pos3*yWeights.z;
mixed3 zdir = cross(xdir, ydir);
xdir *= rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
zdir *= rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
ydir = cross(zdir, xdir);
pos.x = origin.x + xdir.x*localPosition.x + ydir.x*localPosition.y + zdir.x*localPosition.z;
pos.y = origin.y + xdir.y*localPosition.x + ydir.y*localPosition.y + zdir.y*localPosition.z;
pos.z = origin.z + xdir.z*localPosition.x + ydir.z*localPosition.y + zdir.z*localPosition.z;
storePos(posq, posqCorrection, atoms.x, pos);
}
}
inline __device__ real3 loadForce(int index, long long* __restrict__ force) {
......@@ -753,7 +783,8 @@ inline __device__ void addForce(int index, long long* __restrict__ force, real3
extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__ posq, const real4* __restrict__ posqCorrection, long long* __restrict__ force,
const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights) {
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights,
const int4* __restrict__ localCoordsAtoms, const real* __restrict__ localCoordsParams) {
// Two particle average sites.
......@@ -797,6 +828,92 @@ extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__
addForce(atoms.z, force, fp2);
addForce(atoms.w, force, fp3);
}
// Local coordinates sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_LOCAL_COORDS; index += blockDim.x*gridDim.x) {
int4 atoms = localCoordsAtoms[index];
const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed3 pos1 = make_mixed3(pos1_4.x, pos1_4.y, pos1_4.z);
mixed3 pos2 = make_mixed3(pos2_4.x, pos2_4.y, pos2_4.z);
mixed3 pos3 = make_mixed3(pos3_4.x, pos3_4.y, pos3_4.z);
mixed3 originWeights = make_mixed3(params[0], params[1], params[2]);
mixed3 wx = make_mixed3(params[3], params[4], params[5]);
mixed3 wy = make_mixed3(params[6], params[7], params[8]);
mixed3 localPosition = make_mixed3(params[9], params[10], params[11]);
mixed3 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed3 xdir = pos1*wx.x + pos2*wx.y + pos3*wx.z;
mixed3 ydir = pos1*wy.x + pos2*wy.y + pos3*wy.z;
mixed3 zdir = cross(xdir, ydir);
mixed invNormXdir = rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed invNormZdir = rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed3 dx = xdir*invNormXdir;
mixed3 dz = zdir*invNormZdir;
mixed3 dy = cross(dz, dx);
// The derivatives for this case are very complicated. They were computed with SymPy then simplified by hand.
mixed t11 = (wx.x*ydir.x-wy.x*xdir.x)*invNormZdir;
mixed t12 = (wx.x*ydir.y-wy.x*xdir.y)*invNormZdir;
mixed t13 = (wx.x*ydir.z-wy.x*xdir.z)*invNormZdir;
mixed t21 = (wx.y*ydir.x-wy.y*xdir.x)*invNormZdir;
mixed t22 = (wx.y*ydir.y-wy.y*xdir.y)*invNormZdir;
mixed t23 = (wx.y*ydir.z-wy.y*xdir.z)*invNormZdir;
mixed t31 = (wx.z*ydir.x-wy.z*xdir.x)*invNormZdir;
mixed t32 = (wx.z*ydir.y-wy.z*xdir.y)*invNormZdir;
mixed t33 = (wx.z*ydir.z-wy.z*xdir.z)*invNormZdir;
mixed sx1 = t13*dz.y-t12*dz.z;
mixed sy1 = t11*dz.z-t13*dz.x;
mixed sz1 = t12*dz.x-t11*dz.y;
mixed sx2 = t23*dz.y-t22*dz.z;
mixed sy2 = t21*dz.z-t23*dz.x;
mixed sz2 = t22*dz.x-t21*dz.y;
mixed sx3 = t33*dz.y-t32*dz.z;
mixed sy3 = t31*dz.z-t33*dz.x;
mixed sz3 = t32*dz.x-t31*dz.y;
mixed3 wxScaled = wx*invNormXdir;
real3 f = loadForce(atoms.x, force);
mixed3 fp1 = localPosition*f.x;
mixed3 fp2 = localPosition*f.y;
mixed3 fp3 = localPosition*f.z;
real3 f1 = make_real3(0);
real3 f2 = make_real3(0);
real3 f3 = make_real3(0);
f1.x += fp1.x*wxScaled.x*(1-dx.x*dx.x) + fp1.z*(dz.x*sx1 ) + fp1.y*((-dx.x*dy.x )*wxScaled.x + dy.x*sx1 - dx.y*t12 - dx.z*t13) + f.x*originWeights.x;
f1.y += fp1.x*wxScaled.x*( -dx.x*dx.y) + fp1.z*(dz.x*sy1+t13) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.x + dy.x*sy1 + dx.y*t11);
f1.z += fp1.x*wxScaled.x*( -dx.x*dx.z) + fp1.z*(dz.x*sz1-t12) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.x + dy.x*sz1 + dx.z*t11);
f2.x += fp1.x*wxScaled.y*(1-dx.x*dx.x) + fp1.z*(dz.x*sx2 ) + fp1.y*((-dx.x*dy.x )*wxScaled.y + dy.x*sx2 - dx.y*t22 - dx.z*t23) + f.x*originWeights.y;
f2.y += fp1.x*wxScaled.y*( -dx.x*dx.y) + fp1.z*(dz.x*sy2+t23) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.y + dy.x*sy2 + dx.y*t21);
f2.z += fp1.x*wxScaled.y*( -dx.x*dx.z) + fp1.z*(dz.x*sz2-t22) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.y + dy.x*sz2 + dx.z*t21);
f3.x += fp1.x*wxScaled.z*(1-dx.x*dx.x) + fp1.z*(dz.x*sx3 ) + fp1.y*((-dx.x*dy.x )*wxScaled.z + dy.x*sx3 - dx.y*t32 - dx.z*t33) + f.x*originWeights.z;
f3.y += fp1.x*wxScaled.z*( -dx.x*dx.y) + fp1.z*(dz.x*sy3+t33) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.z + dy.x*sy3 + dx.y*t31);
f3.z += fp1.x*wxScaled.z*( -dx.x*dx.z) + fp1.z*(dz.x*sz3-t32) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.z + dy.x*sz3 + dx.z*t31);
f1.x += fp2.x*wxScaled.x*( -dx.y*dx.x) + fp2.z*(dz.y*sx1-t13) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.x - dy.y*sx1 - dx.x*t12);
f1.y += fp2.x*wxScaled.x*(1-dx.y*dx.y) + fp2.z*(dz.y*sy1 ) - fp2.y*(( dx.y*dy.y )*wxScaled.x - dy.y*sy1 + dx.x*t11 + dx.z*t13) + f.y*originWeights.x;
f1.z += fp2.x*wxScaled.x*( -dx.y*dx.z) + fp2.z*(dz.y*sz1+t11) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.x - dy.y*sz1 - dx.z*t12);
f2.x += fp2.x*wxScaled.y*( -dx.y*dx.x) + fp2.z*(dz.y*sx2-t23) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.y - dy.y*sx2 - dx.x*t22);
f2.y += fp2.x*wxScaled.y*(1-dx.y*dx.y) + fp2.z*(dz.y*sy2 ) - fp2.y*(( dx.y*dy.y )*wxScaled.y - dy.y*sy2 + dx.x*t21 + dx.z*t23) + f.y*originWeights.y;
f2.z += fp2.x*wxScaled.y*( -dx.y*dx.z) + fp2.z*(dz.y*sz2+t21) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.y - dy.y*sz2 - dx.z*t22);
f3.x += fp2.x*wxScaled.z*( -dx.y*dx.x) + fp2.z*(dz.y*sx3-t33) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.z - dy.y*sx3 - dx.x*t32);
f3.y += fp2.x*wxScaled.z*(1-dx.y*dx.y) + fp2.z*(dz.y*sy3 ) - fp2.y*(( dx.y*dy.y )*wxScaled.z - dy.y*sy3 + dx.x*t31 + dx.z*t33) + f.y*originWeights.z;
f3.z += fp2.x*wxScaled.z*( -dx.y*dx.z) + fp2.z*(dz.y*sz3+t31) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.z - dy.y*sz3 - dx.z*t32);
f1.x += fp3.x*wxScaled.x*( -dx.z*dx.x) + fp3.z*(dz.z*sx1+t12) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.x + dy.z*sx1 + dx.x*t13);
f1.y += fp3.x*wxScaled.x*( -dx.z*dx.y) + fp3.z*(dz.z*sy1-t11) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.x + dy.z*sy1 + dx.y*t13);
f1.z += fp3.x*wxScaled.x*(1-dx.z*dx.z) + fp3.z*(dz.z*sz1 ) + fp3.y*((-dx.z*dy.z )*wxScaled.x + dy.z*sz1 - dx.x*t11 - dx.y*t12) + f.z*originWeights.x;
f2.x += fp3.x*wxScaled.y*( -dx.z*dx.x) + fp3.z*(dz.z*sx2+t22) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.y + dy.z*sx2 + dx.x*t23);
f2.y += fp3.x*wxScaled.y*( -dx.z*dx.y) + fp3.z*(dz.z*sy2-t21) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.y + dy.z*sy2 + dx.y*t23);
f2.z += fp3.x*wxScaled.y*(1-dx.z*dx.z) + fp3.z*(dz.z*sz2 ) + fp3.y*((-dx.z*dy.z )*wxScaled.y + dy.z*sz2 - dx.x*t21 - dx.y*t22) + f.z*originWeights.y;
f3.x += fp3.x*wxScaled.z*( -dx.z*dx.x) + fp3.z*(dz.z*sx3+t32) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.z + dy.z*sx3 + dx.x*t33);
f3.y += fp3.x*wxScaled.z*( -dx.z*dx.y) + fp3.z*(dz.z*sy3-t31) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.z + dy.z*sy3 + dx.y*t33);
f3.z += fp3.x*wxScaled.z*(1-dx.z*dx.z) + fp3.z*(dz.z*sz3 ) + fp3.y*((-dx.z*dy.z )*wxScaled.z + dy.z*sz3 - dx.x*t31 - dx.y*t32) + f.z*originWeights.z;
addForce(atoms.y, force, f1);
addForce(atoms.z, force, f2);
addForce(atoms.w, force, f3);
}
}
/**
......
platforms/cuda/tests/TestCudaVirtualSites.cpp
View file @ c88213f8
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012 Stanford University and the Authors. *
* Portions copyright (c) 2012-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -213,6 +213,84 @@ void testOutOfPlane() {
}
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
// Set the particles at random positions.
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
context.applyConstraints(0.0001);
// See if the virtual site is positioned correctly.
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
context.applyConstraints(0.0001);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
context.applyConstraints(0.0001);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-3)
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
......@@ -278,6 +356,26 @@ void testConservationLaws() {
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Create a molecule with a LocalCoordinatesSite virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(14, new LocalCoordinatesSite(11, 12, 13, Vec3(0.3, 0.3, 0.4), Vec3(1.0, -0.5, -0.5), Vec3(0, -1.0, 1.0), Vec3(0.2, 0.2, 1.0)));
system.addConstraint(11, 12, 1.0);
system.addConstraint(11, 13, 1.0);
system.addConstraint(12, 13, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+11, j+11, 0, 1, 0);
}
positions.push_back(Vec3(1, 2, 0));
positions.push_back(Vec3(2, 2, 0));
positions.push_back(Vec3(1, 3, 0));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
......@@ -315,7 +413,7 @@ void testConservationLaws() {
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.02);
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.03);
integrator.step(1);
}
}
......@@ -434,6 +532,7 @@ int main(int argc, char* argv[]) {
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testConservationLaws();
testReordering();
}
......
platforms/opencl/include/OpenCLIntegrationUtilities.h
View file @ c88213f8
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009-2013 Stanford University and the Authors. *
* Portions copyright (c) 2009-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -148,6 +148,8 @@ private:
OpenCLArray* vsite3AvgWeights;
OpenCLArray* vsiteOutOfPlaneAtoms;
OpenCLArray* vsiteOutOfPlaneWeights;
OpenCLArray* vsiteLocalCoordsAtoms;
OpenCLArray* vsiteLocalCoordsParams;
int randomPos;
int lastSeed, numVsites;
bool hasInitializedPosConstraintKernels, hasInitializedVelConstraintKernels, ccmaUseDirectBuffer;
......
platforms/opencl/src/OpenCLIntegrationUtilities.cpp
View file @ c88213f8
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009-2013 Stanford University and the Authors. *
* Portions copyright (c) 2009-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -100,7 +100,8 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
ccmaReducedMass(NULL), ccmaAtomConstraints(NULL), ccmaNumAtomConstraints(NULL), ccmaConstraintMatrixColumn(NULL),
ccmaConstraintMatrixValue(NULL), ccmaDelta1(NULL), ccmaDelta2(NULL), ccmaConverged(NULL), ccmaConvergedHostBuffer(NULL),
vsite2AvgAtoms(NULL), vsite2AvgWeights(NULL), vsite3AvgAtoms(NULL), vsite3AvgWeights(NULL),
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false) {
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsAtoms(NULL), vsiteLocalCoordsParams(NULL),
hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false) {
// Create workspace arrays.
if (context.getUseDoublePrecision() || context.getUseMixedPrecision()) {
......@@ -595,6 +596,8 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vector<mm_double4> vsite3AvgWeightVec;
vector<mm_int4> vsiteOutOfPlaneAtomVec;
vector<mm_double4> vsiteOutOfPlaneWeightVec;
vector<mm_int4> vsiteLocalCoordsAtomVec;
vector<cl_double> vsiteLocalCoordsParamVec;
for (int i = 0; i < numAtoms; i++) {
if (system.isVirtualSite(i)) {
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -618,35 +621,65 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsiteOutOfPlaneAtomVec.push_back(mm_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
vsiteOutOfPlaneWeightVec.push_back(mm_double4(site.getWeight12(), site.getWeight13(), site.getWeightCross(), 0.0));
}
else if (dynamic_cast<const LocalCoordinatesSite*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
vsiteLocalCoordsAtomVec.push_back(mm_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
Vec3 origin = site.getOriginWeights();
Vec3 x = site.getXWeights();
Vec3 y = site.getYWeights();
Vec3 pos = site.getLocalPosition();
vsiteLocalCoordsParamVec.push_back(origin[0]);
vsiteLocalCoordsParamVec.push_back(origin[1]);
vsiteLocalCoordsParamVec.push_back(origin[2]);
vsiteLocalCoordsParamVec.push_back(x[0]);
vsiteLocalCoordsParamVec.push_back(x[1]);
vsiteLocalCoordsParamVec.push_back(x[2]);
vsiteLocalCoordsParamVec.push_back(y[0]);
vsiteLocalCoordsParamVec.push_back(y[1]);
vsiteLocalCoordsParamVec.push_back(y[2]);
vsiteLocalCoordsParamVec.push_back(pos[0]);
vsiteLocalCoordsParamVec.push_back(pos[1]);
vsiteLocalCoordsParamVec.push_back(pos[2]);
}
}
}
int num2Avg = vsite2AvgAtomVec.size();
int num3Avg = vsite3AvgAtomVec.size();
int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
int numLocalCoords = vsiteLocalCoordsAtomVec.size();
vsite2AvgAtoms = OpenCLArray::create<mm_int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
vsite3AvgAtoms = OpenCLArray::create<mm_int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
vsiteOutOfPlaneAtoms = OpenCLArray::create<mm_int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
vsiteLocalCoordsAtoms = OpenCLArray::create<mm_int4>(context, max(1, numLocalCoords), "vsiteLocalCoordinatesAtoms");
if (num2Avg > 0)
vsite2AvgAtoms->upload(vsite2AvgAtomVec);
if (num3Avg > 0)
vsite3AvgAtoms->upload(vsite3AvgAtomVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneAtoms->upload(vsiteOutOfPlaneAtomVec);
if (numLocalCoords > 0)
vsiteLocalCoordsAtoms->upload(vsiteLocalCoordsAtomVec);
if (context.getUseDoublePrecision()) {
vsite2AvgWeights = OpenCLArray::create<mm_double2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = OpenCLArray::create<mm_double4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = OpenCLArray::create<mm_double4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = OpenCLArray::create<cl_double>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
if (num2Avg > 0)
vsite2AvgWeights->upload(vsite2AvgWeightVec);
if (num3Avg > 0)
vsite3AvgWeights->upload(vsite3AvgWeightVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
if (numLocalCoords > 0)
vsiteLocalCoordsParams->upload(vsiteLocalCoordsParamVec);
}
else {
vsite2AvgWeights = OpenCLArray::create<mm_float2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = OpenCLArray::create<mm_float4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = OpenCLArray::create<mm_float4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = OpenCLArray::create<float>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
if (num2Avg > 0) {
vector<mm_float2> floatWeights(num2Avg);
for (int i = 0; i < num2Avg; i++)
......@@ -665,6 +698,12 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
floatWeights[i] = mm_float4((float) vsiteOutOfPlaneWeightVec[i].x, (float) vsiteOutOfPlaneWeightVec[i].y, (float) vsiteOutOfPlaneWeightVec[i].z, 0.0f);
vsiteOutOfPlaneWeights->upload(floatWeights);
}
if (numLocalCoords > 0) {
vector<cl_float> floatParams(vsiteLocalCoordsParamVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsParamVec.size(); i++)
floatParams[i] = (cl_float) vsiteLocalCoordsParamVec[i];
vsiteLocalCoordsParams->upload(floatParams);
}
}
// Create the kernels for virtual sites.
......@@ -673,6 +712,7 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
defines["NUM_2_AVERAGE"] = context.intToString(num2Avg);
defines["NUM_3_AVERAGE"] = context.intToString(num3Avg);
defines["NUM_OUT_OF_PLANE"] = context.intToString(numOutOfPlane);
defines["NUM_LOCAL_COORDS"] = context.intToString(numLocalCoords);
cl::Program vsiteProgram = context.createProgram(OpenCLKernelSources::virtualSites, defines);
vsitePositionKernel = cl::Kernel(vsiteProgram, "computeVirtualSites");
int index = 0;
......@@ -685,6 +725,8 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsitePositionKernel.setArg<cl::Buffer>(index++, vsite3AvgWeights->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneAtoms->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneWeights->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsAtoms->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsParams->getDeviceBuffer());
vsiteForceKernel = cl::Kernel(vsiteProgram, "distributeForces");
index = 0;
vsiteForceKernel.setArg<cl::Buffer>(index++, context.getPosq().getDeviceBuffer());
......@@ -697,7 +739,9 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsiteForceKernel.setArg<cl::Buffer>(index++, vsite3AvgWeights->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneAtoms->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneWeights->getDeviceBuffer());
numVsites = num2Avg+num3Avg+numOutOfPlane;
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsAtoms->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsParams->getDeviceBuffer());
numVsites = num2Avg+num3Avg+numOutOfPlane+numLocalCoords;
}
OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
......@@ -751,6 +795,10 @@ OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
delete vsiteOutOfPlaneAtoms;
if (vsiteOutOfPlaneWeights != NULL)
delete vsiteOutOfPlaneWeights;
if (vsiteLocalCoordsAtoms != NULL)
delete vsiteLocalCoordsAtoms;
if (vsiteLocalCoordsParams != NULL)
delete vsiteLocalCoordsParams;
}
void OpenCLIntegrationUtilities::applyConstraints(double tol) {
......
platforms/opencl/src/kernels/virtualSites.cl
View file @ c88213f8
......@@ -32,7 +32,8 @@ __kernel void computeVirtualSites(__global real4* restrict posq,
#endif
__global const int4* restrict avg2Atoms, __global const real2* restrict avg2Weights,
__global const int4* restrict avg3Atoms, __global const real4* restrict avg3Weights,
__global const int4* restrict outOfPlaneAtoms, __global const real4* restrict outOfPlaneWeights) {
__global const int4* restrict outOfPlaneAtoms, __global const real4* restrict outOfPlaneWeights,
__global const int4* restrict localCoordsAtoms, __global const real* restrict localCoordsParams) {
#ifndef USE_MIXED_PRECISION
__global real4* posqCorrection = 0;
#endif
......@@ -76,6 +77,35 @@ __kernel void computeVirtualSites(__global real4* restrict posq,
pos.xyz = pos1.xyz + v12.xyz*weights.x + v13.xyz*weights.y + cross(v12, v13).xyz*weights.z;
storePos(posq, posqCorrection, atoms.x, pos);
}
// Local coordinates sites.
for (int index = get_global_id(0); index < NUM_LOCAL_COORDS; index += get_global_size(0)) {
int4 atoms = localCoordsAtoms[index];
__global const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed4 pos1 = (mixed4) (pos1_4.x, pos1_4.y, pos1_4.z, 0);
mixed4 pos2 = (mixed4) (pos2_4.x, pos2_4.y, pos2_4.z, 0);
mixed4 pos3 = (mixed4) (pos3_4.x, pos3_4.y, pos3_4.z, 0);
mixed4 originWeights = (mixed4) (params[0], params[1], params[2], 0);
mixed4 xWeights = (mixed4) (params[3], params[4], params[5], 0);
mixed4 yWeights = (mixed4) (params[6], params[7], params[8], 0);
mixed4 localPosition = (mixed4) (params[9], params[10], params[11], 0);
mixed4 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed4 xdir = pos1*xWeights.x + pos2*xWeights.y + pos3*xWeights.z;
mixed4 ydir = pos1*yWeights.x + pos2*yWeights.y + pos3*yWeights.z;
mixed4 zdir = cross(xdir, ydir);
xdir *= rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
zdir *= rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
ydir = cross(zdir, xdir);
pos.x = origin.x + xdir.x*localPosition.x + ydir.x*localPosition.y + zdir.x*localPosition.z;
pos.y = origin.y + xdir.y*localPosition.x + ydir.y*localPosition.y + zdir.y*localPosition.z;
pos.z = origin.z + xdir.z*localPosition.x + ydir.z*localPosition.y + zdir.z*localPosition.z;
storePos(posq, posqCorrection, atoms.x, pos);
}
}
/**
......@@ -87,7 +117,8 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
#endif
__global const int4* restrict avg2Atoms, __global const real2* restrict avg2Weights,
__global const int4* restrict avg3Atoms, __global const real4* restrict avg3Weights,
__global const int4* restrict outOfPlaneAtoms, __global const real4* restrict outOfPlaneWeights) {
__global const int4* restrict outOfPlaneAtoms, __global const real4* restrict outOfPlaneWeights,
__global const int4* restrict localCoordsAtoms, __global const real* restrict localCoordsParams) {
#ifndef USE_MIXED_PRECISION
__global real4* posqCorrection = 0;
#endif
......@@ -150,4 +181,90 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
force[atoms.z] = f2;
force[atoms.w] = f3;
}
// Local coordinates sites.
for (int index = get_global_id(0); index < NUM_LOCAL_COORDS; index += get_global_size(0)) {
int4 atoms = localCoordsAtoms[index];
__global const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed4 pos1 = (mixed4) (pos1_4.x, pos1_4.y, pos1_4.z, 0);
mixed4 pos2 = (mixed4) (pos2_4.x, pos2_4.y, pos2_4.z, 0);
mixed4 pos3 = (mixed4) (pos3_4.x, pos3_4.y, pos3_4.z, 0);
mixed4 originWeights = (mixed4) (params[0], params[1], params[2], 0);
mixed4 wx = (mixed4) (params[3], params[4], params[5], 0);
mixed4 wy = (mixed4) (params[6], params[7], params[8], 0);
mixed4 localPosition = (mixed4) (params[9], params[10], params[11], 0);
mixed4 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed4 xdir = pos1*wx.x + pos2*wx.y + pos3*wx.z;
mixed4 ydir = pos1*wy.x + pos2*wy.y + pos3*wy.z;
mixed4 zdir = cross(xdir, ydir);
mixed invNormXdir = rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed invNormZdir = rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed4 dx = xdir*invNormXdir;
mixed4 dz = zdir*invNormZdir;
mixed4 dy = cross(dz, dx);
// The derivatives for this case are very complicated. They were computed with SymPy then simplified by hand.
mixed t11 = (wx.x*ydir.x-wy.x*xdir.x)*invNormZdir;
mixed t12 = (wx.x*ydir.y-wy.x*xdir.y)*invNormZdir;
mixed t13 = (wx.x*ydir.z-wy.x*xdir.z)*invNormZdir;
mixed t21 = (wx.y*ydir.x-wy.y*xdir.x)*invNormZdir;
mixed t22 = (wx.y*ydir.y-wy.y*xdir.y)*invNormZdir;
mixed t23 = (wx.y*ydir.z-wy.y*xdir.z)*invNormZdir;
mixed t31 = (wx.z*ydir.x-wy.z*xdir.x)*invNormZdir;
mixed t32 = (wx.z*ydir.y-wy.z*xdir.y)*invNormZdir;
mixed t33 = (wx.z*ydir.z-wy.z*xdir.z)*invNormZdir;
mixed sx1 = t13*dz.y-t12*dz.z;
mixed sy1 = t11*dz.z-t13*dz.x;
mixed sz1 = t12*dz.x-t11*dz.y;
mixed sx2 = t23*dz.y-t22*dz.z;
mixed sy2 = t21*dz.z-t23*dz.x;
mixed sz2 = t22*dz.x-t21*dz.y;
mixed sx3 = t33*dz.y-t32*dz.z;
mixed sy3 = t31*dz.z-t33*dz.x;
mixed sz3 = t32*dz.x-t31*dz.y;
mixed4 wxScaled = wx*invNormXdir;
real4 f = force[atoms.x];
real4 f1 = force[atoms.y];
real4 f2 = force[atoms.z];
real4 f3 = force[atoms.w];
mixed4 fp1 = localPosition*f.x;
mixed4 fp2 = localPosition*f.y;
mixed4 fp3 = localPosition*f.z;
f1.x += fp1.x*wxScaled.x*(1-dx.x*dx.x) + fp1.z*(dz.x*sx1 ) + fp1.y*((-dx.x*dy.x )*wxScaled.x + dy.x*sx1 - dx.y*t12 - dx.z*t13) + f.x*originWeights.x;
f1.y += fp1.x*wxScaled.x*( -dx.x*dx.y) + fp1.z*(dz.x*sy1+t13) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.x + dy.x*sy1 + dx.y*t11);
f1.z += fp1.x*wxScaled.x*( -dx.x*dx.z) + fp1.z*(dz.x*sz1-t12) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.x + dy.x*sz1 + dx.z*t11);
f2.x += fp1.x*wxScaled.y*(1-dx.x*dx.x) + fp1.z*(dz.x*sx2 ) + fp1.y*((-dx.x*dy.x )*wxScaled.y + dy.x*sx2 - dx.y*t22 - dx.z*t23) + f.x*originWeights.y;
f2.y += fp1.x*wxScaled.y*( -dx.x*dx.y) + fp1.z*(dz.x*sy2+t23) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.y + dy.x*sy2 + dx.y*t21);
f2.z += fp1.x*wxScaled.y*( -dx.x*dx.z) + fp1.z*(dz.x*sz2-t22) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.y + dy.x*sz2 + dx.z*t21);
f3.x += fp1.x*wxScaled.z*(1-dx.x*dx.x) + fp1.z*(dz.x*sx3 ) + fp1.y*((-dx.x*dy.x )*wxScaled.z + dy.x*sx3 - dx.y*t32 - dx.z*t33) + f.x*originWeights.z;
f3.y += fp1.x*wxScaled.z*( -dx.x*dx.y) + fp1.z*(dz.x*sy3+t33) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.z + dy.x*sy3 + dx.y*t31);
f3.z += fp1.x*wxScaled.z*( -dx.x*dx.z) + fp1.z*(dz.x*sz3-t32) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.z + dy.x*sz3 + dx.z*t31);
f1.x += fp2.x*wxScaled.x*( -dx.y*dx.x) + fp2.z*(dz.y*sx1-t13) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.x - dy.y*sx1 - dx.x*t12);
f1.y += fp2.x*wxScaled.x*(1-dx.y*dx.y) + fp2.z*(dz.y*sy1 ) - fp2.y*(( dx.y*dy.y )*wxScaled.x - dy.y*sy1 + dx.x*t11 + dx.z*t13) + f.y*originWeights.x;
f1.z += fp2.x*wxScaled.x*( -dx.y*dx.z) + fp2.z*(dz.y*sz1+t11) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.x - dy.y*sz1 - dx.z*t12);
f2.x += fp2.x*wxScaled.y*( -dx.y*dx.x) + fp2.z*(dz.y*sx2-t23) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.y - dy.y*sx2 - dx.x*t22);
f2.y += fp2.x*wxScaled.y*(1-dx.y*dx.y) + fp2.z*(dz.y*sy2 ) - fp2.y*(( dx.y*dy.y )*wxScaled.y - dy.y*sy2 + dx.x*t21 + dx.z*t23) + f.y*originWeights.y;
f2.z += fp2.x*wxScaled.y*( -dx.y*dx.z) + fp2.z*(dz.y*sz2+t21) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.y - dy.y*sz2 - dx.z*t22);
f3.x += fp2.x*wxScaled.z*( -dx.y*dx.x) + fp2.z*(dz.y*sx3-t33) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.z - dy.y*sx3 - dx.x*t32);
f3.y += fp2.x*wxScaled.z*(1-dx.y*dx.y) + fp2.z*(dz.y*sy3 ) - fp2.y*(( dx.y*dy.y )*wxScaled.z - dy.y*sy3 + dx.x*t31 + dx.z*t33) + f.y*originWeights.z;
f3.z += fp2.x*wxScaled.z*( -dx.y*dx.z) + fp2.z*(dz.y*sz3+t31) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.z - dy.y*sz3 - dx.z*t32);
f1.x += fp3.x*wxScaled.x*( -dx.z*dx.x) + fp3.z*(dz.z*sx1+t12) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.x + dy.z*sx1 + dx.x*t13);
f1.y += fp3.x*wxScaled.x*( -dx.z*dx.y) + fp3.z*(dz.z*sy1-t11) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.x + dy.z*sy1 + dx.y*t13);
f1.z += fp3.x*wxScaled.x*(1-dx.z*dx.z) + fp3.z*(dz.z*sz1 ) + fp3.y*((-dx.z*dy.z )*wxScaled.x + dy.z*sz1 - dx.x*t11 - dx.y*t12) + f.z*originWeights.x;
f2.x += fp3.x*wxScaled.y*( -dx.z*dx.x) + fp3.z*(dz.z*sx2+t22) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.y + dy.z*sx2 + dx.x*t23);
f2.y += fp3.x*wxScaled.y*( -dx.z*dx.y) + fp3.z*(dz.z*sy2-t21) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.y + dy.z*sy2 + dx.y*t23);
f2.z += fp3.x*wxScaled.y*(1-dx.z*dx.z) + fp3.z*(dz.z*sz2 ) + fp3.y*((-dx.z*dy.z )*wxScaled.y + dy.z*sz2 - dx.x*t21 - dx.y*t22) + f.z*originWeights.y;
f3.x += fp3.x*wxScaled.z*( -dx.z*dx.x) + fp3.z*(dz.z*sx3+t32) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.z + dy.z*sx3 + dx.x*t33);
f3.y += fp3.x*wxScaled.z*( -dx.z*dx.y) + fp3.z*(dz.z*sy3-t31) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.z + dy.z*sy3 + dx.y*t33);
f3.z += fp3.x*wxScaled.z*(1-dx.z*dx.z) + fp3.z*(dz.z*sz3 ) + fp3.y*((-dx.z*dy.z )*wxScaled.z + dy.z*sz3 - dx.x*t31 - dx.y*t32) + f.z*originWeights.z;
force[atoms.y] = f1;
force[atoms.z] = f2;
force[atoms.w] = f3;
}
}
platforms/opencl/tests/TestOpenCLVirtualSites.cpp
View file @ c88213f8
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012 Stanford University and the Authors. *
* Portions copyright (c) 2012-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -213,6 +213,84 @@ void testOutOfPlane() {
}
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
// Set the particles at random positions.
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
context.applyConstraints(0.0001);
// See if the virtual site is positioned correctly.
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
context.applyConstraints(0.0001);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
context.applyConstraints(0.0001);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-3)
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
......@@ -278,6 +356,26 @@ void testConservationLaws() {
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Create a molecule with a LocalCoordinatesSite virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(14, new LocalCoordinatesSite(11, 12, 13, Vec3(0.3, 0.3, 0.4), Vec3(1.0, -0.5, -0.5), Vec3(0, -1.0, 1.0), Vec3(0.2, 0.2, 1.0)));
system.addConstraint(11, 12, 1.0);
system.addConstraint(11, 13, 1.0);
system.addConstraint(12, 13, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+11, j+11, 0, 1, 0);
}
positions.push_back(Vec3(1, 2, 0));
positions.push_back(Vec3(2, 2, 0));
positions.push_back(Vec3(1, 3, 0));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
......@@ -315,7 +413,7 @@ void testConservationLaws() {
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.02);
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.03);
integrator.step(1);
}
}
......@@ -434,6 +532,7 @@ int main(int argc, char* argv[]) {
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testConservationLaws();
testReordering();
}
......
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