Skip to content

GitLab

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

Continuing to implement new CUDA platform: CustomBondForce, CustomAngleForce,... 

Continuing to implement new CUDA platform: CustomBondForce, CustomAngleForce, CustomTorsionForce, CustomExternalForce, CustomCompoundBondForce
parent 72d6b6a9
Show whitespace changes
Inline Side-by-side
Showing with 3224 additions and 1543 deletions
platforms/cuda2/src/CudaBondedUtilities.cpp
View file @ eb64fa2f
...@@ -58,6 +58,9 @@ std::string CudaBondedUtilities::addArgument(CUdeviceptr data, const string& typ ...@@ -58,6 +58,9 @@ std::string CudaBondedUtilities::addArgument(CUdeviceptr data, const string& typ
} }
void CudaBondedUtilities::addPrefixCode(const string& source) { void CudaBondedUtilities::addPrefixCode(const string& source) {
for (int i = 0; i < (int) prefixCode.size(); i++)
if (prefixCode[i] == source)
return;
prefixCode.push_back(source); prefixCode.push_back(source);
} }
...@@ -121,16 +124,13 @@ void CudaBondedUtilities::initialize(const System& system) { ...@@ -121,16 +124,13 @@ void CudaBondedUtilities::initialize(const System& system) {
string CudaBondedUtilities::createForceSource(int forceIndex, int numBonds, int numAtoms, int group, const string& computeForce) { string CudaBondedUtilities::createForceSource(int forceIndex, int numBonds, int numAtoms, int group, const string& computeForce) {
maxBonds = max(maxBonds, numBonds); maxBonds = max(maxBonds, numBonds);
string suffix1[] = {""}; string suffix[] = {".x", ".y", ".z", ".w"};
string suffix4[] = {".x", ".y", ".z", ".w"};
string* suffix;
stringstream s; stringstream s;
s<<"if ((groups&"<<(1<<group)<<") != 0)\n"; s<<"if ((groups&"<<(1<<group)<<") != 0)\n";
s<<"for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < "<<numBonds<<"; index += blockDim.x*gridDim.x) {\n"; s<<"for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < "<<numBonds<<"; index += blockDim.x*gridDim.x) {\n";
int startAtom = 0; int startAtom = 0;
for (int i = 0; i < (int) atomIndices[forceIndex].size(); i++) { for (int i = 0; i < (int) atomIndices[forceIndex].size(); i++) {
int indexWidth = atomIndices[forceIndex][i]->getElementSize()/4; int indexWidth = atomIndices[forceIndex][i]->getElementSize()/4;
suffix = (indexWidth == 1 ? suffix1 : suffix4);
string indexType = "uint"+context.intToString(indexWidth); string indexType = "uint"+context.intToString(indexWidth);
s<<" "<<indexType<<" atoms"<<i<<" = atomIndices"<<forceIndex<<"_"<<i<<"[index];\n"; s<<" "<<indexType<<" atoms"<<i<<" = atomIndices"<<forceIndex<<"_"<<i<<"[index];\n";
int atomsToLoad = min(indexWidth, numAtoms-startAtom); int atomsToLoad = min(indexWidth, numAtoms-startAtom);
......
platforms/cuda2/src/CudaKernelFactory.cpp
View file @ eb64fa2f
...@@ -78,20 +78,20 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform ...@@ -78,20 +78,20 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
return new CudaVirtualSitesKernel(name, platform, cu); return new CudaVirtualSitesKernel(name, platform, cu);
if (name == CalcHarmonicBondForceKernel::Name()) if (name == CalcHarmonicBondForceKernel::Name())
return new CudaCalcHarmonicBondForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcHarmonicBondForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomBondForceKernel::Name()) if (name == CalcCustomBondForceKernel::Name())
// return new CudaCalcCustomBondForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomBondForceKernel(name, platform, cu, context.getSystem());
if (name == CalcHarmonicAngleForceKernel::Name()) if (name == CalcHarmonicAngleForceKernel::Name())
return new CudaCalcHarmonicAngleForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcHarmonicAngleForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomAngleForceKernel::Name()) if (name == CalcCustomAngleForceKernel::Name())
// return new CudaCalcCustomAngleForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomAngleForceKernel(name, platform, cu, context.getSystem());
if (name == CalcPeriodicTorsionForceKernel::Name()) if (name == CalcPeriodicTorsionForceKernel::Name())
return new CudaCalcPeriodicTorsionForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcPeriodicTorsionForceKernel(name, platform, cu, context.getSystem());
if (name == CalcRBTorsionForceKernel::Name()) if (name == CalcRBTorsionForceKernel::Name())
return new CudaCalcRBTorsionForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcRBTorsionForceKernel(name, platform, cu, context.getSystem());
if (name == CalcCMAPTorsionForceKernel::Name()) if (name == CalcCMAPTorsionForceKernel::Name())
return new CudaCalcCMAPTorsionForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCMAPTorsionForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomTorsionForceKernel::Name()) if (name == CalcCustomTorsionForceKernel::Name())
// return new CudaCalcCustomTorsionForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomTorsionForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcNonbondedForceKernel::Name()) // if (name == CalcNonbondedForceKernel::Name())
// return new CudaCalcNonbondedForceKernel(name, platform, cu, context.getSystem()); // return new CudaCalcNonbondedForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomNonbondedForceKernel::Name()) // if (name == CalcCustomNonbondedForceKernel::Name())
...@@ -100,12 +100,12 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform ...@@ -100,12 +100,12 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
// return new CudaCalcGBSAOBCForceKernel(name, platform, cu); // return new CudaCalcGBSAOBCForceKernel(name, platform, cu);
// if (name == CalcCustomGBForceKernel::Name()) // if (name == CalcCustomGBForceKernel::Name())
// return new CudaCalcCustomGBForceKernel(name, platform, cu, context.getSystem()); // return new CudaCalcCustomGBForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomExternalForceKernel::Name()) if (name == CalcCustomExternalForceKernel::Name())
// return new CudaCalcCustomExternalForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomExternalForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomHbondForceKernel::Name()) // if (name == CalcCustomHbondForceKernel::Name())
// return new CudaCalcCustomHbondForceKernel(name, platform, cu, context.getSystem()); // return new CudaCalcCustomHbondForceKernel(name, platform, cu, context.getSystem());
// if (name == CalcCustomCompoundBondForceKernel::Name()) if (name == CalcCustomCompoundBondForceKernel::Name())
// return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem());
if (name == IntegrateVerletStepKernel::Name()) if (name == IntegrateVerletStepKernel::Name())
return new CudaIntegrateVerletStepKernel(name, platform, cu); return new CudaIntegrateVerletStepKernel(name, platform, cu);
// if (name == IntegrateLangevinStepKernel::Name()) // if (name == IntegrateLangevinStepKernel::Name())
......
platforms/cuda2/src/CudaKernels.cpp
View file @ eb64fa2f
This source diff could not be displayed because it is too large. You can view the blob instead.
platforms/cuda2/src/CudaKernels.h
View file @ eb64fa2f
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
#include "CudaPlatform.h" #include "CudaPlatform.h"
#include "CudaArray.h" #include "CudaArray.h"
#include "CudaContext.h" #include "CudaContext.h"
//#include "CudaParameterSet.h" #include "CudaParameterSet.h"
#include "CudaSort.h" #include "CudaSort.h"
#include "openmm/kernels.h" #include "openmm/kernels.h"
#include "openmm/System.h" #include "openmm/System.h"
...@@ -257,48 +257,48 @@ private: ...@@ -257,48 +257,48 @@ private:
CudaArray* params; CudaArray* params;
}; };
///** /**
// * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system.
// */ */
//class CudaCalcCustomBondForceKernel : public CalcCustomBondForceKernel { class CudaCalcCustomBondForceKernel : public CalcCustomBondForceKernel {
//public: public:
// CudaCalcCustomBondForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomBondForceKernel(name, platform), CudaCalcCustomBondForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomBondForceKernel(name, platform),
// hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) { hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) {
// } }
// ~CudaCalcCustomBondForceKernel(); ~CudaCalcCustomBondForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the CustomBondForce this kernel will be used for * @param force the CustomBondForce this kernel will be used for
// */ */
// void initialize(const System& system, const CustomBondForce& force); void initialize(const System& system, const CustomBondForce& force);
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Execute the kernel to calculate the forces and/or energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force * @return the potential energy due to the force
// */ */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// /** /**
// * Copy changed parameters over to a context. * Copy changed parameters over to a context.
// * *
// * @param context the context to copy parameters to * @param context the context to copy parameters to
// * @param force the CustomBondForce to copy the parameters from * @param force the CustomBondForce to copy the parameters from
// */ */
// void copyParametersToContext(ContextImpl& context, const CustomBondForce& force); void copyParametersToContext(ContextImpl& context, const CustomBondForce& force);
//private: private:
// int numBonds; int numBonds;
// bool hasInitializedKernel; bool hasInitializedKernel;
// CudaContext& cu; CudaContext& cu;
// System& system; System& system;
// CudaParameterSet* params; CudaParameterSet* params;
// CudaArray<cl_float>* globals; CudaArray* globals;
// std::vector<std::string> globalParamNames; std::vector<std::string> globalParamNames;
// std::vector<cl_float> globalParamValues; std::vector<float> globalParamValues;
//}; };
/** /**
* This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
...@@ -340,48 +340,48 @@ private: ...@@ -340,48 +340,48 @@ private:
CudaArray* params; CudaArray* params;
}; };
///** /**
// * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
// */ */
//class CudaCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel { class CudaCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
//public: public:
// CudaCalcCustomAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomAngleForceKernel(name, platform), CudaCalcCustomAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomAngleForceKernel(name, platform),
// hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) { hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) {
// } }
// ~CudaCalcCustomAngleForceKernel(); ~CudaCalcCustomAngleForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the CustomAngleForce this kernel will be used for * @param force the CustomAngleForce this kernel will be used for
// */ */
// void initialize(const System& system, const CustomAngleForce& force); void initialize(const System& system, const CustomAngleForce& force);
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Execute the kernel to calculate the forces and/or energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force * @return the potential energy due to the force
// */ */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// /** /**
// * Copy changed parameters over to a context. * Copy changed parameters over to a context.
// * *
// * @param context the context to copy parameters to * @param context the context to copy parameters to
// * @param force the CustomAngleForce to copy the parameters from * @param force the CustomAngleForce to copy the parameters from
// */ */
// void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force); void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force);
//private: private:
// int numAngles; int numAngles;
// bool hasInitializedKernel; bool hasInitializedKernel;
// CudaContext& cu; CudaContext& cu;
// System& system; System& system;
// CudaParameterSet* params; CudaParameterSet* params;
// CudaArray<cl_float>* globals; CudaArray* globals;
// std::vector<std::string> globalParamNames; std::vector<std::string> globalParamNames;
// std::vector<cl_float> globalParamValues; std::vector<float> globalParamValues;
//}; };
/** /**
* This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
...@@ -499,49 +499,49 @@ private: ...@@ -499,49 +499,49 @@ private:
CudaArray* torsionMaps; CudaArray* torsionMaps;
}; };
///** /**
// * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
// */ */
//class CudaCalcCustomTorsionForceKernel : public CalcCustomTorsionForceKernel { class CudaCalcCustomTorsionForceKernel : public CalcCustomTorsionForceKernel {
//public: public:
// CudaCalcCustomTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomTorsionForceKernel(name, platform), CudaCalcCustomTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomTorsionForceKernel(name, platform),
// hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) { hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) {
// } }
// ~CudaCalcCustomTorsionForceKernel(); ~CudaCalcCustomTorsionForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the CustomTorsionForce this kernel will be used for * @param force the CustomTorsionForce this kernel will be used for
// */ */
// void initialize(const System& system, const CustomTorsionForce& force); void initialize(const System& system, const CustomTorsionForce& force);
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Execute the kernel to calculate the forces and/or energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force * @return the potential energy due to the force
// */ */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// /** /**
// * Copy changed parameters over to a context. * Copy changed parameters over to a context.
// * *
// * @param context the context to copy parameters to * @param context the context to copy parameters to
// * @param force the CustomTorsionForce to copy the parameters from * @param force the CustomTorsionForce to copy the parameters from
// */ */
// void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force); void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force);
//private: private:
// int numTorsions; int numTorsions;
// bool hasInitializedKernel; bool hasInitializedKernel;
// CudaContext& cu; CudaContext& cu;
// System& system; System& system;
// CudaParameterSet* params; CudaParameterSet* params;
// CudaArray<cl_float>* globals; CudaArray* globals;
// std::vector<std::string> globalParamNames; std::vector<std::string> globalParamNames;
// std::vector<cl_float> globalParamValues; std::vector<float> globalParamValues;
//}; };
//
///** ///**
// * This kernel is invoked by NonbondedForce to calculate the forces acting on the system. // * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
// */ // */
...@@ -616,6 +616,7 @@ private: ...@@ -616,6 +616,7 @@ private:
// CUfunction pmeConvolutionKernel; // CUfunction pmeConvolutionKernel;
// CUfunction pmeInterpolateForceKernel; // CUfunction pmeInterpolateForceKernel;
// std::map<std::string, std::string> pmeDefines; // std::map<std::string, std::string> pmeDefines;
// std::vector<std::pair<int, int> > exceptionAtoms;
// double ewaldSelfEnergy, dispersionCoefficient, alpha; // double ewaldSelfEnergy, dispersionCoefficient, alpha;
// int interpolateForceThreads; // int interpolateForceThreads;
// bool hasCoulomb, hasLJ; // bool hasCoulomb, hasLJ;
...@@ -769,50 +770,50 @@ private: ...@@ -769,50 +770,50 @@ private:
// System& system; // System& system;
// CUfunction pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, gradientChainRuleKernel; // CUfunction pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, gradientChainRuleKernel;
//}; //};
//
///** /**
// * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
// */ */
//class CudaCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel { class CudaCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
//public: public:
// CudaCalcCustomExternalForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomExternalForceKernel(name, platform), CudaCalcCustomExternalForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomExternalForceKernel(name, platform),
// hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) { hasInitializedKernel(false), cu(cu), system(system), params(NULL), globals(NULL) {
// } }
// ~CudaCalcCustomExternalForceKernel(); ~CudaCalcCustomExternalForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the CustomExternalForce this kernel will be used for * @param force the CustomExternalForce this kernel will be used for
// */ */
// void initialize(const System& system, const CustomExternalForce& force); void initialize(const System& system, const CustomExternalForce& force);
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Execute the kernel to calculate the forces and/or energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force * @return the potential energy due to the force
// */ */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// /** /**
// * Copy changed parameters over to a context. * Copy changed parameters over to a context.
// * *
// * @param context the context to copy parameters to * @param context the context to copy parameters to
// * @param force the CustomExternalForce to copy the parameters from * @param force the CustomExternalForce to copy the parameters from
// */ */
// void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force); void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
//private: private:
// int numParticles; int numParticles;
// bool hasInitializedKernel; bool hasInitializedKernel;
// CudaContext& cu; CudaContext& cu;
// System& system; System& system;
// CudaParameterSet* params; CudaParameterSet* params;
// CudaArray<cl_float>* globals; CudaArray* globals;
// std::vector<std::string> globalParamNames; std::vector<std::string> globalParamNames;
// std::vector<cl_float> globalParamValues; std::vector<float> globalParamValues;
//}; };
//
///** ///**
// * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system. // * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
// */ // */
...@@ -867,51 +868,51 @@ private: ...@@ -867,51 +868,51 @@ private:
// System& system; // System& system;
// CUfunction donorKernel, acceptorKernel; // CUfunction donorKernel, acceptorKernel;
//}; //};
//
///** /**
// * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system. * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
// */ */
//class CudaCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel { class CudaCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
//public: public:
// CudaCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomCompoundBondForceKernel(name, platform), CudaCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCustomCompoundBondForceKernel(name, platform),
// cu(cu), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) { cu(cu), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
// } }
// ~CudaCalcCustomCompoundBondForceKernel(); ~CudaCalcCustomCompoundBondForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the CustomCompoundBondForce this kernel will be used for * @param force the CustomCompoundBondForce this kernel will be used for
// */ */
// void initialize(const System& system, const CustomCompoundBondForce& force); void initialize(const System& system, const CustomCompoundBondForce& force);
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Execute the kernel to calculate the forces and/or energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force * @return the potential energy due to the force
// */ */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// /** /**
// * Copy changed parameters over to a context. * Copy changed parameters over to a context.
// * *
// * @param context the context to copy parameters to * @param context the context to copy parameters to
// * @param force the CustomCompoundBondForce to copy the parameters from * @param force the CustomCompoundBondForce to copy the parameters from
// */ */
// void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force); void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
//
//private: private:
// int numBonds; int numBonds;
// CudaContext& cu; CudaContext& cu;
// CudaParameterSet* params; CudaParameterSet* params;
// CudaArray<cl_float>* globals; CudaArray* globals;
// CudaArray<mm_float4>* tabulatedFunctionParams; CudaArray* tabulatedFunctionParams;
// std::vector<std::string> globalParamNames; std::vector<std::string> globalParamNames;
// std::vector<cl_float> globalParamValues; std::vector<float> globalParamValues;
// std::vector<CudaArray<mm_float4>*> tabulatedFunctions; std::vector<CudaArray*> tabulatedFunctions;
// System& system; System& system;
//}; };
/** /**
* This kernel is invoked by VerletIntegrator to take one time step. * This kernel is invoked by VerletIntegrator to take one time step.
......
platforms/cuda2/src/CudaNonbondedUtilities.h 0 → 100644
View file @ eb64fa2f
#ifndef OPENMM_CUDANONBONDEDUTILITIES_H_
#define OPENMM_CUDANONBONDEDUTILITIES_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "CudaContext.h"
#include "openmm/System.h"
#include "CudaExpressionUtilities.h"
#include <sstream>
#include <string>
#include <vector>
namespace OpenMM {
/**
* This class provides a generic interface for calculating nonbonded interactions. It does this in two
* ways. First, it can be used to create Kernels that evaluate nonbonded interactions. Clients
* only need to provide the code for evaluating a single interaction and the list of parameters it depends on.
* A complete kernel is then synthesized using an appropriate algorithm to evaluate all interactions on all
* atoms.
*
* Second, this class itself creates and invokes a single "default" interaction kernel, allowing several
* different forces to be evaluated at once for greater efficiency. Call addInteraction() and addParameter()
* to add interactions to this default kernel.
*
* During each force or energy evaluation, the following sequence of steps takes place:
*
* 1. Data structures (e.g. neighbor lists) are calculated to allow nonbonded interactions to be evaluated
* quickly.
*
* 2. calcForcesAndEnergy() is called on each ForceImpl in the System.
*
* 3. Finally, the default interaction kernel is invoked to calculate all interactions that were added
* to it.
*
* This sequence means that the default interaction kernel may depend on quantities that were calculated
* by ForceImpls during calcForcesAndEnergy().
*/
class OPENMM_EXPORT CudaNonbondedUtilities {
public:
class ParameterInfo;
// CudaNonbondedUtilities(CudaContext& context);
// ~CudaNonbondedUtilities();
// /**
// * Add a nonbonded interaction to be evaluated by the default interaction kernel.
// *
// * @param usesCutoff specifies whether a cutoff should be applied to this interaction
// * @param usesPeriodic specifies whether periodic boundary conditions should be applied to this interaction
// * @param usesExclusions specifies whether this interaction uses exclusions. If this is true, it must have identical exclusions to every other interaction.
// * @param cutoffDistance the cutoff distance for this interaction (ignored if usesCutoff is false)
// * @param exclusionList for each atom, specifies the list of other atoms whose interactions should be excluded
// * @param kernel the code to evaluate the interaction
// * @param forceGroup the force group in which the interaction should be calculated
// */
// void addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList, const std::string& kernel, int forceGroup);
// /**
// * Add a per-atom parameter that the default interaction kernel may depend on.
// */
// void addParameter(const ParameterInfo& parameter);
// /**
// * Add an array (other than a per-atom parameter) that should be passed as an argument to the default interaction kernel.
// */
// void addArgument(const ParameterInfo& parameter);
// /**
// * Specify the list of exclusions that an interaction outside the default kernel will depend on.
// *
// * @param exclusionList for each atom, specifies the list of other atoms whose interactions should be excluded
// */
// void requestExclusions(const std::vector<std::vector<int> >& exclusionList);
// /**
// * Initialize this object in preparation for a simulation.
// */
// void initialize(const System& system);
// /**
// * Get the number of force buffers required for nonbonded forces.
// */
// int getNumForceBuffers() {
// return numForceBuffers;
// }
// /**
// * Get the number of energy buffers required for nonbonded forces.
// */
// int getNumEnergyBuffers() {
// return numForceThreadBlocks*forceThreadBlockSize;
// }
// /**
// * Get whether a cutoff is being used.
// */
// bool getUseCutoff() {
// return useCutoff;
// }
// /**
// * Get whether periodic boundary conditions are being used.
// */
// bool getUsePeriodic() {
// return usePeriodic;
// }
// /**
// * Get whether there is one force buffer per atom block.
// */
// bool getForceBufferPerAtomBlock() {
// return forceBufferPerAtomBlock;
// }
// /**
// * Get the number of work groups used for computing nonbonded forces.
// */
// int getNumForceThreadBlocks() {
// return numForceThreadBlocks;
// }
// /**
// * Get the size of each work group used for computing nonbonded forces.
// */
// int getForceThreadBlockSize() {
// return forceThreadBlockSize;
// }
// /**
// * Get the cutoff distance.
// */
// double getCutoffDistance() {
// return cutoff;
// }
// /**
// * Get whether any interactions have been added.
// */
// bool getHasInteractions() {
// return cutoff != -1.0;
// }
// /**
// * Get the force group in which nonbonded interactions should be computed.
// */
// int getForceGroup() {
// return nonbondedForceGroup;
// }
// /**
// * Prepare to compute interactions. This updates the neighbor list.
// */
// void prepareInteractions();
// /**
// * Compute the nonbonded interactions.
// */
// void computeInteractions();
// /**
// * Check to see if the neighbor list arrays are large enough, and make them bigger if necessary.
// */
// void updateNeighborListSize();
// /**
// * Get the array containing the center of each atom block.
// */
// CudaArray<mm_float4>& getBlockCenters() {
// return *blockCenter;
// }
// /**
// * Get the array containing the dimensions of each atom block.
// */
// CudaArray<mm_float4>& getBlockBoundingBoxes() {
// return *blockBoundingBox;
// }
// /**
// * Get the array whose first element contains the number of tiles with interactions.
// */
// CudaArray<cl_uint>& getInteractionCount() {
// return *interactionCount;
// }
// /**
// * Get the array containing tiles with interactions.
// */
// CudaArray<mm_ushort2>& getInteractingTiles() {
// return *interactingTiles;
// }
// /**
// * Get the array containing flags for tiles with interactions.
// */
// CudaArray<cl_uint>& getInteractionFlags() {
// return *interactionFlags;
// }
// /**
// * Get the array containing exclusion flags.
// */
// CudaArray<cl_uint>& getExclusions() {
// return *exclusions;
// }
// /**
// * Get the array containing the index into the exclusion array for each tile.
// */
// CudaArray<cl_uint>& getExclusionIndices() {
// return *exclusionIndices;
// }
// /**
// * Get the array listing where the exclusion data starts for each row.
// */
// CudaArray<cl_uint>& getExclusionRowIndices() {
// return *exclusionRowIndices;
// }
// /**
// * Get the index of the first tile this context is responsible for processing.
// */
// int getStartTileIndex() const {
// return startTileIndex;
// }
// /**
// * Get the total number of tiles this context is responsible for processing.
// */
// int getNumTiles() const {
// return numTiles;
// }
// /**
// * Set the range of tiles that should be processed by this context.
// */
// void setTileRange(int startTileIndex, int numTiles);
// /**
// * Create a Kernel for evaluating a nonbonded interaction. Cutoffs and periodic boundary conditions
// * are assumed to be the same as those for the default interaction Kernel, since this kernel will use
// * the same neighbor list.
// *
// * @param source the source code for evaluating the force and energy
// * @param params the per-atom parameters this kernel may depend on
// * @param arguments arrays (other than per-atom parameters) that should be passed as arguments to the kernel
// * @param useExclusions specifies whether exclusions are applied to this interaction
// * @param isSymmetric specifies whether the interaction is symmetric
// */
// cl::Kernel createInteractionKernel(const std::string& source, const std::vector<ParameterInfo>& params, const std::vector<ParameterInfo>& arguments, bool useExclusions, bool isSymmetric) const;
private:
// static int findExclusionIndex(int x, int y, const std::vector<cl_uint>& exclusionIndices, const std::vector<cl_uint>& exclusionRowIndices);
// CudaContext& context;
// cl::Kernel forceKernel;
// cl::Kernel findBlockBoundsKernel;
// cl::Kernel findInteractingBlocksKernel;
// cl::Kernel findInteractionsWithinBlocksKernel;
// CudaArray<cl_uint>* exclusions;
// CudaArray<cl_uint>* exclusionIndices;
// CudaArray<cl_uint>* exclusionRowIndices;
// CudaArray<mm_ushort2>* interactingTiles;
// CudaArray<cl_uint>* interactionFlags;
// CudaArray<cl_uint>* interactionCount;
// CudaArray<mm_float4>* blockCenter;
// CudaArray<mm_float4>* blockBoundingBox;
// std::vector<std::vector<int> > atomExclusions;
// std::vector<ParameterInfo> parameters;
// std::vector<ParameterInfo> arguments;
// std::string kernelSource;
// std::map<std::string, std::string> kernelDefines;
// double cutoff;
// bool useCutoff, usePeriodic, forceBufferPerAtomBlock, deviceIsCpu, anyExclusions;
// int numForceBuffers, startTileIndex, numTiles, numForceThreadBlocks, forceThreadBlockSize, nonbondedForceGroup;
};
/**
* This class stores information about a per-atom parameter that may be used in a nonbonded kernel.
*/
class CudaNonbondedUtilities::ParameterInfo {
public:
/**
* Create a ParameterInfo object.
*
* @param name the name of the parameter
* @param type the data type of the parameter's components
* @param numComponents the number of components in the parameter
* @param size the size of the parameter in bytes
* @param memory the memory containing the parameter values
*/
ParameterInfo(const std::string& name, const std::string& componentType, int numComponents, int size, CUdeviceptr memory) :
name(name), componentType(componentType), numComponents(numComponents), size(size), memory(memory) {
if (numComponents == 1)
type = componentType;
else {
std::stringstream s;
s << componentType << numComponents;
type = s.str();
}
}
const std::string& getName() const {
return name;
}
const std::string& getComponentType() const {
return componentType;
}
const std::string& getType() const {
return type;
}
int getNumComponents() const {
return numComponents;
}
int getSize() const {
return size;
}
CUdeviceptr getMemory() const {
return memory;
}
private:
std::string name;
std::string componentType;
std::string type;
int size, numComponents;
CUdeviceptr memory;
};
} // namespace OpenMM
#endif /*OPENMM_CUDANONBONDEDUTILITIES_H_*/
platforms/cuda2/src/CudaParameterSet.cpp 0 → 100644
View file @ eb64fa2f
/* -------------------------------------------------------------------------- *
* 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) 2009-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "CudaParameterSet.h"
#include "openmm/OpenMMException.h"
#include <cmath>
#include <sstream>
using namespace OpenMM;
using namespace std;
#define CHECK_RESULT(result) \
if (result != CUDA_SUCCESS) { \
std::stringstream m; \
m<<errorMessage<<": "<<context.getErrorString(result)<<" ("<<result<<")"; \
throw OpenMMException(m.str());\
}
CudaParameterSet::CudaParameterSet(CudaContext& context, int numParameters, int numObjects, const string& name, bool bufferPerParameter) :
context(context), numParameters(numParameters), numObjects(numObjects), name(name) {
int params = numParameters;
int bufferCount = 0;
int elementSize = 4;
CUdeviceptr pointer;
string errorMessage = "Error creating parameter set "+name;
if (!bufferPerParameter) {
while (params > 2) {
CHECK_RESULT(cuMemAlloc(&pointer, numObjects*elementSize*4));
std::stringstream name;
name << "param" << (++bufferCount);
buffers.push_back(CudaNonbondedUtilities::ParameterInfo(name.str(), "float", 4, elementSize*4, pointer));
params -= 4;
}
if (params > 1) {
CHECK_RESULT(cuMemAlloc(&pointer, numObjects*elementSize*2));
std::stringstream name;
name << "param" << (++bufferCount);
buffers.push_back(CudaNonbondedUtilities::ParameterInfo(name.str(), "float", 2, elementSize*2, pointer));
params -= 2;
}
}
while (params > 0) {
CHECK_RESULT(cuMemAlloc(&pointer, numObjects*elementSize));
std::stringstream name;
name << "param" << (++bufferCount);
buffers.push_back(CudaNonbondedUtilities::ParameterInfo(name.str(), "float", 1, elementSize, pointer));
params--;
}
}
CudaParameterSet::~CudaParameterSet() {
string errorMessage = "Error freeing device memory";
for (int i = 0; i < (int) buffers.size(); i++)
CHECK_RESULT(cuMemFree(buffers[i].getMemory()));
}
void CudaParameterSet::getParameterValues(vector<vector<float> >& values) const {
values.resize(numObjects);
for (int i = 0; i < numObjects; i++)
values[i].resize(numParameters);
int base = 0;
string errorMessage = "Error downloading parameter set "+name;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getType() == "float4") {
vector<float4> data(numObjects);
CHECK_RESULT(cuMemcpyDtoH(&data[0], buffers[i].getMemory(), numObjects*buffers[i].getSize()));
for (int j = 0; j < numObjects; j++) {
values[j][base] = data[j].x;
if (base+1 < numParameters)
values[j][base+1] = data[j].y;
if (base+2 < numParameters)
values[j][base+2] = data[j].z;
if (base+3 < numParameters)
values[j][base+3] = data[j].w;
}
base += 4;
}
else if (buffers[i].getType() == "float2") {
vector<float2> data(numObjects);
CHECK_RESULT(cuMemcpyDtoH(&data[0], buffers[i].getMemory(), numObjects*buffers[i].getSize()));
for (int j = 0; j < numObjects; j++) {
values[j][base] = data[j].x;
if (base+1 < numParameters)
values[j][base+1] = data[j].y;
}
base += 2;
}
else if (buffers[i].getType() == "float") {
vector<float> data(numObjects);
CHECK_RESULT(cuMemcpyDtoH(&data[0], buffers[i].getMemory(), numObjects*buffers[i].getSize()));
for (int j = 0; j < numObjects; j++)
values[j][base] = data[j];
base++;
}
else
throw OpenMMException("Internal error: Unknown buffer type in CudaParameterSet");
}
}
void CudaParameterSet::setParameterValues(const vector<vector<float> >& values) {
int base = 0;
string errorMessage = "Error uploading parameter set "+name;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getType() == "float4") {
vector<float4> data(numObjects);
for (int j = 0; j < numObjects; j++) {
data[j].x = values[j][base];
if (base+1 < numParameters)
data[j].y = values[j][base+1];
if (base+2 < numParameters)
data[j].z = values[j][base+2];
if (base+3 < numParameters)
data[j].w = values[j][base+3];
}
CHECK_RESULT(cuMemcpyHtoD(buffers[i].getMemory(), &data[0], numObjects*buffers[i].getSize()));
base += 4;
}
else if (buffers[i].getType() == "float2") {
vector<float2> data(numObjects);
for (int j = 0; j < numObjects; j++) {
data[j].x = values[j][base];
if (base+1 < numParameters)
data[j].y = values[j][base+1];
}
CHECK_RESULT(cuMemcpyHtoD(buffers[i].getMemory(), &data[0], numObjects*buffers[i].getSize()));
base += 2;
}
else if (buffers[i].getType() == "float") {
vector<float> data(numObjects);
for (int j = 0; j < numObjects; j++)
data[j] = values[j][base];
CHECK_RESULT(cuMemcpyHtoD(buffers[i].getMemory(), &data[0], numObjects*buffers[i].getSize()));
base++;
}
else
throw OpenMMException("Internal error: Unknown buffer type in CudaParameterSet");
}
}
string CudaParameterSet::getParameterSuffix(int index, const std::string& extraSuffix) const {
const string suffixes[] = {".x", ".y", ".z", ".w"};
int buffer = -1;
for (int i = 0; buffer == -1 && i < (int) buffers.size(); i++) {
if (index*sizeof(float) < buffers[i].getSize())
buffer = i;
else
index -= buffers[i].getSize()/sizeof(float);
}
if (buffer == -1)
throw OpenMMException("Internal error: Illegal argument to CudaParameterSet::getParameterSuffix() ("+name+")");
stringstream suffix;
suffix << (buffer+1) << extraSuffix;
if (buffers[buffer].getType() != "float")
suffix << suffixes[index];
return suffix.str();
}
platforms/cuda2/src/CudaParameterSet.h 0 → 100644
View file @ eb64fa2f
#ifndef OPENMM_CUDAPARAMETERSET_H_
#define OPENMM_CUDAPARAMETERSET_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "CudaContext.h"
#include "CudaNonbondedUtilities.h"
namespace OpenMM {
class CudaNonbondedUtilities;
/**
* This class represents a set of floating point parameter values for a set of objects (particles, bonds, etc.).
* It automatically creates an appropriate set of device buffers to hold the parameter values, based
* on the number of parameters required.
*/
class OPENMM_EXPORT CudaParameterSet {
public:
/**
* Create an CudaParameterSet.
*
* @param context the context for which to create the parameter set
* @param numParameters the number of parameters for each object
* @param numObjects the number of objects to store parameter values for
* @param name the name of the parameter set
* @param bufferPerParameter if true, a separate buffer is created for each parameter. If false,
* multiple parameters may be combined into a single buffer.
*/
CudaParameterSet(CudaContext& context, int numParameters, int numObjects, const std::string& name, bool bufferPerParameter=false);
~CudaParameterSet();
/**
* Get the number of parameters.
*/
int getNumParameters() const {
return numParameters;
}
/**
* Get the number of objects.
*/
int getNumObjects() const {
return numObjects;
}
/**
* Get the values of all parameters.
*
* @param values on exit, values[i][j] contains the value of parameter j for object i
*/
void getParameterValues(std::vector<std::vector<float> >& values) const;
/**
* Set the values of all parameters.
*
* @param values values[i][j] contains the value of parameter j for object i
*/
void setParameterValues(const std::vector<std::vector<float> >& values);
/**
* Get a set of CudaNonbondedUtilities::ParameterInfo objects which describe the Buffers
* containing the data.
*/
const std::vector<CudaNonbondedUtilities::ParameterInfo>& getBuffers() const {
return buffers;
}
/**
* Get a suffix to add to variable names when accessing a certain parameter.
*
* @param index the index of the parameter
* @param extraSuffix an extra suffix to add to the variable name
* @return the suffix to append
*/
std::string getParameterSuffix(int index, const std::string& extraSuffix = "") const;
private:
CudaContext& context;
int numParameters;
int numObjects;
std::string name;
std::vector<CudaNonbondedUtilities::ParameterInfo> buffers;
};
} // namespace OpenMM
#endif /*OPENMM_CUDAPARAMETERSET_H_*/
platforms/cuda2/src/kernels/customCompoundBond.cu 0 → 100644
View file @ eb64fa2f
/**
* Convert a real4 to a real3 by removing its last element.
*/
__device__ real3 ccb_trim(real4 v) {
return make_real3(v.x, v.y, v.z);
}
/**
* Compute the difference between two vectors, setting the fourth component to the squared magnitude.
*/
__device__ real4 ccb_delta(real4 vec1, real4 vec2) {
real4 result = make_real4(vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0);
result.w = result.x*result.x + result.y*result.y + result.z*result.z;
return result;
}
/**
* Compute the angle between two vectors. The w component of each vector should contain the squared magnitude.
*/
__device__ real ccb_computeAngle(real4 vec1, real4 vec2) {
real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
real angle;
if (cosine > 0.99f || cosine < -0.99f) {
// We're close to the singularity in acos(), so take the cross product and use asin() instead.
real3 crossProduct = cross(vec1, vec2);
real scale = vec1.w*vec2.w;
angle = ASIN(SQRT(dot(crossProduct, crossProduct)/scale));
if (cosine < 0.0f)
angle = M_PI-angle;
}
else
angle = ACOS(cosine);
return angle;
}
/**
* Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
*/
__device__ real4 ccb_computeCross(real4 vec1, real4 vec2) {
real3 cp = cross(vec1, vec2);
return make_real4(cp.x, cp.y, cp.z, cp.x*cp.x+cp.y*cp.y+cp.z*cp.z);
}
platforms/cuda2/src/kernels/customExternalForce.cu 0 → 100644
View file @ eb64fa2f
COMPUTE_FORCE
real3 force1 = make_real3(-dEdX, -dEdY, -dEdZ);
platforms/cuda2/src/kernels/vectorOps.cu
View file @ eb64fa2f
...@@ -456,6 +456,50 @@ inline __device__ double4 operator*(double a, double4 b) { ...@@ -456,6 +456,50 @@ inline __device__ double4 operator*(double a, double4 b) {
return make_double4(a*b.x, a*b.y, a*b.z, a*b.w); return make_double4(a*b.x, a*b.y, a*b.z, a*b.w);
} }
// Divide a vector by a constant.
inline __device__ int2 operator/(int2 a, int b) {
return make_int2(a.x/b, a.y/b);
}
inline __device__ int3 operator/(int3 a, int b) {
return make_int3(a.x/b, a.y/b, a.z/b);
}
inline __device__ int4 operator/(int4 a, int b) {
return make_int4(a.x/b, a.y/b, a.z/b, a.w/b);
}
inline __device__ float2 operator/(float2 a, float b) {
float scale = 1.0f/b;
return a*scale;
}
inline __device__ float3 operator/(float3 a, float b) {
float scale = 1.0f/b;
return a*scale;
}
inline __device__ float4 operator/(float4 a, float b) {
float scale = 1.0f/b;
return a*scale;
}
inline __device__ double2 operator/(double2 a, double b) {
double scale = 1.0/b;
return a*scale;
}
inline __device__ double3 operator/(double3 a, double b) {
double scale = 1.0/b;
return a*scale;
}
inline __device__ double4 operator/(double4 a, double b) {
double scale = 1.0/b;
return a*scale;
}
// *= operator (multiply vector by constant) // *= operator (multiply vector by constant)
inline __device__ void operator*=(int2& a, int b) { inline __device__ void operator*=(int2& a, int b) {
......
platforms/cuda2/tests/TestCudaCustomAngleForce.cpp 0 → 100644
View file @ eb64fa2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of CustomAngleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomAngleForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testAngles() {
CudaPlatform platform;
// Create a system using a CustomAngleForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomAngleForce* custom = new CustomAngleForce("scale*k*(theta-theta0)^2");
custom->addPerAngleParameter("theta0");
custom->addPerAngleParameter("k");
custom->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 0.8;
custom->addAngle(0, 1, 2, parameters);
parameters[0] = 2.0;
parameters[1] = 0.5;
custom->addAngle(1, 2, 3, parameters);
customSystem.addForce(custom);
// Create an identical system using a HarmonicAngleForce.
System harmonicSystem;
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
HarmonicAngleForce* harmonic = new HarmonicAngleForce();
harmonic->addAngle(0, 1, 2, 1.5, 0.8);
harmonic->addAngle(1, 2, 3, 2.0, 0.5);
harmonicSystem.addForce(harmonic);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(4);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(harmonicSystem, integrator2, platform);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the angle parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
custom->setAngleParameters(0, 0, 1, 2, parameters);
parameters[0] = 2.1;
parameters[1] = 0.6;
custom->setAngleParameters(1, 1, 2, 3, parameters);
custom->updateParametersInContext(c1);
harmonic->setAngleParameters(0, 0, 1, 2, 1.6, 0.9);
harmonic->setAngleParameters(1, 1, 2, 3, 2.1, 0.6);
harmonic->updateParametersInContext(c2);
{
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void testParallelComputation() {
CudaPlatform platform;
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomAngleForce* force = new CustomAngleForce("(theta-1.1)^2");
vector<double> params;
for (int i = 2; i < numParticles; i++)
force->addAngle(i-2, i-1, i, params);
system.addForce(force);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(i, i%2, 0);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
map<string, string> props;
props[CudaPlatform::CudaDeviceIndex()] = deviceIndex+","+deviceIndex;
Context context2(system, integrator2, platform, props);
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
int main() {
try {
testAngles();
// testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda2/tests/TestCudaCustomBondForce.cpp 0 → 100644
View file @ eb64fa2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of CustomBondForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBonds() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomBondForce* forceField = new CustomBondForce("scale*k*(r-r0)^2");
forceField->addPerBondParameter("r0");
forceField->addPerBondParameter("k");
forceField->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 0.8;
forceField->addBond(0, 1, parameters);
parameters[0] = 1.2;
parameters[1] = 0.7;
forceField->addBond(1, 2, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
// Try changing the bond parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
forceField->setBondParameters(0, 0, 1, parameters);
parameters[0] = 1.3;
parameters[1] = 0.8;
forceField->setBondParameters(1, 1, 2, parameters);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.9*0.4, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.8*0.3, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.9*0.4*0.4 + 0.5*0.8*0.3*0.3, state.getPotentialEnergy(), TOL);
}
}
void testManyParameters() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomBondForce* forceField = new CustomBondForce("(a+b+c+d+e+f+g+h+i)*r");
forceField->addPerBondParameter("a");
forceField->addPerBondParameter("b");
forceField->addPerBondParameter("c");
forceField->addPerBondParameter("d");
forceField->addPerBondParameter("e");
forceField->addPerBondParameter("f");
forceField->addPerBondParameter("g");
forceField->addPerBondParameter("h");
forceField->addPerBondParameter("i");
vector<double> parameters(forceField->getNumPerBondParameters());
for (int i = 0; i < parameters.size(); i++)
parameters[i] = i;
forceField->addBond(0, 1, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 2.5, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double f = 1+2+3+4+5+6+7+8;
ASSERT_EQUAL_VEC(Vec3(0, f, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -f, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(f*2.5, state.getPotentialEnergy(), TOL);
}
void testParallelComputation() {
CudaPlatform platform;
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomBondForce* force = new CustomBondForce(("(r-1.1)^2"));
vector<double> params;
for (int i = 1; i < numParticles; i++)
force->addBond(i-1, i, params);
system.addForce(force);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(i, 0, 0);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
map<string, string> props;
props[CudaPlatform::CudaDeviceIndex()] = deviceIndex+","+deviceIndex;
Context context2(system, integrator2, platform, props);
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
int main() {
try {
testBonds();
testManyParameters();
// testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda2/tests/TestCudaCustomCompoundBondForce.cpp 0 → 100644
View file @ eb64fa2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of CustomCompoundBondForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBond() {
CudaPlatform platform;
// Create a system using a CustomCompoundBondForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))");
custom->addPerBondParameter("kb");
custom->addPerBondParameter("ka");
custom->addPerBondParameter("kt");
custom->addPerBondParameter("b0");
custom->addPerBondParameter("a0");
custom->addPerBondParameter("t0");
vector<int> particles(4);
particles[0] = 0;
particles[1] = 1;
particles[2] = 3;
particles[3] = 2;
vector<double> parameters(6);
parameters[0] = 1.5;
parameters[1] = 0.8;
parameters[2] = 0.6;
parameters[3] = 1.1;
parameters[4] = 2.9;
parameters[5] = 1.3;
custom->addBond(particles, parameters);
customSystem.addForce(custom);
// Create an identical system using standard forces.
System standardSystem;
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
HarmonicBondForce* bonds = new HarmonicBondForce();
bonds->addBond(0, 1, 1.1, 1.5);
bonds->addBond(1, 3, 1.1, 1.5);
standardSystem.addForce(bonds);
HarmonicAngleForce* angles = new HarmonicAngleForce();
angles->addAngle(1, 3, 2, 2.9, 0.8);
standardSystem.addForce(angles);
PeriodicTorsionForce* torsions = new PeriodicTorsionForce();
torsions->addTorsion(0, 1, 3, 2, 1, 1.3, 0.6);
standardSystem.addForce(torsions);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(standardSystem, integrator2, platform);
vector<Vec3> positions(4);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the bond parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[3] = 1.3;
custom->setBondParameters(0, particles, parameters);
custom->updateParametersInContext(c1);
bonds->setBondParameters(0, 0, 1, 1.3, 1.6);
bonds->setBondParameters(1, 1, 3, 1.3, 1.6);
bonds->updateParametersInContext(c2);
{
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces();
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void testPositionDependence() {
CudaPlatform platform;
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(2, "scale1*distance(p1,p2)+scale2*x1+2*y2");
custom->addGlobalParameter("scale1", 0.3);
custom->addGlobalParameter("scale2", 0.2);
vector<int> particles(2);
particles[0] = 0;
particles[1] = 1;
vector<double> parameters;
custom->addBond(particles, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(2);
positions[0] = Vec3(0.5, 1, 0);
positions[1] = Vec3(1.5, 1, 0);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.3*1.0+0.2*0.5+2*1, state.getPotentialEnergy(), 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3-0.2, 0, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(-0.3, -2, 0), state.getForces()[1], 1e-5);
}
void testParallelComputation() {
CudaPlatform platform;
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomCompoundBondForce* force = new CustomCompoundBondForce(2, ("(distance(p1,p2)-1.1)^2"));
vector<int> particles(2);
vector<double> params;
for (int i = 1; i < numParticles; i++) {
particles[0] = i-1;
particles[1] = i;
force->addBond(particles, params);
}
system.addForce(force);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(i, 0, 0);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
map<string, string> props;
props[CudaPlatform::CudaDeviceIndex()] = deviceIndex+","+deviceIndex;
Context context2(system, integrator2, platform, props);
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
int main() {
try {
testBond();
testPositionDependence();
// testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda2/tests/TestCudaCustomExternalForce.cpp 0 → 100644
View file @ eb64fa2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of CustomExternalForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testForce() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomExternalForce* forceField = new CustomExternalForce("scale*(x+yscale*(y-y0)^2)");
forceField->addPerParticleParameter("y0");
forceField->addPerParticleParameter("yscale");
forceField->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 0.5;
parameters[1] = 2.0;
forceField->addParticle(0, parameters);
parameters[0] = 1.5;
parameters[1] = 3.0;
forceField->addParticle(2, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 1);
positions[2] = Vec3(1, 0, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.0*2.0*1.5, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.0*1.5*1.5), state.getPotentialEnergy(), TOL);
}
// Try changing the parameters and make sure it's still correct.
parameters[0] = 1.4;
parameters[1] = 3.5;
forceField->setParticleParameters(1, 2, parameters);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.5*2.0*1.4, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.5*1.4*1.4), state.getPotentialEnergy(), TOL);
}
}
void testManyParameters() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomExternalForce* forceField = new CustomExternalForce("xscale*(x-x0)^2+yscale*(y-y0)^2+zscale*(z-z0)^2");
forceField->addPerParticleParameter("x0");
forceField->addPerParticleParameter("y0");
forceField->addPerParticleParameter("z0");
forceField->addPerParticleParameter("xscale");
forceField->addPerParticleParameter("yscale");
forceField->addPerParticleParameter("zscale");
vector<double> parameters(6);
parameters[0] = 1.0;
parameters[1] = 2.0;
parameters[2] = 3.0;
parameters[3] = 0.1;
parameters[4] = 0.2;
parameters[5] = 0.3;
forceField->addParticle(0, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, -1, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(2*0.1*1.0, 2*0.2*3.0, 2*0.3*3.0), forces[0], TOL);
ASSERT_EQUAL_TOL(0.1*1*1 + 0.2*3*3 + 0.3*3*3, state.getPotentialEnergy(), TOL);
}
void testParallelComputation() {
CudaPlatform platform;
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomExternalForce* force = new CustomExternalForce("x^2+y^2+z^2");
vector<double> params;
for (int i = 0; i < numParticles; i++)
force->addParticle(i, params);
system.addForce(force);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
map<string, string> props;
props[CudaPlatform::CudaDeviceIndex()] = deviceIndex+","+deviceIndex;
Context context2(system, integrator2, platform, props);
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
int main() {
try {
testForce();
testManyParameters();
// testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda2/tests/TestCudaCustomTorsionForce.cpp 0 → 100644
View file @ eb64fa2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of CustomTorsionForce.
*/
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomTorsionForce.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testTorsions() {
CudaPlatform platform;
// Create a system using a CustomTorsionForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomTorsionForce* custom = new CustomTorsionForce("k*(1+cos(n*theta-theta0))");
custom->addPerTorsionParameter("theta0");
custom->addPerTorsionParameter("n");
custom->addGlobalParameter("k", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 1;
custom->addTorsion(0, 1, 2, 3, parameters);
parameters[0] = 2.0;
parameters[1] = 2;
custom->addTorsion(1, 2, 3, 4, parameters);
customSystem.addForce(custom);
// Create an identical system using a PeriodicTorsionForce.
System harmonicSystem;
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
periodic->addTorsion(0, 1, 2, 3, 1, 1.5, 0.5);
periodic->addTorsion(1, 2, 3, 4, 2, 2.0, 0.5);
harmonicSystem.addForce(periodic);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(5);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(harmonicSystem, integrator2, platform);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 2;
custom->setTorsionParameters(0, 0, 1, 2, 3, parameters);
parameters[0] = 2.1;
parameters[1] = 3;
custom->setTorsionParameters(1, 1, 2, 3, 4, parameters);
custom->updateParametersInContext(c1);
periodic->setTorsionParameters(0, 0, 1, 2, 3, 2, 1.6, 0.5);
periodic->setTorsionParameters(1, 1, 2, 3, 4, 3, 2.1, 0.5);
periodic->updateParametersInContext(c2);
{
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void testRange() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
CustomTorsionForce* custom = new CustomTorsionForce("theta");
custom->addTorsion(0, 1, 2, 3, vector<double>());
system.addForce(custom);
// Set the atoms in various positions, and verify that the angle is always in the expected range.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(4);
VerletIntegrator integrator(0.01);
double minAngle = 1000;
double maxAngle = -1000;
Context context(system, integrator, platform);
for (int i = 0; i < 100; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
context.setPositions(positions);
double angle = context.getState(State::Energy).getPotentialEnergy();
if (angle < minAngle)
minAngle = angle;
if (angle > maxAngle)
maxAngle = angle;
}
ASSERT(minAngle >= -M_PI);
ASSERT(maxAngle <= M_PI);
}
void testParallelComputation() {
CudaPlatform platform;
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomTorsionForce* force = new CustomTorsionForce("sin(theta-1.1)");
vector<double> params;
for (int i = 3; i < numParticles; i++)
force->addTorsion(i-3, i-2, i-1, i, params);
system.addForce(force);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(i, i%2, i%3);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
map<string, string> props;
props[CudaPlatform::CudaDeviceIndex()] = deviceIndex+","+deviceIndex;
Context context2(system, integrator2, platform, props);
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
int main() {
try {
testTorsions();
testRange();
// testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment