Skip to content

GitLab

Unverified Commit be19e022 authored by Peter Eastman's avatar Peter Eastman Committed by GitHub
Browse files

Converted RPMD plugin to common platform (#3079) 

* Converted RPMD plugin to common platform

* Merged RPMD tests for different platforms

* Try to fix errors on CPU OpenCL
parent 98d81730
Hide whitespace changes
Inline Side-by-side
Showing with 611 additions and 2474 deletions
plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.cpp
View file @ be19e022
......@@ -29,545 +29,19 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of RPMDIntegrator.
*/
#include "CudaTests.h"
#include "TestRpmd.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/Context.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/RPMDIntegrator.h"
#include "openmm/RPMDMonteCarloBarostat.h"
#include "openmm/VirtualSite.h"
#include "SimTKOpenMMUtilities.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
extern "C" void registerRPMDCudaKernelFactories();
using namespace OpenMM;
using namespace std;
extern "C" OPENMM_EXPORT void registerRPMDCudaKernelFactories();
void testFreeParticles() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
integ.setPositions(i, positions);
}
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0);
const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Positions | State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
double rg2 = 0.0;
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*mass*v.dot(v);
for (int m = 0; m < numCopies; m++) {
Vec3 delta = state[k].getPositions()[j]-state[m].getPositions()[j];
rg2 += delta.dot(delta);
}
}
rg[j] += rg2/(2*numCopies*numCopies);
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
double meanRg2 = 0.0;
for (int i = 0; i < numParticles; i++)
meanRg2 += rg[i];
meanRg2 /= numSteps*numParticles;
double expectedRg = hbar/(2*sqrt(mass*BOLTZ*temperature));
ASSERT_USUALLY_EQUAL_TOL(expectedRg, sqrt(meanRg2), 1e-3);
}
void testParaHydrogen() {
const int numParticles = 32;
const int numCopies = 12;
const double temperature = 25.0;
const double mass = 2.0;
const double boxSize = 1.1896;
const int numSteps = 1000;
const int numBins = 200;
const double reference[] = {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 4.932814042206152e-5, 1.244331241336431e-4, 4.052316284060125e-4,
1.544810863683946e-3, 4.376197806690222e-3, 1.025847561714293e-2, 2.286702037465422e-2,
4.371052180263602e-2, 7.518538770734748e-2, 0.122351534531647, 0.185758975626622,
0.266399984652322, 0.363380262153250, 0.473696401293219, 0.595312098494172,
0.726049519422861, 0.862264551954547, 0.991102029379444, 1.1147503922535,
1.23587006992066, 1.33495411932817, 1.42208208736987, 1.49273884004107,
1.54633319690403, 1.58714702233941, 1.60439217751355, 1.61804190608902,
1.60680198476058, 1.58892222973695, 1.56387607986781, 1.52629494593350,
1.48421439018970, 1.43656176771959, 1.38752775598872, 1.33310695719931,
1.28363477223121, 1.23465642750248, 1.18874848666326, 1.14350496170519,
1.10292486009936, 1.06107270157688, 1.02348927970441, 0.989729345271297,
0.959273446941802, 0.932264875865758, 0.908818658748942, 0.890946420768315,
0.869332737718165, 0.856401736350349, 0.842370069917020, 0.834386614237393,
0.826268072171045, 0.821547250199453, 0.818786865315836, 0.819441757028076,
0.819156933383128, 0.822275325148621, 0.828919078023881, 0.837233720599450,
0.846961908186718, 0.855656955481099, 0.864520333201247, 0.876082425547566,
0.886950044046000, 0.900275658318995
};
// Create a box of para-hydrogen.
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize,0,0), Vec3(0,boxSize,0), Vec3(0,0,boxSize));
CustomNonbondedForce* nb = new CustomNonbondedForce("2625.49963*(exp(1.713-1.5671*p-0.00993*p*p)-(12.14/p^6+215.2/p^8-143.1/p^9+4813.9/p^10)*(step(rc-p)*exp(-(rc/p-1)^2)+1-step(rc-p))); p=r/0.05291772108; rc=8.32");
nb->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
nb->setCutoffDistance(boxSize/2);
vector<double> params;
for (int i = 0; i < numParticles; i++)
nb->addParticle(params);
system.addForce(nb);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.0005);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
for (int i = 0; i < numCopies; i++)
integ.setPositions(i, positions);
integ.step(1000);
// Simulate it.
vector<int> counts(numBins, 0);
const double invBoxSize = 1.0/boxSize;
double meanKE = 0.0;
for (int step = 0; step < numSteps; step++) {
integ.step(20);
vector<State> states(numCopies);
for (int i = 0; i < numCopies; i++)
states[i] = integ.getState(i, State::Positions | State::Forces);
// Record the radial distribution function.
const vector<Vec3>& pos = states[0].getPositions();
for (int j = 0; j < numParticles; j++)
for (int k = 0; k < j; k++) {
Vec3 delta = pos[j]-pos[k];
delta[0] -= floor(delta[0]*invBoxSize+0.5)*boxSize;
delta[1] -= floor(delta[1]*invBoxSize+0.5)*boxSize;
delta[2] -= floor(delta[2]*invBoxSize+0.5)*boxSize;
double dist = sqrt(delta.dot(delta));
int bin = (int) (numBins*(dist/boxSize));
counts[bin]++;
}
// Calculate the quantum contribution to the kinetic energy.
vector<Vec3> centroids(numParticles, Vec3());
for (int i = 0; i < numCopies; i++) {
const vector<Vec3>& pos = states[i].getPositions();
for (int j = 0; j < numParticles; j++)
centroids[j] += pos[j];
}
for (int j = 0; j < numParticles; j++)
centroids[j] *= 1.0/numCopies;
double ke = 0.0;
for (int i = 0; i < numCopies; i++) {
const vector<Vec3>& pos = states[i].getPositions();
const vector<Vec3>& f = states[i].getForces();
for (int j = 0; j < numParticles; j++) {
Vec3 delta = centroids[j]-pos[j];
ke += delta.dot(f[j]);
}
}
meanKE += ke/(2*numCopies*numParticles);
}
// Check against expected values.
double scale = (boxSize*boxSize*boxSize)/(numSteps*0.5*numParticles*numParticles);
for (int i = 0; i < numBins/2; i++) {
double r1 = i*boxSize/numBins;
double r2 = (i+1)*boxSize/numBins;
double volume = (4.0/3.0)*M_PI*(r2*r2*r2-r1*r1*r1);
ASSERT_USUALLY_EQUAL_TOL(reference[i], scale*counts[i]/volume, 0.1);
}
meanKE /= numSteps*BOLTZ;
ASSERT_USUALLY_EQUAL_TOL(60.0, 1.5*temperature+meanKE, 0.02);
}
Vec3 calcCM(const vector<Vec3>& values, System& system) {
Vec3 cm;
for (int j = 0; j < system.getNumParticles(); ++j) {
cm[0] += values[j][0]*system.getParticleMass(j);
cm[1] += values[j][1]*system.getParticleMass(j);
cm[2] += values[j][2]*system.getParticleMass(j);
}
return cm;
}
void testCMMotionRemoval() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.addForce(new CMMotionRemover());
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
Vec3 cmPos = calcCM(positions, system);
for (int j = 0; j < numParticles; j++)
positions[j] -= cmPos*(1/(mass*numParticles));
integ.setPositions(i, positions);
}
// Make sure the CMMotionRemover is getting applied.
for (int i = 0; i < 200; ++i) {
integ.step(1);
Vec3 pos;
for (int j = 0; j < numCopies; j++) {
State state = integ.getState(0, State::Positions | State::Velocities);
pos += calcCM(state.getPositions(), system);
}
pos *= 1.0/numCopies;
ASSERT_EQUAL_VEC(Vec3(0,0,0), pos, 0.5);
}
}
void testVirtualSites() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*3;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
}
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[3*index] = pos;
positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
positions[3*index+2] = Vec3();
}
integ.setPositions(copy, positions);
}
// Check the temperature and virtual site locations.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++) {
state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
const vector<Vec3>& pos = state[j].getPositions();
for (int k = 0; k < numMolecules; k++)
ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
}
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testContractions() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
map<int, int> contractions;
contractions[1] = 3;
contractions[2] = 1;
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001, contractions);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testWithoutThermostat() {
const int numParticles = 20;
const int numCopies = 10;
const double temperature = 300.0;
const double mass = 2.0;
// Create a chain of particles.
System system;
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
for (int i = 0; i < numParticles; i++) {
system.addParticle(mass);
if (i > 0)
bonds->addBond(i-1, i, 1.0, 1000.0);
}
RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
integ.setApplyThermostat(false);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<vector<Vec3> > positions(numCopies);
for (int i = 0; i < numCopies; i++) {
positions[i].resize(numParticles);
for (int j = 0; j < numParticles; j++)
positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
integ.setPositions(i, positions[i]);
}
// Simulate it and see if the energy remains constant.
double initialEnergy;
int numSteps = 100;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
double energy = integ.getTotalEnergy();
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
}
}
void testWithBarostat() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
Platform& platform = Platform::getPlatformByName("CUDA");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(100);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
void runPlatformTests() {
testParaHydrogen();
}
int main(int argc, char* argv[]) {
try {
registerRPMDCudaKernelFactories();
if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("Precision", string(argv[1]));
testFreeParticles();
testParaHydrogen();
testCMMotionRemoval();
testVirtualSites();
testContractions();
testWithoutThermostat();
testWithBarostat();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
void setupKernels (int argc, char* argv[]) {
registerRPMDCudaKernelFactories();
platform = dynamic_cast<CudaPlatform&>(Platform::getPlatformByName("CUDA"));
initializeTests(argc, argv);
}
plugins/rpmd/platforms/opencl/CMakeLists.txt
View file @ be19e022
......@@ -12,7 +12,7 @@
# The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
SET(OPENMM_SOURCE_SUBDIRS .)
SET(OPENMM_SOURCE_SUBDIRS . ../common)
# Collect up information about the version of the OpenMM library we're building
......@@ -59,32 +59,25 @@ FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
ENDFOREACH(subdir)
SET(COMMON_KERNELS_CPP ${CMAKE_CURRENT_BINARY_DIR}/../common/src/CommonRpmdKernelSources.cpp)
SET(SOURCE_FILES ${SOURCE_FILES} ${COMMON_KERNELS_CPP})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/../common/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/opencl/include)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/opencl/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_BINARY_DIR}/platforms/opencl/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/common/include)
# Set variables needed for encoding kernel sources into a C++ class
SET(KERNEL_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
SET(KERNEL_SOURCE_CLASS OpenCLRpmdKernelSources)
SET(KERNELS_CPP ${CMAKE_CURRENT_BINARY_DIR}/src/${KERNEL_SOURCE_CLASS}.cpp)
SET(KERNELS_H ${CMAKE_CURRENT_BINARY_DIR}/src/${KERNEL_SOURCE_CLASS}.h)
SET(SOURCE_FILES ${SOURCE_FILES} ${KERNELS_CPP} ${KERNELS_H})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_BINARY_DIR}/platforms/common/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}/../common/src)
# Create the library
INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR})
FILE(GLOB OPENCL_KERNELS ${KERNEL_SOURCE_DIR}/kernels/*.cl)
ADD_CUSTOM_COMMAND(OUTPUT ${KERNELS_CPP} ${KERNELS_H}
COMMAND ${CMAKE_COMMAND}
ARGS -D KERNEL_SOURCE_DIR=${KERNEL_SOURCE_DIR} -D KERNELS_CPP=${KERNELS_CPP} -D KERNELS_H=${KERNELS_H} -D KERNEL_SOURCE_CLASS=${KERNEL_SOURCE_CLASS} -D KERNEL_FILE_EXTENSION=cl -P ${CMAKE_SOURCE_DIR}/cmake_modules/EncodeKernelFiles.cmake
DEPENDS ${OPENCL_KERNELS}
)
SET_SOURCE_FILES_PROPERTIES(${KERNELS_CPP} ${KERNELS_H} PROPERTIES GENERATED TRUE)
SET_SOURCE_FILES_PROPERTIES(${COMMON_KERNELS_CPP} PROPERTIES GENERATED TRUE)
ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
ADD_DEPENDENCIES(${SHARED_TARGET} RpmdCommonKernels)
TARGET_LINK_LIBRARIES(${SHARED_TARGET} ${OPENMM_LIBRARY_NAME} ${OPENCL_LIBRARIES} ${PTHREADS_LIB})
TARGET_LINK_LIBRARIES(${SHARED_TARGET} ${OPENMM_LIBRARY_NAME}OpenCL)
......
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernelFactory.cpp
View file @ be19e022
......@@ -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) 2011 Stanford University and the Authors. *
* Portions copyright (c) 2011-2021 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -27,7 +27,8 @@
#include <exception>
#include "OpenCLRpmdKernelFactory.h"
#include "OpenCLRpmdKernels.h"
#include "CommonRpmdKernels.h"
#include "OpenCLContext.h"
#include "openmm/internal/windowsExportRpmd.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/OpenMMException.h"
......@@ -61,6 +62,6 @@ extern "C" OPENMM_EXPORT void registerRPMDOpenCLKernelFactories() {
KernelImpl* OpenCLRpmdKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
OpenCLContext& cl = *static_cast<OpenCLPlatform::PlatformData*>(context.getPlatformData())->contexts[0];
if (name == IntegrateRPMDStepKernel::Name())
return new OpenCLIntegrateRPMDStepKernel(name, platform, cl);
return new CommonIntegrateRPMDStepKernel(name, platform, cl);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
}
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernelSources.cpp.in deleted 100644 → 0
View file @ 98d81730
/* -------------------------------------------------------------------------- *
* 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) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "OpenCLRpmdKernelSources.h"
using namespace OpenMM;
using namespace std;
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernelSources.h.in deleted 100644 → 0
View file @ 98d81730
#ifndef OPENMM_OPENCLRPMDKERNELSOURCES_H_
#define OPENMM_OPENCLRPMDKERNELSOURCES_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include <string>
namespace OpenMM {
/**
* This class is a central holding place for the source code of OpenCL kernels.
* The CMake build script inserts declarations into it based on the .cl files in the
* kernels subfolder.
*/
class OpenCLRpmdKernelSources {
public:
@KERNEL_FILE_DECLARATIONS@
};
} // namespace OpenMM
#endif /*OPENMM_OPENCLRPMDKERNELSOURCES_H_*/
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp deleted 100644 → 0
View file @ 98d81730
/* -------------------------------------------------------------------------- *
* 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) 2011-2020 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. *
* -------------------------------------------------------------------------- */
#include "OpenCLRpmdKernels.h"
#include "OpenCLRpmdKernelSources.h"
#include "openmm/internal/ContextImpl.h"
#include "OpenCLIntegrationUtilities.h"
#include "OpenCLExpressionUtilities.h"
#include "OpenCLFFT3D.h"
#include "OpenCLNonbondedUtilities.h"
#include "SimTKOpenMMRealType.h"
using namespace OpenMM;
using namespace std;
void OpenCLIntegrateRPMDStepKernel::initialize(const System& system, const RPMDIntegrator& integrator) {
cl.getPlatformData().initializeContexts(system);
numCopies = integrator.getNumCopies();
numParticles = system.getNumParticles();
workgroupSize = numCopies;
if (numCopies != OpenCLFFT3D::findLegalDimension(numCopies))
throw OpenMMException("RPMDIntegrator: the number of copies must be a multiple of powers of 2, 3, and 5.");
int paddedParticles = cl.getPaddedNumAtoms();
int forceElementSize = (cl.getUseDoublePrecision() ? sizeof(mm_double4) : sizeof(mm_float4));
forces.initialize(cl, numCopies*paddedParticles, forceElementSize, "rpmdForces");
bool useDoublePrecision = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision());
int elementSize = (useDoublePrecision ? sizeof(mm_double4) : sizeof(mm_float4));
positions.initialize(cl, numCopies*paddedParticles, elementSize, "rpmdPositions");
velocities.initialize(cl, numCopies*paddedParticles, elementSize, "rpmdVelocities");
cl.getIntegrationUtilities().initRandomNumberGenerator((unsigned int) integrator.getRandomNumberSeed());
// Fill in the posq and velm arrays with safe values to avoid a risk of nans.
if (useDoublePrecision) {
vector<mm_double4> temp(positions.getSize());
for (int i = 0; i < positions.getSize(); i++)
temp[i] = mm_double4(0, 0, 0, 0);
positions.upload(temp);
for (int i = 0; i < velocities.getSize(); i++)
temp[i] = mm_double4(0, 0, 0, 1);
velocities.upload(temp);
}
else {
vector<mm_float4> temp(positions.getSize());
for (int i = 0; i < positions.getSize(); i++)
temp[i] = mm_float4(0, 0, 0, 0);
positions.upload(temp);
for (int i = 0; i < velocities.getSize(); i++)
temp[i] = mm_float4(0, 0, 0, 1);
velocities.upload(temp);
}
// Build a list of contractions.
groupsNotContracted = -1;
const map<int, int>& contractions = integrator.getContractions();
int maxContractedCopies = 0;
for (auto& c : contractions) {
int group = c.first;
int copies = c.second;
if (group < 0 || group > 31)
throw OpenMMException("RPMDIntegrator: Force group must be between 0 and 31");
if (copies < 0 || copies > numCopies)
throw OpenMMException("RPMDIntegrator: Number of copies for contraction cannot be greater than the total number of copies being simulated");
if (copies != OpenCLFFT3D::findLegalDimension(copies))
throw OpenMMException("RPMDIntegrator: Number of copies for contraction must be a multiple of powers of 2, 3, and 5.");
if (copies != numCopies) {
if (groupsByCopies.find(copies) == groupsByCopies.end()) {
groupsByCopies[copies] = 1<<group;
if (copies > maxContractedCopies)
maxContractedCopies = copies;
}
else
groupsByCopies[copies] |= 1<<group;
groupsNotContracted -= 1<<group;
}
}
groupsNotContracted &= integrator.getIntegrationForceGroups();
if (maxContractedCopies > 0) {
contractedForces.initialize(cl, maxContractedCopies*paddedParticles, forceElementSize, "rpmdContractedForces");
contractedPositions.initialize(cl, maxContractedCopies*paddedParticles, elementSize, "rpmdContractedPositions");
}
// Create kernels.
map<string, string> defines;
defines["NUM_ATOMS"] = cl.intToString(cl.getNumAtoms());
defines["PADDED_NUM_ATOMS"] = cl.intToString(cl.getPaddedNumAtoms());
defines["NUM_COPIES"] = cl.intToString(numCopies);
defines["THREAD_BLOCK_SIZE"] = cl.intToString(workgroupSize);
defines["HBAR"] = cl.doubleToString(1.054571628e-34*AVOGADRO/(1000*1e-12));
defines["SCALE"] = cl.doubleToString(1.0/sqrt((double) numCopies));
defines["M_PI"] = cl.doubleToString(M_PI);
map<string, string> replacements;
replacements["FFT_Q_FORWARD"] = createFFT(numCopies, "q", true);
replacements["FFT_Q_BACKWARD"] = createFFT(numCopies, "q", false);
replacements["FFT_V_FORWARD"] = createFFT(numCopies, "v", true);
replacements["FFT_V_BACKWARD"] = createFFT(numCopies, "v", false);
cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLRpmdKernelSources::rpmd, replacements), defines, "");
pileKernel = cl::Kernel(program, "applyPileThermostat");
stepKernel = cl::Kernel(program, "integrateStep");
velocitiesKernel = cl::Kernel(program, "advanceVelocities");
copyToContextKernel = cl::Kernel(program, "copyDataToContext");
copyFromContextKernel = cl::Kernel(program, "copyDataFromContext");
translateKernel = cl::Kernel(program, "applyCellTranslations");
// Create kernels for doing contractions.
for (auto& g : groupsByCopies) {
int copies = g.first;
replacements.clear();
replacements["NUM_CONTRACTED_COPIES"] = cl.intToString(copies);
replacements["POS_SCALE"] = cl.doubleToString(1.0/numCopies);
replacements["FORCE_SCALE"] = cl.doubleToString(1.0/copies);
replacements["FFT_Q_FORWARD"] = createFFT(numCopies, "q", true);
replacements["FFT_Q_BACKWARD"] = createFFT(copies, "q", false);
replacements["FFT_F_FORWARD"] = createFFT(copies, "f", true);
replacements["FFT_F_BACKWARD"] = createFFT(numCopies, "f", false);
program = cl.createProgram(cl.replaceStrings(OpenCLRpmdKernelSources::rpmdContraction, replacements), defines, "");
positionContractionKernels[copies] = cl::Kernel(program, "contractPositions");
forceContractionKernels[copies] = cl::Kernel(program, "contractForces");
}
}
void OpenCLIntegrateRPMDStepKernel::initializeKernels(ContextImpl& context) {
hasInitializedKernel = true;
pileKernel.setArg<cl::Buffer>(0, velocities.getDeviceBuffer());
stepKernel.setArg<cl::Buffer>(0, positions.getDeviceBuffer());
stepKernel.setArg<cl::Buffer>(1, velocities.getDeviceBuffer());
stepKernel.setArg<cl::Buffer>(2, forces.getDeviceBuffer());
velocitiesKernel.setArg<cl::Buffer>(0, velocities.getDeviceBuffer());
velocitiesKernel.setArg<cl::Buffer>(1, forces.getDeviceBuffer());
translateKernel.setArg<cl::Buffer>(0, positions.getDeviceBuffer());
translateKernel.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
translateKernel.setArg<cl::Buffer>(2, cl.getAtomIndexArray().getDeviceBuffer());
copyToContextKernel.setArg<cl::Buffer>(0, velocities.getDeviceBuffer());
copyToContextKernel.setArg<cl::Buffer>(1, cl.getVelm().getDeviceBuffer());
copyToContextKernel.setArg<cl::Buffer>(3, cl.getPosq().getDeviceBuffer());
copyToContextKernel.setArg<cl::Buffer>(4, cl.getAtomIndexArray().getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(0, cl.getForce().getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(3, velocities.getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(6, cl.getAtomIndexArray().getDeviceBuffer());
for (auto& g : groupsByCopies) {
int copies = g.first;
positionContractionKernels[copies].setArg<cl::Buffer>(0, positions.getDeviceBuffer());
positionContractionKernels[copies].setArg<cl::Buffer>(1, contractedPositions.getDeviceBuffer());
forceContractionKernels[copies].setArg<cl::Buffer>(0, forces.getDeviceBuffer());
forceContractionKernels[copies].setArg<cl::Buffer>(1, contractedForces.getDeviceBuffer());
}
}
void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator, bool forcesAreValid) {
OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
if (!hasInitializedKernel)
initializeKernels(context);
// Loop over copies and compute the force on each one.
if (!forcesAreValid)
computeForces(context);
// Apply the PILE-L thermostat.
bool useDoublePrecision = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision());
const double dt = integrator.getStepSize();
pileKernel.setArg<cl_uint>(2, integration.prepareRandomNumbers(numParticles*numCopies));
pileKernel.setArg<cl::Buffer>(1, integration.getRandom().getDeviceBuffer()); // Do this *after* prepareRandomNumbers(), which might rebuild the array.
if (useDoublePrecision) {
pileKernel.setArg<cl_double>(3, dt);
pileKernel.setArg<cl_double>(4, integrator.getTemperature()*BOLTZ);
pileKernel.setArg<cl_double>(5, integrator.getFriction());
stepKernel.setArg<cl_double>(3, dt);
stepKernel.setArg<cl_double>(4, integrator.getTemperature()*BOLTZ);
velocitiesKernel.setArg<cl_double>(2, dt);
}
else {
pileKernel.setArg<cl_float>(3, (cl_float) dt);
pileKernel.setArg<cl_float>(4, (cl_float) (integrator.getTemperature()*BOLTZ));
pileKernel.setArg<cl_float>(5, (cl_float) integrator.getFriction());
stepKernel.setArg<cl_float>(3, (cl_float) dt);
stepKernel.setArg<cl_float>(4, (cl_float) (integrator.getTemperature()*BOLTZ));
velocitiesKernel.setArg<cl_float>(2, (cl_float) dt);
}
if (integrator.getApplyThermostat())
cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
// Update positions and velocities.
cl.executeKernel(stepKernel, numParticles*numCopies, workgroupSize);
// Calculate forces based on the updated positions.
computeForces(context);
// Update velocities.
cl.executeKernel(velocitiesKernel, numParticles*numCopies, workgroupSize);
// Apply the PILE-L thermostat again.
if (integrator.getApplyThermostat()) {
pileKernel.setArg<cl_uint>(2, integration.prepareRandomNumbers(numParticles*numCopies));
cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
}
// Update the time and step count.
cl.setTime(cl.getTime()+dt);
cl.setStepCount(cl.getStepCount()+1);
cl.reorderAtoms();
if (cl.getAtomsWereReordered() && cl.getNonbondedUtilities().getUsePeriodic()) {
// Atoms may have been translated into a different periodic box, so apply
// the same translation to all the beads.
translateKernel.setArg<cl_int>(3, numCopies-1);
cl.executeKernel(translateKernel, cl.getNumAtoms());
}
}
void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
// Compute forces from all groups that didn't have a specified contraction.
copyToContextKernel.setArg<cl::Buffer>(2, positions.getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(1, forces.getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(5, positions.getDeviceBuffer());
for (int i = 0; i < numCopies; i++) {
copyToContextKernel.setArg<cl_int>(5, i);
cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
context.computeVirtualSites();
Vec3 initialBox[3];
context.getPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
context.updateContextState();
Vec3 finalBox[3];
context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2])
throw OpenMMException("Standard barostats cannot be used with RPMDIntegrator. Use RPMDMonteCarloBarostat instead.");
context.calcForcesAndEnergy(true, false, groupsNotContracted);
copyFromContextKernel.setArg<cl_int>(7, i);
cl.executeKernel(copyFromContextKernel, cl.getNumAtoms());
}
// Now loop over contractions and compute forces from them.
if (groupsByCopies.size() > 0) {
copyToContextKernel.setArg<cl::Buffer>(2, contractedPositions.getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(1, contractedForces.getDeviceBuffer());
copyFromContextKernel.setArg<cl::Buffer>(5, contractedPositions.getDeviceBuffer());
for (auto& g : groupsByCopies) {
int copies = g.first;
int groupFlags = g.second;
// Find the contracted positions.
cl.executeKernel(positionContractionKernels[copies], numParticles*numCopies, workgroupSize);
// Compute forces.
for (int i = 0; i < copies; i++) {
copyToContextKernel.setArg<cl_int>(5, i);
cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
context.computeVirtualSites();
context.calcForcesAndEnergy(true, false, groupFlags);
copyFromContextKernel.setArg<cl_int>(7, i);
cl.executeKernel(copyFromContextKernel, cl.getNumAtoms());
}
// Apply the forces to the original copies.
cl.executeKernel(forceContractionKernels[copies], numParticles*numCopies, workgroupSize);
}
}
if (groupsByCopies.size() > 0) {
// Ensure the Context contains the positions from the last copy, since we'll assume that later.
copyToContextKernel.setArg<cl::Buffer>(2, positions.getDeviceBuffer());
copyToContextKernel.setArg<cl_int>(5, numCopies-1);
cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
}
}
double OpenCLIntegrateRPMDStepKernel::computeKineticEnergy(ContextImpl& context, const RPMDIntegrator& integrator) {
return cl.getIntegrationUtilities().computeKineticEnergy(0);
}
void OpenCLIntegrateRPMDStepKernel::setPositions(int copy, const vector<Vec3>& pos) {
if (!positions.isInitialized())
throw OpenMMException("RPMDIntegrator: Cannot set positions before the integrator is added to a Context");
if (pos.size() != numParticles)
throw OpenMMException("RPMDIntegrator: wrong number of values passed to setPositions()");
// Adjust the positions based on the current cell offsets.
const vector<int>& order = cl.getAtomIndex();
mm_double4 periodicBoxSize = cl.getPeriodicBoxSizeDouble();
vector<Vec3> offsetPos(numParticles);
for (int i = 0; i < numParticles; ++i) {
mm_int4 offset = cl.getPosCellOffsets()[i];
offsetPos[order[i]] = pos[order[i]] + Vec3(offset.x*periodicBoxSize.x, offset.y*periodicBoxSize.y, offset.z*periodicBoxSize.z);
}
// Record the positions.
if (cl.getUseDoublePrecision()) {
vector<mm_double4> posq(cl.getPaddedNumAtoms());
cl.getPosq().download(posq);
for (int i = 0; i < numParticles; i++)
posq[i] = mm_double4(offsetPos[i][0], offsetPos[i][1], offsetPos[i][2], posq[i].w);
cl.getQueue().enqueueWriteBuffer(positions.getDeviceBuffer(), CL_TRUE, copy*cl.getPaddedNumAtoms()*sizeof(mm_double4), numParticles*sizeof(mm_double4), &posq[0]);
}
else if (cl.getUseMixedPrecision()) {
vector<mm_float4> posqf(cl.getPaddedNumAtoms());
cl.getPosq().download(posqf);
vector<mm_double4> posq(cl.getPaddedNumAtoms());
for (int i = 0; i < numParticles; i++)
posq[i] = mm_double4(offsetPos[i][0], offsetPos[i][1], offsetPos[i][2], posqf[i].w);
cl.getQueue().enqueueWriteBuffer(positions.getDeviceBuffer(), CL_TRUE, copy*cl.getPaddedNumAtoms()*sizeof(mm_double4), numParticles*sizeof(mm_double4), &posq[0]);
}
else {
vector<mm_float4> posq(cl.getPaddedNumAtoms());
cl.getPosq().download(posq);
for (int i = 0; i < numParticles; i++)
posq[i] = mm_float4((cl_float) offsetPos[i][0], (cl_float) offsetPos[i][1], (cl_float) offsetPos[i][2], posq[i].w);
cl.getQueue().enqueueWriteBuffer(positions.getDeviceBuffer(), CL_TRUE, copy*cl.getPaddedNumAtoms()*sizeof(mm_float4), numParticles*sizeof(mm_float4), &posq[0]);
}
}
void OpenCLIntegrateRPMDStepKernel::setVelocities(int copy, const vector<Vec3>& vel) {
if (!velocities.isInitialized())
throw OpenMMException("RPMDIntegrator: Cannot set velocities before the integrator is added to a Context");
if (vel.size() != numParticles)
throw OpenMMException("RPMDIntegrator: wrong number of values passed to setVelocities()");
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
vector<mm_double4> velm(cl.getPaddedNumAtoms());
cl.getVelm().download(velm);
for (int i = 0; i < numParticles; i++)
velm[i] = mm_double4(vel[i][0], vel[i][1], vel[i][2], velm[i].w);
cl.getQueue().enqueueWriteBuffer(velocities.getDeviceBuffer(), CL_TRUE, copy*cl.getPaddedNumAtoms()*sizeof(mm_double4), numParticles*sizeof(mm_double4), &velm[0]);
}
else {
vector<mm_float4> velm(cl.getPaddedNumAtoms());
cl.getVelm().download(velm);
for (int i = 0; i < numParticles; i++)
velm[i] = mm_float4((cl_float) vel[i][0], (cl_float) vel[i][1], (cl_float) vel[i][2], velm[i].w);
cl.getQueue().enqueueWriteBuffer(velocities.getDeviceBuffer(), CL_TRUE, copy*cl.getPaddedNumAtoms()*sizeof(mm_float4), numParticles*sizeof(mm_float4), &velm[0]);
}
}
void OpenCLIntegrateRPMDStepKernel::copyToContext(int copy, ContextImpl& context) {
if (!hasInitializedKernel)
initializeKernels(context);
copyToContextKernel.setArg<cl::Buffer>(2, positions.getDeviceBuffer());
copyToContextKernel.setArg<cl_int>(5, copy);
cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
}
string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable, bool forward) {
stringstream source;
int stage = 0;
int L = size;
int m = 1;
string sign = (forward ? "1.0f" : "-1.0f");
string multReal = (forward ? "multiplyComplexRealPart" : "multiplyComplexRealPartConj");
string multImag = (forward ? "multiplyComplexImagPart" : "multiplyComplexImagPartConj");
source<<"{\n";
source<<"__local mixed4* real0 = "<<variable<<"real;\n";
source<<"__local mixed4* imag0 = "<<variable<<"imag;\n";
source<<"__local mixed4* real1 = &temp[blockStart];\n";
source<<"__local mixed4* imag1 = &temp[blockStart+get_local_size(0)];\n";
// Factor size, generating an appropriate block of code for each factor.
while (L > 1) {
int input = stage%2;
int output = 1-input;
int radix;
if (L%5 == 0)
radix = 5;
else if (L%4 == 0)
radix = 4;
else if (L%3 == 0)
radix = 3;
else if (L%2 == 0)
radix = 2;
else
throw OpenMMException("Illegal size for FFT: "+cl.intToString(size));
source<<"{\n";
L = L/radix;
source<<"// Pass "<<(stage+1)<<" (radix "<<radix<<")\n";
source<<"if (indexInBlock < "<<(L*m)<<") {\n";
source<<"int i = indexInBlock;\n";
source<<"int j = i/"<<m<<";\n";
if (radix == 5) {
source<<"mixed4 c0r = real"<<input<<"[i];\n";
source<<"mixed4 c0i = imag"<<input<<"[i];\n";
source<<"mixed4 c1r = real"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c1i = imag"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c2r = real"<<input<<"[i+"<<(2*L*m)<<"];\n";
source<<"mixed4 c2i = imag"<<input<<"[i+"<<(2*L*m)<<"];\n";
source<<"mixed4 c3r = real"<<input<<"[i+"<<(3*L*m)<<"];\n";
source<<"mixed4 c3i = imag"<<input<<"[i+"<<(3*L*m)<<"];\n";
source<<"mixed4 c4r = real"<<input<<"[i+"<<(4*L*m)<<"];\n";
source<<"mixed4 c4i = imag"<<input<<"[i+"<<(4*L*m)<<"];\n";
source<<"mixed4 d0r = c1r+c4r;\n";
source<<"mixed4 d0i = c1i+c4i;\n";
source<<"mixed4 d1r = c2r+c3r;\n";
source<<"mixed4 d1i = c2i+c3i;\n";
source<<"mixed4 d2r = "<<cl.doubleToString(sin(0.4*M_PI))<<"*(c1r-c4r);\n";
source<<"mixed4 d2i = "<<cl.doubleToString(sin(0.4*M_PI))<<"*(c1i-c4i);\n";
source<<"mixed4 d3r = "<<cl.doubleToString(sin(0.4*M_PI))<<"*(c2r-c3r);\n";
source<<"mixed4 d3i = "<<cl.doubleToString(sin(0.4*M_PI))<<"*(c2i-c3i);\n";
source<<"mixed4 d4r = d0r+d1r;\n";
source<<"mixed4 d4i = d0i+d1i;\n";
source<<"mixed4 d5r = "<<cl.doubleToString(0.25*sqrt(5.0))<<"*(d0r-d1r);\n";
source<<"mixed4 d5i = "<<cl.doubleToString(0.25*sqrt(5.0))<<"*(d0i-d1i);\n";
source<<"mixed4 d6r = c0r-0.25f*d4r;\n";
source<<"mixed4 d6i = c0i-0.25f*d4i;\n";
source<<"mixed4 d7r = d6r+d5r;\n";
source<<"mixed4 d7i = d6i+d5i;\n";
source<<"mixed4 d8r = d6r-d5r;\n";
source<<"mixed4 d8i = d6i-d5i;\n";
string coeff = cl.doubleToString(sin(0.2*M_PI)/sin(0.4*M_PI));
source<<"mixed4 d9r = "<<sign<<"*(d2i+"<<coeff<<"*d3i);\n";
source<<"mixed4 d9i = "<<sign<<"*(-d2r-"<<coeff<<"*d3r);\n";
source<<"mixed4 d10r = "<<sign<<"*("<<coeff<<"*d2i-d3i);\n";
source<<"mixed4 d10i = "<<sign<<"*(d3r-"<<coeff<<"*d2r);\n";
source<<"real"<<output<<"[i+4*j*"<<m<<"] = c0r+d4r;\n";
source<<"imag"<<output<<"[i+4*j*"<<m<<"] = c0i+d4i;\n";
source<<"real"<<output<<"[i+(4*j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(5*L)<<"], d7r+d9r, d7i+d9i);\n";
source<<"imag"<<output<<"[i+(4*j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(5*L)<<"], d7r+d9r, d7i+d9i);\n";
source<<"real"<<output<<"[i+(4*j+2)*"<<m<<"] = "<<multReal<<"(w[j*"<<(2*size)<<"/"<<(5*L)<<"], d8r+d10r, d8i+d10i);\n";
source<<"imag"<<output<<"[i+(4*j+2)*"<<m<<"] = "<<multImag<<"(w[j*"<<(2*size)<<"/"<<(5*L)<<"], d8r+d10r, d8i+d10i);\n";
source<<"real"<<output<<"[i+(4*j+3)*"<<m<<"] = "<<multReal<<"(w[j*"<<(3*size)<<"/"<<(5*L)<<"], d8r-d10r, d8i-d10i);\n";
source<<"imag"<<output<<"[i+(4*j+3)*"<<m<<"] = "<<multImag<<"(w[j*"<<(3*size)<<"/"<<(5*L)<<"], d8r-d10r, d8i-d10i);\n";
source<<"real"<<output<<"[i+(4*j+4)*"<<m<<"] = "<<multReal<<"(w[j*"<<(4*size)<<"/"<<(5*L)<<"], d7r-d9r, d7i-d9i);\n";
source<<"imag"<<output<<"[i+(4*j+4)*"<<m<<"] = "<<multImag<<"(w[j*"<<(4*size)<<"/"<<(5*L)<<"], d7r-d9r, d7i-d9i);\n";
}
else if (radix == 4) {
source<<"mixed4 c0r = real"<<input<<"[i];\n";
source<<"mixed4 c0i = imag"<<input<<"[i];\n";
source<<"mixed4 c1r = real"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c1i = imag"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c2r = real"<<input<<"[i+"<<(2*L*m)<<"];\n";
source<<"mixed4 c2i = imag"<<input<<"[i+"<<(2*L*m)<<"];\n";
source<<"mixed4 c3r = real"<<input<<"[i+"<<(3*L*m)<<"];\n";
source<<"mixed4 c3i = imag"<<input<<"[i+"<<(3*L*m)<<"];\n";
source<<"mixed4 d0r = c0r+c2r;\n";
source<<"mixed4 d0i = c0i+c2i;\n";
source<<"mixed4 d1r = c0r-c2r;\n";
source<<"mixed4 d1i = c0i-c2i;\n";
source<<"mixed4 d2r = c1r+c3r;\n";
source<<"mixed4 d2i = c1i+c3i;\n";
source<<"mixed4 d3r = "<<sign<<"*(c1i-c3i);\n";
source<<"mixed4 d3i = "<<sign<<"*(c3r-c1r);\n";
source<<"real"<<output<<"[i+3*j*"<<m<<"] = d0r+d2r;\n";
source<<"imag"<<output<<"[i+3*j*"<<m<<"] = d0i+d2i;\n";
source<<"real"<<output<<"[i+(3*j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(4*L)<<"], d1r+d3r, d1i+d3i);\n";
source<<"imag"<<output<<"[i+(3*j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(4*L)<<"], d1r+d3r, d1i+d3i);\n";
source<<"real"<<output<<"[i+(3*j+2)*"<<m<<"] = "<<multReal<<"(w[j*"<<(2*size)<<"/"<<(4*L)<<"], d0r-d2r, d0i-d2i);\n";
source<<"imag"<<output<<"[i+(3*j+2)*"<<m<<"] = "<<multImag<<"(w[j*"<<(2*size)<<"/"<<(4*L)<<"], d0r-d2r, d0i-d2i);\n";
source<<"real"<<output<<"[i+(3*j+3)*"<<m<<"] = "<<multReal<<"(w[j*"<<(3*size)<<"/"<<(4*L)<<"], d1r-d3r, d1i-d3i);\n";
source<<"imag"<<output<<"[i+(3*j+3)*"<<m<<"] = "<<multImag<<"(w[j*"<<(3*size)<<"/"<<(4*L)<<"], d1r-d3r, d1i-d3i);\n";
}
else if (radix == 3) {
source<<"mixed4 c0r = real"<<input<<"[i];\n";
source<<"mixed4 c0i = imag"<<input<<"[i];\n";
source<<"mixed4 c1r = real"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c1i = imag"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c2r = real"<<input<<"[i+"<<(2*L*m)<<"];\n";
source<<"mixed4 c2i = imag"<<input<<"[i+"<<(2*L*m)<<"];\n";
source<<"mixed4 d0r = c1r+c2r;\n";
source<<"mixed4 d0i = c1i+c2i;\n";
source<<"mixed4 d1r = c0r-0.5f*d0r;\n";
source<<"mixed4 d1i = c0i-0.5f*d0i;\n";
source<<"mixed4 d2r = "<<sign<<"*"<<cl.doubleToString(sin(M_PI/3.0))<<"*(c1i-c2i);\n";
source<<"mixed4 d2i = "<<sign<<"*"<<cl.doubleToString(sin(M_PI/3.0))<<"*(c2r-c1r);\n";
source<<"real"<<output<<"[i+2*j*"<<m<<"] = c0r+d0r;\n";
source<<"imag"<<output<<"[i+2*j*"<<m<<"] = c0i+d0i;\n";
source<<"real"<<output<<"[i+(2*j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(3*L)<<"], d1r+d2r, d1i+d2i);\n";
source<<"imag"<<output<<"[i+(2*j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(3*L)<<"], d1r+d2r, d1i+d2i);\n";
source<<"real"<<output<<"[i+(2*j+2)*"<<m<<"] = "<<multReal<<"(w[j*"<<(2*size)<<"/"<<(3*L)<<"], d1r-d2r, d1i-d2i);\n";
source<<"imag"<<output<<"[i+(2*j+2)*"<<m<<"] = "<<multImag<<"(w[j*"<<(2*size)<<"/"<<(3*L)<<"], d1r-d2r, d1i-d2i);\n";
}
else if (radix == 2) {
source<<"mixed4 c0r = real"<<input<<"[i];\n";
source<<"mixed4 c0i = imag"<<input<<"[i];\n";
source<<"mixed4 c1r = real"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"mixed4 c1i = imag"<<input<<"[i+"<<(L*m)<<"];\n";
source<<"real"<<output<<"[i+j*"<<m<<"] = c0r+c1r;\n";
source<<"imag"<<output<<"[i+j*"<<m<<"] = c0i+c1i;\n";
source<<"real"<<output<<"[i+(j+1)*"<<m<<"] = "<<multReal<<"(w[j*"<<size<<"/"<<(2*L)<<"], c0r-c1r, c0i-c1i);\n";
source<<"imag"<<output<<"[i+(j+1)*"<<m<<"] = "<<multImag<<"(w[j*"<<size<<"/"<<(2*L)<<"], c0r-c1r, c0i-c1i);\n";
}
source<<"}\n";
m = m*radix;
source<<"barrier(CLK_LOCAL_MEM_FENCE);\n";
source<<"}\n";
++stage;
}
// Create the kernel.
if (stage%2 == 1) {
source<<"real0[indexInBlock] = real1[indexInBlock];\n";
source<<"imag0[indexInBlock] = imag1[indexInBlock];\n";
source<<"barrier(CLK_LOCAL_MEM_FENCE);\n";
}
source<<"}\n";
return source.str();
}
plugins/rpmd/platforms/opencl/src/kernels/rpmd.cl deleted 100644 → 0
View file @ 98d81730
mixed4 multiplyComplexRealPart(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2r-c1.y*c2i;
}
mixed4 multiplyComplexImagPart(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2i+c1.y*c2r;
}
mixed4 multiplyComplexRealPartConj(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2r+c1.y*c2i;
}
mixed4 multiplyComplexImagPartConj(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2i-c1.y*c2r;
}
/**
* Apply the PILE-L thermostat.
*/
__kernel void applyPileThermostat(__global mixed4* velm, __global float4* random, unsigned int randomIndex,
mixed dt, mixed kT, mixed friction) {
const int numBlocks = get_global_size(0)/NUM_COPIES;
const int blockStart = NUM_COPIES*(get_local_id(0)/NUM_COPIES);
const int indexInBlock = get_local_id(0)-blockStart;
const mixed nkT = NUM_COPIES*kT;
const mixed twown = 2.0f*nkT/HBAR;
const mixed c1_0 = exp(-0.5f*dt*friction);
const mixed c2_0 = sqrt(1.0f-c1_0*c1_0);
__local mixed4 v[2*THREAD_BLOCK_SIZE];
__local mixed4 temp[2*THREAD_BLOCK_SIZE];
__local mixed2 w[NUM_COPIES];
__local mixed4* vreal = &v[blockStart];
__local mixed4* vimag = &v[blockStart+get_local_size(0)];
if (get_local_id(0) < NUM_COPIES)
w[indexInBlock] = (mixed2) (cos(-indexInBlock*2*M_PI/NUM_COPIES), sin(-indexInBlock*2*M_PI/NUM_COPIES));
barrier(CLK_LOCAL_MEM_FENCE);
randomIndex += NUM_COPIES*(get_global_id(0)/NUM_COPIES);
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
mixed4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
mixed invMass = particleVelm.w;
mixed c3_0 = c2_0*sqrt(nkT*invMass);
// Forward FFT.
vreal[indexInBlock] = SCALE*particleVelm;
vimag[indexInBlock] = (mixed4) (0.0f, 0.0f, 0.0f, 0.0f);
barrier(CLK_GLOBAL_MEM_FENCE);
FFT_V_FORWARD
// Apply the thermostat.
if (indexInBlock == 0) {
// Apply a local Langevin thermostat to the centroid mode.
vreal[0].xyz = vreal[0].xyz*c1_0 + c3_0*convert_mixed4(random[randomIndex]).xyz;
}
else {
// Use critical damping white noise for the remaining modes.
int k = (indexInBlock <= NUM_COPIES/2 ? indexInBlock : NUM_COPIES-indexInBlock);
const bool isCenter = (NUM_COPIES%2 == 0 && k == NUM_COPIES/2);
const mixed wk = twown*sin(k*M_PI/NUM_COPIES);
const mixed c1 = exp(-wk*dt);
const mixed c2 = sqrt((1.0f-c1*c1)/2.0f) * (isCenter ? sqrt(2.0f) : 1.0f);
const mixed c3 = c2*sqrt(nkT*invMass);
mixed4 rand1 = c3*convert_mixed4(random[randomIndex+k]);
mixed4 rand2 = (isCenter ? 0.0f : c3*convert_mixed4(random[randomIndex+NUM_COPIES-k]));
vreal[indexInBlock].xyz = c1*vreal[indexInBlock].xyz + rand1.xyz;
vimag[indexInBlock].xyz = c1*vimag[indexInBlock].xyz + (indexInBlock < NUM_COPIES/2 ? rand2.xyz : -rand2.xyz);
}
barrier(CLK_GLOBAL_MEM_FENCE);
// Inverse FFT.
FFT_V_BACKWARD
if (invMass != 0)
velm[particle+indexInBlock*PADDED_NUM_ATOMS].xyz = SCALE*vreal[indexInBlock].xyz;
randomIndex += get_global_size(0);
}
}
/**
* Advance the positions and velocities.
*/
__kernel void integrateStep(__global mixed4* posq, __global mixed4* velm, __global real4* force, mixed dt, mixed kT) {
const int numBlocks = get_global_size(0)/NUM_COPIES;
const int blockStart = NUM_COPIES*(get_local_id(0)/NUM_COPIES);
const int indexInBlock = get_local_id(0)-blockStart;
const mixed nkT = NUM_COPIES*kT;
const mixed twown = 2.0f*nkT/HBAR;
__local mixed4 q[2*THREAD_BLOCK_SIZE];
__local mixed4 v[2*THREAD_BLOCK_SIZE];
__local mixed4 temp[2*THREAD_BLOCK_SIZE];
__local mixed2 w[NUM_COPIES];
// Update velocities.
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
int index = particle+indexInBlock*PADDED_NUM_ATOMS;
mixed4 particleVelm = velm[index];
particleVelm.xyz += convert_mixed4(force[index]).xyz*(0.5f*dt*particleVelm.w);
if (particleVelm.w != 0)
velm[index] = particleVelm;
}
// Evolve the free ring polymer by transforming to the frequency domain.
__local mixed4* qreal = &q[blockStart];
__local mixed4* qimag = &q[blockStart+get_local_size(0)];
__local mixed4* vreal = &v[blockStart];
__local mixed4* vimag = &v[blockStart+get_local_size(0)];
if (get_local_id(0) < NUM_COPIES)
w[indexInBlock] = (mixed2) (cos(-indexInBlock*2*M_PI/NUM_COPIES), sin(-indexInBlock*2*M_PI/NUM_COPIES));
barrier(CLK_LOCAL_MEM_FENCE);
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
mixed4 particlePosq = posq[particle+indexInBlock*PADDED_NUM_ATOMS];
mixed4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
// Forward FFT.
qreal[indexInBlock] = SCALE*particlePosq;
qimag[indexInBlock] = (mixed4) (0.0f, 0.0f, 0.0f, 0.0f);
vreal[indexInBlock] = SCALE*particleVelm;
vimag[indexInBlock] = (mixed4) (0.0f, 0.0f, 0.0f, 0.0f);
barrier(CLK_GLOBAL_MEM_FENCE);
FFT_Q_FORWARD
FFT_V_FORWARD
// Apply the thermostat.
if (indexInBlock == 0) {
qreal[0].xyz += vreal[0].xyz*dt;
qimag[0].xyz += vimag[0].xyz*dt;
}
else {
const mixed wk = twown*sin(indexInBlock*M_PI/NUM_COPIES);
const mixed wt = wk*dt;
const mixed coswt = cos(wt);
const mixed sinwt = sin(wt);
const mixed4 vprimereal = vreal[indexInBlock]*coswt - qreal[indexInBlock]*(wk*sinwt); // Advance velocity from t to t+dt
const mixed4 vprimeimag = vimag[indexInBlock]*coswt - qimag[indexInBlock]*(wk*sinwt);
qreal[indexInBlock] = vreal[indexInBlock]*(sinwt/wk) + qreal[indexInBlock]*coswt; // Advance position from t to t+dt
qimag[indexInBlock] = vimag[indexInBlock]*(sinwt/wk) + qimag[indexInBlock]*coswt;
vreal[indexInBlock] = vprimereal;
vimag[indexInBlock] = vprimeimag;
}
barrier(CLK_GLOBAL_MEM_FENCE);
// Inverse FFT.
FFT_Q_BACKWARD
FFT_V_BACKWARD
if (particleVelm.w != 0) {
posq[particle+indexInBlock*PADDED_NUM_ATOMS].xyz = SCALE*qreal[indexInBlock].xyz;
velm[particle+indexInBlock*PADDED_NUM_ATOMS].xyz = SCALE*vreal[indexInBlock].xyz;
}
}
}
/**
* Advance the velocities by a half step.
*/
__kernel void advanceVelocities(__global mixed4* velm, __global real4* force, mixed dt) {
const int numBlocks = get_global_size(0)/NUM_COPIES;
const int blockStart = NUM_COPIES*(get_local_id(0)/NUM_COPIES);
const int indexInBlock = get_local_id(0)-blockStart;
// Update velocities.
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
int index = particle+indexInBlock*PADDED_NUM_ATOMS;
mixed4 particleVelm = velm[index];
particleVelm.xyz += convert_mixed4(force[index]).xyz*(0.5f*dt*particleVelm.w);
if (particleVelm.w != 0)
velm[index] = particleVelm;
}
}
/**
* Copy a set of positions and velocities from the integrator's arrays to the context.
*/
__kernel void copyDataToContext(__global mixed4* srcVel, __global mixed4* dstVel, __global mixed4* srcPos,
__global real4* dstPos, __global int* order, int copy) {
const int base = copy*PADDED_NUM_ATOMS;
for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {
int index = base+order[particle];
dstVel[particle] = srcVel[index];
dstPos[particle].xyz = convert_real4(srcPos[index]).xyz;
}
}
/**
* Copy a set of positions, velocities, and forces from the context to the integrator's arrays.
*/
__kernel void copyDataFromContext(__global real4* srcForce, __global real4* dstForce, __global mixed4* srcVel,
__global mixed4* dstVel, __global real4* srcPos, __global mixed4* dstPos, __global int* order, int copy) {
const int base = copy*PADDED_NUM_ATOMS;
for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {
int index = base+order[particle];
dstForce[index] = srcForce[particle];
dstVel[index] = srcVel[particle];
dstPos[index].xyz = convert_mixed4(srcPos[particle]).xyz;
}
}
/**
* Atom positions in one copy have been modified. Apply the same offsets to all the other copies.
*/
__kernel void applyCellTranslations(__global mixed4* posq, __global real4* movedPos, __global int* order, int movedCopy) {
for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {
int index = order[particle];
mixed4 delta = convert_mixed4(movedPos[particle])-posq[movedCopy*PADDED_NUM_ATOMS+index];
for (int copy = 0; copy < NUM_COPIES; copy++)
posq[copy*PADDED_NUM_ATOMS+index].xyz += delta.xyz;
}
}
plugins/rpmd/platforms/opencl/src/kernels/rpmdContraction.cl deleted 100644 → 0
View file @ 98d81730
mixed4 multiplyComplexRealPart(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2r-c1.y*c2i;
}
mixed4 multiplyComplexImagPart(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2i+c1.y*c2r;
}
mixed4 multiplyComplexRealPartConj(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2r+c1.y*c2i;
}
mixed4 multiplyComplexImagPartConj(mixed2 c1, mixed4 c2r, mixed4 c2i) {
return c1.x*c2i-c1.y*c2r;
}
/**
* Compute the contracted positions
*/
__kernel void contractPositions(__global mixed4* posq, __global mixed4* contracted) {
const int numBlocks = get_global_size(0)/NUM_COPIES;
const int blockStart = NUM_COPIES*(get_local_id(0)/NUM_COPIES);
const int indexInBlock = get_local_id(0)-blockStart;
__local mixed4 q[2*THREAD_BLOCK_SIZE];
__local mixed4 temp[2*THREAD_BLOCK_SIZE];
__local mixed2 w1[NUM_COPIES];
__local mixed2 w2[NUM_CONTRACTED_COPIES];
__local mixed4* qreal = &q[blockStart];
__local mixed4* qimag = &q[blockStart+get_local_size(0)];
__local mixed4* tempreal = &temp[blockStart];
__local mixed4* tempimag = &temp[blockStart+get_local_size(0)];
if (get_local_id(0) < NUM_COPIES)
w1[indexInBlock] = (mixed2) (cos(-indexInBlock*2*M_PI/NUM_COPIES), sin(-indexInBlock*2*M_PI/NUM_COPIES));
if (get_local_id(0) < NUM_CONTRACTED_COPIES)
w2[indexInBlock] = (mixed2) (cos(-indexInBlock*2*M_PI/NUM_CONTRACTED_COPIES), sin(-indexInBlock*2*M_PI/NUM_CONTRACTED_COPIES));
barrier(CLK_LOCAL_MEM_FENCE);
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
// Load the particle position.
mixed4 particlePosq = convert_mixed4(posq[particle+indexInBlock*PADDED_NUM_ATOMS]);
qreal[indexInBlock] = particlePosq;
qimag[indexInBlock] = (mixed4) (0.0f, 0.0f, 0.0f, 0.0f);
// Forward FFT.
barrier(CLK_LOCAL_MEM_FENCE);
__local mixed2* w = w1;
FFT_Q_FORWARD
if (NUM_CONTRACTED_COPIES > 1) {
// Compress the data to remove high frequencies.
int start = (NUM_CONTRACTED_COPIES+1)/2;
tempreal[indexInBlock] = qreal[indexInBlock];
tempimag[indexInBlock] = qimag[indexInBlock];
barrier(CLK_LOCAL_MEM_FENCE);
if (indexInBlock < NUM_CONTRACTED_COPIES) {
qreal[indexInBlock] = tempreal[indexInBlock < start ? indexInBlock : indexInBlock+(NUM_COPIES-NUM_CONTRACTED_COPIES)];
qimag[indexInBlock] = tempimag[indexInBlock < start ? indexInBlock : indexInBlock+(NUM_COPIES-NUM_CONTRACTED_COPIES)];
}
barrier(CLK_LOCAL_MEM_FENCE);
w = w2;
FFT_Q_BACKWARD
}
// Store results.
if (indexInBlock < NUM_CONTRACTED_COPIES)
contracted[particle+indexInBlock*PADDED_NUM_ATOMS] = (mixed4) (POS_SCALE*qreal[indexInBlock].x, POS_SCALE*qreal[indexInBlock].y, POS_SCALE*qreal[indexInBlock].z, particlePosq.w);
}
}
/**
* Apply the contracted forces to all copies.
*/
__kernel void contractForces(__global real4* force, __global real4* contracted) {
const int numBlocks = get_global_size(0)/NUM_COPIES;
const int blockStart = NUM_COPIES*(get_local_id(0)/NUM_COPIES);
const int indexInBlock = get_local_id(0)-blockStart;
__local mixed4 f[2*THREAD_BLOCK_SIZE];
__local mixed4 temp[2*THREAD_BLOCK_SIZE];
__local mixed2 w1[NUM_COPIES];
__local mixed2 w2[NUM_CONTRACTED_COPIES];
__local mixed4* freal = &f[blockStart];
__local mixed4* fimag = &f[blockStart+get_local_size(0)];
__local mixed4* tempreal = &temp[blockStart];
__local mixed4* tempimag = &temp[blockStart+get_local_size(0)];
if (get_local_id(0) < NUM_COPIES)
w1[indexInBlock] = (mixed2) (cos(-indexInBlock*2*M_PI/NUM_COPIES), sin(-indexInBlock*2*M_PI/NUM_COPIES));
if (get_local_id(0) < NUM_CONTRACTED_COPIES)
w2[indexInBlock] = (mixed2) (cos(-indexInBlock*2*M_PI/NUM_CONTRACTED_COPIES), sin(-indexInBlock*2*M_PI/NUM_CONTRACTED_COPIES));
barrier(CLK_LOCAL_MEM_FENCE);
for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
// Load the force.
int index = particle+indexInBlock*PADDED_NUM_ATOMS;
if (indexInBlock < NUM_CONTRACTED_COPIES) {
freal[indexInBlock] = convert_mixed4(contracted[index]);
fimag[indexInBlock] = (mixed4) (0.0f, 0.0f, 0.0f, 0.0f);
}
barrier(CLK_LOCAL_MEM_FENCE);
// Forward FFT.
__local mixed2* w = w2;
if (NUM_CONTRACTED_COPIES > 1) {
FFT_F_FORWARD
}
// Set the high frequency components to 0.
int start = (NUM_CONTRACTED_COPIES+1)/2;
int end = NUM_COPIES-NUM_CONTRACTED_COPIES+start;
tempreal[indexInBlock] = freal[indexInBlock];
tempimag[indexInBlock] = fimag[indexInBlock];
barrier(CLK_LOCAL_MEM_FENCE);
if (indexInBlock >= start) {
freal[indexInBlock] = (indexInBlock < end ? (mixed4) (0.0f, 0.0f, 0.0f, 0.0f) : tempreal[indexInBlock-(NUM_COPIES-NUM_CONTRACTED_COPIES)]);
fimag[indexInBlock] = (indexInBlock < end ? (mixed4) (0.0f, 0.0f, 0.0f, 0.0f) : tempimag[indexInBlock-(NUM_COPIES-NUM_CONTRACTED_COPIES)]);
}
barrier(CLK_LOCAL_MEM_FENCE);
w = w1;
FFT_F_BACKWARD
// Store results.
force[index] += convert_real4(FORCE_SCALE*freal[indexInBlock]);
}
}
plugins/rpmd/platforms/opencl/tests/CMakeLists.txt
View file @ be19e022
......@@ -5,6 +5,8 @@
ENABLE_TESTING()
INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR})
INCLUDE_DIRECTORIES(${OPENMM_DIR}/plugins/rpmd/tests)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/opencl/tests)
# Automatically create tests using files named "Test*.cpp"
FILE(GLOB TEST_PROGS "*Test*.cpp")
......
plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.cpp
View file @ be19e022
......@@ -29,546 +29,19 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of RPMDIntegrator.
*/
#include "OpenCLTests.h"
#include "TestRpmd.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/Context.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/RPMDIntegrator.h"
#include "openmm/RPMDMonteCarloBarostat.h"
#include "openmm/VirtualSite.h"
#include "SimTKOpenMMUtilities.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
extern "C" void registerRPMDOpenCLKernelFactories();
using namespace OpenMM;
using namespace std;
extern "C" OPENMM_EXPORT void registerRPMDOpenCLKernelFactories();
void testFreeParticles() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
integ.setPositions(i, positions);
}
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0);
const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Positions | State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
double rg2 = 0.0;
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*mass*v.dot(v);
for (int m = 0; m < numCopies; m++) {
Vec3 delta = state[k].getPositions()[j]-state[m].getPositions()[j];
rg2 += delta.dot(delta);
}
}
rg[j] += rg2/(2*numCopies*numCopies);
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
double meanRg2 = 0.0;
for (int i = 0; i < numParticles; i++)
meanRg2 += rg[i];
meanRg2 /= numSteps*numParticles;
double expectedRg = hbar/(2*sqrt(mass*BOLTZ*temperature));
ASSERT_USUALLY_EQUAL_TOL(expectedRg, sqrt(meanRg2), 1e-3);
}
void testParaHydrogen() {
const int numParticles = 32;
const int numCopies = 12;
const double temperature = 25.0;
const double mass = 2.0;
const double boxSize = 1.1896;
const int numSteps = 1000;
const int numBins = 200;
const double reference[] = {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 4.932814042206152e-5, 1.244331241336431e-4, 4.052316284060125e-4,
1.544810863683946e-3, 4.376197806690222e-3, 1.025847561714293e-2, 2.286702037465422e-2,
4.371052180263602e-2, 7.518538770734748e-2, 0.122351534531647, 0.185758975626622,
0.266399984652322, 0.363380262153250, 0.473696401293219, 0.595312098494172,
0.726049519422861, 0.862264551954547, 0.991102029379444, 1.1147503922535,
1.23587006992066, 1.33495411932817, 1.42208208736987, 1.49273884004107,
1.54633319690403, 1.58714702233941, 1.60439217751355, 1.61804190608902,
1.60680198476058, 1.58892222973695, 1.56387607986781, 1.52629494593350,
1.48421439018970, 1.43656176771959, 1.38752775598872, 1.33310695719931,
1.28363477223121, 1.23465642750248, 1.18874848666326, 1.14350496170519,
1.10292486009936, 1.06107270157688, 1.02348927970441, 0.989729345271297,
0.959273446941802, 0.932264875865758, 0.908818658748942, 0.890946420768315,
0.869332737718165, 0.856401736350349, 0.842370069917020, 0.834386614237393,
0.826268072171045, 0.821547250199453, 0.818786865315836, 0.819441757028076,
0.819156933383128, 0.822275325148621, 0.828919078023881, 0.837233720599450,
0.846961908186718, 0.855656955481099, 0.864520333201247, 0.876082425547566,
0.886950044046000, 0.900275658318995
};
// Create a box of para-hydrogen.
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize,0,0), Vec3(0,boxSize,0), Vec3(0,0,boxSize));
CustomNonbondedForce* nb = new CustomNonbondedForce("2625.49963*(exp(1.713-1.5671*p-0.00993*p*p)-(12.14/p^6+215.2/p^8-143.1/p^9+4813.9/p^10)*(step(rc-p)*exp(-(rc/p-1)^2)+1-step(rc-p))); p=r/0.05291772108; rc=8.32");
nb->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
nb->setCutoffDistance(boxSize/2);
vector<double> params;
for (int i = 0; i < numParticles; i++)
nb->addParticle(params);
system.addForce(nb);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.0005);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
for (int i = 0; i < numCopies; i++)
integ.setPositions(i, positions);
integ.step(1000);
// Simulate it.
vector<int> counts(numBins, 0);
const double invBoxSize = 1.0/boxSize;
double meanKE = 0.0;
const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int step = 0; step < numSteps; step++) {
integ.step(20);
vector<State> states(numCopies);
for (int i = 0; i < numCopies; i++)
states[i] = integ.getState(i, State::Positions | State::Forces);
// Record the radial distribution function.
const vector<Vec3>& pos = states[0].getPositions();
for (int j = 0; j < numParticles; j++)
for (int k = 0; k < j; k++) {
Vec3 delta = pos[j]-pos[k];
delta[0] -= floor(delta[0]*invBoxSize+0.5)*boxSize;
delta[1] -= floor(delta[1]*invBoxSize+0.5)*boxSize;
delta[2] -= floor(delta[2]*invBoxSize+0.5)*boxSize;
double dist = sqrt(delta.dot(delta));
int bin = (int) (numBins*(dist/boxSize));
counts[bin]++;
}
// Calculate the quantum contribution to the kinetic energy.
vector<Vec3> centroids(numParticles, Vec3());
for (int i = 0; i < numCopies; i++) {
const vector<Vec3>& pos = states[i].getPositions();
for (int j = 0; j < numParticles; j++)
centroids[j] += pos[j];
}
for (int j = 0; j < numParticles; j++)
centroids[j] *= 1.0/numCopies;
double ke = 0.0;
for (int i = 0; i < numCopies; i++) {
const vector<Vec3>& pos = states[i].getPositions();
const vector<Vec3>& f = states[i].getForces();
for (int j = 0; j < numParticles; j++) {
Vec3 delta = centroids[j]-pos[j];
ke += delta.dot(f[j]);
}
}
meanKE += ke/(2*numCopies*numParticles);
}
// Check against expected values.
double scale = (boxSize*boxSize*boxSize)/(numSteps*0.5*numParticles*numParticles);
for (int i = 0; i < numBins/2; i++) {
double r1 = i*boxSize/numBins;
double r2 = (i+1)*boxSize/numBins;
double volume = (4.0/3.0)*M_PI*(r2*r2*r2-r1*r1*r1);
ASSERT_USUALLY_EQUAL_TOL(reference[i], scale*counts[i]/volume, 0.1);
}
meanKE /= numSteps*BOLTZ;
ASSERT_USUALLY_EQUAL_TOL(60.0, 1.5*temperature+meanKE, 0.02);
}
Vec3 calcCM(const vector<Vec3>& values, System& system) {
Vec3 cm;
for (int j = 0; j < system.getNumParticles(); ++j) {
cm[0] += values[j][0]*system.getParticleMass(j);
cm[1] += values[j][1]*system.getParticleMass(j);
cm[2] += values[j][2]*system.getParticleMass(j);
}
return cm;
}
void testCMMotionRemoval() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.addForce(new CMMotionRemover());
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
Vec3 cmPos = calcCM(positions, system);
for (int j = 0; j < numParticles; j++)
positions[j] -= cmPos*(1/(mass*numParticles));
integ.setPositions(i, positions);
}
// Make sure the CMMotionRemover is getting applied.
for (int i = 0; i < 200; ++i) {
integ.step(1);
Vec3 pos;
for (int j = 0; j < numCopies; j++) {
State state = integ.getState(0, State::Positions | State::Velocities);
pos += calcCM(state.getPositions(), system);
}
pos *= 1.0/numCopies;
ASSERT_EQUAL_VEC(Vec3(0,0,0), pos, 0.5);
}
}
void testVirtualSites() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*3;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
}
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[3*index] = pos;
positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
positions[3*index+2] = Vec3();
}
integ.setPositions(copy, positions);
}
// Check the temperature and virtual site locations.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++) {
state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
const vector<Vec3>& pos = state[j].getPositions();
for (int k = 0; k < numMolecules; k++)
ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
}
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testContractions() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
map<int, int> contractions;
contractions[1] = 3;
contractions[2] = 1;
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001, contractions);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testWithoutThermostat() {
const int numParticles = 20;
const int numCopies = 10;
const double temperature = 300.0;
const double mass = 2.0;
// Create a chain of particles.
System system;
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
for (int i = 0; i < numParticles; i++) {
system.addParticle(mass);
if (i > 0)
bonds->addBond(i-1, i, 1.0, 1000.0);
}
RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
integ.setApplyThermostat(false);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<vector<Vec3> > positions(numCopies);
for (int i = 0; i < numCopies; i++) {
positions[i].resize(numParticles);
for (int j = 0; j < numParticles; j++)
positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
integ.setPositions(i, positions[i]);
}
// Simulate it and see if the energy remains constant.
double initialEnergy;
int numSteps = 100;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
double energy = integ.getTotalEnergy();
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
}
}
void testWithBarostat() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(100);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
void runPlatformTests() {
testParaHydrogen();
}
int main(int argc, char* argv[]) {
try {
registerRPMDOpenCLKernelFactories();
if (argc > 1)
Platform::getPlatformByName("OpenCL").setPropertyDefaultValue("Precision", string(argv[1]));
testFreeParticles();
testParaHydrogen();
testCMMotionRemoval();
testVirtualSites();
testContractions();
testWithoutThermostat();
testWithBarostat();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
void setupKernels (int argc, char* argv[]) {
registerRPMDOpenCLKernelFactories();
platform = dynamic_cast<OpenCLPlatform&>(Platform::getPlatformByName("OpenCL"));
initializeTests(argc, argv);
}
plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.h
View file @ be19e022
......@@ -73,11 +73,11 @@ public:
*/
double computeKineticEnergy(ContextImpl& context, const RPMDIntegrator& integrator);
/**
* Get the positions of all particles in one copy of the system.
* Set the positions of all particles in one copy of the system.
*/
void setPositions(int copy, const std::vector<Vec3>& positions);
/**
* Get the velocities of all particles in one copy of the system.
* Set the velocities of all particles in one copy of the system.
*/
void setVelocities(int copy, const std::vector<Vec3>& velocities);
/**
......
plugins/rpmd/platforms/reference/tests/CMakeLists.txt
View file @ be19e022
......@@ -2,9 +2,8 @@
# Testing
#
ENABLE_TESTING()
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/include)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/openmmapi/include/openmm)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/tests)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/plugins/rpmd/tests)
SET(SHARED_OPENMM_RPMD_TARGET OpenMMRPMD)
......
plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.cpp
View file @ be19e022
......@@ -29,426 +29,17 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of RPMDIntegrator.
*/
#include "ReferenceTests.h"
#include "TestRpmd.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/RPMDIntegrator.h"
#include "openmm/RPMDMonteCarloBarostat.h"
#include "openmm/VirtualSite.h"
#include "SimTKOpenMMUtilities.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
extern "C" void registerRpmdReferenceKernelFactories();
using namespace OpenMM;
using namespace std;
extern "C" OPENMM_EXPORT void registerRpmdReferenceKernelFactories();
void runPlatformTests() {}
void testFreeParticles() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
integ.setPositions(i, positions);
}
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0);
const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Positions | State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
double rg2 = 0.0;
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*mass*v.dot(v);
for (int m = 0; m < numCopies; m++) {
Vec3 delta = state[k].getPositions()[j]-state[m].getPositions()[j];
rg2 += delta.dot(delta);
}
}
rg[j] += rg2/(2*numCopies*numCopies);
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
double meanRg2 = 0.0;
for (int i = 0; i < numParticles; i++)
meanRg2 += rg[i];
meanRg2 /= numSteps*numParticles;
double expectedRg = hbar/(2*sqrt(mass*BOLTZ*temperature));
ASSERT_USUALLY_EQUAL_TOL(expectedRg, sqrt(meanRg2), 1e-3);
}
Vec3 calcCM(const vector<Vec3>& values, System& system) {
Vec3 cm;
for (int j = 0; j < system.getNumParticles(); ++j) {
cm[0] += values[j][0]*system.getParticleMass(j);
cm[1] += values[j][1]*system.getParticleMass(j);
cm[2] += values[j][2]*system.getParticleMass(j);
}
return cm;
}
void testCMMotionRemoval() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.addForce(new CMMotionRemover());
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
Vec3 cmPos = calcCM(positions, system);
for (int j = 0; j < numParticles; j++)
positions[j] -= cmPos*(1/(mass*numParticles));
integ.setPositions(i, positions);
}
// Make sure the CMMotionRemover is getting applied.
for (int i = 0; i < 200; ++i) {
integ.step(1);
Vec3 pos;
for (int j = 0; j < numCopies; j++) {
State state = integ.getState(0, State::Positions | State::Velocities);
pos += calcCM(state.getPositions(), system);
}
pos *= 1.0/numCopies;
ASSERT_EQUAL_VEC(Vec3(0,0,0), pos, 0.5);
}
}
void testVirtualSites() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*3;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
}
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[3*index] = pos;
positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
positions[3*index+2] = Vec3();
}
integ.setPositions(copy, positions);
}
// Check the temperature and virtual site locations.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++) {
state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
const vector<Vec3>& pos = state[j].getPositions();
for (int k = 0; k < numMolecules; k++)
ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
}
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testContractions() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
map<int, int> contractions;
contractions[1] = 3;
contractions[2] = 1;
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001, contractions);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testWithoutThermostat() {
const int numParticles = 20;
const int numCopies = 10;
const double temperature = 300.0;
const double mass = 2.0;
// Create a chain of particles.
System system;
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
for (int i = 0; i < numParticles; i++) {
system.addParticle(mass);
if (i > 0)
bonds->addBond(i-1, i, 1.0, 1000.0);
}
RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
integ.setApplyThermostat(false);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<vector<Vec3> > positions(numCopies);
for (int i = 0; i < numCopies; i++) {
positions[i].resize(numParticles);
for (int j = 0; j < numParticles; j++)
positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
integ.setPositions(i, positions[i]);
}
// Simulate it and see if the energy remains constant.
double initialEnergy;
int numSteps = 100;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
double energy = integ.getTotalEnergy();
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
}
}
void testWithBarostat() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 5;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
Platform& platform = Platform::getPlatformByName("Reference");
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(100);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
int main() {
try {
registerRpmdReferenceKernelFactories();
testFreeParticles();
testCMMotionRemoval();
testVirtualSites();
testContractions();
testWithoutThermostat();
testWithBarostat();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
void setupKernels (int argc, char* argv[]) {
registerRpmdReferenceKernelFactories();
platform = dynamic_cast<ReferencePlatform&>(Platform::getPlatformByName("Reference"));
initializeTests(argc, argv);
}
plugins/rpmd/tests/TestRpmd.h 0 → 100644
View file @ be19e022
/* -------------------------------------------------------------------------- *
* 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) 2011-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of RPMDIntegrator.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/Context.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/RPMDIntegrator.h"
#include "openmm/RPMDMonteCarloBarostat.h"
#include "openmm/VirtualSite.h"
#include "SimTKOpenMMUtilities.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testFreeParticles() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
integ.setPositions(i, positions);
}
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0);
const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Positions | State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
double rg2 = 0.0;
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*mass*v.dot(v);
for (int m = 0; m < numCopies; m++) {
Vec3 delta = state[k].getPositions()[j]-state[m].getPositions()[j];
rg2 += delta.dot(delta);
}
}
rg[j] += rg2/(2*numCopies*numCopies);
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
double meanRg2 = 0.0;
for (int i = 0; i < numParticles; i++)
meanRg2 += rg[i];
meanRg2 /= numSteps*numParticles;
double expectedRg = hbar/(2*sqrt(mass*BOLTZ*temperature));
ASSERT_USUALLY_EQUAL_TOL(expectedRg, sqrt(meanRg2), 1e-3);
}
void testParaHydrogen() {
const int numParticles = 32;
const int numCopies = 12;
const double temperature = 25.0;
const double mass = 2.0;
const double boxSize = 1.1896;
const int numSteps = 1000;
const int numBins = 200;
const double reference[] = {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 4.932814042206152e-5, 1.244331241336431e-4, 4.052316284060125e-4,
1.544810863683946e-3, 4.376197806690222e-3, 1.025847561714293e-2, 2.286702037465422e-2,
4.371052180263602e-2, 7.518538770734748e-2, 0.122351534531647, 0.185758975626622,
0.266399984652322, 0.363380262153250, 0.473696401293219, 0.595312098494172,
0.726049519422861, 0.862264551954547, 0.991102029379444, 1.1147503922535,
1.23587006992066, 1.33495411932817, 1.42208208736987, 1.49273884004107,
1.54633319690403, 1.58714702233941, 1.60439217751355, 1.61804190608902,
1.60680198476058, 1.58892222973695, 1.56387607986781, 1.52629494593350,
1.48421439018970, 1.43656176771959, 1.38752775598872, 1.33310695719931,
1.28363477223121, 1.23465642750248, 1.18874848666326, 1.14350496170519,
1.10292486009936, 1.06107270157688, 1.02348927970441, 0.989729345271297,
0.959273446941802, 0.932264875865758, 0.908818658748942, 0.890946420768315,
0.869332737718165, 0.856401736350349, 0.842370069917020, 0.834386614237393,
0.826268072171045, 0.821547250199453, 0.818786865315836, 0.819441757028076,
0.819156933383128, 0.822275325148621, 0.828919078023881, 0.837233720599450,
0.846961908186718, 0.855656955481099, 0.864520333201247, 0.876082425547566,
0.886950044046000, 0.900275658318995
};
// Create a box of para-hydrogen.
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize,0,0), Vec3(0,boxSize,0), Vec3(0,0,boxSize));
CustomNonbondedForce* nb = new CustomNonbondedForce("2625.49963*(exp(1.713-1.5671*p-0.00993*p*p)-(12.14/p^6+215.2/p^8-143.1/p^9+4813.9/p^10)*(step(rc-p)*exp(-(rc/p-1)^2)+1-step(rc-p))); p=r/0.05291772108; rc=8.32");
nb->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
nb->setCutoffDistance(boxSize/2);
vector<double> params;
for (int i = 0; i < numParticles; i++)
nb->addParticle(params);
system.addForce(nb);
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.0005);
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
for (int i = 0; i < numCopies; i++)
integ.setPositions(i, positions);
integ.step(1000);
// Simulate it.
vector<int> counts(numBins, 0);
const double invBoxSize = 1.0/boxSize;
double meanKE = 0.0;
for (int step = 0; step < numSteps; step++) {
integ.step(20);
vector<State> states(numCopies);
for (int i = 0; i < numCopies; i++)
states[i] = integ.getState(i, State::Positions | State::Forces);
// Record the radial distribution function.
const vector<Vec3>& pos = states[0].getPositions();
for (int j = 0; j < numParticles; j++)
for (int k = 0; k < j; k++) {
Vec3 delta = pos[j]-pos[k];
delta[0] -= floor(delta[0]*invBoxSize+0.5)*boxSize;
delta[1] -= floor(delta[1]*invBoxSize+0.5)*boxSize;
delta[2] -= floor(delta[2]*invBoxSize+0.5)*boxSize;
double dist = sqrt(delta.dot(delta));
int bin = (int) (numBins*(dist/boxSize));
counts[bin]++;
}
// Calculate the quantum contribution to the kinetic energy.
vector<Vec3> centroids(numParticles, Vec3());
for (int i = 0; i < numCopies; i++) {
const vector<Vec3>& pos = states[i].getPositions();
for (int j = 0; j < numParticles; j++)
centroids[j] += pos[j];
}
for (int j = 0; j < numParticles; j++)
centroids[j] *= 1.0/numCopies;
double ke = 0.0;
for (int i = 0; i < numCopies; i++) {
const vector<Vec3>& pos = states[i].getPositions();
const vector<Vec3>& f = states[i].getForces();
for (int j = 0; j < numParticles; j++) {
Vec3 delta = centroids[j]-pos[j];
ke += delta.dot(f[j]);
}
}
meanKE += ke/(2*numCopies*numParticles);
}
// Check against expected values.
double scale = (boxSize*boxSize*boxSize)/(numSteps*0.5*numParticles*numParticles);
for (int i = 0; i < numBins/2; i++) {
double r1 = i*boxSize/numBins;
double r2 = (i+1)*boxSize/numBins;
double volume = (4.0/3.0)*M_PI*(r2*r2*r2-r1*r1*r1);
ASSERT_USUALLY_EQUAL_TOL(reference[i], scale*counts[i]/volume, 0.1);
}
meanKE /= numSteps*BOLTZ;
ASSERT_USUALLY_EQUAL_TOL(60.0, 1.5*temperature+meanKE, 0.02);
}
Vec3 calcCM(const vector<Vec3>& values, System& system) {
Vec3 cm;
for (int j = 0; j < system.getNumParticles(); ++j) {
cm[0] += values[j][0]*system.getParticleMass(j);
cm[1] += values[j][1]*system.getParticleMass(j);
cm[2] += values[j][2]*system.getParticleMass(j);
}
return cm;
}
void testCMMotionRemoval() {
const int numParticles = 100;
const int numCopies = 30;
const double temperature = 300.0;
const double mass = 1.0;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(mass);
system.addForce(new CMMotionRemover());
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numCopies; i++)
{
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
Vec3 cmPos = calcCM(positions, system);
for (int j = 0; j < numParticles; j++)
positions[j] -= cmPos*(1/(mass*numParticles));
integ.setPositions(i, positions);
}
// Make sure the CMMotionRemover is getting applied.
for (int i = 0; i < 200; ++i) {
integ.step(1);
Vec3 pos;
for (int j = 0; j < numCopies; j++) {
State state = integ.getState(0, State::Positions | State::Velocities);
pos += calcCM(state.getPositions(), system);
}
pos *= 1.0/numCopies;
ASSERT_EQUAL_VEC(Vec3(0,0,0), pos, 0.5);
}
}
void testVirtualSites() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*3;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addParticle(0.1, 0.2, 0.2);
nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
}
RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[3*index] = pos;
positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
positions[3*index+2] = Vec3();
}
integ.setPositions(copy, positions);
}
// Check the temperature and virtual site locations.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++) {
state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
const vector<Vec3>& pos = state[j].getPositions();
for (int k = 0; k < numMolecules; k++)
ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
}
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testContractions() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
map<int, int> contractions;
contractions[1] = 3;
contractions[2] = 1;
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001, contractions);
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(1000);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void testWithoutThermostat() {
const int numParticles = 20;
const int numCopies = 10;
const double temperature = 300.0;
const double mass = 2.0;
// Create a chain of particles.
System system;
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
for (int i = 0; i < numParticles; i++) {
system.addParticle(mass);
if (i > 0)
bonds->addBond(i-1, i, 1.0, 1000.0);
}
RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
integ.setApplyThermostat(false);
Context context(system, integ, platform);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<vector<Vec3> > positions(numCopies);
for (int i = 0; i < numCopies; i++) {
positions[i].resize(numParticles);
for (int j = 0; j < numParticles; j++)
positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
integ.setPositions(i, positions[i]);
}
// Simulate it and see if the energy remains constant.
double initialEnergy;
int numSteps = 100;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
double energy = integ.getTotalEnergy();
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
}
}
void testWithBarostat() {
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const int numParticles = numMolecules*2;
const int numCopies = 10;
const double spacing = 2.0;
const double cutoff = 3.0;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setForceGroup(1);
nonbonded->setReciprocalSpaceForceGroup(2);
system.addForce(nonbonded);
system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.2, 0.2, 0.2);
nonbonded->addParticle(0.2, 0.2, 0.2);
nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
}
RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
Context context(system, integ, platform);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
int index = k+gridSize*(j+gridSize*i);
positions[2*index] = pos;
positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
}
integ.setPositions(copy, positions);
}
// Check the temperature.
const int numSteps = 1000;
integ.step(100);
vector<double> ke(numCopies, 0.0);
for (int i = 0; i < numSteps; i++) {
integ.step(1);
vector<State> state(numCopies);
for (int j = 0; j < numCopies; j++)
state[j] = integ.getState(j, State::Velocities, true);
for (int j = 0; j < numParticles; j++) {
for (int k = 0; k < numCopies; k++) {
Vec3 v = state[k].getVelocities()[j];
ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
}
}
}
double meanKE = 0.0;
for (int i = 0; i < numCopies; i++)
meanKE += ke[i];
meanKE /= numSteps*numCopies;
double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
}
void setupKernels(int argc, char* argv[]);
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
setupKernels(argc, argv);
testFreeParticles();
testCMMotionRemoval();
testVirtualSites();
testContractions();
testWithoutThermostat();
testWithBarostat();
runPlatformTests();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::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