Skip to content

GitLab

Commit 543d51d7 authored by peastman's avatar peastman
Browse files

Parallelized the computation of torsions with CPU platform

parent 4f339321
Hide whitespace changes
Inline Side-by-side
Showing with 796 additions and 0 deletions
platforms/cpu/include/CpuBondForce.h 0 → 100644
View file @ 543d51d7
#ifndef OPENMM_CPUBONDFORCE_H_
#define OPENMM_CPUBONDFORCE_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) 2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceBondIxn.h"
#include "windowsExportCpu.h"
#include "openmm/internal/ThreadPool.h"
#include <list>
#include <set>
#include <vector>
namespace OpenMM {
/**
* This class parallelizes the calculation of bonded forces.
*/
class OPENMM_EXPORT_CPU CpuBondForce {
public:
class ComputeForceTask;
CpuBondForce();
/**
* Analyze the set of bonds and decide which to compute with each thread.
*/
void initialize(int numAtoms, int numBonds, int numAtomsPerBond, int** bondAtoms, ThreadPool& threads);
/**
* Compute the forces from all bonds.
*/
void calculateForce(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** parameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn);
/**
* This routine contains the code executed by each thread.
*/
void threadComputeForce(ThreadPool& threads, int threadIndex, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** parameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn);
private:
bool canAssignBond(int bond, int thread, std::vector<int>& atomThread);
void assignBond(int bond, int thread, std::vector<int>& atomThread, std::vector<int>& bondThread, std::vector<std::set<int> >& atomBonds, std::list<int>& candidateBonds);
int numBonds, numAtomsPerBond;
int** bondAtoms;
ThreadPool* threads;
std::vector<std::vector<int> > threadBonds;
std::vector<int> extraBonds;
};
} // namespace OpenMM
#endif /*OPENMM_CPUBONDFORCE_H_*/
platforms/cpu/include/CpuKernels.h
View file @ 543d51d7
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
#include "CpuBondForce.h"
#include "CpuGBSAOBCForce.h" #include "CpuGBSAOBCForce.h"
#include "CpuLangevinDynamics.h" #include "CpuLangevinDynamics.h"
#include "CpuNeighborList.h" #include "CpuNeighborList.h"
...@@ -85,6 +86,86 @@ private: ...@@ -85,6 +86,86 @@ private:
Kernel referenceKernel; Kernel referenceKernel;
}; };
/**
* This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class CpuCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
public:
CpuCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
CalcPeriodicTorsionForceKernel(name, platform), data(data), torsionIndexArray(NULL), torsionParamArray(NULL) {
}
~CpuCalcPeriodicTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the PeriodicTorsionForce this kernel will be used for
*/
void initialize(const System& system, const PeriodicTorsionForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the PeriodicTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
private:
CpuPlatform::PlatformData& data;
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
CpuBondForce bondForce;
};
/**
* This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class CpuCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
public:
CpuCalcRBTorsionForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
CalcRBTorsionForceKernel(name, platform), data(data), torsionIndexArray(NULL), torsionParamArray(NULL) {
}
~CpuCalcRBTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the RBTorsionForce this kernel will be used for
*/
void initialize(const System& system, const RBTorsionForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the RBTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
private:
CpuPlatform::PlatformData& data;
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
CpuBondForce bondForce;
};
/** /**
* 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.
*/ */
......
platforms/cpu/src/CpuBondForce.cpp 0 → 100644
View file @ 543d51d7
/* -------------------------------------------------------------------------- *
* 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) 2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "CpuBondForce.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
using namespace std;
class CpuBondForce::ComputeForceTask : public ThreadPool::Task {
public:
ComputeForceTask(CpuBondForce& owner, vector<RealVec>& atomCoordinates, RealOpenMM** parameters, vector<RealVec>& forces,
vector<RealOpenMM>& threadEnergy, RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn) : owner(owner), atomCoordinates(atomCoordinates),
parameters(parameters), forces(forces), threadEnergy(threadEnergy), totalEnergy(totalEnergy), referenceBondIxn(referenceBondIxn) {
}
void execute(ThreadPool& threads, int threadIndex) {
RealOpenMM* energy = (totalEnergy == NULL ? NULL : &threadEnergy[threadIndex]);
owner.threadComputeForce(threads, threadIndex, atomCoordinates, parameters, forces, energy, referenceBondIxn);
}
CpuBondForce& owner;
vector<RealVec>& atomCoordinates;
RealOpenMM** parameters;
vector<RealVec>& forces;
vector<RealOpenMM>& threadEnergy;
RealOpenMM* totalEnergy;
ReferenceBondIxn& referenceBondIxn;
};
CpuBondForce::CpuBondForce() {
}
void CpuBondForce::initialize(int numAtoms, int numBonds, int numAtomsPerBond, int** bondAtoms, ThreadPool& threads) {
this->numBonds = numBonds;
this->numAtomsPerBond = numAtomsPerBond;
this->bondAtoms = bondAtoms;
this->threads = &threads;
int numThreads = threads.getNumThreads();
int targetBondsPerThread = numBonds/numThreads;
// Record the bonds that include each atom.
vector<set<int> > atomBonds(numAtoms);
for (int bond = 0; bond < numBonds; bond++) {
for (int i = 0; i < numAtomsPerBond; i++)
atomBonds[bondAtoms[bond][i]].insert(bond);
}
// Divide bonds into groups.
vector<int> atomThread(numAtoms, -1);
vector<int> bondThread(numBonds, -1);
threadBonds.resize(numThreads);
int numProcessed = 0;
int thread = 0;
list<int> candidateBonds;
while (thread < numThreads) {
// Find the next unassigned bond.
while (numProcessed < numBonds && bondThread[numProcessed] != -1)
numProcessed++;
if (numProcessed == numBonds)
break; // We've gone through the whole list of bonds.
// See if this bond can be assigned to this thread.
if (!canAssignBond(numProcessed, thread, atomThread)) {
numProcessed++;
continue;
}
// Assign this bond to the thread.
assignBond(numProcessed++, thread, atomThread, bondThread, atomBonds, candidateBonds);
// Assign additional bonds that have been identified as involving atoms assigned to this thread.
while (!candidateBonds.empty() && threadBonds[thread].size() < targetBondsPerThread) {
int bond = *candidateBonds.begin();
if (bondThread[bond] == -1 && canAssignBond(bond, thread, atomThread))
assignBond(bond, thread, atomThread, bondThread, atomBonds, candidateBonds);
candidateBonds.pop_front();
}
// If we have assigned enough bonds to this thread, move on to the next one.
if (threadBonds[thread].size() >= targetBondsPerThread) {
candidateBonds.clear();
thread++;
}
}
// Look through the remaining bonds and see whether any of them can be assigned.
candidateBonds.clear();
for (int bond = 0; bond < numBonds; bond++) {
if (bondThread[bond] == -1) {
// See whether this bond can be assigned to a thread.
bool canAssign = true;
int assignment = -1;
for (int i = 0; i < numAtomsPerBond; i++) {
int thread = atomThread[bondAtoms[bond][i]];
if (thread == -1 || thread == assignment)
continue;
if (assignment == -1)
assignment = thread;
else {
canAssign = false;
break;
}
}
if (canAssign) {
// Assign this bond to a thread.
if (assignment == -1)
assignment = numThreads-1;
assignBond(bond, assignment, atomThread, bondThread, atomBonds, candidateBonds);
}
else {
// Add it to the list of "extra" bonds.
extraBonds.push_back(bond);
}
}
}
}
bool CpuBondForce::canAssignBond(int bond, int thread, vector<int>& atomThread) {
for (int i = 0; i < numAtomsPerBond; i++) {
int atom = bondAtoms[bond][i];
if (atomThread[atom] != -1 && atomThread[atom] != thread)
return false;
}
return true;
}
void CpuBondForce::assignBond(int bond, int thread, vector<int>& atomThread, vector<int>& bondThread, vector<set<int> >& atomBonds, list<int>& candidateBonds) {
// Assign the bond to a thread.
bondThread[bond] = thread;
threadBonds[thread].push_back(bond);
// Mark every atom in this bond as also belonging to the thread, and add all of their
// bonds to the list of candidates.
for (int i = 0; i < numAtomsPerBond; i++) {
int& atom = atomThread[bondAtoms[bond][i]];
if (atom == thread)
continue;
if (atom != -1)
throw OpenMMException("CpuBondForce: Internal error: atoms assigned to threads incorrectly");
atom = thread;
for (set<int>::const_iterator iter = atomBonds[atom].begin(); iter != atomBonds[atom].end(); ++iter)
candidateBonds.push_back(*iter);
}
}
void CpuBondForce::calculateForce(vector<RealVec>& atomCoordinates, RealOpenMM** parameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn) {
// Have the worker threads compute their forces.
vector<RealOpenMM> threadEnergy(threads->getNumThreads(), 0);
ComputeForceTask task(*this, atomCoordinates, parameters, forces, threadEnergy, totalEnergy, referenceBondIxn);
threads->execute(task);
threads->waitForThreads();
// Compute any "extra" bonds.
for (int i = 0; i < extraBonds.size(); i++) {
int bond = extraBonds[i];
referenceBondIxn.calculateBondIxn(bondAtoms[bond], atomCoordinates, parameters[bond], forces, totalEnergy);
}
// Compute the total energy.
if (totalEnergy != NULL)
for (int i = 0; i < threads->getNumThreads(); i++)
*totalEnergy += threadEnergy[i];
}
void CpuBondForce::threadComputeForce(ThreadPool& threads, int threadIndex, vector<RealVec>& atomCoordinates, RealOpenMM** parameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn) {
vector<int>& bonds = threadBonds[threadIndex];
int numBonds = bonds.size();
for (int i = 0; i < numBonds; i++) {
int bond = bonds[i];
referenceBondIxn.calculateBondIxn(bondAtoms[bond], atomCoordinates, parameters[bond], forces, totalEnergy);
}
}
\ No newline at end of file
platforms/cpu/src/CpuKernelFactory.cpp
View file @ 543d51d7
...@@ -41,6 +41,10 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform& ...@@ -41,6 +41,10 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
CpuPlatform::PlatformData& data = CpuPlatform::getPlatformData(context); CpuPlatform::PlatformData& data = CpuPlatform::getPlatformData(context);
if (name == CalcForcesAndEnergyKernel::Name()) if (name == CalcForcesAndEnergyKernel::Name())
return new CpuCalcForcesAndEnergyKernel(name, platform, data, context); return new CpuCalcForcesAndEnergyKernel(name, platform, data, context);
if (name == CalcPeriodicTorsionForceKernel::Name())
return new CpuCalcPeriodicTorsionForceKernel(name, platform, data);
if (name == CalcRBTorsionForceKernel::Name())
return new CpuCalcRBTorsionForceKernel(name, platform, data);
if (name == CalcNonbondedForceKernel::Name()) if (name == CalcNonbondedForceKernel::Name())
return new CpuCalcNonbondedForceKernel(name, platform, data); return new CpuCalcNonbondedForceKernel(name, platform, data);
if (name == CalcGBSAOBCForceKernel::Name()) if (name == CalcGBSAOBCForceKernel::Name())
......
platforms/cpu/src/CpuKernels.cpp
View file @ 543d51d7
...@@ -35,6 +35,8 @@ ...@@ -35,6 +35,8 @@
#include "ReferenceKernelFactory.h" #include "ReferenceKernelFactory.h"
#include "ReferenceKernels.h" #include "ReferenceKernels.h"
#include "ReferenceLJCoulomb14.h" #include "ReferenceLJCoulomb14.h"
#include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h"
#include "openmm/Context.h" #include "openmm/Context.h"
#include "openmm/OpenMMException.h" #include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
...@@ -183,6 +185,134 @@ double CpuCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, boo ...@@ -183,6 +185,134 @@ double CpuCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, boo
return referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().finishComputation(context, includeForce, includeEnergy, groups); return referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().finishComputation(context, includeForce, includeEnergy, groups);
} }
CpuCalcPeriodicTorsionForceKernel::~CpuCalcPeriodicTorsionForceKernel() {
if (torsionIndexArray != NULL) {
for (int i = 0; i < numTorsions; i++) {
delete[] torsionIndexArray[i];
delete[] torsionParamArray[i];
}
delete[] torsionIndexArray;
delete[] torsionParamArray;
}
}
void CpuCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionIndexArray = new int*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionIndexArray[i] = new int[4];
torsionParamArray = new RealOpenMM*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionParamArray[i] = new RealOpenMM[3];
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
torsionIndexArray[i][0] = particle1;
torsionIndexArray[i][1] = particle2;
torsionIndexArray[i][2] = particle3;
torsionIndexArray[i][3] = particle4;
torsionParamArray[i][0] = (RealOpenMM) k;
torsionParamArray[i][1] = (RealOpenMM) phase;
torsionParamArray[i][2] = (RealOpenMM) periodicity;
}
bondForce.initialize(system.getNumParticles(), numTorsions, 4, torsionIndexArray, data.threads);
}
double CpuCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
ReferenceProperDihedralBond periodicTorsionBond;
bondForce.calculateForce(posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
return energy;
}
void CpuCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
torsionParamArray[i][0] = (RealOpenMM) k;
torsionParamArray[i][1] = (RealOpenMM) phase;
torsionParamArray[i][2] = (RealOpenMM) periodicity;
}
}
CpuCalcRBTorsionForceKernel::~CpuCalcRBTorsionForceKernel() {
if (torsionIndexArray != NULL) {
for (int i = 0; i < numTorsions; i++) {
delete[] torsionIndexArray[i];
delete[] torsionParamArray[i];
}
delete[] torsionIndexArray;
delete[] torsionParamArray;
}
}
void CpuCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionIndexArray = new int*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionIndexArray[i] = new int[4];
torsionParamArray = new RealOpenMM*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionParamArray[i] = new RealOpenMM[6];
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
torsionIndexArray[i][0] = particle1;
torsionIndexArray[i][1] = particle2;
torsionIndexArray[i][2] = particle3;
torsionIndexArray[i][3] = particle4;
torsionParamArray[i][0] = (RealOpenMM) c0;
torsionParamArray[i][1] = (RealOpenMM) c1;
torsionParamArray[i][2] = (RealOpenMM) c2;
torsionParamArray[i][3] = (RealOpenMM) c3;
torsionParamArray[i][4] = (RealOpenMM) c4;
torsionParamArray[i][5] = (RealOpenMM) c5;
}
bondForce.initialize(system.getNumParticles(), numTorsions, 4, torsionIndexArray, data.threads);
}
double CpuCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
ReferenceRbDihedralBond rbTorsionBond;
bondForce.calculateForce(posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
return energy;
}
void CpuCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& context, const RBTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
torsionParamArray[i][0] = (RealOpenMM) c0;
torsionParamArray[i][1] = (RealOpenMM) c1;
torsionParamArray[i][2] = (RealOpenMM) c2;
torsionParamArray[i][3] = (RealOpenMM) c3;
torsionParamArray[i][4] = (RealOpenMM) c4;
torsionParamArray[i][5] = (RealOpenMM) c5;
}
}
class CpuCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO { class CpuCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
public: public:
PmeIO(float* posq, float* force, int numParticles) : posq(posq), force(force), numParticles(numParticles) { PmeIO(float* posq, float* force, int numParticles) : posq(posq), force(force), numParticles(numParticles) {
......
platforms/cpu/src/CpuPlatform.cpp
View file @ 543d51d7
...@@ -51,6 +51,8 @@ map<ContextImpl*, CpuPlatform::PlatformData*> CpuPlatform::contextData; ...@@ -51,6 +51,8 @@ map<ContextImpl*, CpuPlatform::PlatformData*> CpuPlatform::contextData;
CpuPlatform::CpuPlatform() { CpuPlatform::CpuPlatform() {
CpuKernelFactory* factory = new CpuKernelFactory(); CpuKernelFactory* factory = new CpuKernelFactory();
registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory); registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory); registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory); registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory); registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
......
platforms/cpu/tests/TestCpuPeriodicTorsionForce.cpp 0 → 100644
View file @ 543d51d7
/* -------------------------------------------------------------------------- *
* 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-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CPU implementation of PeriodicTorsionForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CpuPlatform platform;
const double TOL = 1e-5;
void testPeriodicTorsions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta = (3*PI_M/2)-(PI_M/3.2);
double torque = -3*1.3*std::sin(dtheta);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
}
}
void testParallelComputation() {
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
PeriodicTorsionForce* force = new PeriodicTorsionForce();
for (int i = 3; i < numParticles; i++)
force->addTorsion(i-3, i-2, i-1, i, 2, 1.1, i);
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);
ReferencePlatform reference;
Context context1(system, integrator1, reference);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
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(int argc, char* argv[]) {
try {
testPeriodicTorsions();
testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cpu/tests/TestCpuRBTorsionForce.cpp 0 → 100644
View file @ 543d51d7
/* -------------------------------------------------------------------------- *
* 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-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of RBTorsionForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/RBTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CpuPlatform platform;
const double TOL = 1e-5;
void testRBTorsions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
RBTorsionForce* forceField = new RBTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.1*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.1*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 0.11, 0.22, 0.33, 0.44, 0.55, 0.66);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.11*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.11*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testParallelComputation() {
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
RBTorsionForce* force = new RBTorsionForce();
for (int i = 3; i < numParticles; i++)
force->addTorsion(i-3, i-2, i-1, i, 2, 0.1*i, 0.5*i, i, 1, 1);
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);
ReferencePlatform reference;
Context context1(system, integrator1, reference);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
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(int argc, char* argv[]) {
try {
testRBTorsions();
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