Skip to content

GitLab

Commit 6e842d8c authored by Peter Eastman's avatar Peter Eastman
Browse files

Continuing to implement new CUDA platform: AndersenThermostat, MonteCarloBarostat, CMMotionRemover

parent c962c2dd
Hide whitespace changes
Inline Side-by-side
Showing with 1038 additions and 239 deletions
platforms/cuda2/src/CudaKernelFactory.cpp
View file @ 6e842d8c
......@@ -118,13 +118,13 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
return new CudaIntegrateVariableLangevinStepKernel(name, platform, cu);
// if (name == IntegrateCustomStepKernel::Name())
// return new CudaIntegrateCustomStepKernel(name, platform, cu);
// if (name == ApplyAndersenThermostatKernel::Name())
// return new CudaApplyAndersenThermostatKernel(name, platform, cu);
// if (name == ApplyMonteCarloBarostatKernel::Name())
// return new CudaApplyMonteCarloBarostatKernel(name, platform, cu);
if (name == ApplyAndersenThermostatKernel::Name())
return new CudaApplyAndersenThermostatKernel(name, platform, cu);
if (name == ApplyMonteCarloBarostatKernel::Name())
return new CudaApplyMonteCarloBarostatKernel(name, platform, cu);
if (name == CalcKineticEnergyKernel::Name())
return new CudaCalcKineticEnergyKernel(name, platform, cu);
// if (name == RemoveCMMotionKernel::Name())
// return new CudaRemoveCMMotionKernel(name, platform, cu);
if (name == RemoveCMMotionKernel::Name())
return new CudaRemoveCMMotionKernel(name, platform, cu);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
}
platforms/cuda2/src/CudaKernels.cpp
View file @ 6e842d8c
......@@ -4961,106 +4961,109 @@ double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
// localPerDofValues[3*i+j][variable] = (float) values[order[i]][j];
// deviceValuesAreCurrent = false;
//}
//
//CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
// cu.setAsCurrent();
// if (atomGroups != NULL)
// delete atomGroups;
//}
//
//void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
// cu.setAsCurrent();
// randomSeed = thermostat.getRandomNumberSeed();
// map<string, string> defines;
// defines["NUM_ATOMS"] = cu.intToString(cu.getNumAtoms());
// CUmodule module = cu.createModule(CudaKernelSources::andersenThermostat, defines);
// kernel = cu.getKernel(module, "applyAndersenThermostat");
// cu.getIntegrationUtilities().initRandomNumberGenerator(randomSeed);
//
// // Create the arrays with the group definitions.
//
// vector<vector<int> > groups = AndersenThermostatImpl::calcParticleGroups(system);
// atomGroups = new CudaArray<int>(cu, cu.getNumAtoms(), "atomGroups");
// vector<int> atoms(atomGroups->getSize());
// for (int i = 0; i < (int) groups.size(); i++) {
// for (int j = 0; j < (int) groups[i].size(); j++)
// atoms[groups[i][j]] = i;
// }
// atomGroups->upload(atoms);
//}
//
//void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
// if (!hasInitializedKernels) {
// hasInitializedKernels = true;
// kernel.setArg<cu::Buffer>(2, cu.getVelm().getDevicePointer());
// kernel.setArg<cu::Buffer>(3, cu.getIntegrationUtilities().getStepSize().getDevicePointer());
// kernel.setArg<cu::Buffer>(4, cu.getIntegrationUtilities().getRandom().getDevicePointer());
// kernel.setArg<cu::Buffer>(6, atomGroups->getDevicePointer());
// }
// kernel.setArg<cl_float>(0, (cl_float) context.getParameter(AndersenThermostat::CollisionFrequency()));
// kernel.setArg<cl_float>(1, (cl_float) (BOLTZ*context.getParameter(AndersenThermostat::Temperature())));
// kernel.setArg<cl_uint>(5, cu.getIntegrationUtilities().prepareRandomNumbers(cu.getPaddedNumAtoms()));
// cu.executeKernel(kernel, cu.getNumAtoms());
//}
//
//CudaApplyMonteCarloBarostatKernel::~CudaApplyMonteCarloBarostatKernel() {
// cu.setAsCurrent();
// if (savedPositions != NULL)
// delete savedPositions;
// if (moleculeAtoms != NULL)
// delete moleculeAtoms;
// if (moleculeStartIndex != NULL)
// delete moleculeStartIndex;
//}
//
//void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
// cu.setAsCurrent();
// savedPositions = new CudaArray<mm_float4>(cu, cu.getPaddedNumAtoms(), "savedPositions");
// CUmodule module = cu.createModule(CudaKernelSources::monteCarloBarostat);
// kernel = cu.getKernel(module, "scalePositions");
//}
//
//void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
// if (!hasInitializedKernels) {
// hasInitializedKernels = true;
//
// // Create the arrays with the molecule definitions.
//
// vector<vector<int> > molecules = context.getMolecules();
// numMolecules = molecules.size();
// moleculeAtoms = new CudaArray<int>(cu, cu.getNumAtoms(), "moleculeAtoms");
// moleculeStartIndex = new CudaArray<int>(cu, numMolecules+1, "moleculeStartIndex");
// vector<int> atoms(moleculeAtoms->getSize());
// vector<int> startIndex(moleculeStartIndex->getSize());
// int index = 0;
// for (int i = 0; i < numMolecules; i++) {
// startIndex[i] = index;
// for (int j = 0; j < (int) molecules[i].size(); j++)
// atoms[index++] = molecules[i][j];
// }
// startIndex[numMolecules] = index;
// moleculeAtoms->upload(atoms);
// moleculeStartIndex->upload(startIndex);
//
// // Initialize the kernel arguments.
//
// kernel.setArg<cl_int>(1, numMolecules);
// kernel.setArg<cu::Buffer>(4, cu.getPosq().getDevicePointer());
// kernel.setArg<cu::Buffer>(5, moleculeAtoms->getDevicePointer());
// kernel.setArg<cu::Buffer>(6, moleculeStartIndex->getDevicePointer());
// }
// cu.getQueue().enqueueCopyBuffer(cu.getPosq().getDevicePointer(), savedPositions->getDevicePointer(), 0, 0, cu.getPosq().getSize()*sizeof(mm_float4));
// kernel.setArg<cl_float>(0, (cl_float) scale);
// kernel.setArg<mm_float4>(2, cu.getPeriodicBoxSize());
// kernel.setArg<mm_float4>(3, cu.getInvPeriodicBoxSize());
// cu.executeKernel(kernel, cu.getNumAtoms());
// for (int i = 0; i < (int) cu.getPosCellOffsets().size(); i++)
// cu.getPosCellOffsets()[i] = mm_int4(0, 0, 0, 0);
//}
//
//void CudaApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
// cu.getQueue().enqueueCopyBuffer(savedPositions->getDevicePointer(), cu.getPosq().getDevicePointer(), 0, 0, cu.getPosq().getSize()*sizeof(mm_float4));
//}
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
cu.setAsCurrent();
if (atomGroups != NULL)
delete atomGroups;
}
void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
cu.setAsCurrent();
randomSeed = thermostat.getRandomNumberSeed();
map<string, string> defines;
defines["NUM_ATOMS"] = cu.intToString(cu.getNumAtoms());
CUmodule module = cu.createModule(CudaKernelSources::andersenThermostat, defines);
kernel = cu.getKernel(module, "applyAndersenThermostat");
cu.getIntegrationUtilities().initRandomNumberGenerator(randomSeed);
// Create the arrays with the group definitions.
vector<vector<int> > groups = AndersenThermostatImpl::calcParticleGroups(system);
atomGroups = CudaArray::create<int>(cu, cu.getNumAtoms(), "atomGroups");
vector<int> atoms(atomGroups->getSize());
for (int i = 0; i < (int) groups.size(); i++) {
for (int j = 0; j < (int) groups[i].size(); j++)
atoms[groups[i][j]] = i;
}
atomGroups->upload(atoms);
}
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
float frequency = (float) context.getParameter(AndersenThermostat::CollisionFrequency());
float kT = (float) (BOLTZ*context.getParameter(AndersenThermostat::Temperature()));
int randomIndex = cu.getIntegrationUtilities().prepareRandomNumbers(cu.getPaddedNumAtoms());
void* args[] = {&frequency, &kT, &cu.getVelm().getDevicePointer(), &cu.getIntegrationUtilities().getStepSize().getDevicePointer(),
&cu.getIntegrationUtilities().getRandom().getDevicePointer(), &randomIndex, &atomGroups->getDevicePointer()};
cu.executeKernel(kernel, args, cu.getNumAtoms());
}
CudaApplyMonteCarloBarostatKernel::~CudaApplyMonteCarloBarostatKernel() {
cu.setAsCurrent();
if (savedPositions != NULL)
delete savedPositions;
if (moleculeAtoms != NULL)
delete moleculeAtoms;
if (moleculeStartIndex != NULL)
delete moleculeStartIndex;
}
void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
cu.setAsCurrent();
savedPositions = new CudaArray(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "savedPositions");
CUmodule module = cu.createModule(CudaKernelSources::monteCarloBarostat);
kernel = cu.getKernel(module, "scalePositions");
}
void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
if (!hasInitializedKernels) {
hasInitializedKernels = true;
// Create the arrays with the molecule definitions.
vector<vector<int> > molecules = context.getMolecules();
numMolecules = molecules.size();
moleculeAtoms = CudaArray::create<int>(cu, cu.getNumAtoms(), "moleculeAtoms");
moleculeStartIndex = CudaArray::create<int>(cu, numMolecules+1, "moleculeStartIndex");
vector<int> atoms(moleculeAtoms->getSize());
vector<int> startIndex(moleculeStartIndex->getSize());
int index = 0;
for (int i = 0; i < numMolecules; i++) {
startIndex[i] = index;
for (int j = 0; j < (int) molecules[i].size(); j++)
atoms[index++] = molecules[i][j];
}
startIndex[numMolecules] = index;
moleculeAtoms->upload(atoms);
moleculeStartIndex->upload(startIndex);
// Initialize the kernel arguments.
}
int bytesToCopy = cu.getPosq().getSize()*(cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4));
CUresult result = cuMemcpyDtoD(savedPositions->getDevicePointer(), cu.getPosq().getDevicePointer(), bytesToCopy);
if (result != CUDA_SUCCESS) {
std::stringstream m;
m<<"Error saving positions for MC barostat: "<<cu.getErrorString(result)<<" ("<<result<<")";
throw OpenMMException(m.str());
}
float scalef = (float) scale;
void* args[] = {&scalef, &numMolecules, cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(),
&cu.getPosq().getDevicePointer(), &moleculeAtoms->getDevicePointer(), &moleculeStartIndex->getDevicePointer()};
cu.executeKernel(kernel, args, cu.getNumAtoms());
for (int i = 0; i < (int) cu.getPosCellOffsets().size(); i++)
cu.getPosCellOffsets()[i] = make_int4(0, 0, 0, 0);
}
void CudaApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
int bytesToCopy = cu.getPosq().getSize()*(cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4));
CUresult result = cuMemcpyDtoD(cu.getPosq().getDevicePointer(), savedPositions->getDevicePointer(), bytesToCopy);
if (result != CUDA_SUCCESS) {
std::stringstream m;
m<<"Error restoring positions for MC barostat: "<<cu.getErrorString(result)<<" ("<<result<<")";
throw OpenMMException(m.str());
}
}
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
cu.setAsCurrent();
......@@ -5095,36 +5098,30 @@ double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
return 0.5*energy;
}
//CudaRemoveCMMotionKernel::~CudaRemoveCMMotionKernel() {
// cu.setAsCurrent();
// if (cmMomentum != NULL)
// delete cmMomentum;
//}
//
//void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
// cu.setAsCurrent();
// frequency = force.getFrequency();
// int numAtoms = cu.getNumAtoms();
// cmMomentum = new CudaArray<mm_float4>(cu, (numAtoms+CudaContext::ThreadBlockSize-1)/CudaContext::ThreadBlockSize, "cmMomentum");
// double totalMass = 0.0;
// for (int i = 0; i < numAtoms; i++)
// totalMass += system.getParticleMass(i);
// map<string, string> defines;
// defines["INVERSE_TOTAL_MASS"] = cu.doubleToString(1.0/totalMass);
// CUmodule module = cu.createModule(CudaKernelSources::removeCM, defines);
// kernel1 = cu.getKernel(module, "calcCenterOfMassMomentum");
// kernel1.setArg<cl_int>(0, numAtoms);
// kernel1.setArg<cu::Buffer>(1, cu.getVelm().getDevicePointer());
// kernel1.setArg<cu::Buffer>(2, cmMomentum->getDevicePointer());
// kernel1.setArg(3, CudaContext::ThreadBlockSize*sizeof(mm_float4), NULL);
// kernel2 = cu.getKernel(module, "removeCenterOfMassMomentum");
// kernel2.setArg<cl_int>(0, numAtoms);
// kernel2.setArg<cu::Buffer>(1, cu.getVelm().getDevicePointer());
// kernel2.setArg<cu::Buffer>(2, cmMomentum->getDevicePointer());
// kernel2.setArg(3, CudaContext::ThreadBlockSize*sizeof(mm_float4), NULL);
//}
//
//void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
// cu.executeKernel(kernel1, cu.getNumAtoms());
// cu.executeKernel(kernel2, cu.getNumAtoms());
//}
CudaRemoveCMMotionKernel::~CudaRemoveCMMotionKernel() {
cu.setAsCurrent();
if (cmMomentum != NULL)
delete cmMomentum;
}
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
cu.setAsCurrent();
frequency = force.getFrequency();
int numAtoms = cu.getNumAtoms();
cmMomentum = CudaArray::create<float4>(cu, (numAtoms+CudaContext::ThreadBlockSize-1)/CudaContext::ThreadBlockSize, "cmMomentum");
double totalMass = 0.0;
for (int i = 0; i < numAtoms; i++)
totalMass += system.getParticleMass(i);
map<string, string> defines;
defines["INVERSE_TOTAL_MASS"] = cu.doubleToString(1.0/totalMass);
CUmodule module = cu.createModule(CudaKernelSources::removeCM, defines);
kernel1 = cu.getKernel(module, "calcCenterOfMassMomentum");
kernel2 = cu.getKernel(module, "removeCenterOfMassMomentum");
}
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
int numAtoms = cu.getNumAtoms();
void* args[] = {&numAtoms, &cu.getVelm().getDevicePointer(), &cmMomentum->getDevicePointer()};
cu.executeKernel(kernel1, args, cu.getNumAtoms(), cu.ThreadBlockSize, cu.ThreadBlockSize*sizeof(float4));
cu.executeKernel(kernel2, args, cu.getNumAtoms(), cu.ThreadBlockSize, cu.ThreadBlockSize*sizeof(float4));
}
platforms/cuda2/src/CudaKernels.h
View file @ 6e842d8c
......@@ -1150,79 +1150,78 @@ private:
// std::vector<int> requiredUniform;
// std::vector<std::string> parameterNames;
//};
//
///**
// * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
// */
//class CudaApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
//public:
// CudaApplyAndersenThermostatKernel(std::string name, const Platform& platform, CudaContext& cu) : ApplyAndersenThermostatKernel(name, platform), cu(cu),
// hasInitializedKernels(false), atomGroups(NULL) {
// }
// ~CudaApplyAndersenThermostatKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param thermostat the AndersenThermostat this kernel will be used for
// */
// void initialize(const System& system, const AndersenThermostat& thermostat);
// /**
// * Execute the kernel.
// *
// * @param context the context in which to execute this kernel
// */
// void execute(ContextImpl& context);
//private:
// CudaContext& cu;
// bool hasInitializedKernels;
// int randomSeed;
// CudaArray<cl_int>* atomGroups;
// CUfunction kernel;
//};
//
///**
// * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
// */
//class CudaApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
//public:
// CudaApplyMonteCarloBarostatKernel(std::string name, const Platform& platform, CudaContext& cu) : ApplyMonteCarloBarostatKernel(name, platform), cu(cu),
// hasInitializedKernels(false), savedPositions(NULL), moleculeAtoms(NULL), moleculeStartIndex(NULL) {
// }
// ~CudaApplyMonteCarloBarostatKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param barostat the MonteCarloBarostat this kernel will be used for
// */
// void initialize(const System& system, const MonteCarloBarostat& barostat);
// /**
// * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
// * This is called BEFORE the periodic box size is modified. It should begin by translating each particle
// * or cluster into the first periodic box, so that coordinates will still be correct after the box size
// * is changed.
// *
// * @param context the context in which to execute this kernel
// * @param scale the scale factor by which to multiply particle positions
// */
// void scaleCoordinates(ContextImpl& context, double scale);
// /**
// * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
// * scaleCoordinates() was last called.
// *
// * @param context the context in which to execute this kernel
// */
// void restoreCoordinates(ContextImpl& context);
//private:
// CudaContext& cu;
// bool hasInitializedKernels;
// int numMolecules;
// CudaArray<mm_float4>* savedPositions;
// CudaArray<cl_int>* moleculeAtoms;
// CudaArray<cl_int>* moleculeStartIndex;
// CUfunction kernel;
//};
/**
* This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
*/
class CudaApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
public:
CudaApplyAndersenThermostatKernel(std::string name, const Platform& platform, CudaContext& cu) : ApplyAndersenThermostatKernel(name, platform), cu(cu),
atomGroups(NULL) {
}
~CudaApplyAndersenThermostatKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param thermostat the AndersenThermostat this kernel will be used for
*/
void initialize(const System& system, const AndersenThermostat& thermostat);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
*/
void execute(ContextImpl& context);
private:
CudaContext& cu;
int randomSeed;
CudaArray* atomGroups;
CUfunction kernel;
};
/**
* This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
*/
class CudaApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
public:
CudaApplyMonteCarloBarostatKernel(std::string name, const Platform& platform, CudaContext& cu) : ApplyMonteCarloBarostatKernel(name, platform), cu(cu),
hasInitializedKernels(false), savedPositions(NULL), moleculeAtoms(NULL), moleculeStartIndex(NULL) {
}
~CudaApplyMonteCarloBarostatKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param barostat the MonteCarloBarostat this kernel will be used for
*/
void initialize(const System& system, const MonteCarloBarostat& barostat);
/**
* Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
* This is called BEFORE the periodic box size is modified. It should begin by translating each particle
* or cluster into the first periodic box, so that coordinates will still be correct after the box size
* is changed.
*
* @param context the context in which to execute this kernel
* @param scale the scale factor by which to multiply particle positions
*/
void scaleCoordinates(ContextImpl& context, double scale);
/**
* Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
* scaleCoordinates() was last called.
*
* @param context the context in which to execute this kernel
*/
void restoreCoordinates(ContextImpl& context);
private:
CudaContext& cu;
bool hasInitializedKernels;
int numMolecules;
CudaArray* savedPositions;
CudaArray* moleculeAtoms;
CudaArray* moleculeStartIndex;
CUfunction kernel;
};
/**
* This kernel is invoked to calculate the kinetic energy of the system.
......@@ -1248,33 +1247,33 @@ private:
std::vector<double> masses;
};
///**
// * This kernel is invoked to remove center of mass motion from the system.
// */
//class CudaRemoveCMMotionKernel : public RemoveCMMotionKernel {
//public:
// CudaRemoveCMMotionKernel(std::string name, const Platform& platform, CudaContext& cu) : RemoveCMMotionKernel(name, platform), cu(cu), cmMomentum(NULL) {
// }
// ~CudaRemoveCMMotionKernel();
// /**
// * Initialize the kernel, setting up the particle masses.
// *
// * @param system the System this kernel will be applied to
// * @param force the CMMotionRemover this kernel will be used for
// */
// void initialize(const System& system, const CMMotionRemover& force);
// /**
// * Execute the kernel.
// *
// * @param context the context in which to execute this kernel
// */
// void execute(ContextImpl& context);
//private:
// CudaContext& cu;
// int frequency;
// CudaArray<mm_float4>* cmMomentum;
// CUfunction kernel1, kernel2;
//};
/**
* This kernel is invoked to remove center of mass motion from the system.
*/
class CudaRemoveCMMotionKernel : public RemoveCMMotionKernel {
public:
CudaRemoveCMMotionKernel(std::string name, const Platform& platform, CudaContext& cu) : RemoveCMMotionKernel(name, platform), cu(cu), cmMomentum(NULL) {
}
~CudaRemoveCMMotionKernel();
/**
* Initialize the kernel, setting up the particle masses.
*
* @param system the System this kernel will be applied to
* @param force the CMMotionRemover this kernel will be used for
*/
void initialize(const System& system, const CMMotionRemover& force);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
*/
void execute(ContextImpl& context);
private:
CudaContext& cu;
int frequency;
CudaArray* cmMomentum;
CUfunction kernel1, kernel2;
};
} // namespace OpenMM
......
platforms/cuda2/src/kernels/andersenThermostat.cu 0 → 100644
View file @ 6e842d8c
/**
* Apply the Andersen thermostat to adjust particle velocities.
*/
extern "C" __global__ void applyAndersenThermostat(float collisionFrequency, float kT, real4* velm, const real2* __restrict__ stepSize, const float4* __restrict__ random,
unsigned int randomIndex, const int* __restrict__ atomGroups) {
float collisionProbability = 1.0f-expf(-collisionFrequency*stepSize[0].y);
float randomRange = erff(collisionProbability/sqrtf(2.0f));
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
real4 velocity = velm[index];
float4 selectRand = random[randomIndex+atomGroups[index]];
float4 velRand = random[randomIndex+index];
real scale = (selectRand.w > -randomRange && selectRand.w < randomRange ? 0 : 1);
real add = (1-scale)*SQRT(kT*velocity.w);
velocity.x = scale*velocity.x + add*velRand.x;
velocity.y = scale*velocity.y + add*velRand.y;
velocity.z = scale*velocity.z + add*velRand.z;
velm[index] = velocity;
}
}
platforms/cuda2/src/kernels/monteCarloBarostat.cu 0 → 100644
View file @ 6e842d8c
/**
* Scale the particle positions.
*/
extern "C" __global__ void scalePositions(float scale, int numMolecules, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4* __restrict__ posq,
const int* __restrict__ moleculeAtoms, const int* __restrict__ moleculeStartIndex) {
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numMolecules; index += blockDim.x*gridDim.x) {
int first = moleculeStartIndex[index];
int last = moleculeStartIndex[index+1];
int numAtoms = last-first;
// Find the center of each molecule.
real3 center = make_real3(0, 0, 0);
for (int atom = first; atom < last; atom++) {
real4 pos = posq[moleculeAtoms[atom]];
center.x += pos.x;
center.y += pos.y;
center.z += pos.z;
}
real invNumAtoms = RECIP(numAtoms);
center.x *= invNumAtoms;
center.y *= invNumAtoms;
center.z *= invNumAtoms;
// Move it into the first periodic box.
int xcell = (int) floor(center.x*invPeriodicBoxSize.x);
int ycell = (int) floor(center.y*invPeriodicBoxSize.y);
int zcell = (int) floor(center.z*invPeriodicBoxSize.z);
real3 delta = make_real3(xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z);
center.x -= delta.x;
center.y -= delta.y;
center.z -= delta.z;
// Now scale the position of the molecule center.
delta.x = center.x*(scale-1)-delta.x;
delta.y = center.y*(scale-1)-delta.y;
delta.z = center.z*(scale-1)-delta.z;
for (int atom = first; atom < last; atom++) {
real4 pos = posq[moleculeAtoms[atom]];
pos.x += delta.x;
pos.y += delta.y;
pos.z += delta.z;
posq[moleculeAtoms[atom]] = pos;
}
}
}
platforms/cuda2/src/kernels/removeCM.cu 0 → 100644
View file @ 6e842d8c
/**
* Calculate the center of mass momentum.
*/
extern "C" __global__ void calcCenterOfMassMomentum(int numAtoms, const real4* __restrict__ velm, float4* __restrict__ cmMomentum) {
extern __shared__ volatile float3 temp[];
float3 cm = make_float3(0, 0, 0);
for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
real4 velocity = velm[index];
if (velocity.w != 0.0) {
real mass = RECIP(velocity.w);
cm.x += (float) velocity.x*mass;
cm.y += (float) velocity.y*mass;
cm.z += (float) velocity.z*mass;
}
}
// Sum the threads in this group.
int thread = threadIdx.x;
temp[thread].x = cm.x;
temp[thread].y = cm.y;
temp[thread].z = cm.z;
__syncthreads();
if (thread < 32) {
temp[thread].x += temp[thread+32].x;
temp[thread].y += temp[thread+32].y;
temp[thread].z += temp[thread+32].z;
if (thread < 16) {
temp[thread].x += temp[thread+16].x;
temp[thread].y += temp[thread+16].y;
temp[thread].z += temp[thread+16].z;
}
if (thread < 8) {
temp[thread].x += temp[thread+8].x;
temp[thread].y += temp[thread+8].y;
temp[thread].z += temp[thread+8].z;
}
if (thread < 4) {
temp[thread].x += temp[thread+4].x;
temp[thread].y += temp[thread+4].y;
temp[thread].z += temp[thread+4].z;
}
if (thread < 2) {
temp[thread].x += temp[thread+2].x;
temp[thread].y += temp[thread+2].y;
temp[thread].z += temp[thread+2].z;
}
}
if (thread == 0) {
float3 sum = make_float3(temp[thread].x+temp[thread+1].x, temp[thread].y+temp[thread+1].y, temp[thread].z+temp[thread+1].z);
cmMomentum[blockIdx.x] = make_float4(sum.x, sum.y, sum.z, 0.0f);
}
}
/**
* Remove center of mass motion.
*/
extern "C" __global__ void removeCenterOfMassMomentum(unsigned int numAtoms, real4* __restrict__ velm, const float4* __restrict__ cmMomentum) {
// First sum all of the momenta that were calculated by individual groups.
extern volatile float3 temp[];
float3 cm = make_float3(0, 0, 0);
for (unsigned int index = threadIdx.x; index < gridDim.x; index += blockDim.x) {
float4 momentum = cmMomentum[index];
cm.x += momentum.x;
cm.y += momentum.y;
cm.z += momentum.z;
}
int thread = threadIdx.x;
temp[thread].x = cm.x;
temp[thread].y = cm.y;
temp[thread].z = cm.z;
__syncthreads();
if (thread < 32) {
temp[thread].x += temp[thread+32].x;
temp[thread].y += temp[thread+32].y;
temp[thread].z += temp[thread+32].z;
if (thread < 16) {
temp[thread].x += temp[thread+16].x;
temp[thread].y += temp[thread+16].y;
temp[thread].z += temp[thread+16].z;
}
if (thread < 8) {
temp[thread].x += temp[thread+8].x;
temp[thread].y += temp[thread+8].y;
temp[thread].z += temp[thread+8].z;
}
if (thread < 4) {
temp[thread].x += temp[thread+4].x;
temp[thread].y += temp[thread+4].y;
temp[thread].z += temp[thread+4].z;
}
if (thread < 2) {
temp[thread].x += temp[thread+2].x;
temp[thread].y += temp[thread+2].y;
temp[thread].z += temp[thread+2].z;
}
}
__syncthreads();
cm = make_float3(INVERSE_TOTAL_MASS*(temp[0].x+temp[1].x), INVERSE_TOTAL_MASS*(temp[0].y+temp[1].y), INVERSE_TOTAL_MASS*(temp[0].z+temp[1].z));
// Now remove the center of mass velocity from each atom.
for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
real4 velocity = velm[index];
velocity.x -= cm.x;
velocity.y -= cm.y;
velocity.z -= cm.z;
velm[index] = velocity;
}
}
platforms/cuda2/tests/TestCudaAndersenThermostat.cpp 0 → 100644
View file @ 6e842d8c
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of AndersenThermostat.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/AndersenThermostat.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
const int numSteps = 10000;
CudaPlatform platform;
System system;
VerletIntegrator integrator(0.005);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
system.addForce(thermstat);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < numSteps; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= numSteps;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, ke, 6/std::sqrt((double) numSteps));
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
const int numSteps = 10000;
CudaPlatform platform;
System system;
VerletIntegrator integrator(0.005);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
system.addConstraint(0, 1, 1);
system.addConstraint(1, 2, 1);
system.addConstraint(2, 3, 1);
system.addConstraint(3, 0, 1);
system.addConstraint(4, 5, 1);
system.addConstraint(5, 6, 1);
system.addConstraint(6, 7, 1);
system.addConstraint(7, 4, 1);
AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
system.addForce(thermstat);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(1, 1, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(1, 0, 1);
positions[5] = Vec3(1, 1, 1);
positions[6] = Vec3(0, 1, 1);
positions[7] = Vec3(0, 0, 1);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < numSteps; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= numSteps;
double expected = 0.5*(numParticles*3-system.getNumConstraints())*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, ke, 6/std::sqrt((double) numSteps));
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
CudaPlatform platform;
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
system.addForce(thermostat);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
thermostat->setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
thermostat->setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
int main() {
try {
testTemperature();
testConstraints();
testRandomSeed();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda2/tests/TestCudaCMMotionRemover.cpp 0 → 100644
View file @ 6e842d8c
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of AndersenThermostat.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
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 testMotionRemoval(Integrator& integrator) {
const int numParticles = 8;
CudaPlatform platform;
System system;
HarmonicBondForce* bonds = new HarmonicBondForce();
bonds->addBond(2, 3, 2.0, 0.5);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i+1);
nonbonded->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(nonbonded);
CMMotionRemover* remover = new CMMotionRemover();
system.addForce(remover);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Now run it for a while and see if the center of mass remains fixed.
Vec3 cmPos = calcCM(context.getState(State::Positions).getPositions(), system);
for (int i = 0; i < 1000; ++i) {
integrator.step(1);
State state = context.getState(State::Positions | State::Velocities);
Vec3 pos = calcCM(state.getPositions(), system);
ASSERT_EQUAL_VEC(cmPos, pos, 1e-2);
Vec3 vel = calcCM(state.getVelocities(), system);
if (i > 0) {
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), vel, 1e-2);
}
}
}
int main() {
try {
LangevinIntegrator langevin(0.0, 1e-5, 0.01);
testMotionRemoval(langevin);
VerletIntegrator verlet(0.01);
testMotionRemoval(verlet);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda2/tests/TestCudaMonteCarloBarostat.cpp 0 → 100644
View file @ 6e842d8c
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of MonteCarloBarostat.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/MonteCarloBarostat.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testChangingBoxSize() {
CudaPlatform platform;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
system.addParticle(1.0);
NonbondedForce* nb = new NonbondedForce();
nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nb->setCutoffDistance(2.0);
nb->addParticle(1, 0.5, 0.5);
system.addForce(nb);
LangevinIntegrator integrator(300.0, 1.0, 0.01);
Context context(system, integrator, platform);
Vec3 x, y, z;
context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
// Shrinking the box too small should produce an exception.
context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
bool ok = true;
try {
context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
ok = false;
}
catch (exception& ex) {
}
ASSERT(ok);
}
void testIdealGas() {
const int numParticles = 64;
const int frequency = 10;
const int steps = 1000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25);
const double temp[] = {300.0, 600.0, 1000.0};
const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
CudaPlatform platform;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
}
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
system.addForce(barostat);
// Test it for three different temperatures.
for (int i = 0; i < 3; i++) {
barostat->setTemperature(temp[i]);
LangevinIntegrator integrator(temp[i], 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
integrator.step(frequency);
}
volume /= steps;
double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
}
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
const double pressure = 1.5;
CudaPlatform platform;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(8, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 8));
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, 1);
system.addForce(barostat);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
barostat->setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
barostat->setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
void testWater() {
const int gridSize = 8;
const int numMolecules = gridSize*gridSize*gridSize;
const int frequency = 10;
const int steps = 400;
const double temp = 273.15;
const double pressure = 3;
const double spacing = 0.32;
const double angle = 109.47*M_PI/180;
const double dOH = 0.1;
const double dHH = dOH*2*std::sin(0.5*angle);
// Create a box of SPC water molecules.
CudaPlatform platform;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setUseDispersionCorrection(true);
vector<Vec3> positions;
Vec3 offset1(dOH, 0, 0);
Vec3 offset2(dOH*std::cos(angle), dOH*std::sin(angle), 0);
for (int i = 0; i < gridSize; ++i) {
for (int j = 0; j < gridSize; ++j) {
for (int k = 0; k < gridSize; ++k) {
int firstParticle = system.getNumParticles();
system.addParticle(16.0);
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.82, 0.316557, 0.650194);
nonbonded->addParticle(0.41, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
Vec3 pos = Vec3(spacing*i, spacing*j, spacing*k);
positions.push_back(pos);
positions.push_back(pos+offset1);
positions.push_back(pos+offset2);
system.addConstraint(firstParticle, firstParticle+1, dOH);
system.addConstraint(firstParticle, firstParticle+2, dOH);
system.addConstraint(firstParticle+1, firstParticle+2, dHH);
nonbonded->addException(firstParticle, firstParticle+1, 0, 1, 0);
nonbonded->addException(firstParticle, firstParticle+2, 0, 1, 0);
nonbonded->addException(firstParticle+1, firstParticle+2, 0, 1, 0);
}
}
}
system.addForce(nonbonded);
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
system.addForce(barostat);
// Simulate it and see if the density matches the expected value (1 g/mL).
LangevinIntegrator integrator(temp, 1.0, 0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
integrator.step(2000);
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
double density = numMolecules*18/(AVOGADRO*volume*1e-21);
ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.02);
}
int main() {
try {
testChangingBoxSize();
testIdealGas();
testRandomSeed();
testWater();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment