Skip to content

GitLab

Commit 4d827cd9 authored by peastman's avatar peastman
Browse files

More optimizations to CPU CustomGBForce

parent 468a8a5c
Show whitespace changes
Inline Side-by-side
Showing with 137 additions and 135 deletions
platforms/cpu/include/CpuCustomGBForce.h
View file @ 4d827cd9
......@@ -25,36 +25,36 @@
#ifndef OPENMM_CPU_CUSTOM_GB_FORCE_H__
#define OPENMM_CPU_CUSTOM_GB_FORCE_H__
#include "ReferenceNeighborList.h"
#include "CompiledExpressionSet.h"
#include "CpuNeighborList.h"
#include "lepton/CompiledExpression.h"
#include "openmm/CustomGBForce.h"
#include <map>
#include <set>
#include <vector>
namespace OpenMM {
class CpuCustomGBForce {
private:
bool cutoff;
bool periodic;
const OpenMM::NeighborList* neighborList;
const CpuNeighborList* neighborList;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
OpenMM::CompiledExpressionSet expressionSet;
CompiledExpressionSet expressionSet;
std::vector<Lepton::CompiledExpression> valueExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
std::vector<std::string> valueNames;
std::vector<int> valueIndex;
std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
std::vector<CustomGBForce::ComputationType> valueTypes;
std::vector<Lepton::CompiledExpression> energyExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
std::vector<std::string> paramNames;
std::vector<int> paramIndex;
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
std::vector<std::string> particleParamNames;
std::vector<std::string> particleValueNames;
std::vector<CustomGBForce::ComputationType> energyTypes;
std::vector<int> particleParamIndex;
std::vector<int> particleValueIndex;
int xindex, yindex, zindex, rindex;
......@@ -66,12 +66,11 @@ private:
* @param numAtoms number of atoms
* @param atomCoordinates atom coordinates
* @param values the vector to store computed values into
* @param globalParameters the values of global parameters
* @param atomParameters atomParameters[atomIndex][paramterIndex]
*/
void calculateSingleParticleValue(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters);
void calculateSingleParticleValue(int index, int numAtoms, std::vector<RealVec>& atomCoordinates, std::vector<std::vector<RealOpenMM> >& values,
RealOpenMM** atomParameters);
/**
* Calculate a computed value that is based on particle pairs
......@@ -81,14 +80,12 @@ private:
* @param atomCoordinates atom coordinates
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param values the vector to store computed values into
* @param globalParameters the values of global parameters
* @param exclusions exclusions[i] is the set of excluded indices for atom i
* @param useExclusions specifies whether to use exclusions
*/
void calculateParticlePairValue(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
void calculateParticlePairValue(int index, int numAtoms, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions, bool useExclusions);
/**
......@@ -99,12 +96,10 @@ private:
* @param atom2 the index of the second atom in the pair
* @param atomCoordinates atom coordinates
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param globalParameters the values of global parameters
* @param values the vector to store computed values into
*/
void calculateOnePairValue(int index, int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
void calculateOnePairValue(int index, int atom1, int atom2, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
std::vector<std::vector<RealOpenMM> >& values);
/**
......@@ -114,15 +109,14 @@ private:
* @param numAtoms number of atoms
* @param atomCoordinates atom coordinates
* @param values the vector containing computed values
* @param globalParameters the values of global parameters
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param forces forces on atoms are added to this
* @param totalEnergy the energy contribution is added to this
* @param dEdV the derivative of energy with respect to computed values is stored in this
*/
void calculateSingleParticleEnergyTerm(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
void calculateSingleParticleEnergyTerm(int index, int numAtoms, std::vector<RealVec>& atomCoordinates, const std::vector<std::vector<RealOpenMM> >& values,
RealOpenMM** atomParameters, std::vector<RealVec>& forces,
RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
/**
......@@ -133,7 +127,6 @@ private:
* @param atomCoordinates atom coordinates
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param values the vector containing computed values
* @param globalParameters the values of global parameters
* @param exclusions exclusions[i] is the set of excluded indices for atom i
* @param useExclusions specifies whether to use exclusions
* @param forces forces on atoms are added to this
......@@ -141,11 +134,10 @@ private:
* @param dEdV the derivative of energy with respect to computed values is stored in this
*/
void calculateParticlePairEnergyTerm(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
void calculateParticlePairEnergyTerm(int index, int numAtoms, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions, bool useExclusions,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
std::vector<RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
/**
* Evaluate a single atom pair as part of calculating an energy term
......@@ -155,17 +147,15 @@ private:
* @param atom2 the index of the second atom in the pair
* @param atomCoordinates atom coordinates
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param globalParameters the values of global parameters
* @param values the vector containing computed values
* @param forces forces on atoms are added to this
* @param totalEnergy the energy contribution is added to this
* @param dEdV the derivative of energy with respect to computed values is stored in this
*/
void calculateOnePairEnergyTerm(int index, int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
void calculateOnePairEnergyTerm(int index, int atom1, int atom2, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
std::vector<RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
/**
* Apply the chain rule to compute forces on atoms
......@@ -174,17 +164,15 @@ private:
* @param atomCoordinates atom coordinates
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param values the vector containing computed values
* @param globalParameters the values of global parameters
* @param exclusions exclusions[i] is the set of excluded indices for atom i
* @param forces forces on atoms are added to this
* @param dEdV the derivative of energy with respect to computed values is stored in this
*/
void calculateChainRuleForces(int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
void calculateChainRuleForces(int numAtoms, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions,
std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV);
std::vector<RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV);
/**
* Evaluate a single atom pair as part of applying the chain rule
......@@ -193,17 +181,15 @@ private:
* @param atom2 the index of the second atom in the pair
* @param atomCoordinates atom coordinates
* @param atomParameters atomParameters[atomIndex][paramterIndex]
* @param globalParameters the values of global parameters
* @param values the vector containing computed values
* @param forces forces on atoms are added to this
* @param dEdV the derivative of energy with respect to computed values is stored in this
* @param isExcluded specifies whether this is an excluded pair
*/
void calculateOnePairChainRule(int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
void calculateOnePairChainRule(int atom1, int atom2, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV,
std::vector<RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV,
bool isExcluded);
public:
......@@ -216,11 +202,11 @@ public:
const std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions,
const std::vector<std::string>& valueNames,
const std::vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
const std::vector<CustomGBForce::ComputationType>& valueTypes,
const std::vector<Lepton::CompiledExpression>& energyExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions,
const std::vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
const std::vector<CustomGBForce::ComputationType>& energyTypes,
const std::vector<std::string>& parameterNames);
~CpuCustomGBForce();
......@@ -232,7 +218,7 @@ public:
* @param neighbors the neighbor list to use
*/
void setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors);
void setUseCutoff(RealOpenMM distance, const CpuNeighborList& neighbors);
/**
* Set the force to use periodic boundary conditions. This requires that a cutoff has
......@@ -242,7 +228,7 @@ public:
* @param boxSize the X, Y, and Z widths of the periodic box
*/
void setPeriodic(OpenMM::RealVec& boxSize);
void setPeriodic(RealVec& boxSize);
/**
* Calculate custom GB ixn
......@@ -256,8 +242,10 @@ public:
* @param totalEnergy total energy
*/
void calculateIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
std::map<std::string, double>& globalParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy);
void calculateIxn(int numberOfAtoms, std::vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
std::map<std::string, double>& globalParameters, std::vector<RealVec>& forces, RealOpenMM* totalEnergy);
};
} // namespace OpenMM
#endif // OPENMM_CPU_CUSTOM_GB_FORCE_H__
platforms/cpu/include/CpuKernels.h
View file @ 4d827cd9
......@@ -342,18 +342,13 @@ private:
bool isPeriodic;
RealOpenMM **particleParamArray;
RealOpenMM nonbondedCutoff;
CpuCustomGBForce* ixn;
std::vector<std::set<int> > exclusions;
std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
std::vector<Lepton::CompiledExpression> valueExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
std::vector<Lepton::CompiledExpression> energyExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
NonbondedMethod nonbondedMethod;
NeighborList* neighborList;
CpuNeighborList* neighborList;
};
/**
......
platforms/cpu/src/CpuCustomGBForce.cpp
View file @ 4d827cd9
......@@ -31,22 +31,18 @@
#include "ReferenceForce.h"
#include "CpuCustomGBForce.h"
using std::map;
using std::set;
using std::string;
using std::stringstream;
using std::vector;
using OpenMM::RealVec;
using namespace OpenMM;
using namespace std;
CpuCustomGBForce::CpuCustomGBForce(const vector<Lepton::CompiledExpression>& valueExpressions,
const vector<vector<Lepton::CompiledExpression> > valueDerivExpressions,
const vector<vector<Lepton::CompiledExpression> > valueGradientExpressions,
const vector<string>& valueNames,
const vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
const vector<CustomGBForce::ComputationType>& valueTypes,
const vector<Lepton::CompiledExpression>& energyExpressions,
const vector<vector<Lepton::CompiledExpression> > energyDerivExpressions,
const vector<vector<Lepton::CompiledExpression> > energyGradientExpressions,
const vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
const vector<CustomGBForce::ComputationType>& energyTypes,
const vector<string>& parameterNames) :
cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions),
valueNames(valueNames), valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions),
......@@ -76,7 +72,6 @@ CpuCustomGBForce::CpuCustomGBForce(const vector<Lepton::CompiledExpression>& val
for (int j = 1; j < 3; j++) {
stringstream name;
name << paramNames[i] << j;
particleParamNames.push_back(name.str());
particleParamIndex.push_back(expressionSet.getVariableIndex(name.str()));
}
}
......@@ -85,7 +80,6 @@ CpuCustomGBForce::CpuCustomGBForce(const vector<Lepton::CompiledExpression>& val
for (int j = 1; j < 3; j++) {
stringstream name;
name << valueNames[i] << j;
particleValueNames.push_back(name.str());
particleValueIndex.push_back(expressionSet.getVariableIndex(name.str()));
}
}
......@@ -94,11 +88,11 @@ CpuCustomGBForce::CpuCustomGBForce(const vector<Lepton::CompiledExpression>& val
CpuCustomGBForce::~CpuCustomGBForce() {
}
void CpuCustomGBForce::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors) {
void CpuCustomGBForce::setUseCutoff(RealOpenMM distance, const CpuNeighborList& neighbors) {
cutoff = true;
cutoffDistance = distance;
neighborList = &neighbors;
}
}
void CpuCustomGBForce::setPeriodic(RealVec& boxSize) {
......@@ -124,33 +118,33 @@ void CpuCustomGBForce::calculateIxn(int numberOfAtoms, vector<RealVec>& atomCoor
int numValues = valueTypes.size();
vector<vector<RealOpenMM> > values(numValues);
for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
if (valueTypes[valueIndex] == OpenMM::CustomGBForce::SingleParticle)
calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters);
else if (valueTypes[valueIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true);
if (valueTypes[valueIndex] == CustomGBForce::SingleParticle)
calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, values, atomParameters);
else if (valueTypes[valueIndex] == CustomGBForce::ParticlePair)
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, true);
else
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false);
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, false);
}
// Now calculate the energy and its derivatives.
vector<vector<RealOpenMM> > dEdV(numValues, vector<RealOpenMM>(numberOfAtoms, (RealOpenMM) 0));
for (int termIndex = 0; termIndex < (int) energyExpressions.size(); termIndex++) {
if (energyTypes[termIndex] == OpenMM::CustomGBForce::SingleParticle)
calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters, forces, totalEnergy, dEdV);
else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy, dEdV);
if (energyTypes[termIndex] == CustomGBForce::SingleParticle)
calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, values, atomParameters, forces, totalEnergy, dEdV);
else if (energyTypes[termIndex] == CustomGBForce::ParticlePair)
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, true, forces, totalEnergy, dEdV);
else
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy, dEdV);
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, false, forces, totalEnergy, dEdV);
}
// Apply the chain rule to evaluate forces.
calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, forces, dEdV);
calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, forces, dEdV);
}
void CpuCustomGBForce::calculateSingleParticleValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters) {
RealOpenMM** atomParameters) {
values[index].resize(numAtoms);
for (int i = 0; i < numAtoms; i++) {
expressionSet.setVariable(xindex, atomCoordinates[i][0]);
......@@ -165,19 +159,29 @@ void CpuCustomGBForce::calculateSingleParticleValue(int index, int numAtoms, vec
}
void CpuCustomGBForce::calculateParticlePairValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions) {
vector<vector<RealOpenMM> >& values, const vector<set<int> >& exclusions, bool useExclusions) {
values[index].resize(numAtoms);
for (int i = 0; i < numAtoms; i++)
values[index][i] = (RealOpenMM) 0.0;
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
for (int i = 0; i < (int) neighbors.size(); i++) {
int first = neighbors[i];
for (int k = 0; k < 4; k++) {
if ((blockExclusions[i] & (1<<k)) == 0) {
int second = blockAtom[k];
if (useExclusions && exclusions[first].find(second) != exclusions[first].end())
continue;
calculateOnePairValue(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, first, second, atomCoordinates, atomParameters, values);
calculateOnePairValue(index, second, first, atomCoordinates, atomParameters, values);
}
}
}
}
}
else {
......@@ -187,15 +191,15 @@ void CpuCustomGBForce::calculateParticlePairValue(int index, int numAtoms, vecto
for (int j = i+1; j < numAtoms; j++) {
if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue;
calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, j, i, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, values);
calculateOnePairValue(index, j, i, atomCoordinates, atomParameters, values);
}
}
}
}
void CpuCustomGBForce::calculateOnePairValue(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, vector<vector<RealOpenMM> >& values) {
vector<vector<RealOpenMM> >& values) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
......@@ -217,7 +221,7 @@ void CpuCustomGBForce::calculateOnePairValue(int index, int atom1, int atom2, ve
}
void CpuCustomGBForce::calculateSingleParticleEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates, const vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM* totalEnergy,
RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) {
for (int i = 0; i < numAtoms; i++) {
expressionSet.setVariable(xindex, atomCoordinates[i][0]);
......@@ -238,16 +242,26 @@ void CpuCustomGBForce::calculateSingleParticleEnergyTerm(int index, int numAtoms
}
void CpuCustomGBForce::calculateParticlePairEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions,
const vector<vector<RealOpenMM> >& values, const vector<set<int> >& exclusions, bool useExclusions,
vector<RealVec>& forces, RealOpenMM* totalEnergy, vector<vector<RealOpenMM> >& dEdV) {
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
for (int i = 0; i < (int) neighbors.size(); i++) {
int first = neighbors[i];
for (int k = 0; k < 4; k++) {
if ((blockExclusions[i] & (1<<k)) == 0) {
int second = blockAtom[k];
if (useExclusions && exclusions[first].find(second) != exclusions[first].end())
continue;
calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
calculateOnePairEnergyTerm(index, first, second, atomCoordinates, atomParameters, values, forces, totalEnergy, dEdV);
}
}
}
}
}
else {
......@@ -257,14 +271,14 @@ void CpuCustomGBForce::calculateParticlePairEnergyTerm(int index, int numAtoms,
for (int j = i+1; j < numAtoms; j++) {
if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue;
calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, values, forces, totalEnergy, dEdV);
}
}
}
}
void CpuCustomGBForce::calculateOnePairEnergyTerm(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces, RealOpenMM* totalEnergy,
const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) {
// Compute the displacement.
......@@ -306,16 +320,25 @@ void CpuCustomGBForce::calculateOnePairEnergyTerm(int index, int atom1, int atom
}
void CpuCustomGBForce::calculateChainRuleForces(int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters,
const vector<set<int> >& exclusions, vector<RealVec>& forces, vector<vector<RealOpenMM> >& dEdV) {
const vector<vector<RealOpenMM> >& values, const vector<set<int> >& exclusions, vector<RealVec>& forces, vector<vector<RealOpenMM> >& dEdV) {
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
bool isExcluded = (exclusions[pair.first].find(pair.second) != exclusions[pair.first].end());
calculateOnePairChainRule(pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
for (int i = 0; i < (int) neighbors.size(); i++) {
int first = neighbors[i];
for (int k = 0; k < 4; k++) {
if ((blockExclusions[i] & (1<<k)) == 0) {
int second = blockAtom[k];
bool isExcluded = (exclusions[first].find(second) != exclusions[first].end());
calculateOnePairChainRule(first, second, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(second, first, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
}
}
}
}
}
else {
......@@ -324,8 +347,8 @@ void CpuCustomGBForce::calculateChainRuleForces(int numAtoms, vector<RealVec>& a
for (int i = 0; i < numAtoms; i++) {
for (int j = i+1; j < numAtoms; j++) {
bool isExcluded = (exclusions[i].find(j) != exclusions[i].end());
calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
}
}
}
......@@ -360,8 +383,7 @@ void CpuCustomGBForce::calculateChainRuleForces(int numAtoms, vector<RealVec>& a
}
void CpuCustomGBForce::calculateOnePairChainRule(int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces,
vector<vector<RealOpenMM> >& dEdV, bool isExcluded) {
const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces, vector<vector<RealOpenMM> >& dEdV, bool isExcluded) {
// Compute the displacement.
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
......@@ -391,7 +413,7 @@ void CpuCustomGBForce::calculateOnePairChainRule(int atom1, int atom2, vector<Re
deltaR[2] *= rinv;
vector<RealOpenMM> dVdR1(valueDerivExpressions.size(), 0.0);
vector<RealOpenMM> dVdR2(valueDerivExpressions.size(), 0.0);
if (!isExcluded || valueTypes[0] != OpenMM::CustomGBForce::ParticlePair) {
if (!isExcluded || valueTypes[0] != CustomGBForce::ParticlePair) {
dVdR1[0] = (RealOpenMM) valueDerivExpressions[0][0].evaluate();
dVdR2[0] = -dVdR1[0];
for (int i = 0; i < 3; i++) {
......
platforms/cpu/src/CpuCustomNonbondedForce.cpp
View file @ 4d827cd9
......@@ -182,9 +182,7 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
if (blockIndex >= neighborList->getNumBlocks())
break;
int blockAtom[4];
for (int i = 0; i < 4; i++)
blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
for (int i = 0; i < (int) neighbors.size(); i++) {
......
platforms/cpu/src/CpuKernels.cpp
View file @ 4d827cd9
......@@ -843,6 +843,8 @@ CpuCalcCustomGBForceKernel::~CpuCalcCustomGBForceKernel() {
}
if (neighborList != NULL)
delete neighborList;
if (ixn != NULL)
delete ixn;
}
void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
......@@ -891,7 +893,7 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
if (nonbondedMethod == NoCutoff)
neighborList = NULL;
else
neighborList = new NeighborList();
neighborList = new CpuNeighborList(4);
// Create custom functions for the tabulated functions.
......@@ -901,8 +903,10 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
// Parse the expressions for computed values.
valueDerivExpressions.resize(force.getNumComputedValues());
valueGradientExpressions.resize(force.getNumComputedValues());
vector<vector<Lepton::CompiledExpression> > valueDerivExpressions(force.getNumComputedValues());
vector<vector<Lepton::CompiledExpression> > valueGradientExpressions(force.getNumComputedValues());
vector<Lepton::CompiledExpression> valueExpressions;
vector<Lepton::CompiledExpression> energyExpressions;
for (int i = 0; i < force.getNumComputedValues(); i++) {
string name, expression;
CustomGBForce::ComputationType type;
......@@ -924,8 +928,8 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
// Parse the expressions for energy terms.
energyDerivExpressions.resize(force.getNumEnergyTerms());
energyGradientExpressions.resize(force.getNumEnergyTerms());
vector<vector<Lepton::CompiledExpression> > energyDerivExpressions(force.getNumEnergyTerms());
vector<vector<Lepton::CompiledExpression> > energyGradientExpressions(force.getNumEnergyTerms());
for (int i = 0; i < force.getNumEnergyTerms(); i++) {
string expression;
CustomGBForce::ComputationType type;
......@@ -953,25 +957,28 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
delete iter->second;
ixn = new CpuCustomGBForce(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
data.isPeriodic = (force.getNonbondedMethod() == CustomGBForce::CutoffPeriodic);
}
double CpuCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
CpuCustomGBForce ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (periodic)
ixn.setPeriodic(extractBoxSize(context));
RealVec& box = extractBoxSize(context);
float floatBoxSize[3] = {(float) box[0], (float) box[1], (float) box[2]};
if (data.isPeriodic)
ixn->setPeriodic(extractBoxSize(context));
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
ixn.setUseCutoff(nonbondedCutoff, *neighborList);
vector<set<int> > noExclusions(numParticles);
neighborList->computeNeighborList(numParticles, data.posq, exclusions, floatBoxSize, data.isPeriodic, nonbondedCutoff, data.threads);
ixn->setUseCutoff(nonbondedCutoff, *neighborList);
}
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
ixn->calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
return energy;
}
......
platforms/cpu/src/CpuNonbondedForceVec4.cpp
View file @ 4d827cd9
......@@ -48,13 +48,11 @@ CpuNonbondedForceVec4::CpuNonbondedForceVec4() {
void CpuNonbondedForceVec4::calculateBlockIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
// Load the positions and parameters of the atoms in the block.
int blockAtom[4];
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
fvec4 blockAtomPosq[4];
fvec4 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
for (int i = 0; i < 4; i++) {
blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
for (int i = 0; i < 4; i++)
blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
}
fvec4 blockAtomX = fvec4(blockAtomPosq[0][0], blockAtomPosq[1][0], blockAtomPosq[2][0], blockAtomPosq[3][0]);
fvec4 blockAtomY = fvec4(blockAtomPosq[0][1], blockAtomPosq[1][1], blockAtomPosq[2][1], blockAtomPosq[3][1]);
fvec4 blockAtomZ = fvec4(blockAtomPosq[0][2], blockAtomPosq[1][2], blockAtomPosq[2][2], blockAtomPosq[3][2]);
......@@ -159,13 +157,11 @@ void CpuNonbondedForceVec4::calculateBlockIxn(int blockIndex, float* forces, dou
void CpuNonbondedForceVec4::calculateBlockEwaldIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
// Load the positions and parameters of the atoms in the block.
int blockAtom[4];
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
fvec4 blockAtomPosq[4];
fvec4 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
for (int i = 0; i < 4; i++) {
blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
for (int i = 0; i < 4; i++)
blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
}
fvec4 blockAtomX = fvec4(blockAtomPosq[0][0], blockAtomPosq[1][0], blockAtomPosq[2][0], blockAtomPosq[3][0]);
fvec4 blockAtomY = fvec4(blockAtomPosq[0][1], blockAtomPosq[1][1], blockAtomPosq[2][1], blockAtomPosq[3][1]);
fvec4 blockAtomZ = fvec4(blockAtomPosq[0][2], blockAtomPosq[1][2], blockAtomPosq[2][2], blockAtomPosq[3][2]);
......
platforms/cpu/src/CpuNonbondedForceVec8.cpp
View file @ 4d827cd9
......@@ -74,14 +74,12 @@ CpuNonbondedForceVec8::CpuNonbondedForceVec8() {
void CpuNonbondedForceVec8::calculateBlockIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
// Load the positions and parameters of the atoms in the block.
int blockAtom[8];
const int* blockAtom = &neighborList->getSortedAtoms()[8*blockIndex];
fvec4 blockAtomPosq[8];
fvec8 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
fvec8 blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge;
for (int i = 0; i < 8; i++) {
blockAtom[i] = neighborList->getSortedAtoms()[8*blockIndex+i];
for (int i = 0; i < 8; i++)
blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
}
transpose(blockAtomPosq[0], blockAtomPosq[1], blockAtomPosq[2], blockAtomPosq[3], blockAtomPosq[4], blockAtomPosq[5], blockAtomPosq[6], blockAtomPosq[7], blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge);
blockAtomCharge *= ONE_4PI_EPS0;
fvec8 blockAtomSigma(atomParameters[blockAtom[0]].first, atomParameters[blockAtom[1]].first, atomParameters[blockAtom[2]].first, atomParameters[blockAtom[3]].first, atomParameters[blockAtom[4]].first, atomParameters[blockAtom[5]].first, atomParameters[blockAtom[6]].first, atomParameters[blockAtom[7]].first);
......@@ -184,14 +182,12 @@ void CpuNonbondedForceVec8::calculateBlockIxn(int blockIndex, float* forces, dou
void CpuNonbondedForceVec8::calculateBlockEwaldIxn(int blockIndex, float* forces, double* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
// Load the positions and parameters of the atoms in the block.
int blockAtom[8];
const int* blockAtom = &neighborList->getSortedAtoms()[8*blockIndex];
fvec4 blockAtomPosq[8];
fvec8 blockAtomForceX(0.0f), blockAtomForceY(0.0f), blockAtomForceZ(0.0f);
fvec8 blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge;
for (int i = 0; i < 8; i++) {
blockAtom[i] = neighborList->getSortedAtoms()[8*blockIndex+i];
for (int i = 0; i < 8; i++)
blockAtomPosq[i] = fvec4(posq+4*blockAtom[i]);
}
transpose(blockAtomPosq[0], blockAtomPosq[1], blockAtomPosq[2], blockAtomPosq[3], blockAtomPosq[4], blockAtomPosq[5], blockAtomPosq[6], blockAtomPosq[7], blockAtomX, blockAtomY, blockAtomZ, blockAtomCharge);
blockAtomCharge *= ONE_4PI_EPS0;
fvec8 blockAtomSigma(atomParameters[blockAtom[0]].first, atomParameters[blockAtom[1]].first, atomParameters[blockAtom[2]].first, atomParameters[blockAtom[3]].first, atomParameters[blockAtom[4]].first, atomParameters[blockAtom[5]].first, atomParameters[blockAtom[6]].first, atomParameters[blockAtom[7]].first);
......
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