Skip to content

GitLab

"plugins/rpmd/platforms/cuda/src/CudaRpmdKernelSources.h.in" did not exist on "a402046652cab8ba297aa423e4cb57c904525144"
Unverified Commit f6431a42 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2351 from andysim/nhc 

WIP: Nosé-Hoover Thermostat
parents 2fc77d89 44dca26d
Hide whitespace changes
Inline Side-by-side
Showing with 2107 additions and 6 deletions
platforms/opencl/src/OpenCLIntegrationUtilities.cpp
View file @ f6431a42
......@@ -848,6 +848,24 @@ int OpenCLIntegrationUtilities::prepareRandomNumbers(int numValues) {
}
void OpenCLIntegrationUtilities::createCheckpoint(ostream& stream) {
size_t numChains = noseHooverChainState.size();
bool useDouble = context.getUseDoublePrecision() || context.getUseMixedPrecision();
stream.write((char*) &numChains, sizeof(size_t));
for (auto &chainState: noseHooverChainState){
int chainID = chainState.first;
size_t chainLength = chainState.second.getSize();
stream.write((char*) &chainID, sizeof(int));
stream.write((char*) &chainLength, sizeof(size_t));
if (useDouble) {
vector<mm_double2> stateVec;
chainState.second.download(stateVec);
stream.write((char*) stateVec.data(), sizeof(mm_double2)*chainLength);
} else {
vector<mm_float2> stateVec;
chainState.second.download(stateVec);
stream.write((char*) stateVec.data(), sizeof(mm_float2)*chainLength);
}
}
if (!random.isInitialized())
return;
stream.write((char*) &randomPos, sizeof(int));
......@@ -860,6 +878,28 @@ void OpenCLIntegrationUtilities::createCheckpoint(ostream& stream) {
}
void OpenCLIntegrationUtilities::loadCheckpoint(istream& stream) {
size_t numChains, chainLength;
bool useDouble = context.getUseDoublePrecision() || context.getUseMixedPrecision();
stream.read((char*) &numChains, sizeof(size_t));
noseHooverChainState.clear();
for (size_t i=0; i<numChains; i++){
int chainID;
stream.read((char*) &chainID, sizeof(int));
stream.read((char*) &chainLength, sizeof(size_t));
if (useDouble) {
noseHooverChainState[chainID] = OpenCLArray();
noseHooverChainState[chainID].initialize<mm_double2>(context, chainLength, "chainState" + std::to_string(chainID));
std::vector<mm_double2> stateVec(chainLength);
stream.read((char*) &stateVec[0], sizeof(mm_double2)*chainLength);
noseHooverChainState[chainID].upload(stateVec);
} else {
noseHooverChainState[chainID] = OpenCLArray();
noseHooverChainState[chainID].initialize<mm_float2>(context, chainLength, "chainState" + std::to_string(chainID));
std::vector<mm_float2> stateVec(chainLength);
stream.read((char*) &stateVec[0], sizeof(mm_float2)*chainLength);
noseHooverChainState[chainID].upload(stateVec);
}
}
if (!random.isInitialized())
return;
stream.read((char*) &randomPos, sizeof(int));
......@@ -918,3 +958,7 @@ double OpenCLIntegrationUtilities::computeKineticEnergy(double timeShift) {
posDelta.copyTo(context.getVelm());
return 0.5*energy;
}
std::map<int, OpenCLArray>& OpenCLIntegrationUtilities::getNoseHooverChainState(){
return noseHooverChainState;
};
platforms/opencl/src/OpenCLKernelFactory.cpp
View file @ f6431a42
......@@ -130,6 +130,10 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
return new OpenCLIntegrateCustomStepKernel(name, platform, cl);
if (name == ApplyAndersenThermostatKernel::Name())
return new OpenCLApplyAndersenThermostatKernel(name, platform, cl);
if (name == NoseHooverChainKernel::Name())
return new OpenCLNoseHooverChainKernel(name, platform, cl);
if (name == IntegrateVelocityVerletStepKernel::Name())
return new OpenCLIntegrateVelocityVerletStepKernel(name, platform, cl);
if (name == ApplyMonteCarloBarostatKernel::Name())
return new OpenCLApplyMonteCarloBarostatKernel(name, platform, cl);
if (name == RemoveCMMotionKernel::Name())
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ f6431a42
......@@ -53,6 +53,7 @@
#include "SimTKOpenMMUtilities.h"
#include "jama_eig.h"
#include <algorithm>
#include <assert.h>
#include <cmath>
#include <iterator>
#include <set>
......@@ -413,7 +414,7 @@ void OpenCLUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, co
}
void OpenCLUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
int version = 2;
int version = 3;
stream.write((char*) &version, sizeof(int));
int precision = (cl.getUseDoublePrecision() ? 2 : cl.getUseMixedPrecision() ? 1 : 0);
stream.write((char*) &precision, sizeof(int));
......@@ -444,7 +445,7 @@ void OpenCLUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream
void OpenCLUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& stream) {
int version;
stream.read((char*) &version, sizeof(int));
if (version != 2)
if (version != 3)
throw OpenMMException("Checkpoint was created with a different version of OpenMM");
int precision;
stream.read((char*) &precision, sizeof(int));
......@@ -7391,6 +7392,155 @@ double OpenCLIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& contex
return cl.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
}
void OpenCLIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
cl.getPlatformData().initializeContexts(system);
map<string, string> defines;
defines["BOLTZ"] = cl.doubleToString(BOLTZ);
cl::Program program = cl.createProgram(OpenCLKernelSources::velocityVerlet, defines, "");
kernel1 = cl::Kernel(program, "integrateVelocityVerletPart1");
kernel2 = cl::Kernel(program, "integrateVelocityVerletPart2");
kernel3 = cl::Kernel(program, "integrateVelocityVerletPart3");
kernelHardWall = cl::Kernel(program, "integrateVelocityVerletHardWall");
prevMaxPairDistance = (cl_float) -1.0;
maxPairDistanceBuffer.initialize<cl_float>(cl, 1, "maxPairDistanceBuffer");
}
void OpenCLIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
int paddedNumAtoms = cl.getPaddedNumAtoms();
double dt = integrator.getStepSize();
cl.getIntegrationUtilities().setNextStepSize(dt);
if( !forcesAreValid ) context.calcForcesAndEnergy(true, false);
const auto& atomList = integrator.getAllThermostatedIndividualParticles();
const auto& pairList = integrator.getAllThermostatedPairs();
int numAtoms = atomList.size();
int numPairs = pairList.size();
int numParticles = numAtoms + 2*numPairs;
float maxPairDistance = integrator.getMaximumPairDistance();
// Make sure atom and pair metadata is uploaded and has the correct dimensions
if (prevMaxPairDistance != maxPairDistance) {
std::vector<float> tmp(1, maxPairDistance);
maxPairDistanceBuffer.upload(tmp);
prevMaxPairDistance = maxPairDistance;
}
if (numAtoms !=0 && (!atomListBuffer.isInitialized() || atomListBuffer.getSize() != numAtoms)) {
if (atomListBuffer.isInitialized()) {
atomListBuffer.resize(atomList.size());
} else {
atomListBuffer.initialize<cl_int>(cl, atomList.size(), "atomListBuffer");
}
atomListBuffer.upload(atomList);
}
if (numPairs !=0 && (!pairListBuffer.isInitialized() || pairListBuffer.getSize() != numPairs)) {
if (pairListBuffer.isInitialized()) {
pairListBuffer.resize(pairList.size());
pairTemperatureBuffer.resize(pairList.size());
} else {
pairListBuffer.initialize<mm_int2>(cl, pairList.size(), "pairListBuffer");
pairTemperatureBuffer.initialize<cl_float>(cl, pairList.size(), "pairTemperatureBuffer");
}
std::vector<mm_int2> tmp;
std::vector<float> tmp2;
for(const auto &pair : pairList) {
tmp.push_back(mm_int2(std::get<0>(pair), std::get<1>(pair)));
tmp2.push_back(std::get<2>(pair));
}
pairListBuffer.upload(tmp);
pairTemperatureBuffer.upload(tmp2);
}
//// Call the first integration kernel.
kernel1.setArg<cl_int>(0, numAtoms);
kernel1.setArg<cl_int>(1, numPairs);
kernel1.setArg<cl_int>(2, paddedNumAtoms);
kernel1.setArg<cl::Buffer>(3, cl.getIntegrationUtilities().getStepSize().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
setPosqCorrectionArg(cl, kernel1, 5);
kernel1.setArg<cl::Buffer>(6, cl.getVelm().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(7, cl.getForce().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(8, integration.getPosDelta().getDeviceBuffer());
if (numAtoms > 0) {
kernel1.setArg<cl::Buffer>(9, atomListBuffer.getDeviceBuffer());
} else {
kernel1.setArg<void*>(9, NULL);
}
if (numPairs > 0) {
kernel1.setArg<cl::Buffer>(10, pairListBuffer.getDeviceBuffer());
} else {
kernel1.setArg<void*>(10, NULL);
}
cl.executeKernel(kernel1, std::max(numAtoms, numPairs));
//// Apply constraints.
integration.applyConstraints(integrator.getConstraintTolerance());
//// Call the second integration kernel.
kernel2.setArg<cl_int>(0, numParticles);
kernel2.setArg<cl::Buffer>(1, cl.getIntegrationUtilities().getStepSize().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
setPosqCorrectionArg(cl, kernel2, 3);
kernel2.setArg<cl::Buffer>(4, cl.getVelm().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(5, integration.getPosDelta().getDeviceBuffer());
cl.executeKernel(kernel2, numParticles);
if (numPairs > 0) {
//// Enforce hard wall constraint
kernelHardWall.setArg<cl_int>(0, numPairs);
kernelHardWall.setArg<cl::Buffer>(1, maxPairDistanceBuffer.getDeviceBuffer());
kernelHardWall.setArg<cl::Buffer>(2, cl.getIntegrationUtilities().getStepSize().getDeviceBuffer());
kernelHardWall.setArg<cl::Buffer>(3, cl.getPosq().getDeviceBuffer());
setPosqCorrectionArg(cl, kernelHardWall, 4);
kernelHardWall.setArg<cl::Buffer>(5, cl.getVelm().getDeviceBuffer());
kernelHardWall.setArg<cl::Buffer>(6, pairListBuffer.getDeviceBuffer());
kernelHardWall.setArg<cl::Buffer>(7, pairTemperatureBuffer.getDeviceBuffer());
cl.executeKernel(kernelHardWall, numPairs);
}
integration.computeVirtualSites();
//// Update forces
context.calcForcesAndEnergy(true, false);
forcesAreValid = true;
//// Call the third integration kernel.
kernel3.setArg<cl_int>(0, numAtoms);
kernel3.setArg<cl_int>(1, numPairs);
kernel3.setArg<cl_int>(2, paddedNumAtoms);
kernel3.setArg<cl::Buffer>(3, cl.getIntegrationUtilities().getStepSize().getDeviceBuffer());
kernel3.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
setPosqCorrectionArg(cl, kernel3, 5);
kernel3.setArg<cl::Buffer>(6, cl.getVelm().getDeviceBuffer());
kernel3.setArg<cl::Buffer>(7, cl.getForce().getDeviceBuffer());
kernel3.setArg<cl::Buffer>(8, integration.getPosDelta().getDeviceBuffer());
if (numAtoms > 0) {
kernel3.setArg<cl::Buffer>(9, atomListBuffer.getDeviceBuffer());
} else {
kernel3.setArg<void*>(9, NULL);
}
if (numPairs > 0) {
kernel3.setArg<cl::Buffer>(10, pairListBuffer.getDeviceBuffer());
} else {
kernel3.setArg<void*>(10, NULL);
}
cl.executeKernel(kernel3, std::max(numAtoms, numPairs));
integration.applyVelocityConstraints(integrator.getConstraintTolerance());
//// Update the time and step count.
cl.setTime(cl.getTime()+dt);
cl.setStepCount(cl.getStepCount()+1);
cl.reorderAtoms();
}
double OpenCLIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
return cl.getIntegrationUtilities().computeKineticEnergy(0);
}
void OpenCLIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
cl.getPlatformData().initializeContexts(system);
cl.getIntegrationUtilities().initRandomNumberGenerator(integrator.getRandomNumberSeed());
......@@ -8776,6 +8926,388 @@ void OpenCLApplyAndersenThermostatKernel::execute(ContextImpl& context) {
kernel.setArg<cl_uint>(5, cl.getIntegrationUtilities().prepareRandomNumbers(cl.getPaddedNumAtoms()));
cl.executeKernel(kernel, cl.getNumAtoms());
}
void OpenCLNoseHooverChainKernel::initialize() {
bool useDouble = cl.getUseDoublePrecision() || cl.getUseMixedPrecision();
map<string, string> defines;
defines["BEGIN_YS_LOOP"] = "const real arr[1] = {1.0}; for(int i=0;i<1;++i) { const real ys = arr[i];";
defines["END_YS_LOOP"] = "}";
cl::Program program = cl.createProgram(OpenCLKernelSources::noseHooverChain, defines);
propagateKernels[1] = cl::Kernel(program, "propagateNoseHooverChain");
defines["BEGIN_YS_LOOP"] = "const real arr[3] = {0.828981543588751, -0.657963087177502, 0.828981543588751}; for(int i=0;i<3;++i) { const real ys = arr[i];";
program = cl.createProgram(OpenCLKernelSources::noseHooverChain, defines);
propagateKernels[3] = cl::Kernel(program, "propagateNoseHooverChain");
defines["BEGIN_YS_LOOP"] = "const real arr[5] = {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065}; for(int i=0;i<5;++i) { const real ys = arr[i];";
program = cl.createProgram(OpenCLKernelSources::noseHooverChain, defines);
propagateKernels[5] = cl::Kernel(program, "propagateNoseHooverChain");
program = cl.createProgram(OpenCLKernelSources::noseHooverChain, defines);
reduceEnergyKernel = cl::Kernel(program, "reduceEnergyPair");
computeHeatBathEnergyKernel = cl::Kernel(program, "computeHeatBathEnergy");
computeAtomsKineticEnergyKernel = cl::Kernel(program, "computeAtomsKineticEnergy");
computePairsKineticEnergyKernel = cl::Kernel(program, "computePairsKineticEnergy");
scaleAtomsVelocitiesKernel = cl::Kernel(program, "scaleAtomsVelocities");
scalePairsVelocitiesKernel = cl::Kernel(program, "scalePairsVelocities");
int energyBufferSize = cl.getEnergyBuffer().getSize();
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
energyBuffer.initialize<mm_double2>(cl, energyBufferSize, "energyBuffer");
} else {
energyBuffer.initialize<mm_float2>(cl, energyBufferSize, "energyBuffer");
}
}
std::pair<double, double> OpenCLNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergies, double timeStep) {
bool useDouble = cl.getUseDoublePrecision() || cl.getUseMixedPrecision();
int chainID = nhc.getChainID();
int nAtoms = nhc.getThermostatedAtoms().size();
int nPairs = nhc.getThermostatedPairs().size();
int chainLength = nhc.getChainLength();
int numYS = nhc.getNumYoshidaSuzukiTimeSteps();
int numMTS = nhc.getNumMultiTimeSteps();
int numDOFs = nhc.getNumDegreesOfFreedom();
double temperature = nhc.getTemperature();
double frequency = nhc.getCollisionFrequency();
double relativeTemperature = nhc.getRelativeTemperature();
double relativeFrequency = nhc.getRelativeCollisionFrequency();
if (numYS != 1 && numYS != 3 && numYS != 5) {
throw OpenMMException("Number of Yoshida Suzuki time steps has to be 1, 3, or 5.");
}
auto & chainState = cl.getIntegrationUtilities().getNoseHooverChainState();
if (!scaleFactorBuffer.isInitialized() ||scaleFactorBuffer.getSize() == 0) {
if(useDouble){
std::vector<mm_double2> zeros{{0,0}};
if (scaleFactorBuffer.isInitialized()) {
scaleFactorBuffer.resize(1);
} else {
scaleFactorBuffer.initialize<mm_double2>(cl, 1, "scaleFactorBuffer");
}
scaleFactorBuffer.upload(zeros);
} else {
std::vector<mm_float2> zeros{{0,0}};
if (scaleFactorBuffer.isInitialized()) {
scaleFactorBuffer.resize(1);
} else {
scaleFactorBuffer.initialize<mm_float2>(cl, 1, "scaleFactorBuffer");
}
scaleFactorBuffer.upload(zeros);
}
}
if (!chainForces.isInitialized() || !chainMasses.isInitialized() ){
if(useDouble){
std::vector<cl_double> zeros(chainLength,0);
if (chainForces.isInitialized()) {
chainMasses.resize(chainLength);
chainForces.resize(chainLength);
} else {
chainMasses.initialize<cl_double>(cl, chainLength, "chainMasses");
chainForces.initialize<cl_double>(cl, chainLength, "chainForces");
}
chainMasses.upload(zeros);
chainForces.upload(zeros);
} else {
std::vector<cl_float> zeros(chainLength,0);
if (chainForces.isInitialized()) {
chainMasses.resize(chainLength);
chainForces.resize(chainLength);
} else {
chainMasses.initialize<cl_float>(cl, chainLength, "chainMasses");
chainForces.initialize<cl_float>(cl, chainLength, "chainForces");
}
chainMasses.upload(zeros);
chainForces.upload(zeros);
}
}
if (chainForces.getSize() < chainLength) chainMasses.resize(chainLength);
if (chainMasses.getSize() < chainLength) chainMasses.resize(chainLength);
float timeStepFloat = (float) timeStep;
// N.B. We ignore the incoming kineticEnergy and grab it from the device buffer instead
if (nAtoms) {
if (!chainState.count(2*chainID)) chainState[2*chainID] = OpenCLArray();
if (chainState.at(2*chainID).getSize() != chainLength) {
// We need to upload the OpenCL array
if(useDouble){
if (chainState.at(2*chainID).isInitialized()) {
chainState.at(2*chainID).resize(chainLength);
} else {
chainState.at(2*chainID).initialize<mm_double2>(cl, chainLength, "chainState" + std::to_string(2*chainID));
}
std::vector<mm_double2> zeros(chainLength, mm_double2(0.0, 0.0));
chainState.at(2*chainID).upload(zeros.data());
} else {
if (chainState.at(2*chainID).isInitialized()) {
chainState.at(2*chainID).resize(chainLength);
} else {
chainState.at(2*chainID).initialize<mm_float2>(cl, chainLength, "chainState" + std::to_string(2*chainID));
}
std::vector<mm_float2> zeros(chainLength, mm_float2(0.0f, 0.0f));
chainState.at(2*chainID).upload(zeros.data());
}
}
int chainType = 0;
double kT = BOLTZ * temperature;
float kTfloat = (float) kT;
float frequencyFloat = (float) frequency;
propagateKernels[numYS].setArg<cl::Buffer>(0, chainState[2*chainID].getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(1, kineticEnergyBuffer.getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(2, scaleFactorBuffer.getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(3, chainMasses.getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(4, chainForces.getDeviceBuffer());
propagateKernels[numYS].setArg<cl_int>(5, chainType);
propagateKernels[numYS].setArg<cl_int>(6, chainLength);
propagateKernels[numYS].setArg<cl_int>(7, numMTS);
propagateKernels[numYS].setArg<cl_int>(8, numDOFs);
propagateKernels[numYS].setArg<cl_float>(9, timeStepFloat);
if (useDouble)
propagateKernels[numYS].setArg<cl_double>(10, kT);
else
propagateKernels[numYS].setArg<cl_float>(10, kTfloat);
propagateKernels[numYS].setArg<cl_float>(11, frequencyFloat);
cl.executeKernel(propagateKernels[numYS], 1, 1);
}
if (nPairs) {
if (!chainState.count(2*chainID+1)) chainState[2*chainID+1] = OpenCLArray();
if (chainState.at(2*chainID+1).getSize() != chainLength) {
// We need to upload the OpenCL array
if(useDouble){
if (chainState.at(2*chainID+1).isInitialized()) {
chainState.at(2*chainID+1).resize(chainLength);
} else {
chainState.at(2*chainID+1).initialize<mm_double2>(cl, chainLength, "chainState" + std::to_string(2*chainID+1));
}
std::vector<mm_double2> zeros(chainLength, mm_double2(0.0, 0.0));
chainState.at(2*chainID+1).upload(zeros.data());
} else {
if (chainState.at(2*chainID+1).isInitialized()) {
chainState.at(2*chainID+1).resize(chainLength);
} else {
chainState.at(2*chainID+1).initialize<mm_float2>(cl, chainLength, "chainState" + std::to_string(2*chainID+1));
}
std::vector<mm_float2> zeros(chainLength, mm_float2(0.0f, 0.0f));
chainState.at(2*chainID+1).upload(zeros.data());
}
}
int chainType = 1;
double kT = BOLTZ * relativeTemperature;
int ndf = 3*nPairs;
float kTfloat = (float) kT;
float frequencyFloat = (float) relativeFrequency;
propagateKernels[numYS].setArg<cl::Buffer>(0, chainState[2*chainID+1].getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(1, kineticEnergyBuffer.getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(2, scaleFactorBuffer.getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(3, chainMasses.getDeviceBuffer());
propagateKernels[numYS].setArg<cl::Buffer>(4, chainForces.getDeviceBuffer());
propagateKernels[numYS].setArg<cl_int>(5, chainType);
propagateKernels[numYS].setArg<cl_int>(6, chainLength);
propagateKernels[numYS].setArg<cl_int>(7, numMTS);
propagateKernels[numYS].setArg<cl_int>(8, ndf);
propagateKernels[numYS].setArg<cl_float>(9, timeStepFloat);
if (useDouble)
propagateKernels[numYS].setArg<cl_double>(10, kT);
else
propagateKernels[numYS].setArg<cl_float>(10, kTfloat);
propagateKernels[numYS].setArg<cl_float>(11, frequencyFloat);
cl.executeKernel(propagateKernels[numYS], 1, 1);
}
return {0, 0};
}
double OpenCLNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
bool useDouble = cl.getUseDoublePrecision() || cl.getUseMixedPrecision();
int chainID = nhc.getChainID();
int chainLength = nhc.getChainLength();
auto & chainState = cl.getIntegrationUtilities().getNoseHooverChainState();
bool absChainIsValid = chainState.count(2*chainID) != 0 &&
chainState[2*chainID].isInitialized() &&
chainState[2*chainID].getSize() == chainLength;
bool relChainIsValid = chainState.count(2*chainID+1) != 0 &&
chainState[2*chainID+1].isInitialized() &&
chainState[2*chainID+1].getSize() == chainLength;
if (!absChainIsValid && !relChainIsValid) return 0.0;
if (!heatBathEnergy.isInitialized() || heatBathEnergy.getSize() == 0) {
if(useDouble){
std::vector<cl_double> one(1);
heatBathEnergy.initialize<cl_double>(cl, 1, "heatBathEnergy");
heatBathEnergy.upload(one);
} else {
std::vector<cl_float> one(1);
heatBathEnergy.initialize<cl_float>(cl, 1, "heatBathEnergy");
heatBathEnergy.upload(one);
}
}
cl.clearBuffer(heatBathEnergy);
if (absChainIsValid) {
int numDOFs = nhc.getNumDegreesOfFreedom();
double temperature = nhc.getTemperature();
double frequency = nhc.getCollisionFrequency();
double kT = BOLTZ * temperature;
float kTfloat = (float) kT;
float frequencyFloat = (float) frequency;
computeHeatBathEnergyKernel.setArg<cl::Buffer>(0, heatBathEnergy.getDeviceBuffer());
computeHeatBathEnergyKernel.setArg<cl_int>(1, chainLength);
computeHeatBathEnergyKernel.setArg<cl_int>(2, numDOFs);
if (useDouble)
computeHeatBathEnergyKernel.setArg<cl_double>(3, kT);
else
computeHeatBathEnergyKernel.setArg<cl_float>(3, kTfloat);
computeHeatBathEnergyKernel.setArg<cl_float>(4, frequencyFloat);
computeHeatBathEnergyKernel.setArg<cl::Buffer>(5, chainState[2*chainID].getDeviceBuffer());
cl.executeKernel(computeHeatBathEnergyKernel, 1, 1);
}
if (relChainIsValid) {
int numDOFs = 3 * nhc.getThermostatedPairs().size();
double temperature = nhc.getRelativeTemperature();
double frequency = nhc.getRelativeCollisionFrequency();
double kT = BOLTZ * temperature;
float kTfloat = (float) kT;
float frequencyFloat = (float) frequency;
computeHeatBathEnergyKernel.setArg<cl::Buffer>(0, heatBathEnergy.getDeviceBuffer());
computeHeatBathEnergyKernel.setArg<cl_int>(1, chainLength);
computeHeatBathEnergyKernel.setArg<cl_int>(2, numDOFs);
if (useDouble)
computeHeatBathEnergyKernel.setArg<cl_double>(3, kT);
else
computeHeatBathEnergyKernel.setArg<cl_float>(3, kTfloat);
computeHeatBathEnergyKernel.setArg<cl_float>(4, frequencyFloat);
computeHeatBathEnergyKernel.setArg<cl::Buffer>(5, chainState[2*chainID+1].getDeviceBuffer());
cl.executeKernel(computeHeatBathEnergyKernel, 1, 1);
}
void * pinnedBuffer = cl.getPinnedBuffer();
heatBathEnergy.download(pinnedBuffer);
if (useDouble){
return *((double*) pinnedBuffer);
} else {
return *((float*) pinnedBuffer);
}
}
std::pair<double, double> OpenCLNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &nhc, bool downloadValue) {
bool useDouble = cl.getUseDoublePrecision() || cl.getUseMixedPrecision();
int chainID = nhc.getChainID();
const auto & nhcAtoms = nhc.getThermostatedAtoms();
const auto & nhcPairs = nhc.getThermostatedPairs();
auto nAtoms = nhcAtoms.size();
auto nPairs = nhcPairs.size();
if (nAtoms) {
if (!atomlists.count(chainID)) {
// We need to upload the OpenCL array
atomlists[chainID] = OpenCLArray();
atomlists[chainID].initialize<int>(cl, nAtoms, "atomlist" + std::to_string(chainID));
atomlists[chainID].upload(nhcAtoms);
}
if (atomlists[chainID].getSize() != nAtoms) {
throw OpenMMException("Number of atoms changed. Cannot be handled by the same Nose-Hoover thermostat.");
}
}
if (nPairs) {
if (!pairlists.count(chainID)) {
// We need to upload the OpenCL array
pairlists[chainID] = OpenCLArray();
pairlists[chainID].initialize<mm_int2>(cl, nPairs, "pairlist" + std::to_string(chainID));
std::vector<mm_int2> int2vec;
for(const auto &p : nhcPairs) int2vec.push_back(mm_int2(p.first, p.second));
pairlists[chainID].upload(int2vec);
}
if (pairlists[chainID].getSize() != nPairs) {
throw OpenMMException("Number of thermostated pairs changed. Cannot be handled by the same Nose-Hoover thermostat.");
}
}
if (!kineticEnergyBuffer.isInitialized() || kineticEnergyBuffer.getSize() == 0) {
if(useDouble){
std::vector<mm_double2> zeros{{0,0}};
kineticEnergyBuffer.initialize<mm_double2>(cl, 1, "kineticEnergyBuffer");
kineticEnergyBuffer.upload(zeros);
} else {
std::vector<mm_float2> zeros{{0,0}};
kineticEnergyBuffer.initialize<mm_float2>(cl, 1, "kineticEnergyBuffer");
kineticEnergyBuffer.upload(zeros);
}
}
cl.clearBuffer(cl.getEnergyBuffer());
if (nAtoms) {
computeAtomsKineticEnergyKernel.setArg<cl::Buffer>(0, energyBuffer.getDeviceBuffer());
computeAtomsKineticEnergyKernel.setArg<cl_int>(1, nAtoms);
computeAtomsKineticEnergyKernel.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
computeAtomsKineticEnergyKernel.setArg<cl::Buffer>(3, atomlists[chainID].getDeviceBuffer());
cl.executeKernel(computeAtomsKineticEnergyKernel, nAtoms);
}
if (nPairs) {
computePairsKineticEnergyKernel.setArg<cl::Buffer>(0, energyBuffer.getDeviceBuffer());
computePairsKineticEnergyKernel.setArg<cl_int>(1, nPairs);
computePairsKineticEnergyKernel.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
computePairsKineticEnergyKernel.setArg<cl::Buffer>(3, pairlists[chainID].getDeviceBuffer());
cl.executeKernel(computePairsKineticEnergyKernel, nPairs);
}
int bufferSize = energyBuffer.getSize();
int workGroupSize = cl.getDevice().getInfo<CL_DEVICE_MAX_WORK_GROUP_SIZE>();
if (workGroupSize > 512)
workGroupSize = 512;
reduceEnergyKernel.setArg<cl::Buffer>(0, energyBuffer.getDeviceBuffer());
reduceEnergyKernel.setArg<cl::Buffer>(1, kineticEnergyBuffer.getDeviceBuffer());
reduceEnergyKernel.setArg<cl_int>(2, bufferSize);
reduceEnergyKernel.setArg<cl_int>(3, workGroupSize);
reduceEnergyKernel.setArg(4, workGroupSize*energyBuffer.getElementSize(), NULL);
cl.executeKernel(reduceEnergyKernel, workGroupSize, workGroupSize);
std::pair<double, double> KEs = {0, 0};
if (downloadValue) {
if (useDouble) {
mm_double2 tmp;
kineticEnergyBuffer.download(&tmp);
KEs.first = tmp.x;
KEs.second = tmp.y;
} else {
mm_float2 tmp;
kineticEnergyBuffer.download(&tmp);
KEs.first = tmp.x;
KEs.second = tmp.y;
}
}
return KEs;
}
void OpenCLNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> scaleFactor) {
// For now we assume that the atoms and pairs info is valid, because compute{Atoms|Pairs}KineticEnergy must have been
// called before this kernel. If that ever ceases to be true, some sanity checks are needed here.
int chainID = nhc.getChainID();
auto nAtoms = nhc.getThermostatedAtoms().size();
auto nPairs = nhc.getThermostatedPairs().size();
if(nAtoms) {
scaleAtomsVelocitiesKernel.setArg<cl::Buffer>(0, scaleFactorBuffer.getDeviceBuffer());
scaleAtomsVelocitiesKernel.setArg<cl_int>(1, nAtoms);
scaleAtomsVelocitiesKernel.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
scaleAtomsVelocitiesKernel.setArg<cl::Buffer>(3, atomlists[chainID].getDeviceBuffer());
cl.executeKernel(scaleAtomsVelocitiesKernel, nAtoms);
}
if(nPairs) {
scalePairsVelocitiesKernel.setArg<cl::Buffer>(0, scaleFactorBuffer.getDeviceBuffer());
scalePairsVelocitiesKernel.setArg<cl_int>(1, nPairs);
scalePairsVelocitiesKernel.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
scalePairsVelocitiesKernel.setArg<cl::Buffer>(3, pairlists[chainID].getDeviceBuffer());
cl.executeKernel(scalePairsVelocitiesKernel, nPairs);
}
}
void OpenCLApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& thermostat) {
savedPositions.initialize(cl, cl.getPaddedNumAtoms(), cl.getUseDoublePrecision() ? sizeof(mm_double4) : sizeof(mm_float4), "savedPositions");
......
platforms/opencl/src/OpenCLPlatform.cpp
View file @ f6431a42
......@@ -87,6 +87,7 @@ OpenCLPlatform::OpenCLPlatform() {
registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateVelocityVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateBAOABStepKernel::Name(), factory);
registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
......@@ -94,6 +95,7 @@ OpenCLPlatform::OpenCLPlatform() {
registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
registerKernelFactory(NoseHooverChainKernel::Name(), factory);
registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
platformProperties.push_back(OpenCLDeviceIndex());
......
platforms/opencl/src/kernels/noseHooverChain.cl 0 → 100644
View file @ f6431a42
//#include <initializer_list>
__kernel void propagateNoseHooverChain(__global mixed2* restrict chainData, __global const mixed2 * restrict energySum, __global mixed2* restrict scaleFactor,
__global mixed* restrict chainMasses, __global mixed* restrict chainForces,
int chainType, int chainLength, int numMTS, int numDOFs, float timeStep,
mixed kT, float frequency){
const mixed kineticEnergy = chainType == 0 ? energySum[0].x : energySum[0].y;
mixed scale = 1;
if(kineticEnergy < 1e-8) return;
for (int bead = 0; bead < chainLength; ++bead) chainMasses[bead] = kT / (frequency * frequency);
chainMasses[0] *= numDOFs;
mixed KE2 = 2.0f * kineticEnergy;
mixed timeOverMTS = timeStep / numMTS;
chainForces[0] = (KE2 - numDOFs * kT) / chainMasses[0];
for (int bead = 0; bead < chainLength - 1; ++bead) {
chainForces[bead + 1] = (chainMasses[bead] * chainData[bead].y * chainData[bead].y - kT) / chainMasses[bead + 1];
}
for (int mts = 0; mts < numMTS; ++mts) {
BEGIN_YS_LOOP
mixed wdt = ys * timeOverMTS;
chainData[chainLength-1].y += 0.25f * wdt * chainForces[chainLength-1];
for (int bead = chainLength - 2; bead >= 0; --bead) {
mixed aa = EXP(-0.125f * wdt * chainData[bead + 1].y);
chainData[bead].y = aa * (chainData[bead].y * aa + 0.25f * wdt * chainForces[bead]);
}
// update particle velocities
mixed aa = EXP(-0.5f * wdt * chainData[0].y);
scale *= aa;
// update the thermostat positions
for (int bead = 0; bead < chainLength; ++bead) {
chainData[bead].x += 0.5f * chainData[bead].y * wdt;
}
// update the forces
chainForces[0] = (scale * scale * KE2 - numDOFs * kT) / chainMasses[0];
// update thermostat velocities
for (int bead = 0; bead < chainLength - 1; ++bead) {
mixed aa = EXP(-0.125f * wdt * chainData[bead + 1].y);
chainData[bead].y = aa * (aa * chainData[bead].y + 0.25f * wdt * chainForces[bead]);
chainForces[bead + 1] = (chainMasses[bead] * chainData[bead].y * chainData[bead].y - kT) / chainMasses[bead + 1];
}
chainData[chainLength-1].y += 0.25f * wdt * chainForces[chainLength-1];
END_YS_LOOP
} // MTS loop
if (chainType == 0) {
scaleFactor[0].x = scale;
} else {
scaleFactor[0].y = scale;
}
}
/**
* Compute total (potential + kinetic) energy of the Nose-Hoover beads
*/
__kernel void computeHeatBathEnergy(__global mixed* restrict heatBathEnergy, int chainLength, int numDOFs,
mixed kT, float frequency, __global const mixed2* restrict chainData){
// Note that this is always incremented; make sure it's zeroed properly before the first call
for(int i = 0; i < chainLength; ++i) {
mixed prefac = i ? 1 : numDOFs;
mixed mass = prefac * kT / (frequency * frequency);
mixed velocity = chainData[i].y;
// The kinetic energy of this bead
heatBathEnergy[0] += 0.5f * mass * velocity * velocity;
// The potential energy of this bead
mixed position = chainData[i].x;
heatBathEnergy[0] += prefac * kT * position;
}
}
__kernel void computeAtomsKineticEnergy(__global mixed2 * restrict energyBuffer, int numAtoms,
__global const mixed4* restrict velm, __global const int *restrict atoms){
mixed2 energy = (mixed2) (0,0);
//energy = 1; return;
int index = get_global_id(0);
while (index < numAtoms){
int atom = atoms[index];
mixed4 v = velm[atom];
mixed mass = v.w == 0 ? 0 : 1 / v.w;
energy.x += 0.5f * mass * (v.x*v.x + v.y*v.y + v.z*v.z);
index += get_global_size(0);
}
energyBuffer[get_global_id(0)] = energy;
}
__kernel void computePairsKineticEnergy(__global mixed2 * restrict energyBuffer, int numPairs,
__global const mixed4* restrict velm, __global const int2 *restrict pairs){
mixed2 energy = (mixed2) (0,0);
int index = get_global_id(0);
while (index < numPairs){
int2 pair = pairs[index];
int atom1 = pair.x;
int atom2 = pair.y;
mixed4 v1 = velm[atom1];
mixed4 v2 = velm[atom2];
mixed m1 = v1.w == 0 ? 0 : 1 / v1.w;
mixed m2 = v2.w == 0 ? 0 : 1 / v2.w;
mixed4 cv;
cv.x = (m1*v1.x + m2*v2.x) / (m1 + m2);
cv.y = (m1*v1.y + m2*v2.y) / (m1 + m2);
cv.z = (m1*v1.z + m2*v2.z) / (m1 + m2);
mixed4 rv;
rv.x = v2.x - v1.x;
rv.y = v2.y - v1.y;
rv.z = v2.z - v1.z;
energy.x += 0.5f * (m1 + m2) * (cv.x*cv.x + cv.y*cv.y + cv.z*cv.z);
energy.y += 0.5f * (m1 * m2 / (m1 + m2)) * (rv.x*rv.x + rv.y*rv.y + rv.z*rv.z);
index += get_global_size(0);
}
// The atoms version of this has been called already, so accumulate instead of assigning here
energyBuffer[get_global_id(0)].xy += energy.xy;
}
__kernel void scaleAtomsVelocities(__global mixed2* restrict scaleFactor, int numAtoms,
__global mixed4* restrict velm, __global const int *restrict atoms){
const mixed scale = scaleFactor[0].x;
int index = get_global_id(0);
while (index < numAtoms){
int atom = atoms[index];
velm[atom].x *= scale;
velm[atom].y *= scale;
velm[atom].z *= scale;
index += get_global_size(0);
}
}
__kernel void scalePairsVelocities(__global mixed2 * restrict scaleFactor, int numPairs,
__global mixed4* restrict velm, __global const int2 *restrict pairs){
int index = get_global_id(0);
while (index < numPairs){
int atom1 = pairs[index].x;
int atom2 = pairs[index].y;
mixed m1 = velm[atom1].w == 0 ? 0 : 1 / velm[atom1].w;
mixed m2 = velm[atom2].w == 0 ? 0 : 1 / velm[atom2].w;
mixed4 cv;
cv.xyz = (m1*velm[atom1].xyz + m2*velm[atom2].xyz) / (m1 + m2);
mixed4 rv;
rv.xyz = velm[atom2].xyz - velm[atom1].xyz;
velm[atom1].x = scaleFactor[0].x * cv.x - scaleFactor[0].y * rv.x * m2 / (m1 + m2);
velm[atom1].y = scaleFactor[0].x * cv.y - scaleFactor[0].y * rv.y * m2 / (m1 + m2);
velm[atom1].z = scaleFactor[0].x * cv.z - scaleFactor[0].y * rv.z * m2 / (m1 + m2);
velm[atom2].x = scaleFactor[0].x * cv.x + scaleFactor[0].y * rv.x * m1 / (m1 + m2);
velm[atom2].y = scaleFactor[0].x * cv.y + scaleFactor[0].y * rv.y * m1 / (m1 + m2);
velm[atom2].z = scaleFactor[0].x * cv.z + scaleFactor[0].y * rv.z * m1 / (m1 + m2);
index += get_global_size(0);
}
}
/**
* Sum the energy buffer containing a pair of energies stored as mixed2. This is copied from utilities.cu with small modifications
*/
__kernel void reduceEnergyPair(__global const mixed2* restrict energyBuffer, __global mixed2* restrict result, int bufferSize, int workGroupSize, __local mixed2* restrict tempBuffer) {
const unsigned int thread = get_local_id(0);
mixed2 sum = (mixed2) (0,0);
for (unsigned int index = thread; index < bufferSize; index += get_local_size(0)) {
sum.xy += energyBuffer[index].xy;
}
tempBuffer[thread].xy = sum.xy;
for (int i = 1; i < workGroupSize; i *= 2) {
barrier(CLK_LOCAL_MEM_FENCE);
if (thread%(i*2) == 0 && thread+i < workGroupSize) {
tempBuffer[thread].xy += tempBuffer[thread+i].xy;
}
}
if (thread == 0) {
*result = tempBuffer[0];
}
}
platforms/opencl/src/kernels/velocityVerlet.cl 0 → 100644
View file @ f6431a42
/**
* Perform the first step of Velocity Verlet integration.
*/
__kernel void integrateVelocityVerletPart1(int numAtoms, int numPairs, int paddedNumAtoms, __global const mixed2* restrict dt, __global const real4* restrict posq,
__global const real4* restrict posqCorrection, __global mixed4* restrict velm, __global const real4* restrict force, __global mixed4* restrict posDelta,
__global const int* restrict atomList, __global const int2* restrict pairList) {
const mixed2 stepSize = dt[0];
const mixed dtPos = stepSize.y;
const mixed dtVel = 0.5f*(stepSize.x+stepSize.y);
int index = get_global_id(0);
while (index < numAtoms) {
int atom = atomList[index];
mixed4 velocity = velm[atom];
if (velocity.w != 0.0) {
#ifdef USE_MIXED_PRECISION
real4 pos1 = posq[atom];
real4 pos2 = posqCorrection[atom];
mixed4 pos = (mixed4) (pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
#else
real4 pos = posq[atom];
#endif
velocity.x += 0.5f * dtVel*force[atom].x*velocity.w;
velocity.y += 0.5f * dtVel*force[atom].y*velocity.w;
velocity.z += 0.5f * dtVel*force[atom].z*velocity.w;
pos.x = velocity.x*dtPos;
pos.y = velocity.y*dtPos;
pos.z = velocity.z*dtPos;
posDelta[atom] = pos;
velm[atom] = velocity;
}
index += get_global_size(0);
}
index = get_global_id(0);
while (index < numPairs){
int atom1 = pairList[index].x;
int atom2 = pairList[index].y;
mixed4 v1 = velm[atom1];
mixed4 v2 = velm[atom2];
mixed m1 = v1.w == 0.0f ? 0.0f : 1.0f / v1.w;
mixed m2 = v2.w == 0.0f ? 0.0f : 1.0f / v2.w;
mixed mass1fract = m1 / (m1 + m2);
mixed mass2fract = m2 / (m1 + m2);
mixed invRedMass = (m1 * m2 != 0.0f) ? (m1 + m2)/(m1 * m2) : 0.0f;
mixed invTotMass = (m1 + m2 != 0.0f) ? 1.0f /(m1 + m2) : 0.0f;
mixed3 comVel;
comVel.x= v1.x*mass1fract + v2.x*mass2fract;
comVel.y= v1.y*mass1fract + v2.y*mass2fract;
comVel.z= v1.z*mass1fract + v2.z*mass2fract;
mixed3 relVel;
relVel.x= v2.x - v1.x;
relVel.y= v2.y - v1.y;
relVel.z= v2.z - v1.z;
mixed3 comFrc;
comFrc.x = force[atom1].x + force[atom2].x;
comFrc.y = force[atom1].y + force[atom2].y;
comFrc.z = force[atom1].z + force[atom2].z;
mixed3 relFrc;
relFrc.x = mass1fract*force[atom2].x - mass2fract*force[atom1].x;
relFrc.y = mass1fract*force[atom2].y - mass2fract*force[atom1].y;
relFrc.z = mass1fract*force[atom2].z - mass2fract*force[atom1].z;
comVel.x += comFrc.x * 0.5f * dtVel * invTotMass;
comVel.y += comFrc.y * 0.5f * dtVel * invTotMass;
comVel.z += comFrc.z * 0.5f * dtVel * invTotMass;
relVel.x += relFrc.x * 0.5f * dtVel * invRedMass;
relVel.y += relFrc.y * 0.5f * dtVel * invRedMass;
relVel.z += relFrc.z * 0.5f * dtVel * invRedMass;
#ifdef USE_MIXED_PRECISION
real4 posv1 = posq[atom1];
real4 posv2 = posq[atom2];
real4 posc1 = posqCorrection[atom1];
real4 posc2 = posqCorrection[atom2];
mixed4 pos1 = (mixed4) (posv1.x+(mixed)posc1.x, posv1.y+(mixed)posc1.y, posv1.z+(mixed)posc1.z, posv1.w);
mixed4 pos2 = (mixed4) (posv2.x+(mixed)posc2.x, posv2.y+(mixed)posc2.y, posv2.z+(mixed)posc2.z, posv2.w);
#else
real4 pos1 = posq[atom1];
real4 pos2 = posq[atom2];
#endif
if (v1.w != 0.0f) {
v1.x = comVel.x - relVel.x*mass2fract;
v1.y = comVel.y - relVel.y*mass2fract;
v1.z = comVel.z - relVel.z*mass2fract;
pos1.x = v1.x*dtPos;
pos1.y = v1.y*dtPos;
pos1.z = v1.z*dtPos;
posDelta[atom1] = pos1;
velm[atom1] = v1;
}
if (v2.w != 0.0f) {
v2.x = comVel.x + relVel.x*mass1fract;
v2.y = comVel.y + relVel.y*mass1fract;
v2.z = comVel.z + relVel.z*mass1fract;
pos2.x = v2.x*dtPos;
pos2.y = v2.y*dtPos;
pos2.z = v2.z*dtPos;
posDelta[atom2] = pos2;
velm[atom2] = v2;
}
index += get_global_size(0);
}
}
/**
* Perform the second step of Velocity Verlet integration.
*/
__kernel void integrateVelocityVerletPart2(int numAtoms, __global mixed2* restrict dt, __global real4* restrict posq,
__global real4* restrict posqCorrection, __global mixed4* restrict velm, __global const mixed4* restrict posDelta) {
mixed2 stepSize = dt[0];
int index = get_global_id(0);
if (index == 0)
dt[0].x = stepSize.y;
while(index < numAtoms) {
mixed4 velocity = velm[index];
if (velocity.w != 0.0) {
#ifdef USE_MIXED_PRECISION
real4 pos1 = posq[index];
real4 pos2 = posqCorrection[index];
mixed4 pos = (mixed4) (pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
#else
real4 pos = posq[index];
#endif
mixed4 delta = posDelta[index];
pos.xyz += delta.xyz;
#ifdef USE_MIXED_PRECISION
posq[index] = (real4) ((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
posqCorrection[index] = (real4) (pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
#else
posq[index] = pos;
#endif
}
index += get_global_size(0);
}
}
/**
* Perform the third step of Velocity Verlet integration.
*/
__kernel void integrateVelocityVerletPart3(int numAtoms, int numPairs, int paddedNumAtoms, __global mixed2* restrict dt, __global real4* restrict posq,
__global real4* restrict posqCorrection, __global mixed4* restrict velm, __global const real4* restrict force, __global const mixed4* restrict posDelta,
__global const int* restrict atomList, __global const int2* __restrict__ pairList) {
mixed2 stepSize = dt[0];
#ifndef SUPPORTS_DOUBLE_PRECISION
double oneOverDt = 1.0/stepSize.y;
#else
float oneOverDt = 1.0f/stepSize.y;
float correction = (1.0f-oneOverDt*stepSize.y)/stepSize.y;
#endif
const mixed dtVel = 0.5f*(stepSize.x+stepSize.y);
int index = get_global_id(0);
if (index == 0)
dt[0].x = stepSize.y;
while(index < numAtoms) {
int atom = atomList[index];
mixed4 velocity = velm[atom];
if (velocity.w != 0.0) {
mixed4 deltaXconstrained = posDelta[atom];
velocity.x += 0.5f * dtVel*force[atom].x*velocity.w + (deltaXconstrained.x - velocity.x*stepSize.y)*oneOverDt;
velocity.y += 0.5f * dtVel*force[atom].y*velocity.w + (deltaXconstrained.y - velocity.y*stepSize.y)*oneOverDt;
velocity.z += 0.5f * dtVel*force[atom].z*velocity.w + (deltaXconstrained.z - velocity.z*stepSize.y)*oneOverDt;
#ifdef SUPPORTS_DOUBLE_PRECISION
velocity.x += (deltaXconstrained.x - velocity.x*stepSize.y)*correction;
velocity.y += (deltaXconstrained.y - velocity.y*stepSize.y)*correction;
velocity.z += (deltaXconstrained.z - velocity.z*stepSize.y)*correction;
#endif
velm[atom] = velocity;
}
index += get_global_size(0);
}
index = get_global_id(0);
while(index < numPairs) {
int atom1 = pairList[index].x;
int atom2 = pairList[index].y;
mixed4 v1 = velm[atom1];
mixed4 v2 = velm[atom2];
mixed m1 = v1.w == 0.0f ? 0.0f : 1.0f / v1.w;
mixed m2 = v2.w == 0.0f ? 0.0f : 1.0f / v2.w;
mixed mass1fract = m1 / (m1 + m2);
mixed mass2fract = m2 / (m1 + m2);
mixed invRedMass = (m1 * m2 != 0.0f) ? (m1 + m2)/(m1 * m2) : 0.0f;
mixed invTotMass = (m1 + m2 != 0.0f) ? 1.0f /(m1 + m2) : 0.0f;
mixed3 comVel;
comVel.x= v1.x*mass1fract + v2.x*mass2fract;
comVel.y= v1.y*mass1fract + v2.y*mass2fract;
comVel.z= v1.z*mass1fract + v2.z*mass2fract;
mixed3 relVel;
relVel.x= v2.x - v1.x;
relVel.y= v2.y - v1.y;
relVel.z= v2.z - v1.z;
mixed3 comFrc;
comFrc.x = force[atom1].x + force[atom2].x;
comFrc.y = force[atom1].y + force[atom2].y;
comFrc.z = force[atom1].z + force[atom2].z;
mixed3 relFrc;
relFrc.x = mass1fract*force[atom2].x - mass2fract*force[atom1].x;
relFrc.y = mass1fract*force[atom2].y - mass2fract*force[atom1].y;
relFrc.z = mass1fract*force[atom2].z - mass2fract*force[atom1].z;
comVel.x += comFrc.x * 0.5f * dtVel * invTotMass;
comVel.y += comFrc.y * 0.5f * dtVel * invTotMass;
comVel.z += comFrc.z * 0.5f * dtVel * invTotMass;
relVel.x += relFrc.x * 0.5f * dtVel * invRedMass;
relVel.y += relFrc.y * 0.5f * dtVel * invRedMass;
relVel.z += relFrc.z * 0.5f * dtVel * invRedMass;
if (v1.w != 0.0f) {
mixed4 deltaXconstrained = posDelta[atom1];
v1.x = comVel.x - relVel.x*mass2fract + (deltaXconstrained.x - v1.x*stepSize.y)*oneOverDt;
v1.y = comVel.y - relVel.y*mass2fract + (deltaXconstrained.y - v1.y*stepSize.y)*oneOverDt;
v1.z = comVel.z - relVel.z*mass2fract + (deltaXconstrained.z - v1.z*stepSize.y)*oneOverDt;
#ifdef SUPPORTS_DOUBLE_PRECISION
v1.x += (deltaXconstrained.x - v1.x*stepSize.y)*correction;
v1.y += (deltaXconstrained.y - v1.y*stepSize.y)*correction;
v1.z += (deltaXconstrained.z - v1.z*stepSize.y)*correction;
#endif
velm[atom1] = v1;
}
if (v2.w != 0.0f) {
mixed4 deltaXconstrained = posDelta[atom2];
v2.x = comVel.x + relVel.x*mass1fract + (deltaXconstrained.x - v2.x*stepSize.y)*oneOverDt;
v2.y = comVel.y + relVel.y*mass1fract + (deltaXconstrained.y - v2.y*stepSize.y)*oneOverDt;
v2.z = comVel.z + relVel.z*mass1fract + (deltaXconstrained.z - v2.z*stepSize.y)*oneOverDt;
#ifdef SUPPORTS_DOUBLE_PRECISION
v2.x += (deltaXconstrained.x - v2.x*stepSize.y)*correction;
v2.y += (deltaXconstrained.y - v2.y*stepSize.y)*correction;
v2.z += (deltaXconstrained.z - v2.z*stepSize.y)*correction;
#endif
velm[atom2] = v2;
}
index += get_global_size(0);
}
}
__kernel void integrateVelocityVerletHardWall(int numPairs, __global const float* restrict maxPairDistance,
__global mixed2* restrict dt, __global real4* restrict posq,
__global real4* restrict posqCorrection, __global mixed4* restrict velm,
__global const int2* restrict pairList, __global const float* __restrict__ pairTemperature) {
mixed dtPos = dt[0].y;
mixed maxDelta = (mixed) maxPairDistance[0];
if (maxDelta > 0){
int index = get_global_id(0);
while(index < numPairs) {
const mixed hardWallScale = sqrt( ((mixed) pairTemperature[index]) * ((mixed) BOLTZ));
int2 atom = (int2) (pairList[index].x, pairList[index].y);
#ifdef USE_MIXED_PRECISION
real4 posv1 = posq[atom.x];
real4 posc1 = posqCorrection[atom.x];
mixed4 pos1 = (mixed4) (posv1.x+(mixed)posc1.x, posv1.y+(mixed)posc1.y, posv1.z+(mixed)posc1.z, posv1.w);
real4 posv2 = posq[atom.y];
real4 posc2 = posqCorrection[atom.y];
mixed4 pos2 = (mixed4) (posv2.x+(mixed)posc2.x, posv2.y+(mixed)posc2.y, posv2.z+(mixed)posc2.z, posv2.w);
#else
real4 pos1 = posq[atom.x];
real4 pos2 = posq[atom.y];
#endif
mixed3 delta = (mixed3) (
(mixed) (pos1.x - pos2.x),
(mixed) (pos1.y - pos2.y),
(mixed) (pos1.z - pos2.z)
);
mixed r = sqrt(delta.x*delta.x + delta.y*delta.y + delta.z*delta.z);
mixed rInv = 1/r;
if (rInv*maxDelta < 1.0) {
// The constraint has been violated, so make the inter-particle distance "bounce"
// off the hard wall.
mixed3 bondDir = (mixed3) (delta.x * rInv, delta.y * rInv, delta.z * rInv);
mixed3 vel1 = (mixed3) (velm[atom.x].x, velm[atom.x].y, velm[atom.x].z);
mixed3 vel2 = (mixed3) (velm[atom.y].x, velm[atom.y].y, velm[atom.y].z);
mixed m1 = velm[atom.x].w != 0.0 ? 1.0/velm[atom.x].w : 0.0;
mixed m2 = velm[atom.y].w != 0.0 ? 1.0/velm[atom.y].w : 0.0;
mixed invTotMass = (m1 + m2 != 0.0) ? 1.0 /(m1 + m2) : 0.0;
mixed deltaR = r-maxDelta;
mixed deltaT = dtPos;
mixed dt = dtPos;
mixed dotvr1 = vel1.x*bondDir.x + vel1.y*bondDir.y + vel1.z*bondDir.z;
mixed3 vb1 = (mixed3) (bondDir.x*dotvr1, bondDir.y*dotvr1, bondDir.z*dotvr1);
mixed3 vp1 = (mixed3) (vel1.x-vb1.x, vel1.y-vb1.y, vel1.z-vb1.z);
if (m2 == 0) {
// The parent particle is massless, so move only the Drude particle.
if (dotvr1 != 0.0)
deltaT = deltaR/fabs(dotvr1);
if (deltaT > dtPos)
deltaT = dtPos;
dotvr1 = -dotvr1*hardWallScale/(fabs(dotvr1)*sqrt(m1));
mixed dr = -deltaR + deltaT*dotvr1;
pos1.x += bondDir.x*dr;
pos1.y += bondDir.y*dr;
pos1.z += bondDir.z*dr;
velm[atom.x] = (mixed4) (vp1.x + bondDir.x*dotvr1, vp1.y + bondDir.y*dotvr1, vp1.z + bondDir.z*dotvr1, velm[atom.x].w);
#ifdef USE_MIXED_PRECISION
posq[atom.x] = (real4) ((real) pos1.x, (real) pos1.y, (real) pos1.z, (real) pos1.w);
posqCorrection[atom.x] = (real4) (pos1.x-(real) pos1.x, pos1.y-(real) pos1.y, pos1.z-(real) pos1.z, 0);
#else
posq[atom.x] = pos1;
#endif
}
else {
// Move both particles.
mixed dotvr2 = vel2.x*bondDir.x + vel2.y*bondDir.y + vel2.z*bondDir.z;
mixed3 vb2 = (mixed3) (bondDir.x*dotvr2, bondDir.y*dotvr2, bondDir.z*dotvr2);
mixed3 vp2 = (mixed3) (vel2.x-vb2.x, vel2.y-vb2.y, vel2.z-vb2.z);
mixed vbCMass = (m1*dotvr1 + m2*dotvr2)*invTotMass;
dotvr1 -= vbCMass;
dotvr2 -= vbCMass;
if (dotvr1 != dotvr2)
deltaT = deltaR/fabs(dotvr1-dotvr2);
if (deltaT > dt)
deltaT = dt;
mixed vBond = hardWallScale/sqrt(m1);
dotvr1 = -dotvr1*vBond*m2*invTotMass/fabs(dotvr1);
dotvr2 = -dotvr2*vBond*m1*invTotMass/fabs(dotvr2);
mixed dr1 = -deltaR*m2*invTotMass + deltaT*dotvr1;
mixed dr2 = deltaR*m1*invTotMass + deltaT*dotvr2;
dotvr1 += vbCMass;
dotvr2 += vbCMass;
pos1.x += bondDir.x*dr1;
pos1.y += bondDir.y*dr1;
pos1.z += bondDir.z*dr1;
pos2.x += bondDir.x*dr2;
pos2.y += bondDir.y*dr2;
pos2.z += bondDir.z*dr2;
velm[atom.x] = (mixed4) (vp1.x + bondDir.x*dotvr1, vp1.y + bondDir.y*dotvr1, vp1.z + bondDir.z*dotvr1, velm[atom.x].w);
velm[atom.y] = (mixed4) (vp2.x + bondDir.x*dotvr2, vp2.y + bondDir.y*dotvr2, vp2.z + bondDir.z*dotvr2, velm[atom.y].w);
#ifdef USE_MIXED_PRECISION
posq[atom.x] = (real4) ((real) pos1.x, (real) pos1.y, (real) pos1.z, (real) pos1.w);
posq[atom.y] = (real4) ((real) pos2.x, (real) pos2.y, (real) pos2.z, (real) pos2.w);
posqCorrection[atom.x] = (real4) (pos1.x-(real) pos1.x, pos1.y-(real) pos1.y, pos1.z-(real) pos1.z, 0);
posqCorrection[atom.y] = (real4) (pos2.x-(real) pos2.x, pos2.y-(real) pos2.y, pos2.z-(real) pos2.z, 0);
#else
posq[atom.x] = pos1;
posq[atom.y] = pos2;
#endif
}
}
index += get_global_size(0);
}
}
}
platforms/opencl/tests/TestOpenCLNoseHooverIntegrator.cpp 0 → 100644
View file @ f6431a42
/* -------------------------------------------------------------------------- *
* 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) 2019 Stanford University and the Authors. *
* Authors: Andreas Krämer and Andrew C. Simmmonett *
* 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 "OpenCLTests.h"
#include "TestNoseHooverIntegrator.h"
void runPlatformTests() {
}
platforms/opencl/tests/TestOpenCLNoseHooverThermostat.cpp 0 → 100644
View file @ f6431a42
/* -------------------------------------------------------------------------- *
* 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) 2019 Stanford University and the Authors. *
* Authors: Andreas Krämer and Andrew C. Simmonett *
* 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 "OpenCLTests.h"
#include "TestNoseHooverThermostat.h"
void runPlatformTests() {
}
platforms/reference/include/ReferenceConstraints.h
View file @ f6431a42
......@@ -61,7 +61,7 @@ public:
* Apply the constraint algorithm to velocities.
*
* @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify
* @param velocities the velocities to modify
* @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/
......
platforms/reference/include/ReferenceKernels.h
View file @ f6431a42
......@@ -59,7 +59,9 @@ class ReferenceGayBerneForce;
class ReferenceBrownianDynamics;
class ReferenceStochasticDynamics;
class ReferenceConstraintAlgorithm;
class ReferenceNoseHooverChain;
class ReferenceMonteCarloBarostat;
class ReferenceVelocityVerletDynamics;
class ReferenceVariableStochasticDynamics;
class ReferenceVariableVerletDynamics;
class ReferenceVerletDynamics;
......@@ -1132,6 +1134,45 @@ private:
ReferenceVerletDynamics* dynamics;
std::vector<double> masses;
double prevStepSize;
}
;
/**
* This kernel is invoked by NoseHooverIntegrator to take one time step.
*/
class ReferenceIntegrateVelocityVerletStepKernel : public IntegrateVelocityVerletStepKernel {
public:
ReferenceIntegrateVelocityVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVelocityVerletStepKernel(name, platform),
data(data), dynamics(0) {
}
~ReferenceIntegrateVelocityVerletStepKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the NoseHooverIntegrator this kernel will be used for
*/
void initialize(const System& system, const NoseHooverIntegrator& integrator);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the VerletIntegrator this kernel is being used for
* @param forcesAreValid a reference to the parent integrator's boolean for keeping
* track of the validity of the current forces.
*/
void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the NoseHooverIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
private:
ReferencePlatform::PlatformData& data;
ReferenceVelocityVerletDynamics* dynamics;
std::vector<double> masses;
double prevStepSize;
};
/**
......@@ -1429,6 +1470,62 @@ private:
std::vector<double> masses;
};
/**
* This kernel is invoked by NoseHooverChain at the start of each time step to adjust the thermostat
* and update the associated particle velocities.
*/
class ReferenceNoseHooverChainKernel : public NoseHooverChainKernel {
public:
ReferenceNoseHooverChainKernel(std::string name, const Platform& platform) : NoseHooverChainKernel(name, platform), chainPropagator(0) {
}
~ReferenceNoseHooverChainKernel();
/**
* Initialize the kernel.
*/
virtual void initialize();
/**
* Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the object describing the chain to be propagated.
* @param kineticEnergy the {absolute, relative} kinetic energies of the particles being thermostated by this chain.
* @param timeStep the time step used by the integrator.
* @return the velocity scale factors to apply to the {absolute, relative} motion of particles associated with this heat bath.
*/
virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergy, double timeStep);
/**
* Execute the kernal that computes the total (kinetic + potential) heat bath energy.
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the chain whose energy is to be determined.
* @return the total heat bath energy.
*/
virtual double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc);
/**
* Execute the kernel that computes the kinetic energy for a subset of atoms,
* or the relative kinetic energy of Drude particles with respect to their parent atoms
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the chain whose energy is to be determined.
* @param downloadValue whether the computed value should be downloaded and returned.
*
*/
virtual std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
/**
* Execute the kernel that scales the velocities of particles associated with a nose hoover chain
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the chain whose energy is to be determined.
* @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
*/
virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor);
private:
ReferenceNoseHooverChain* chainPropagator;
};
/**
* This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
*/
......
platforms/reference/include/ReferenceNoseHooverChain.h 0 → 100644
View file @ f6431a42
/* Portions copyright (c) 2019 Stanford University and Simbios.
* Contributors: Andreas Krämer and Andrew C. Simmonett
*
* 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.
*/
#ifndef __ReferenceNoseHooverChain_H__
#define __ReferenceNoseHooverChain_H__
#include "openmm/Vec3.h"
#include <vector>
namespace OpenMM {
using std::vector;
class ReferenceNoseHooverChain {
private:
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceNoseHooverChain();
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceNoseHooverChain();
/**---------------------------------------------------------------------------------------
Propagate the Nose-Hoover chain a half timestep and find the appropriate velocity scaling
@param kineticEnergy the instantaneous kinetic energy of the particles being thermostated
@param chainVelocities the velocities of the chain's beads in nm / ps
@param chainPositions the positions of the chains's beads in nm
@param numDOFs the number of degrees of freedom in the system that this chain thermostats
@param temperature thermostat temperature in Kelvin
@param collisionFrequency collision frequency for each atom in ps^-1
@param timeStep full integration step size in ps (this only propagates half way)
@param numMTS number of multi timestep increments used in the Trotter expansion
@param YSWeights vector of weights used in the Yoshida-Suzuki multi-timestepping.
--------------------------------------------------------------------------------------- */
double propagate(double kineticEnergy, vector<double>& chainVelocities,
vector<double>& chainPositions, int numDOFs,
double temperature, double collisionFrequency, double timeStep,
int numMTS, const vector<double>& YSWeights) const;
};
} // namespace OpenMM
#endif // __ReferenceNoseHooverChain_H__
platforms/reference/include/ReferencePlatform.h
View file @ f6431a42
......@@ -72,6 +72,8 @@ public:
Vec3* periodicBoxVectors;
ReferenceConstraints* constraints;
std::map<std::string, double>* energyParameterDerivatives;
std::vector<std::vector<double>>* noseHooverPositions;
std::vector<std::vector<double>>* noseHooverVelocities;
};
} // namespace OpenMM
......
platforms/reference/include/ReferenceSETTLEAlgorithm.h
View file @ f6431a42
......@@ -75,7 +75,7 @@ public:
* Apply the constraint algorithm to velocities.
*
* @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify
* @param velocities the velocities to modify
* @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/
......
platforms/reference/include/ReferenceVelocityVerletDynamics.h 0 → 100644
View file @ f6431a42
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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.
*/
#ifndef __ReferenceVelocityVerletDynamics_H__
#define __ReferenceVelocityVerletDynamics_H__
#include "ReferenceDynamics.h"
namespace OpenMM {
class ContextImpl;
class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
private:
std::vector<OpenMM::Vec3> xPrime;
std::vector<double> inverseMasses;
public:
/**---------------------------------------------------------------------------------------
Constructor
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param friction friction coefficient
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceVelocityVerletDynamics(int numberOfAtoms, double deltaT);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceVelocityVerletDynamics();
/**---------------------------------------------------------------------------------------
Update
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param tolerance the constraint tolerance
@param forcesAreValid whether the forces are valid (duh!)
@param allAtoms a list of all atoms not involved in a Drude-like pair
@param allPairs a list of all Drude-like pairs, and their KT values, in the system
@param maxPairDistance the maximum separation allowed for a Drude-like pair
--------------------------------------------------------------------------------------- */
void update(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
};
} // namespace OpenMM
#endif // __ReferenceVelocityVerletDynamics_H__
platforms/reference/src/ReferenceKernelFactory.cpp
View file @ f6431a42
......@@ -88,6 +88,10 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcGayBerneForceKernel(name, platform);
if (name == IntegrateVerletStepKernel::Name())
return new ReferenceIntegrateVerletStepKernel(name, platform, data);
if (name == IntegrateVelocityVerletStepKernel::Name())
return new ReferenceIntegrateVelocityVerletStepKernel(name, platform, data);
if (name == NoseHooverChainKernel::Name())
return new ReferenceNoseHooverChainKernel(name, platform);
if (name == IntegrateLangevinStepKernel::Name())
return new ReferenceIntegrateLangevinStepKernel(name, platform, data);
if (name == IntegrateBAOABStepKernel::Name())
......
platforms/reference/src/ReferenceKernels.cpp
View file @ f6431a42
......@@ -56,6 +56,7 @@
#include "ReferenceLJCoulomb14.h"
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceMonteCarloBarostat.h"
#include "ReferenceNoseHooverChain.h"
#include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h"
#include "ReferenceRMSDForce.h"
......@@ -63,6 +64,7 @@
#include "ReferenceTabulatedFunction.h"
#include "ReferenceVariableStochasticDynamics.h"
#include "ReferenceVariableVerletDynamics.h"
#include "ReferenceVelocityVerletDynamics.h"
#include "ReferenceVerletDynamics.h"
#include "ReferenceVirtualSites.h"
#include "openmm/CMMotionRemover.h"
......@@ -125,6 +127,15 @@ static map<string, double>& extractEnergyParameterDerivatives(ContextImpl& conte
return *data->energyParameterDerivatives;
}
static vector<vector<double> >& extractNoseHooverPositions(ContextImpl& context) {
ReferencePlatform::PlatformData *data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<vector<double> >*) data->noseHooverPositions);
}
static vector<vector<double> >& extractNoseHooverVelocities(ContextImpl& context) {
ReferencePlatform::PlatformData *data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<vector<double> >*) data->noseHooverVelocities);
}
/**
* Make sure an expression doesn't use any undefined variables.
*/
......@@ -277,7 +288,7 @@ void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context,
}
void ReferenceUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
int version = 2;
int version = 3;
stream.write((char*) &version, sizeof(int));
stream.write((char*) &data.time, sizeof(data.time));
vector<Vec3>& posData = extractPositions(context);
......@@ -286,13 +297,27 @@ void ReferenceUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostr
stream.write((char*) &velData[0], sizeof(Vec3)*velData.size());
Vec3* vectors = extractBoxVectors(context);
stream.write((char*) vectors, 3*sizeof(Vec3));
auto& allNoseHooverPositions = extractNoseHooverPositions(context);
auto& allNoseHooverVelocities = extractNoseHooverVelocities(context);
size_t numChains = allNoseHooverPositions.size();
assert(numChains == allNoseHooverVelocities.size());
stream.write((char*) &numChains, sizeof(size_t));
for (size_t i=0; i<numChains; i++){
auto & noseHooverPositions = allNoseHooverPositions.at(i);
auto & noseHooverVelocities = allNoseHooverVelocities.at(i);
size_t numBeads = noseHooverPositions.size();
assert(numBeads == noseHooverVelocities.size());
stream.write((char*) &numBeads, sizeof(size_t));
stream.write((char*) noseHooverPositions.data(), sizeof(double)*numBeads);
stream.write((char*) noseHooverVelocities.data(), sizeof(double)*numBeads);
}
SimTKOpenMMUtilities::createCheckpoint(stream);
}
void ReferenceUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& stream) {
int version;
stream.read((char*) &version, sizeof(int));
if (version != 2)
if (version != 3)
throw OpenMMException("Checkpoint was created with a different version of OpenMM");
stream.read((char*) &data.time, sizeof(data.time));
vector<Vec3>& posData = extractPositions(context);
......@@ -301,6 +326,21 @@ void ReferenceUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istrea
stream.read((char*) &velData[0], sizeof(Vec3)*velData.size());
Vec3* vectors = extractBoxVectors(context);
stream.read((char*) vectors, 3*sizeof(Vec3));
size_t numChains, numBeads;
auto& allNoseHooverPositions = extractNoseHooverPositions(context);
auto& allNoseHooverVelocities = extractNoseHooverVelocities(context);
stream.read((char*) &numChains, sizeof(size_t));
allNoseHooverPositions.clear();
allNoseHooverVelocities.clear();
for (size_t i=0; i<numChains; i++){
stream.read((char*) &numBeads, sizeof(size_t));
std::vector<double> noseHooverPositions(numBeads);
std::vector<double> noseHooverVelocities(numBeads);
stream.read((char*) &noseHooverPositions[0], sizeof(double)*numBeads);
stream.read((char*) &noseHooverVelocities[0], sizeof(double)*numBeads);
allNoseHooverPositions.push_back(noseHooverPositions);
allNoseHooverVelocities.push_back(noseHooverVelocities);
}
SimTKOpenMMUtilities::loadCheckpoint(stream);
}
......@@ -2106,6 +2146,41 @@ double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& con
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
}
ReferenceIntegrateVelocityVerletStepKernel::~ReferenceIntegrateVelocityVerletStepKernel() {
if (dynamics)
delete dynamics;
}
void ReferenceIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
int numParticles = system.getNumParticles();
masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
}
void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
double stepSize = integrator.getStepSize();
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& velData = extractVelocities(context);
vector<Vec3>& forceData = extractForces(context);
if (dynamics == 0 || stepSize != prevStepSize) {
// Recreate the computation objects with the new parameters.
if (dynamics)
delete dynamics;
dynamics = new ReferenceVelocityVerletDynamics(context.getSystem().getNumParticles(), stepSize);
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevStepSize = stepSize;
}
dynamics->update(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
data.time += stepSize;
data.stepCount++;
}
double ReferenceIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0);
}
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
if (dynamics)
delete dynamics;
......@@ -2430,6 +2505,195 @@ void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
context.getIntegrator().getStepSize());
}
ReferenceNoseHooverChainKernel::~ReferenceNoseHooverChainKernel() {
if (chainPropagator)
delete chainPropagator;
}
void ReferenceNoseHooverChainKernel::initialize() {
this->chainPropagator = new ReferenceNoseHooverChain();
//SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
}
std::pair<double, double> ReferenceNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergy, double timeStep) {
double absKE = kineticEnergy.first;
double relKE = kineticEnergy.second;
if (absKE < 1e-8) return {1.0, 1.0}; // (catches the problem of zero velocities in the first dynamics step, where we have nothing to scale)
// Get the variables describing the NHC
int chainLength = nhc.getChainLength();
int chainID = nhc.getChainID();
int numDOFs = nhc.getNumDegreesOfFreedom();
int numMTS = nhc.getNumMultiTimeSteps();
// Get the state of the NHC from the context
auto& allChainPositions = extractNoseHooverPositions(context);
auto& allChainVelocities = extractNoseHooverVelocities(context);
int nAtoms = nhc.getThermostatedAtoms().size();
double absScale = 0;
if (nAtoms) {
if (allChainPositions.size() < 2*chainID+1){
allChainPositions.resize(2*chainID+1);
}
if (allChainVelocities.size() < 2*chainID+1){
allChainVelocities.resize(2*chainID+1);
}
auto& chainPositions = allChainPositions.at(2*chainID);
auto& chainVelocities = allChainVelocities.at(2*chainID);
if (chainPositions.size() < chainLength){
chainPositions.resize(chainLength, 0);
}
if (chainVelocities.size() < chainLength){
chainVelocities.resize(chainLength, 0);
}
double temperature = nhc.getTemperature();
double collisionFrequency = nhc.getCollisionFrequency();
absScale = chainPropagator->propagate(absKE, chainVelocities, chainPositions, numDOFs,
temperature, collisionFrequency, timeStep,
numMTS, nhc.getYoshidaSuzukiWeights());
}
double relScale = 0;
int nPairs = nhc.getThermostatedPairs().size();
if (nPairs) {
if (allChainPositions.size() < 2*chainID+2){
allChainPositions.resize(2*chainID+2);
}
if (allChainVelocities.size() < 2*chainID+2){
allChainVelocities.resize(2*chainID+2);
}
auto& chainPositions = allChainPositions.at(2*chainID+1);
auto& chainVelocities = allChainVelocities.at(2*chainID+1);
if (chainPositions.size() < chainLength){
chainPositions.resize(chainLength, 0);
}
if (chainVelocities.size() < chainLength){
chainVelocities.resize(chainLength, 0);
}
double temperature = nhc.getRelativeTemperature();
double collisionFrequency = nhc.getRelativeCollisionFrequency();
relScale = chainPropagator->propagate(relKE, chainVelocities, chainPositions, 3*nPairs,
temperature, collisionFrequency, timeStep,
numMTS, nhc.getYoshidaSuzukiWeights());
}
return {absScale, relScale};
}
double ReferenceNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
double potentialEnergy = 0;
double kineticEnergy = 0;
int chainLength = nhc.getChainLength();
int chainID = nhc.getChainID();
int nAtoms = nhc.getThermostatedAtoms().size();
int nPairs = nhc.getThermostatedPairs().size();
auto& nhcPositions = extractNoseHooverPositions(context);
auto& nhcVelocities = extractNoseHooverVelocities(context);
if (nAtoms) {
double temperature = nhc.getTemperature();
double collisionFrequency = nhc.getCollisionFrequency();
double kT = temperature * BOLTZ;
int numDOFs = nhc.getNumDegreesOfFreedom();
for(int i = 0; i < chainLength; ++i) {
double prefac = i ? 1 : numDOFs;
double mass = prefac * kT / (collisionFrequency * collisionFrequency);
double velocity = nhcVelocities[2*chainID][i];
// The kinetic energy of this bead
kineticEnergy += 0.5 * mass * velocity * velocity;
// The potential energy of this bead
double position = nhcPositions[2*chainID][i];
potentialEnergy += prefac * kT * position;
}
}
if (nPairs) {
double temperature = nhc.getRelativeTemperature();
double collisionFrequency = nhc.getRelativeCollisionFrequency();
double kT = temperature * BOLTZ;
int numDOFs = 3 * nPairs;
for(int i = 0; i < chainLength; ++i) {
double prefac = i ? 1 : numDOFs;
double mass = prefac * kT / (collisionFrequency * collisionFrequency);
double velocity = nhcVelocities[2*chainID+1][i];
// The kinetic energy of this bead
kineticEnergy += 0.5 * mass * velocity * velocity;
// The potential energy of this bead
double position = nhcPositions[2*chainID+1][i];
potentialEnergy += prefac * kT * position;
}
}
return kineticEnergy + potentialEnergy;
}
std::pair<double, double> ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) {
const std::vector<int>& atomsList = noseHooverChain.getThermostatedAtoms();
const std::vector<std::pair<int,int>>& pairsList = noseHooverChain.getThermostatedPairs();
std::vector<Vec3>& velocities = extractVelocities(context);
const System& system = context.getSystem();
int numParticles = system.getNumParticles();
std::vector<double> masses(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
double comKE = 0;
double relKE = 0;
// kinetic energy of individual atoms
for (const auto &m: atomsList){
comKE += 0.5 * masses[m] * velocities[m].dot(velocities[m]);
}
// center of mass kinetic energy of pairs
for (const auto &p: pairsList){
double m1 = masses[p.first];
double m2 = masses[p.second];
Vec3 v1 = velocities[p.first];
Vec3 v2 = velocities[p.second];
double invMass = 1.0 / (m1 + m2);
double redMass = m1 * m2 * invMass;
double fracM1 = m1 * invMass;
double fracM2 = m2 * invMass;
Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
Vec3 relVelocity = v2 - v1;
comKE += 0.5 * (m1 + m2) * comVelocity.dot(comVelocity);
relKE += 0.5 * redMass * relVelocity.dot(relVelocity);
}
// We ignore the downloadValue argument here and always return the correct value
return {comKE, relKE};
}
void ReferenceNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactors) {
const auto& atoms = noseHooverChain.getThermostatedAtoms();
const auto& pairs = noseHooverChain.getThermostatedPairs();
std::vector<Vec3>& velocities = extractVelocities(context);
double absScale = scaleFactors.first;
double relScale = scaleFactors.second;
const System& system = context.getSystem();
int numParticles = system.getNumParticles();
std::vector<double> masses(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
// scale absolute velocities
for (const auto &a: atoms){
velocities[a] *= absScale;
}
// scale relative velocities and absolute center of mass velocities for each pair
for (const auto &p: pairs){
int p1 = p.first;
int p2 = p.second;
double m1 = masses[p.first];
double m2 = masses[p.second];
Vec3 v1 = velocities[p.first];
Vec3 v2 = velocities[p.second];
double invMass = 1.0 / (m1 + m2);
double fracM1 = m1 * invMass;
double fracM2 = m2 * invMass;
Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
Vec3 relVelocity = v2 - v1;
velocities[p1] = absScale * comVelocity - relScale * relVelocity * fracM2;
velocities[p2] = absScale * comVelocity + relScale * relVelocity * fracM1;
}
}
ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
if (barostat)
delete barostat;
......
platforms/reference/src/ReferencePlatform.cpp
View file @ f6431a42
......@@ -66,6 +66,8 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateVelocityVerletStepKernel::Name(), factory);
registerKernelFactory(NoseHooverChainKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateBAOABStepKernel::Name(), factory);
registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
......@@ -102,6 +104,8 @@ ReferencePlatform::PlatformData::PlatformData(const System& system) : time(0.0),
periodicBoxVectors = new Vec3[3];
constraints = new ReferenceConstraints(system);
energyParameterDerivatives = new map<string, double>();
noseHooverPositions = new vector<vector<double> >();
noseHooverVelocities = new vector<vector<double> >();
}
ReferencePlatform::PlatformData::~PlatformData() {
......@@ -112,4 +116,6 @@ ReferencePlatform::PlatformData::~PlatformData() {
delete[] periodicBoxVectors;
delete constraints;
delete energyParameterDerivatives;
delete noseHooverPositions;
delete noseHooverVelocities;
}
platforms/reference/src/SimTKReference/ReferenceNoseHooverChain.cpp 0 → 100644
View file @ f6431a42
/* Portions copyright (c) 2008-2010 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* 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 <cmath>
#include <string.h>
#include <sstream>
#include <exception>
#include "SimTKOpenMMUtilities.h"
#include "ReferenceNoseHooverChain.h"
using std::vector;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceNoseHooverChain::ReferenceNoseHooverChain() {
}
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
ReferenceNoseHooverChain::~ReferenceNoseHooverChain() {
}
double ReferenceNoseHooverChain::propagate(double kineticEnergy, vector<double>& chainVelocities,
vector<double>& chainPositions, int numDOFs,
double temperature, double collisionFrequency, double timeStep,
int numMTS, const vector<double>& YSWeights) const {
double scale = 1;
const double kT = BOLTZ * temperature;
const size_t chainLength = chainPositions.size();
std::vector<double> chainForces(chainLength, 0);
std::vector<double> chainMasses(chainLength, kT/(collisionFrequency*collisionFrequency));
chainMasses[0] *= numDOFs;
double KE2 = 2 * kineticEnergy;
chainForces[0] = (KE2 - numDOFs * kT) / chainMasses[0];
for (int bead = 0; bead < chainLength - 1; ++bead) {
chainForces[bead + 1] = (chainMasses[bead] * chainVelocities[bead] * chainVelocities[bead] - kT) / chainMasses[bead + 1];
}
for (int mts = 0; mts < numMTS; ++mts) {
for (const auto &ys : YSWeights) {
double wdt = ys * timeStep / numMTS;
chainVelocities.back() += 0.25 * wdt * chainForces.back();
for (int bead = chainLength - 2; bead >= 0; --bead) {
double aa = exp(-0.125 * wdt * chainVelocities[bead + 1]);
chainVelocities[bead] = aa * (chainVelocities[bead] * aa + 0.25 * wdt * chainForces[bead]);
}
// update particle velocities
double aa = exp(-0.5 * wdt * chainVelocities[0]);
scale *= aa;
// update the thermostat positions
for (int bead = 0; bead < chainLength; ++bead) {
chainPositions[bead] += 0.5 * chainVelocities[bead] * wdt;
}
// update the forces
chainForces[0] = (scale * scale * KE2 - numDOFs * kT) / chainMasses[0];
// update thermostat velocities
for (int bead = 0; bead < chainLength - 1; ++bead) {
double aa = exp(-0.125 * wdt * chainVelocities[bead + 1]);
chainVelocities[bead] = aa * (aa * chainVelocities[bead] + 0.25 * wdt * chainForces[bead]);
chainForces[bead + 1] = (chainMasses[bead] * chainVelocities[bead] * chainVelocities[bead] - kT) / chainMasses[bead + 1];
}
chainVelocities[chainLength-1] += 0.25 * wdt * chainForces.back();
} // YS loop
} // MTS loop
return scale;
}
platforms/reference/src/SimTKReference/ReferenceVelocityVerletDynamics.cpp 0 → 100644
View file @ f6431a42
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* 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 <cstring>
#include <iostream>
#include <sstream>
#include "openmm/OpenMMException.h"
#include "SimTKOpenMMUtilities.h"
#include "openmm/internal/ContextImpl.h"
#include "ReferenceVelocityVerletDynamics.h"
#include "ReferenceVirtualSites.h"
#include <cstdio>
using std::vector;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
ReferenceVelocityVerletDynamics constructor
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param friction friction coefficient
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceVelocityVerletDynamics::ReferenceVelocityVerletDynamics(int numberOfAtoms, double deltaT) :
ReferenceDynamics(numberOfAtoms, deltaT, 0.0) {
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
}
/**---------------------------------------------------------------------------------------
ReferenceVelocityVerletDynamics destructor
--------------------------------------------------------------------------------------- */
ReferenceVelocityVerletDynamics::~ReferenceVelocityVerletDynamics() {
}
/**---------------------------------------------------------------------------------------
Update -- driver routine for performing Velocity Verlet dynamics update of coordinates
and velocities
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param atomList list of all atoms not involved in a Drude-like pair
@param pairList list of all Drude-like pairs
@param maxPairDistance the maximum separation of any Drude-like pairs
--------------------------------------------------------------------------------------- */
void ReferenceVelocityVerletDynamics::update(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & atomList, const std::vector<std::tuple<int, int, double>> &pairList,
double maxPairDistance) {
// first-time-through initialization
if (!forcesAreValid) context.calcForcesAndEnergy(true, false);
int numberOfAtoms = system.getNumParticles();
if (getTimeStep() == 0) {
// invert masses
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] == 0.0)
inverseMasses[ii] = 0.0;
else
inverseMasses[ii] = 1.0/masses[ii];
}
}
//// Perform the integration.
// Regular atoms
for (const auto &atom : atomList) {
if (masses[atom] != 0.0) {
velocities[atom] += 0.5 * inverseMasses[atom]*forces[atom]*getDeltaT();
xPrime[atom] = atomCoordinates[atom];
atomCoordinates[atom] += velocities[atom]*getDeltaT();
}
}
// Connected particles
for (const auto &pair : pairList) {
const auto &atom1 = std::get<0>(pair);
const auto &atom2 = std::get<1>(pair);
double m1 = system.getParticleMass(atom1);
double m2 = system.getParticleMass(atom2);
double mass1fract = m1 / (m1 + m2);
double mass2fract = m2 / (m1 + m2);
double invRedMass = (m1 * m2 != 0.0) ? (m1 + m2)/(m1 * m2) : 0.0;
double invTotMass = (m1 + m2 != 0.0) ? 1.0 /(m1 + m2) : 0.0;
Vec3 comVel = velocities[atom1]*mass1fract + velocities[atom2]*mass2fract;
Vec3 relVel = velocities[atom2] - velocities[atom1];
Vec3 comForce = forces[atom1] + forces[atom2];
Vec3 relForce = mass1fract*forces[atom2] - mass2fract*forces[atom1];
comVel += 0.5 * comForce * getDeltaT() * invTotMass;
relVel += 0.5 * relForce * getDeltaT() * invRedMass;
if (m1 != 0.0) {
velocities[atom1] = comVel - relVel*mass2fract;
xPrime[atom1] = atomCoordinates[atom1];
atomCoordinates[atom1] += velocities[atom1]*getDeltaT();
}
if (m2 != 0.0) {
velocities[atom2] = comVel + relVel*mass1fract;
xPrime[atom2] = atomCoordinates[atom2];
atomCoordinates[atom2] += velocities[atom2]*getDeltaT();
}
}
//
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(xPrime, atomCoordinates, inverseMasses, tolerance);
// Apply hard wall constraints.
if (maxPairDistance > 0) {
for (const auto & pair : pairList) {
const int atom1 = std::get<0>(pair);
const int atom2 = std::get<1>(pair);
const double hardWallScale = sqrt(std::get<2>(pair)*BOLTZ);
Vec3 delta = atomCoordinates[atom1]-atomCoordinates[atom2];
double r = sqrt(delta.dot(delta));
double rInv = 1/r;
if (rInv*maxPairDistance < 1.0) {
// The constraint has been violated, so make the inter-particle distance "bounce"
// off the hard wall.
//if (rInv*maxPairDistance < 0.5)
// throw OpenMMException("Drude particle moved too far beyond hard wall constraint");
// TODO: Review this - I commented it out to make the NoseHooverThermostat test work
Vec3 bondDir = delta*rInv;
Vec3 vel1 = velocities[atom1];
Vec3 vel2 = velocities[atom2];
double m1 = masses[atom1];
double m2 = masses[atom2];
double invTotMass = (m1 + m2 != 0.0) ? 1.0 /(m1 + m2) : 0.0;
double deltaR = r-maxPairDistance;
double deltaT = getDeltaT();
double dt = getDeltaT();
double dotvr1 = vel1.dot(bondDir);
Vec3 vb1 = bondDir*dotvr1;
Vec3 vp1 = vel1-vb1;
if (m2 == 0) {
// The parent particle is massless, so move only the Drude particle.
if (dotvr1 != 0.0)
deltaT = deltaR/std::abs(dotvr1);
if (deltaT > getDeltaT())
deltaT = getDeltaT();
dotvr1 = -dotvr1*hardWallScale/(std::abs(dotvr1)*sqrt(m1));
double dr = -deltaR + deltaT*dotvr1;
atomCoordinates[atom1] += bondDir*dr;
velocities[atom1] = vp1 + bondDir*dotvr1;
}
else {
// Move both particles.
double dotvr2 = vel2.dot(bondDir);
Vec3 vb2 = bondDir*dotvr2;
Vec3 vp2 = vel2-vb2;
double vbCMass = (m1*dotvr1 + m2*dotvr2)*invTotMass;
dotvr1 -= vbCMass;
dotvr2 -= vbCMass;
if (dotvr1 != dotvr2)
deltaT = deltaR/std::abs(dotvr1-dotvr2);
if (deltaT > dt)
deltaT = dt;
double vBond = hardWallScale/sqrt(m1);
dotvr1 = -dotvr1*vBond*m2*invTotMass/std::abs(dotvr1);
dotvr2 = -dotvr2*vBond*m1*invTotMass/std::abs(dotvr2);
double dr1 = -deltaR*m2*invTotMass + deltaT*dotvr1;
double dr2 = deltaR*m1*invTotMass + deltaT*dotvr2;
dotvr1 += vbCMass;
dotvr2 += vbCMass;
atomCoordinates[atom1] += bondDir*dr1;
atomCoordinates[atom2] += bondDir*dr2;
velocities[atom1] = vp1 + bondDir*dotvr1;
velocities[atom2] = vp2 + bondDir*dotvr2;
}
}
}
} /* end of hard wall constraint part */
ReferenceVirtualSites::computePositions(system, atomCoordinates);
context.calcForcesAndEnergy(true, false);
forcesAreValid = true;
for (int i = 0; i < numberOfAtoms; ++i) {
if (masses[i] != 0.0)
for (int j = 0; j < 3; ++j) {
xPrime[i][j] += velocities[i][j]*getDeltaT();
}
}
// Update the positions and velocities.
// Regular atoms
for (const auto &atom : atomList) {
if (masses[atom] != 0.0) {
velocities[atom] += 0.5 * inverseMasses[atom]*forces[atom]*getDeltaT() + (atomCoordinates[atom] - xPrime[atom])/getDeltaT();
}
}
// Connected particles
for (const auto &pair : pairList) {
const auto &atom1 = std::get<0>(pair);
const auto &atom2 = std::get<1>(pair);
double m1 = system.getParticleMass(atom1);
double m2 = system.getParticleMass(atom2);
double mass1fract = m1 / (m1 + m2);
double mass2fract = m2 / (m1 + m2);
double invRedMass = (m1 * m2 != 0.0) ? (m1 + m2)/(m1 * m2) : 0.0;
double invTotMass = (m1 + m2 != 0.0) ? 1.0 /(m1 + m2) : 0.0;
Vec3 comVel = velocities[atom1]*mass1fract + velocities[atom2]*mass2fract;
Vec3 relVel = velocities[atom2] - velocities[atom1];
Vec3 comForce = forces[atom1] + forces[atom2];
Vec3 relForce = mass1fract*forces[atom2] - mass2fract*forces[atom1];
comVel += 0.5 * comForce * getDeltaT() * invTotMass;
relVel += 0.5 * relForce * getDeltaT() * invRedMass;
if (m1 != 0.0) {
velocities[atom1] = comVel - relVel*mass2fract + (atomCoordinates[atom1] - xPrime[atom1])/getDeltaT();
}
if (m2 != 0.0) {
velocities[atom2] = comVel + relVel*mass1fract + (atomCoordinates[atom2] - xPrime[atom2])/getDeltaT();
}
}
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
incrementTimeStep();
}
platforms/reference/tests/TestReferenceNoseHooverIntegrator.cpp 0 → 100644
View file @ f6431a42
/* -------------------------------------------------------------------------- *
* 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) 2015 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 "ReferenceTests.h"
#include "TestNoseHooverIntegrator.h"
void runPlatformTests() {
}
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