Skip to content

GitLab

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

Continuing to implement new CUDA platform: virtual sites

parent 6e842d8c
Hide whitespace changes
Inline Side-by-side
Showing with 641 additions and 42 deletions
platforms/cuda2/src/CudaIntegrationUtilities.cpp
View file @ bf1f6f32
......@@ -569,11 +569,11 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
// Build the list of virtual sites.
vector<int4> vsite2AvgAtomVec;
vector<float2> vsite2AvgWeightVec;
vector<double2> vsite2AvgWeightVec;
vector<int4> vsite3AvgAtomVec;
vector<float4> vsite3AvgWeightVec;
vector<double4> vsite3AvgWeightVec;
vector<int4> vsiteOutOfPlaneAtomVec;
vector<float4> vsiteOutOfPlaneWeightVec;
vector<double4> vsiteOutOfPlaneWeightVec;
for (int i = 0; i < numAtoms; i++) {
if (system.isVirtualSite(i)) {
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -581,21 +581,21 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
const TwoParticleAverageSite& site = dynamic_cast<const TwoParticleAverageSite&>(system.getVirtualSite(i));
vsite2AvgAtomVec.push_back(make_int4(i, site.getParticle(0), site.getParticle(1), 0));
vsite2AvgWeightVec.push_back(make_float2((float) site.getWeight(0), (float) site.getWeight(1)));
vsite2AvgWeightVec.push_back(make_double2(site.getWeight(0), site.getWeight(1)));
}
else if (dynamic_cast<const ThreeParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
// A three particle average.
const ThreeParticleAverageSite& site = dynamic_cast<const ThreeParticleAverageSite&>(system.getVirtualSite(i));
vsite3AvgAtomVec.push_back(make_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
vsite3AvgWeightVec.push_back(make_float4((float) site.getWeight(0), (float) site.getWeight(1), (float) site.getWeight(2), 0.0f));
vsite3AvgWeightVec.push_back(make_double4(site.getWeight(0), site.getWeight(1), site.getWeight(2), 0.0));
}
else if (dynamic_cast<const OutOfPlaneSite*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
const OutOfPlaneSite& site = dynamic_cast<const OutOfPlaneSite&>(system.getVirtualSite(i));
vsiteOutOfPlaneAtomVec.push_back(make_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
vsiteOutOfPlaneWeightVec.push_back(make_float4((float) site.getWeight12(), (float) site.getWeight13(), (float) site.getWeightCross(), 0.0f));
vsiteOutOfPlaneWeightVec.push_back(make_double4(site.getWeight12(), site.getWeight13(), site.getWeightCross(), 0.0));
}
}
}
......@@ -603,48 +603,59 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
int num3Avg = vsite3AvgAtomVec.size();
int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
vsite2AvgAtoms = CudaArray::create<int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
vsite2AvgWeights = CudaArray::create<float2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgAtoms = CudaArray::create<int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
vsite3AvgWeights = CudaArray::create<float4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneAtoms = CudaArray::create<int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
vsiteOutOfPlaneWeights = CudaArray::create<float4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
if (num2Avg > 0) {
if (num2Avg > 0)
vsite2AvgAtoms->upload(vsite2AvgAtomVec);
vsite2AvgWeights->upload(vsite2AvgWeightVec);
}
if (num3Avg > 0) {
if (num3Avg > 0)
vsite3AvgAtoms->upload(vsite3AvgAtomVec);
vsite3AvgWeights->upload(vsite3AvgWeightVec);
}
if (numOutOfPlane > 0) {
if (numOutOfPlane > 0)
vsiteOutOfPlaneAtoms->upload(vsiteOutOfPlaneAtomVec);
vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
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");
if (num2Avg > 0)
vsite2AvgWeights->upload(vsite2AvgWeightVec);
if (num3Avg > 0)
vsite3AvgWeights->upload(vsite3AvgWeightVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
}
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");
if (num2Avg > 0) {
vector<float2> floatWeights(num2Avg);
for (int i = 0; i < num2Avg; i++)
floatWeights[i] = make_float2((float) vsite2AvgWeightVec[i].x, (float) vsite2AvgWeightVec[i].y);
vsite2AvgWeights->upload(floatWeights);
}
if (num3Avg > 0) {
vector<float4> floatWeights(num3Avg);
for (int i = 0; i < num3Avg; i++)
floatWeights[i] = make_float4((float) vsite3AvgWeightVec[i].x, (float) vsite3AvgWeightVec[i].y, (float) vsite3AvgWeightVec[i].z, 0.0f);
vsite3AvgWeights->upload(floatWeights);
}
if (numOutOfPlane > 0) {
vector<float4> floatWeights(numOutOfPlane);
for (int i = 0; i < numOutOfPlane; i++)
floatWeights[i] = make_float4((float) vsiteOutOfPlaneWeightVec[i].x, (float) vsiteOutOfPlaneWeightVec[i].y, (float) vsiteOutOfPlaneWeightVec[i].z, 0.0f);
vsiteOutOfPlaneWeights->upload(floatWeights);
}
}
// Create the kernels for virtual sites.
map<string, string> defines;
defines["NUM_2_AVERAGE"] = context.intToString(num2Avg);
defines["NUM_3_AVERAGE"] = context.intToString(num3Avg);
defines["NUM_OUT_OF_PLANE"] = context.intToString(numOutOfPlane);
// CUmodule vsiteModule = context.createModule(CudaKernelSources::virtualSites, defines);
// vsitePositionKernel = context.getKernel(vsiteModule, "computeVirtualSites");
// vsitePositionKernel.setArg<cl::Buffer>(0, context.getPosq().getDeviceBuffer());
// vsitePositionKernel.setArg<cl::Buffer>(1, vsite2AvgAtoms->getDeviceBuffer());
// vsitePositionKernel.setArg<cl::Buffer>(2, vsite2AvgWeights->getDeviceBuffer());
// vsitePositionKernel.setArg<cl::Buffer>(3, vsite3AvgAtoms->getDeviceBuffer());
// vsitePositionKernel.setArg<cl::Buffer>(4, vsite3AvgWeights->getDeviceBuffer());
// vsitePositionKernel.setArg<cl::Buffer>(5, vsiteOutOfPlaneAtoms->getDeviceBuffer());
// vsitePositionKernel.setArg<cl::Buffer>(6, vsiteOutOfPlaneWeights->getDeviceBuffer());
// vsiteForceKernel = context.getKernel(vsiteModule, "distributeForces");
// vsiteForceKernel.setArg<cl::Buffer>(0, context.getPosq().getDeviceBuffer());
// // Skip argument 1: the force array hasn't been created yet.
// vsiteForceKernel.setArg<cl::Buffer>(2, vsite2AvgAtoms->getDeviceBuffer());
// vsiteForceKernel.setArg<cl::Buffer>(3, vsite2AvgWeights->getDeviceBuffer());
// vsiteForceKernel.setArg<cl::Buffer>(4, vsite3AvgAtoms->getDeviceBuffer());
// vsiteForceKernel.setArg<cl::Buffer>(5, vsite3AvgWeights->getDeviceBuffer());
// vsiteForceKernel.setArg<cl::Buffer>(6, vsiteOutOfPlaneAtoms->getDeviceBuffer());
// vsiteForceKernel.setArg<cl::Buffer>(7, vsiteOutOfPlaneWeights->getDeviceBuffer());
defines["PADDED_NUM_ATOMS"] = context.intToString(context.getPaddedNumAtoms());
CUmodule vsiteModule = context.createModule(CudaKernelSources::vectorOps+CudaKernelSources::virtualSites, defines);
vsitePositionKernel = context.getKernel(vsiteModule, "computeVirtualSites");
vsiteForceKernel = context.getKernel(vsiteModule, "distributeForces");
numVsites = num2Avg+num3Avg+numOutOfPlane;
}
......@@ -777,15 +788,22 @@ void CudaIntegrationUtilities::applyConstraints(bool constrainVelocities, double
}
void CudaIntegrationUtilities::computeVirtualSites() {
// if (numVsites > 0)
// context.executeKernel(vsitePositionKernel, numVsites);
if (numVsites > 0) {
void* args[] = {&context.getPosq().getDevicePointer(), &vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
&vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer()};
context.executeKernel(vsitePositionKernel, args, numVsites);
}
}
void CudaIntegrationUtilities::distributeForcesFromVirtualSites() {
// if (numVsites > 0) {
// vsiteForceKernel.setArg<cl::Buffer>(1, context.getForce().getDeviceBuffer());
// context.executeKernel(vsiteForceKernel, numVsites);
// }
if (numVsites > 0) {
void* args[] = {&context.getPosq().getDevicePointer(), &context.getForce().getDevicePointer(),
&vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
&vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer()};
context.executeKernel(vsiteForceKernel, args, numVsites);
}
}
void CudaIntegrationUtilities::initRandomNumberGenerator(unsigned int randomNumberSeed) {
......
platforms/cuda2/src/kernels/virtualSites.cu 0 → 100644
View file @ bf1f6f32
/**
* Compute the positions of virtual sites
*/
extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights) {
// Two particle average sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_2_AVERAGE; index += blockDim.x*gridDim.x) {
int4 atoms = avg2Atoms[index];
real2 weights = avg2Weights[index];
real4 pos = posq[atoms.x];
real4 pos1 = posq[atoms.y];
real4 pos2 = posq[atoms.z];
pos.x = pos1.x*weights.x + pos2.x*weights.y;
pos.y = pos1.y*weights.x + pos2.y*weights.y;
pos.z = pos1.z*weights.x + pos2.z*weights.y;
posq[atoms.x] = pos;
}
// Three particle average sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_3_AVERAGE; index += blockDim.x*gridDim.x) {
int4 atoms = avg3Atoms[index];
real4 weights = avg3Weights[index];
real4 pos = posq[atoms.x];
real4 pos1 = posq[atoms.y];
real4 pos2 = posq[atoms.z];
real4 pos3 = posq[atoms.w];
pos.x = pos1.x*weights.x + pos2.x*weights.y + pos3.x*weights.z;
pos.y = pos1.y*weights.x + pos2.y*weights.y + pos3.y*weights.z;
pos.z = pos1.z*weights.x + pos2.z*weights.y + pos3.z*weights.z;
posq[atoms.x] = pos;
}
// Out of plane sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_OUT_OF_PLANE; index += blockDim.x*gridDim.x) {
int4 atoms = outOfPlaneAtoms[index];
real4 weights = outOfPlaneWeights[index];
real4 pos = posq[atoms.x];
real4 pos1 = posq[atoms.y];
real4 pos2 = posq[atoms.z];
real4 pos3 = posq[atoms.w];
real4 v12 = pos2-pos1;
real4 v13 = pos3-pos1;
real3 cr = cross(v12, v13);
pos.x = pos1.x + v12.x*weights.x + v13.x*weights.y + cr.x*weights.z;
pos.y = pos1.y + v12.y*weights.x + v13.y*weights.y + cr.y*weights.z;
pos.z = pos1.z + v12.z*weights.x + v13.z*weights.y + cr.z*weights.z;
posq[atoms.x] = pos;
}
}
inline __device__ real3 loadForce(int index, long long* __restrict__ force) {
real scale = RECIP((real) 0xFFFFFFFF);
return make_real3(scale*force[index], scale*force[index+PADDED_NUM_ATOMS], scale*force[index+PADDED_NUM_ATOMS*2]);
}
inline __device__ void storeForce(int index, long long* __restrict__ force, real3 value) {
force[index] = (long long) (value.x*0xFFFFFFFF);
force[index+PADDED_NUM_ATOMS] = (long long) (value.y*0xFFFFFFFF);
force[index+2*PADDED_NUM_ATOMS] = (long long) (value.z*0xFFFFFFFF);
}
/**
* Distribute forces from virtual sites to the atoms they are based on.
*/
extern "C" __global__ void distributeForces(const real4* __restrict__ posq, 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) {
// Two particle average sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_2_AVERAGE; index += blockDim.x*gridDim.x) {
int4 atoms = avg2Atoms[index];
real2 weights = avg2Weights[index];
real3 f = loadForce(atoms.x, force);
real3 f1 = loadForce(atoms.y, force);
real3 f2 = loadForce(atoms.z, force);
f1 += f*weights.x;
f2 += f*weights.y;
storeForce(atoms.y, force, f1);
storeForce(atoms.z, force, f2);
}
// Three particle average sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_3_AVERAGE; index += blockDim.x*gridDim.x) {
int4 atoms = avg3Atoms[index];
real4 weights = avg3Weights[index];
real3 f = loadForce(atoms.x, force);
real3 f1 = loadForce(atoms.y, force);
real3 f2 = loadForce(atoms.z, force);
real3 f3 = loadForce(atoms.w, force);
f1 += f*weights.x;
f2 += f*weights.y;
f3 += f*weights.z;
storeForce(atoms.y, force, f1);
storeForce(atoms.z, force, f2);
storeForce(atoms.w, force, f3);
}
// Out of plane sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_OUT_OF_PLANE; index += blockDim.x*gridDim.x) {
int4 atoms = outOfPlaneAtoms[index];
real4 weights = outOfPlaneWeights[index];
real4 pos1 = posq[atoms.y];
real4 pos2 = posq[atoms.z];
real4 pos3 = posq[atoms.w];
real4 v12 = pos2-pos1;
real4 v13 = pos3-pos1;
real3 f = loadForce(atoms.x, force);
real3 f1 = loadForce(atoms.y, force);
real3 f2 = loadForce(atoms.z, force);
real3 f3 = loadForce(atoms.w, force);
real3 fp2 = make_real3(weights.x*f.x - weights.z*v13.z*f.y + weights.z*v13.y*f.z,
weights.z*v13.z*f.x + weights.x*f.y - weights.z*v13.x*f.z,
-weights.z*v13.y*f.x + weights.z*v13.x*f.y + weights.x*f.z);
real3 fp3 = make_real3(weights.y*f.x + weights.z*v12.z*f.y - weights.z*v12.y*f.z,
-weights.z*v12.z*f.x + weights.y*f.y + weights.z*v12.x*f.z,
weights.z*v12.y*f.x - weights.z*v12.x*f.y + weights.y*f.z);
f1 += f-fp2-fp3;
f2 += fp2;
f3 += fp3;
storeForce(atoms.y, force, f1);
storeForce(atoms.z, force, f2);
storeForce(atoms.w, force, f3);
}
}
platforms/cuda2/tests/TestCudaVirtualSites.cpp 0 → 100644
View file @ bf1f6f32
/* -------------------------------------------------------------------------- *
* 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) 2012 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 CUDA implementation of virtual sites.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VirtualSite.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
/**
* Check that massless particles are handled correctly.
*/
void testMasslessParticle() {
CudaPlatform platform;
System system;
system.addParticle(0.0);
system.addParticle(1.0);
CustomBondForce* bonds = new CustomBondForce("-1/r");
system.addForce(bonds);
vector<double> params;
bonds->addBond(0, 1, params);
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
vector<Vec3> velocities(2);
velocities[0] = Vec3(0, 0, 0);
velocities[1] = Vec3(0, 1, 0);
context.setVelocities(velocities);
// The second particle should move in a circular orbit around the first one.
// Compare it to the analytical solution.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities | State::Forces);
double time = state.getTime();
ASSERT_EQUAL_VEC(Vec3(), state.getPositions()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(), state.getVelocities()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
integrator.step(1);
}
}
/**
* Test a TwoParticleAverageSite virtual site.
*/
void testTwoParticleAverage() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.8, 0.2));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-4);
ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-4);
integrator.step(1);
}
}
/**
* Test a ThreeParticleAverageSite virtual site.
*/
void testThreeParticleAverage() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new ThreeParticleAverageSite(0, 1, 2, 0.2, 0.3, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-5);
integrator.step(1);
}
}
/**
* Test an OutOfPlaneSite virtual site.
*/
void testOutOfPlane() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new OutOfPlaneSite(0, 1, 2, 0.3, 0.4, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 v12 = pos[1]-pos[0];
Vec3 v13 = pos[2]-pos[0];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-5);
const vector<Vec3>& f = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-5);
Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-5);
integrator.step(1);
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
*/
void testConservationLaws() {
CudaPlatform platform;
System system;
NonbondedForce* forceField = new NonbondedForce();
system.addForce(forceField);
vector<Vec3> positions;
// Create a linear molecule with a TwoParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.4, 0.6));
system.addConstraint(0, 1, 2.0);
for (int i = 0; i < 3; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i, j, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(2, 0, 0));
positions.push_back(Vec3());
// Create a planar molecule with a ThreeParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(6, new ThreeParticleAverageSite(3, 4, 5, 0.3, 0.5, 0.2));
system.addConstraint(3, 4, 1.0);
system.addConstraint(3, 5, 1.0);
system.addConstraint(4, 5, 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+3, j+3, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 1));
positions.push_back(Vec3(1, 0, 1));
positions.push_back(Vec3(0, 1, 1));
positions.push_back(Vec3());
// Create a tetrahedral molecule with an OutOfPlane virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(10, new OutOfPlaneSite(7, 8, 9, 0.3, 0.5, 0.2));
system.addConstraint(7, 8, 1.0);
system.addConstraint(7, 9, 1.0);
system.addConstraint(8, 9, 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+7, j+7, 0, 1, 0);
}
positions.push_back(Vec3(1, 0, -1));
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
int numParticles = system.getNumParticles();
double initialEnergy;
Vec3 initialMomentum, initialAngularMomentum;
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
const vector<Vec3>& pos = state.getPositions();
const vector<Vec3>& vel = state.getVelocities();
const vector<Vec3>& f = state.getForces();
double energy = state.getPotentialEnergy();
for (int j = 0; j < numParticles; j++) {
Vec3 v = vel[j] + f[j]*0.5*integrator.getStepSize();
energy += 0.5*system.getParticleMass(j)*v.dot(v);
}
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
Vec3 momentum;
for (int j = 0; j < numParticles; j++)
momentum += vel[j]*system.getParticleMass(j);
if (i == 0)
initialMomentum = momentum;
else
ASSERT_EQUAL_VEC(initialMomentum, momentum, 0.02);
Vec3 angularMomentum;
for (int j = 0; j < numParticles; j++)
angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.02);
integrator.step(1);
}
}
/**
* Make sure that atom reordering respects virtual sites.
*/
void testReordering() {
const double cutoff = 2.0;
const double boxSize = 20.0;
CudaPlatform platform;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
nonbonded->setCutoffDistance(cutoff);
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
// Create linear molecules with TwoParticleAverage virtual sites.
for (int i = 0; i < 50; i++) {
int start = system.getNumParticles();
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(start+2, new TwoParticleAverageSite(start, start+1, 0.4, 0.6));
system.addConstraint(start, start+1, 2.0);
for (int i = 0; i < 3; i++) {
nonbonded->addParticle(0, 0.2, 1);
for (int j = 0; j < i; j++)
nonbonded->addException(start+i, start+j, 0, 1, 0);
}
Vec3 pos(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions.push_back(pos);
positions.push_back(pos+Vec3(2, 0, 0));
positions.push_back(Vec3());
}
// Create planar molecules with ThreeParticleAverage virtual sites.
for (int i = 0; i < 50; i++) {
int start = system.getNumParticles();
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(start+3, new ThreeParticleAverageSite(start, start+1, start+2, 0.3, 0.5, 0.2));
system.addConstraint(start, start+1, 1.0);
system.addConstraint(start, start+2, 1.0);
system.addConstraint(start+1, start+2, sqrt(2.0));
for (int i = 0; i < 4; i++) {
nonbonded->addParticle(0, 0.2, 1);
for (int j = 0; j < i; j++)
nonbonded->addException(start+i, start+j, 0, 1, 0);
}
Vec3 pos(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions.push_back(pos);
positions.push_back(pos+Vec3(1, 0, 0));
positions.push_back(pos+Vec3(0, 1, 0));
positions.push_back(Vec3());
}
// Create tetrahedral molecules with OutOfPlane virtual sites.
for (int i = 0; i < 50; i++) {
int start = system.getNumParticles();
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(start+3, new OutOfPlaneSite(start, start+1, start+2, 0.3, 0.5, 0.2));
system.addConstraint(start, start+1, 1.0);
system.addConstraint(start, start+2, 1.0);
system.addConstraint(start+1, start+2, sqrt(2.0));
for (int i = 0; i < 4; i++) {
nonbonded->addParticle(0, 0.2, 1);
for (int j = 0; j < i; j++)
nonbonded->addException(start+i, start+j, 0, 1, 0);
}
Vec3 pos(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions.push_back(pos);
positions.push_back(pos+Vec3(1, 0, 0));
positions.push_back(pos+Vec3(0, 1, 0));
positions.push_back(Vec3());
}
// Simulate it and check conservation laws.
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions);
const vector<Vec3>& pos = state.getPositions();
for (int j = 0; j < 150; j += 3)
ASSERT_EQUAL_VEC(pos[j]*0.4+pos[j+1]*0.6, pos[j+2], 1e-5);
for (int j = 150; j < 350; j += 4)
ASSERT_EQUAL_VEC(pos[j]*0.3+pos[j+1]*0.5+pos[j+2]*0.2, pos[j+3], 1e-5);
for (int j = 350; j < 550; j += 4) {
Vec3 v12 = pos[j+1]-pos[j];
Vec3 v13 = pos[j+2]-pos[j];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[j]+v12*0.3+v13*0.5+cross*0.2, pos[j+3], 1e-5);
}
integrator.step(1);
}
}
int main() {
try {
testMasslessParticle();
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testConservationLaws();
testReordering();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment