Skip to content

GitLab

Commit 7fb10336 authored by peastman's avatar peastman
Browse files

Cleaned up lots of formatting to be more consistent with the rest of OpenMM

parent 1d3ffd7b
Show whitespace changes
Inline Side-by-side
Showing with 219 additions and 219 deletions
platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp
View file @ 7fb10336
...@@ -71,7 +71,7 @@ ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesP ...@@ -71,7 +71,7 @@ ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesP
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomCompoundBondIxn::~ReferenceCustomCompoundBondIxn( ){ ReferenceCustomCompoundBondIxn::~ReferenceCustomCompoundBondIxn() {
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -92,7 +92,7 @@ void ReferenceCustomCompoundBondIxn::calculatePairIxn(vector<RealVec>& atomCoord ...@@ -92,7 +92,7 @@ void ReferenceCustomCompoundBondIxn::calculatePairIxn(vector<RealVec>& atomCoord
map<string, double> variables = globalParameters; map<string, double> variables = globalParameters;
int numBonds = bondAtoms.size(); int numBonds = bondAtoms.size();
for (int bond = 0; bond < numBonds; bond++){ for (int bond = 0; bond < numBonds; bond++) {
for (int j = 0; j < (int) bondParamNames.size(); j++) for (int j = 0; j < (int) bondParamNames.size(); j++)
variables[bondParamNames[j]] = bondParameters[bond][j]; variables[bondParamNames[j]] = bondParameters[bond][j];
calculateOneIxn(bond, atomCoordinates, variables, forces, totalEnergy); calculateOneIxn(bond, atomCoordinates, variables, forces, totalEnergy);
......
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp
View file @ 7fb10336
...@@ -96,7 +96,7 @@ ReferenceCustomDynamics::~ReferenceCustomDynamics() { ...@@ -96,7 +96,7 @@ ReferenceCustomDynamics::~ReferenceCustomDynamics() {
void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates, void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses, vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses,
map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid, RealOpenMM tolerance){ map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid, RealOpenMM tolerance) {
int numSteps = stepType.size(); int numSteps = stepType.size();
globals.insert(context.getParameters().begin(), context.getParameters().end()); globals.insert(context.getParameters().begin(), context.getParameters().end());
oldPos = atomCoordinates; oldPos = atomCoordinates;
......
platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
View file @ 7fb10336
...@@ -78,7 +78,7 @@ ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::CompiledExp ...@@ -78,7 +78,7 @@ ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::CompiledExp
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomExternalIxn::~ReferenceCustomExternalIxn( ){ ReferenceCustomExternalIxn::~ReferenceCustomExternalIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -100,11 +100,11 @@ ReferenceCustomExternalIxn::~ReferenceCustomExternalIxn( ){ ...@@ -100,11 +100,11 @@ ReferenceCustomExternalIxn::~ReferenceCustomExternalIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomExternalIxn::calculateForce( int atomIndex, void ReferenceCustomExternalIxn::calculateForce(int atomIndex,
vector<RealVec>& atomCoordinates, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, RealOpenMM* parameters,
vector<RealVec>& forces, vector<RealVec>& forces,
RealOpenMM* energy ) const { RealOpenMM* energy) const {
static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn"; static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
......
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @ 7fb10336
...@@ -85,7 +85,7 @@ ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgra ...@@ -85,7 +85,7 @@ ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgra
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){ ReferenceCustomGBIxn::~ReferenceCustomGBIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -104,7 +104,7 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){ ...@@ -104,7 +104,7 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomGBIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors ) { void ReferenceCustomGBIxn::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors) {
cutoff = true; cutoff = true;
cutoffDistance = distance; cutoffDistance = distance;
...@@ -202,8 +202,8 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v ...@@ -202,8 +202,8 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v
else { else {
// Perform an O(N^2) loop over all atom pairs. // Perform an O(N^2) loop over all atom pairs.
for (int i = 0; i < numAtoms; i++){ for (int i = 0; i < numAtoms; i++) {
for (int j = i+1; j < numAtoms; j++ ){ for (int j = i+1; j < numAtoms; j++) {
if (useExclusions && exclusions[i].find(j) != exclusions[i].end()) if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue; continue;
calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, globalParameters, values); calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, globalParameters, values);
...@@ -274,8 +274,8 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto ...@@ -274,8 +274,8 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto
else { else {
// Perform an O(N^2) loop over all atom pairs. // Perform an O(N^2) loop over all atom pairs.
for (int i = 0; i < numAtoms; i++){ for (int i = 0; i < numAtoms; i++) {
for (int j = i+1; j < numAtoms; j++ ){ for (int j = i+1; j < numAtoms; j++) {
if (useExclusions && exclusions[i].find(j) != exclusions[i].end()) if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue; continue;
calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV); calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
...@@ -343,8 +343,8 @@ void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec ...@@ -343,8 +343,8 @@ void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec
else { else {
// Perform an O(N^2) loop over all atom pairs. // Perform an O(N^2) loop over all atom pairs.
for (int i = 0; i < numAtoms; i++){ for (int i = 0; i < numAtoms; i++) {
for (int j = i+1; j < numAtoms; j++ ){ for (int j = i+1; j < numAtoms; j++) {
bool isExcluded = (exclusions[i].find(j) != exclusions[i].end()); bool isExcluded = (exclusions[i].find(j) != exclusions[i].end());
calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded); calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded); calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
......
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
View file @ 7fb10336
...@@ -64,7 +64,7 @@ ReferenceCustomHbondIxn::ReferenceCustomHbondIxn(const vector<vector<int> >& don ...@@ -64,7 +64,7 @@ ReferenceCustomHbondIxn::ReferenceCustomHbondIxn(const vector<vector<int> >& don
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomHbondIxn::~ReferenceCustomHbondIxn( ){ ReferenceCustomHbondIxn::~ReferenceCustomHbondIxn() {
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -128,7 +128,7 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, ...@@ -128,7 +128,7 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates,
int numDonors = donorAtoms.size(); int numDonors = donorAtoms.size();
int numAcceptors = acceptorAtoms.size(); int numAcceptors = acceptorAtoms.size();
for( int donor = 0; donor < numDonors; donor++ ){ for (int donor = 0; donor < numDonors; donor++) {
// Initialize per-donor parameters. // Initialize per-donor parameters.
for (int j = 0; j < (int) donorParamNames.size(); j++) for (int j = 0; j < (int) donorParamNames.size(); j++)
...@@ -136,7 +136,7 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, ...@@ -136,7 +136,7 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates,
// loop over atom pairs // loop over atom pairs
for( int acceptor = 0; acceptor < numAcceptors; acceptor++ ){ for (int acceptor = 0; acceptor < numAcceptors; acceptor++) {
if (exclusions[donor].find(acceptor) == exclusions[donor].end()) { if (exclusions[donor].find(acceptor) == exclusions[donor].end()) {
for (int j = 0; j < (int) acceptorParamNames.size(); j++) for (int j = 0; j < (int) acceptorParamNames.size(); j++)
variables[acceptorParamNames[j]] = acceptorParameters[acceptor][j]; variables[acceptorParamNames[j]] = acceptorParameters[acceptor][j];
......
platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
View file @ 7fb10336
...@@ -103,7 +103,7 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP ...@@ -103,7 +103,7 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP
CustomManyParticleForceImpl::buildFilterArrays(force, numTypes, particleTypes, orderIndex, particleOrder); CustomManyParticleForceImpl::buildFilterArrays(force, numTypes, particleTypes, orderIndex, particleOrder);
} }
ReferenceCustomManyParticleIxn::~ReferenceCustomManyParticleIxn( ){ ReferenceCustomManyParticleIxn::~ReferenceCustomManyParticleIxn() {
} }
void ReferenceCustomManyParticleIxn::calculateIxn(vector<RealVec>& atomCoordinates, RealOpenMM** particleParameters, void ReferenceCustomManyParticleIxn::calculateIxn(vector<RealVec>& atomCoordinates, RealOpenMM** particleParameters,
......
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @ 7fb10336
...@@ -72,7 +72,7 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledE ...@@ -72,7 +72,7 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledE
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -91,7 +91,7 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){ ...@@ -91,7 +91,7 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors ) { void ReferenceCustomNonbondedIxn::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors) {
cutoff = true; cutoff = true;
cutoffDistance = distance; cutoffDistance = distance;
...@@ -119,7 +119,7 @@ void ReferenceCustomNonbondedIxn::setInteractionGroups(const vector<pair<set<int ...@@ -119,7 +119,7 @@ void ReferenceCustomNonbondedIxn::setInteractionGroups(const vector<pair<set<int
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance ) { void ReferenceCustomNonbondedIxn::setUseSwitchingFunction(RealOpenMM distance) {
useSwitch = true; useSwitch = true;
switchingDistance = distance; switchingDistance = distance;
} }
...@@ -165,10 +165,10 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance ) ...@@ -165,10 +165,10 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates, void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<set<int> >& exclusions, RealOpenMM** atomParameters, vector<set<int> >& exclusions,
RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces, RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) { RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) {
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) { for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(energyExpression, iter->first), iter->second); ReferenceForce::setVariable(ReferenceForce::getVariablePointer(energyExpression, iter->first), iter->second);
...@@ -244,8 +244,8 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re ...@@ -244,8 +244,8 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates, vector<RealVec>& forces, void ReferenceCustomNonbondedIxn::calculateOneIxn(int ii, int jj, vector<RealVec>& atomCoordinates, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) { RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -267,9 +267,9 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe ...@@ -267,9 +267,9 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex]; RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic) if (periodic)
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR ); ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR);
else else
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR ); ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex]; RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (cutoff && r >= cutoffDistance) if (cutoff && r >= cutoffDistance)
return; return;
...@@ -289,7 +289,7 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe ...@@ -289,7 +289,7 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
energy *= switchValue; energy *= switchValue;
} }
} }
for( int kk = 0; kk < 3; kk++ ){ for (int kk = 0; kk < 3; kk++) {
RealOpenMM force = -dEdR*deltaR[kk]; RealOpenMM force = -dEdR*deltaR[kk];
forces[ii][kk] += force; forces[ii][kk] += force;
forces[jj][kk] -= force; forces[jj][kk] -= force;
...@@ -297,10 +297,10 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe ...@@ -297,10 +297,10 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
// accumulate energies // accumulate energies
if( totalEnergy || energyByAtom ) { if (totalEnergy || energyByAtom) {
if( totalEnergy ) if (totalEnergy)
*totalEnergy += energy; *totalEnergy += energy;
if( energyByAtom ){ if (energyByAtom) {
energyByAtom[ii] += energy; energyByAtom[ii] += energy;
energyByAtom[jj] += energy; energyByAtom[jj] += energy;
} }
......
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp
View file @ 7fb10336
...@@ -61,7 +61,7 @@ ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpre ...@@ -61,7 +61,7 @@ ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpre
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn( ){ ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -83,11 +83,11 @@ ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn( ){ ...@@ -83,11 +83,11 @@ ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices, void ReferenceCustomTorsionIxn::calculateBondIxn(int* atomIndices,
vector<RealVec>& atomCoordinates, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, RealOpenMM* parameters,
vector<RealVec>& forces, vector<RealVec>& forces,
RealOpenMM* totalEnergy ) const { RealOpenMM* totalEnergy) const {
static const std::string methodName = "\nReferenceCustomTorsionIxn::calculateTorsionIxn"; static const std::string methodName = "\nReferenceCustomTorsionIxn::calculateTorsionIxn";
...@@ -135,20 +135,20 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices, ...@@ -135,20 +135,20 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
RealOpenMM internalF[4][3]; RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4]; RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3( crossProduct[0], crossProduct[0] ); RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex]; RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1; forceFactors[0] = (-dEdAngle*normBC)/normCross1;
RealOpenMM normCross2 = DOT3( crossProduct[1], crossProduct[1] ); RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2; forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3( deltaR[0], deltaR[1] ); forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index]; forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3( deltaR[2], deltaR[1] ); forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index]; forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii]; internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii]; internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
...@@ -161,7 +161,7 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices, ...@@ -161,7 +161,7 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
// accumulate forces // accumulate forces
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forces[atomAIndex][ii] += internalF[0][ii]; forces[atomAIndex][ii] += internalF[0][ii];
forces[atomBIndex][ii] -= internalF[1][ii]; forces[atomBIndex][ii] -= internalF[1][ii];
forces[atomCIndex][ii] -= internalF[2][ii]; forces[atomCIndex][ii] -= internalF[2][ii];
......
platforms/reference/src/SimTKReference/ReferenceDynamics.cpp
View file @ 7fb10336
...@@ -45,7 +45,7 @@ using namespace OpenMM; ...@@ -45,7 +45,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceDynamics::ReferenceDynamics( int numberOfAtoms, RealOpenMM deltaT, RealOpenMM temperature ) : ReferenceDynamics::ReferenceDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM temperature) :
_numberOfAtoms(numberOfAtoms), _deltaT(deltaT), _temperature(temperature) { _numberOfAtoms(numberOfAtoms), _deltaT(deltaT), _temperature(temperature) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -68,7 +68,7 @@ ReferenceDynamics::ReferenceDynamics( int numberOfAtoms, RealOpenMM deltaT, Rea ...@@ -68,7 +68,7 @@ ReferenceDynamics::ReferenceDynamics( int numberOfAtoms, RealOpenMM deltaT, Rea
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceDynamics::~ReferenceDynamics( ){ ReferenceDynamics::~ReferenceDynamics() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -76,7 +76,7 @@ ReferenceDynamics::~ReferenceDynamics( ){ ...@@ -76,7 +76,7 @@ ReferenceDynamics::~ReferenceDynamics( ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
if( _ownReferenceConstraint ){ if (_ownReferenceConstraint) {
delete _referenceConstraint; delete _referenceConstraint;
} }
} }
...@@ -89,7 +89,7 @@ ReferenceDynamics::~ReferenceDynamics( ){ ...@@ -89,7 +89,7 @@ ReferenceDynamics::~ReferenceDynamics( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
int ReferenceDynamics::getNumberOfAtoms( void ) const { int ReferenceDynamics::getNumberOfAtoms() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -108,7 +108,7 @@ int ReferenceDynamics::getNumberOfAtoms( void ) const { ...@@ -108,7 +108,7 @@ int ReferenceDynamics::getNumberOfAtoms( void ) const {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
int ReferenceDynamics::getTimeStep( void ) const { int ReferenceDynamics::getTimeStep() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -127,7 +127,7 @@ int ReferenceDynamics::getTimeStep( void ) const { ...@@ -127,7 +127,7 @@ int ReferenceDynamics::getTimeStep( void ) const {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
int ReferenceDynamics::incrementTimeStep( void ){ int ReferenceDynamics::incrementTimeStep() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -146,7 +146,7 @@ int ReferenceDynamics::incrementTimeStep( void ){ ...@@ -146,7 +146,7 @@ int ReferenceDynamics::incrementTimeStep( void ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM ReferenceDynamics::getDeltaT( void ) const { RealOpenMM ReferenceDynamics::getDeltaT() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -163,7 +163,7 @@ RealOpenMM ReferenceDynamics::getDeltaT( void ) const { ...@@ -163,7 +163,7 @@ RealOpenMM ReferenceDynamics::getDeltaT( void ) const {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceDynamics::setDeltaT( RealOpenMM deltaT ) { void ReferenceDynamics::setDeltaT(RealOpenMM deltaT) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -182,7 +182,7 @@ void ReferenceDynamics::setDeltaT( RealOpenMM deltaT ) { ...@@ -182,7 +182,7 @@ void ReferenceDynamics::setDeltaT( RealOpenMM deltaT ) {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM ReferenceDynamics::getTemperature( void ) const { RealOpenMM ReferenceDynamics::getTemperature() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -201,7 +201,7 @@ RealOpenMM ReferenceDynamics::getTemperature( void ) const { ...@@ -201,7 +201,7 @@ RealOpenMM ReferenceDynamics::getTemperature( void ) const {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceConstraintAlgorithm* ReferenceDynamics::getReferenceConstraintAlgorithm( void ) const { ReferenceConstraintAlgorithm* ReferenceDynamics::getReferenceConstraintAlgorithm() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -220,7 +220,7 @@ ReferenceConstraintAlgorithm* ReferenceDynamics::getReferenceConstraintAlgorithm ...@@ -220,7 +220,7 @@ ReferenceConstraintAlgorithm* ReferenceDynamics::getReferenceConstraintAlgorithm
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgorithm* referenceConstraint ){ void ReferenceDynamics::setReferenceConstraintAlgorithm(ReferenceConstraintAlgorithm* referenceConstraint) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -230,7 +230,7 @@ void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgo ...@@ -230,7 +230,7 @@ void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgo
// delete if own // delete if own
if( _referenceConstraint && _ownReferenceConstraint ){ if (_referenceConstraint && _ownReferenceConstraint) {
delete _referenceConstraint; delete _referenceConstraint;
} }
......
platforms/reference/src/SimTKReference/ReferenceForce.cpp
View file @ 7fb10336
...@@ -39,7 +39,7 @@ using namespace OpenMM; ...@@ -39,7 +39,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceForce::ReferenceForce( ){ ReferenceForce::ReferenceForce() {
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -48,7 +48,7 @@ ReferenceForce::ReferenceForce( ){ ...@@ -48,7 +48,7 @@ ReferenceForce::ReferenceForce( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceForce::~ReferenceForce( ){ ReferenceForce::~ReferenceForce() {
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -64,28 +64,28 @@ RealOpenMM ReferenceForce::periodicDifference(RealOpenMM val1, RealOpenMM val2, ...@@ -64,28 +64,28 @@ RealOpenMM ReferenceForce::periodicDifference(RealOpenMM val1, RealOpenMM val2,
} }
void ReferenceForce::getDeltaR( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ, void ReferenceForce::getDeltaR(const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
RealOpenMM* deltaR ){ RealOpenMM* deltaR) {
deltaR[XIndex] = atomCoordinatesJ[0] - atomCoordinatesI[0]; deltaR[XIndex] = atomCoordinatesJ[0] - atomCoordinatesI[0];
deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1]; deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1];
deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2]; deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2];
deltaR[R2Index] = DOT3( deltaR, deltaR ); deltaR[R2Index] = DOT3(deltaR, deltaR);
deltaR[RIndex] = (RealOpenMM) SQRT( deltaR[R2Index] ); deltaR[RIndex] = (RealOpenMM) SQRT(deltaR[R2Index]);
} }
void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ, void ReferenceForce::getDeltaRPeriodic(const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
const RealOpenMM* boxSize, RealOpenMM* deltaR ){ const RealOpenMM* boxSize, RealOpenMM* deltaR) {
deltaR[XIndex] = periodicDifference(atomCoordinatesJ[0], atomCoordinatesI[0], boxSize[0]); deltaR[XIndex] = periodicDifference(atomCoordinatesJ[0], atomCoordinatesI[0], boxSize[0]);
deltaR[YIndex] = periodicDifference(atomCoordinatesJ[1], atomCoordinatesI[1], boxSize[1]); deltaR[YIndex] = periodicDifference(atomCoordinatesJ[1], atomCoordinatesI[1], boxSize[1]);
deltaR[ZIndex] = periodicDifference(atomCoordinatesJ[2], atomCoordinatesI[2], boxSize[2]); deltaR[ZIndex] = periodicDifference(atomCoordinatesJ[2], atomCoordinatesI[2], boxSize[2]);
deltaR[R2Index] = DOT3( deltaR, deltaR ); deltaR[R2Index] = DOT3(deltaR, deltaR);
deltaR[RIndex] = (RealOpenMM) SQRT( deltaR[R2Index] ); deltaR[RIndex] = (RealOpenMM) SQRT(deltaR[R2Index]);
} }
void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ, void ReferenceForce::getDeltaRPeriodic(const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
const RealVec* boxVectors, RealOpenMM* deltaR ){ const RealVec* boxVectors, RealOpenMM* deltaR) {
RealVec diff = atomCoordinatesJ-atomCoordinatesI; RealVec diff = atomCoordinatesJ-atomCoordinatesI;
diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5); diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5);
diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5); diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5);
......
platforms/reference/src/SimTKReference/ReferenceHarmonicBondIxn.cpp
View file @ 7fb10336
...@@ -39,7 +39,7 @@ using namespace OpenMM; ...@@ -39,7 +39,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn( ){ ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -55,7 +55,7 @@ ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn( ){ ...@@ -55,7 +55,7 @@ ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn( ){ ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -78,11 +78,11 @@ ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn( ){ ...@@ -78,11 +78,11 @@ ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceHarmonicBondIxn::calculateBondIxn( int* atomIndices, void ReferenceHarmonicBondIxn::calculateBondIxn(int* atomIndices,
vector<RealVec>& atomCoordinates, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, RealOpenMM* parameters,
vector<RealVec>& forces, vector<RealVec>& forces,
RealOpenMM* totalEnergy ) const { RealOpenMM* totalEnergy) const {
static const std::string methodName = "\nReferenceHarmonicBondIxn::calculateBondIxn"; static const std::string methodName = "\nReferenceHarmonicBondIxn::calculateBondIxn";
...@@ -100,7 +100,7 @@ void ReferenceHarmonicBondIxn::calculateBondIxn( int* atomIndices, ...@@ -100,7 +100,7 @@ void ReferenceHarmonicBondIxn::calculateBondIxn( int* atomIndices,
int atomAIndex = atomIndices[0]; int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1]; int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR( atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR ); ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);
// deltaIdeal = r - r_0 // deltaIdeal = r - r_0
......
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp
View file @ 7fb10336
...@@ -39,7 +39,7 @@ using namespace OpenMM; ...@@ -39,7 +39,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceLJCoulomb14::ReferenceLJCoulomb14( ) { ReferenceLJCoulomb14::ReferenceLJCoulomb14() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -55,7 +55,7 @@ ReferenceLJCoulomb14::ReferenceLJCoulomb14( ) { ...@@ -55,7 +55,7 @@ ReferenceLJCoulomb14::ReferenceLJCoulomb14( ) {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){ ReferenceLJCoulomb14::~ReferenceLJCoulomb14() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -80,9 +80,9 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){ ...@@ -80,9 +80,9 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>& atomCoordinates, void ReferenceLJCoulomb14::calculateBondIxn(int* atomIndices, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, vector<RealVec>& forces, RealOpenMM* parameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy ) const { RealOpenMM* totalEnergy) const {
static const std::string methodName = "\nReferenceLJCoulomb14::calculateBondIxn"; static const std::string methodName = "\nReferenceLJCoulomb14::calculateBondIxn";
...@@ -112,7 +112,7 @@ void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>& ...@@ -112,7 +112,7 @@ void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>&
int atomAIndex = atomIndices[0]; int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1]; int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] ); ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index]; RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]); RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
...@@ -120,13 +120,13 @@ void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>& ...@@ -120,13 +120,13 @@ void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>&
sig2 *= sig2; sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2; RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM dEdR = parameters[1]*( twelve*sig6 - six )*sig6; RealOpenMM dEdR = parameters[1]*(twelve*sig6 - six)*sig6;
dEdR += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR); dEdR += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR);
dEdR *= inverseR*inverseR; dEdR *= inverseR*inverseR;
// accumulate forces // accumulate forces
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
RealOpenMM force = dEdR*deltaR[0][ii]; RealOpenMM force = dEdR*deltaR[0][ii];
forces[atomAIndex][ii] += force; forces[atomAIndex][ii] += force;
forces[atomBIndex][ii] -= force; forces[atomBIndex][ii] -= force;
...@@ -135,5 +135,5 @@ void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>& ...@@ -135,5 +135,5 @@ void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, vector<RealVec>&
// accumulate energies // accumulate energies
if (totalEnergy != NULL) if (totalEnergy != NULL)
*totalEnergy += parameters[1]*( sig6 - one )*sig6 + (ONE_4PI_EPS0*parameters[2]*inverseR); *totalEnergy += parameters[1]*(sig6 - one)*sig6 + (ONE_4PI_EPS0*parameters[2]*inverseR);
} }
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ 7fb10336
...@@ -48,7 +48,7 @@ using namespace OpenMM; ...@@ -48,7 +48,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false) { ReferenceLJCoulombIxn::ReferenceLJCoulombIxn() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -64,7 +64,7 @@ ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), useSwitch(false ...@@ -64,7 +64,7 @@ ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), useSwitch(false
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -84,7 +84,7 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){ ...@@ -84,7 +84,7 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceLJCoulombIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) { void ReferenceLJCoulombIxn::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric) {
cutoff = true; cutoff = true;
cutoffDistance = distance; cutoffDistance = distance;
...@@ -101,7 +101,7 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){ ...@@ -101,7 +101,7 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceLJCoulombIxn::setUseSwitchingFunction( RealOpenMM distance ) { void ReferenceLJCoulombIxn::setUseSwitchingFunction(RealOpenMM distance) {
useSwitch = true; useSwitch = true;
switchingDistance = distance; switchingDistance = distance;
} }
...@@ -210,16 +210,16 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -210,16 +210,16 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
// ************************************************************************************** // **************************************************************************************
if (includeReciprocal) { if (includeReciprocal) {
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){ for (int atomID = 0; atomID < numberOfAtoms; atomID++) {
RealOpenMM selfEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI); RealOpenMM selfEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI);
totalSelfEwaldEnergy -= selfEwaldEnergy; totalSelfEwaldEnergy -= selfEwaldEnergy;
if( energyByAtom ){ if (energyByAtom) {
energyByAtom[atomID] -= selfEwaldEnergy; energyByAtom[atomID] -= selfEwaldEnergy;
} }
} }
} }
if( totalEnergy ){ if (totalEnergy) {
*totalEnergy += totalSelfEwaldEnergy; *totalEnergy += totalSelfEwaldEnergy;
} }
...@@ -238,11 +238,11 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -238,11 +238,11 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
charges[i] = atomParameters[i][QIndex]; charges[i] = atomParameters[i][QIndex];
pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxVectors,&recipEnergy); pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxVectors,&recipEnergy);
if( totalEnergy ) if (totalEnergy)
*totalEnergy += recipEnergy; *totalEnergy += recipEnergy;
if( energyByAtom ) if (energyByAtom)
for(int n = 0; n < numberOfAtoms; n++) for (int n = 0; n < numberOfAtoms; n++)
energyByAtom[n] += recipEnergy; energyByAtom[n] += recipEnergy;
pme_destroy(pmedata); pme_destroy(pmedata);
...@@ -267,16 +267,16 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -267,16 +267,16 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
if (kmax < 1) if (kmax < 1)
throw OpenMMException("kmax for Ewald summation < 1"); throw OpenMMException("kmax for Ewald summation < 1");
for(int i = 0; (i < numberOfAtoms); i++) { for (int i = 0; (i < numberOfAtoms); i++) {
for(int m = 0; (m < 3); m++) for (int m = 0; (m < 3); m++)
EIR(0, i, m) = d_complex(1,0); EIR(0, i, m) = d_complex(1,0);
for(int m=0; (m<3); m++) for (int m=0; (m<3); m++)
EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]), EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
sin(atomCoordinates[i][m]*recipBoxSize[m])); sin(atomCoordinates[i][m]*recipBoxSize[m]));
for(int j=2; (j<kmax); j++) for (int j=2; (j<kmax); j++)
for(int m=0; (m<3); m++) for (int m=0; (m<3); m++)
EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m); EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
} }
...@@ -285,40 +285,40 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -285,40 +285,40 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
int lowry = 0; int lowry = 0;
int lowrz = 1; int lowrz = 1;
for(int rx = 0; rx < numRx; rx++) { for (int rx = 0; rx < numRx; rx++) {
RealOpenMM kx = rx * recipBoxSize[0]; RealOpenMM kx = rx * recipBoxSize[0];
for(int ry = lowry; ry < numRy; ry++) { for (int ry = lowry; ry < numRy; ry++) {
RealOpenMM ky = ry * recipBoxSize[1]; RealOpenMM ky = ry * recipBoxSize[1];
if(ry >= 0) { if (ry >= 0) {
for(int n = 0; n < numberOfAtoms; n++) for (int n = 0; n < numberOfAtoms; n++)
tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1); tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
} }
else { else {
for(int n = 0; n < numberOfAtoms; n++) for (int n = 0; n < numberOfAtoms; n++)
tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1)); tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
} }
for (int rz = lowrz; rz < numRz; rz++) { for (int rz = lowrz; rz < numRz; rz++) {
if( rz >= 0) { if (rz >= 0) {
for( int n = 0; n < numberOfAtoms; n++) for (int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2)); tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
} }
else { else {
for( int n = 0; n < numberOfAtoms; n++) for (int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2))); tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
} }
RealOpenMM cs = 0.0f; RealOpenMM cs = 0.0f;
RealOpenMM ss = 0.0f; RealOpenMM ss = 0.0f;
for( int n = 0; n < numberOfAtoms; n++) { for (int n = 0; n < numberOfAtoms; n++) {
cs += tab_qxyz[n].real(); cs += tab_qxyz[n].real();
ss += tab_qxyz[n].imag(); ss += tab_qxyz[n].imag();
} }
...@@ -327,21 +327,21 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -327,21 +327,21 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
RealOpenMM k2 = kx * kx + ky * ky + kz * kz; RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
RealOpenMM ak = exp(k2*factorEwald) / k2; RealOpenMM ak = exp(k2*factorEwald) / k2;
for(int n = 0; n < numberOfAtoms; n++) { for (int n = 0; n < numberOfAtoms; n++) {
RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real()); RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
forces[n][0] += 2 * recipCoeff * force * kx ; forces[n][0] += 2 * recipCoeff * force * kx ;
forces[n][1] += 2 * recipCoeff * force * ky ; forces[n][1] += 2 * recipCoeff * force * ky ;
forces[n][2] += 2 * recipCoeff * force * kz ; forces[n][2] += 2 * recipCoeff * force * kz ;
} }
recipEnergy = recipCoeff * ak * ( cs * cs + ss * ss); recipEnergy = recipCoeff * ak * (cs * cs + ss * ss);
totalRecipEnergy += recipEnergy; totalRecipEnergy += recipEnergy;
if( totalEnergy ) if (totalEnergy)
*totalEnergy += recipEnergy; *totalEnergy += recipEnergy;
if( energyByAtom ) if (energyByAtom)
for(int n = 0; n < numberOfAtoms; n++) for (int n = 0; n < numberOfAtoms; n++)
energyByAtom[n] += recipEnergy; energyByAtom[n] += recipEnergy;
lowrz = 1 - numRz; lowrz = 1 - numRz;
...@@ -366,7 +366,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -366,7 +366,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
int jj = pair.second; int jj = pair.second;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex]; RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0] ); ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
RealOpenMM r = deltaR[0][ReferenceForce::RIndex]; RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]); RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM switchValue = 1, switchDeriv = 0; RealOpenMM switchValue = 1, switchDeriv = 0;
...@@ -379,14 +379,14 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -379,14 +379,14 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR); RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
dEdR = (RealOpenMM) (dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI )); dEdR = (RealOpenMM) (dEdR * (erfc(alphaR) + 2 * alphaR * exp (- alphaR * alphaR) / SQRT_PI));
RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex]; RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex];
RealOpenMM sig2 = inverseR*sig; RealOpenMM sig2 = inverseR*sig;
sig2 *= sig2; sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2; RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex]; RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
dEdR += switchValue*eps*( twelve*sig6 - six )*sig6*inverseR*inverseR; dEdR += switchValue*eps*(twelve*sig6 - six)*sig6*inverseR*inverseR;
vdwEnergy = eps*(sig6-one)*sig6; vdwEnergy = eps*(sig6-one)*sig6;
if (useSwitch) { if (useSwitch) {
dEdR -= vdwEnergy*switchDeriv*inverseR; dEdR -= vdwEnergy*switchDeriv*inverseR;
...@@ -395,7 +395,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -395,7 +395,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
// accumulate forces // accumulate forces
for( int kk = 0; kk < 3; kk++ ){ for (int kk = 0; kk < 3; kk++) {
RealOpenMM force = dEdR*deltaR[0][kk]; RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force; forces[ii][kk] += force;
forces[jj][kk] -= force; forces[jj][kk] -= force;
...@@ -408,14 +408,14 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -408,14 +408,14 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
totalVdwEnergy += vdwEnergy; totalVdwEnergy += vdwEnergy;
totalRealSpaceEwaldEnergy += realSpaceEwaldEnergy; totalRealSpaceEwaldEnergy += realSpaceEwaldEnergy;
if( energyByAtom ){ if (energyByAtom) {
energyByAtom[ii] += realSpaceEwaldEnergy + vdwEnergy; energyByAtom[ii] += realSpaceEwaldEnergy + vdwEnergy;
energyByAtom[jj] += realSpaceEwaldEnergy + vdwEnergy; energyByAtom[jj] += realSpaceEwaldEnergy + vdwEnergy;
} }
} }
if( totalEnergy ) if (totalEnergy)
*totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy; *totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy;
// Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum. // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
...@@ -428,17 +428,17 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -428,17 +428,17 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
int jj = *iter; int jj = *iter;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex]; RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] ); ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
RealOpenMM r = deltaR[0][ReferenceForce::RIndex]; RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]); RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r; RealOpenMM alphaR = alphaEwald * r;
if (erf(alphaR) > 1e-6) { if (erf(alphaR) > 1e-6) {
RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR); RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
dEdR = (RealOpenMM) (dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI )); dEdR = (RealOpenMM) (dEdR * (erf(alphaR) - 2 * alphaR * exp (- alphaR * alphaR) / SQRT_PI));
// accumulate forces // accumulate forces
for( int kk = 0; kk < 3; kk++ ){ for (int kk = 0; kk < 3; kk++) {
RealOpenMM force = dEdR*deltaR[0][kk]; RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] -= force; forces[ii][kk] -= force;
forces[jj][kk] += force; forces[jj][kk] += force;
...@@ -449,7 +449,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -449,7 +449,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
realSpaceEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR)); realSpaceEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR));
totalExclusionEnergy += realSpaceEwaldEnergy; totalExclusionEnergy += realSpaceEwaldEnergy;
if( energyByAtom ){ if (energyByAtom) {
energyByAtom[ii] -= realSpaceEwaldEnergy; energyByAtom[ii] -= realSpaceEwaldEnergy;
energyByAtom[jj] -= realSpaceEwaldEnergy; energyByAtom[jj] -= realSpaceEwaldEnergy;
} }
...@@ -457,7 +457,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec> ...@@ -457,7 +457,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
} }
} }
if( totalEnergy ) if (totalEnergy)
*totalEnergy -= totalExclusionEnergy; *totalEnergy -= totalExclusionEnergy;
} }
...@@ -499,10 +499,10 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>& ...@@ -499,10 +499,10 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>&
} }
} }
else { else {
for( int ii = 0; ii < numberOfAtoms; ii++ ){ for (int ii = 0; ii < numberOfAtoms; ii++) {
// loop over atom pairs // loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ) for (int jj = ii+1; jj < numberOfAtoms; jj++)
if (exclusions[jj].find(ii) == exclusions[jj].end()) if (exclusions[jj].find(ii) == exclusions[jj].end())
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy); calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
} }
...@@ -523,9 +523,9 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>& ...@@ -523,9 +523,9 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>&
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates, void ReferenceLJCoulombIxn::calculateOneIxn(int ii, int jj, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM** atomParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const { RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -552,9 +552,9 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at ...@@ -552,9 +552,9 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
// get deltaR, R2, and R between 2 atoms // get deltaR, R2, and R between 2 atoms
if (periodic) if (periodic)
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0] ); ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
else else
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] ); ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index]; RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]); RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
...@@ -573,7 +573,7 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at ...@@ -573,7 +573,7 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
RealOpenMM sig6 = sig2*sig2*sig2; RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex]; RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
RealOpenMM dEdR = switchValue*eps*( twelve*sig6 - six )*sig6; RealOpenMM dEdR = switchValue*eps*(twelve*sig6 - six)*sig6;
if (cutoff) if (cutoff)
dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2)); dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2));
else else
...@@ -591,7 +591,7 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at ...@@ -591,7 +591,7 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
// accumulate forces // accumulate forces
for( int kk = 0; kk < 3; kk++ ){ for (int kk = 0; kk < 3; kk++) {
RealOpenMM force = dEdR*deltaR[0][kk]; RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force; forces[ii][kk] += force;
forces[jj][kk] -= force; forces[jj][kk] -= force;
...@@ -599,9 +599,9 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at ...@@ -599,9 +599,9 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
// accumulate energies // accumulate energies
if( totalEnergy ) if (totalEnergy)
*totalEnergy += energy; *totalEnergy += energy;
if( energyByAtom ){ if (energyByAtom) {
energyByAtom[ii] += energy; energyByAtom[ii] += energy;
energyByAtom[jj] += energy; energyByAtom[jj] += energy;
} }
......
platforms/reference/src/SimTKReference/ReferenceLincsAlgorithm.cpp
View file @ 7fb10336
...@@ -43,9 +43,9 @@ using namespace OpenMM; ...@@ -43,9 +43,9 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceLincsAlgorithm::ReferenceLincsAlgorithm( int numberOfConstraints, ReferenceLincsAlgorithm::ReferenceLincsAlgorithm(int numberOfConstraints,
int** atomIndices, int** atomIndices,
RealOpenMM* distance ){ RealOpenMM* distance) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -69,7 +69,7 @@ ReferenceLincsAlgorithm::ReferenceLincsAlgorithm( int numberOfConstraints, ...@@ -69,7 +69,7 @@ ReferenceLincsAlgorithm::ReferenceLincsAlgorithm( int numberOfConstraints,
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
int ReferenceLincsAlgorithm::getNumberOfConstraints( void ) const { int ReferenceLincsAlgorithm::getNumberOfConstraints() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -88,7 +88,7 @@ int ReferenceLincsAlgorithm::getNumberOfConstraints( void ) const { ...@@ -88,7 +88,7 @@ int ReferenceLincsAlgorithm::getNumberOfConstraints( void ) const {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
int ReferenceLincsAlgorithm::getNumTerms( void ) const { int ReferenceLincsAlgorithm::getNumTerms() const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -105,7 +105,7 @@ int ReferenceLincsAlgorithm::getNumTerms( void ) const { ...@@ -105,7 +105,7 @@ int ReferenceLincsAlgorithm::getNumTerms( void ) const {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceLincsAlgorithm::setNumTerms( int terms ){ void ReferenceLincsAlgorithm::setNumTerms(int terms) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -218,9 +218,9 @@ void ReferenceLincsAlgorithm::updateAtomPositions(int numberOfAtoms, vector<Real ...@@ -218,9 +218,9 @@ void ReferenceLincsAlgorithm::updateAtomPositions(int numberOfAtoms, vector<Real
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
int ReferenceLincsAlgorithm::apply( int numberOfAtoms, vector<RealVec>& atomCoordinates, int ReferenceLincsAlgorithm::apply(int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& atomCoordinatesP, vector<RealVec>& atomCoordinatesP,
vector<RealOpenMM>& inverseMasses ){ vector<RealOpenMM>& inverseMasses) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -235,7 +235,7 @@ int ReferenceLincsAlgorithm::apply( int numberOfAtoms, vector<RealVec>& atomCoor ...@@ -235,7 +235,7 @@ int ReferenceLincsAlgorithm::apply( int numberOfAtoms, vector<RealVec>& atomCoor
if (_numberOfConstraints == 0) if (_numberOfConstraints == 0)
return SimTKOpenMMCommon::DefaultReturn; return SimTKOpenMMCommon::DefaultReturn;
if( !_hasInitialized ) if (!_hasInitialized)
initialize(numberOfAtoms, inverseMasses); initialize(numberOfAtoms, inverseMasses);
// Calculate the direction of each constraint, along with the initial RHS and solution vectors. // Calculate the direction of each constraint, along with the initial RHS and solution vectors.
......
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
View file @ 7fb10336
...@@ -49,7 +49,7 @@ ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vec ...@@ -49,7 +49,7 @@ ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vec
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceMonteCarloBarostat::~ReferenceMonteCarloBarostat( ) { ReferenceMonteCarloBarostat::~ReferenceMonteCarloBarostat() {
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp
View file @ 7fb10336
...@@ -45,12 +45,12 @@ void OPENMM_EXPORT computeNeighborListNaive( ...@@ -45,12 +45,12 @@ void OPENMM_EXPORT computeNeighborListNaive(
for (AtomIndex atomJ = atomI + 1; atomJ < (AtomIndex) nAtoms; ++atomJ) for (AtomIndex atomJ = atomI + 1; atomJ < (AtomIndex) nAtoms; ++atomJ)
{ {
double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxVectors, usePeriodic); double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxVectors, usePeriodic);
if ( (pairDistanceSquared <= maxDistanceSquared) && (pairDistanceSquared >= minDistanceSquared)) if ((pairDistanceSquared <= maxDistanceSquared) && (pairDistanceSquared >= minDistanceSquared))
if (exclusions[atomI].find(atomJ) == exclusions[atomI].end()) if (exclusions[atomI].find(atomJ) == exclusions[atomI].end())
{ {
neighborList.push_back( AtomPair(atomI, atomJ) ); neighborList.push_back(AtomPair(atomI, atomJ));
if (reportSymmetricPairs) if (reportSymmetricPairs)
neighborList.push_back( AtomPair(atomI, atomJ) ); neighborList.push_back(AtomPair(atomI, atomJ));
} }
} }
} }
...@@ -199,9 +199,9 @@ public: ...@@ -199,9 +199,9 @@ public:
// Ignore exclusions. // Ignore exclusions.
if (exclusions[atomI].find(atomJ) != exclusions[atomI].end()) continue; if (exclusions[atomI].find(atomJ) != exclusions[atomI].end()) continue;
neighbors.push_back( AtomPair(atomI, atomJ) ); neighbors.push_back(AtomPair(atomI, atomJ));
if (reportSymmetricPairs) if (reportSymmetricPairs)
neighbors.push_back( AtomPair(atomJ, atomI) ); neighbors.push_back(AtomPair(atomJ, atomI));
} }
} }
} }
......
platforms/reference/src/SimTKReference/ReferencePME.cpp
View file @ 7fb10336
...@@ -77,7 +77,7 @@ struct pme ...@@ -77,7 +77,7 @@ struct pme
* If particle i has coordinates { 0.543 , 6.235 , -0.73 }, we will get: * If particle i has coordinates { 0.543 , 6.235 , -0.73 }, we will get:
* *
* particleindex[i] = { 5 , 62 , 92 } (-0.73 + 10 = 9.27, we always apply PBC for grid calculations!) * particleindex[i] = { 5 , 62 , 92 } (-0.73 + 10 = 9.27, we always apply PBC for grid calculations!)
* particlefraction[i] = { 0.43 , 0.35 , 0.7 } ( this is the fraction of the cell length where the atom is) * particlefraction[i] = { 0.43 , 0.35 , 0.7 } (this is the fraction of the cell length where the atom is)
* *
* (The reason for precaculating / storing these is that it gets a bit more complex for triclinic cells :-) * (The reason for precaculating / storing these is that it gets a bit more complex for triclinic cells :-)
* *
...@@ -108,7 +108,7 @@ pme_calculate_bsplines_moduli(pme_t pme) ...@@ -108,7 +108,7 @@ pme_calculate_bsplines_moduli(pme_t pme)
RealOpenMM sc,ss,arg; RealOpenMM sc,ss,arg;
nmax = 0; nmax = 0;
for(d=0;d<3;d++) for (d=0;d<3;d++)
{ {
nmax = (pme->ngrid[d] > nmax) ? pme->ngrid[d] : nmax; nmax = (pme->ngrid[d] > nmax) ? pme->ngrid[d] : nmax;
pme->bsplines_moduli[d] = (RealOpenMM *) malloc(sizeof(RealOpenMM)*pme->ngrid[d]); pme->bsplines_moduli[d] = (RealOpenMM *) malloc(sizeof(RealOpenMM)*pme->ngrid[d]);
...@@ -125,11 +125,11 @@ pme_calculate_bsplines_moduli(pme_t pme) ...@@ -125,11 +125,11 @@ pme_calculate_bsplines_moduli(pme_t pme)
data[1]=0; data[1]=0;
data[0]=1; data[0]=1;
for(k=3;k<order;k++) for (k=3;k<order;k++)
{ {
div=(RealOpenMM) (1.0/(k-1.0)); div=(RealOpenMM) (1.0/(k-1.0));
data[k-1]=0; data[k-1]=0;
for(l=1;l<(k-1);l++) for (l=1;l<(k-1);l++)
{ {
data[k-l-1]=div*(l*data[k-l-2]+(k-l)*data[k-l-1]); data[k-l-1]=div*(l*data[k-l-2]+(k-l)*data[k-l-1]);
} }
...@@ -138,7 +138,7 @@ pme_calculate_bsplines_moduli(pme_t pme) ...@@ -138,7 +138,7 @@ pme_calculate_bsplines_moduli(pme_t pme)
/* differentiate */ /* differentiate */
ddata[0]=-data[0]; ddata[0]=-data[0];
for(k=1;k<order;k++) for (k=1;k<order;k++)
{ {
ddata[k]=data[k-1]-data[k]; ddata[k]=data[k-1]-data[k];
} }
...@@ -146,29 +146,29 @@ pme_calculate_bsplines_moduli(pme_t pme) ...@@ -146,29 +146,29 @@ pme_calculate_bsplines_moduli(pme_t pme)
div=(RealOpenMM) (1.0/(order-1)); div=(RealOpenMM) (1.0/(order-1));
data[order-1]=0; data[order-1]=0;
for(l=1;l<(order-1);l++) for (l=1;l<(order-1);l++)
{ {
data[order-l-1]=div*(l*data[order-l-2]+(order-l)*data[order-l-1]); data[order-l-1]=div*(l*data[order-l-2]+(order-l)*data[order-l-1]);
} }
data[0]=div*data[0]; data[0]=div*data[0];
for(i=0;i<nmax;i++) for (i=0;i<nmax;i++)
{ {
bsplines_data[i]=0; bsplines_data[i]=0;
} }
for(i=1;i<=order;i++) for (i=1;i<=order;i++)
{ {
bsplines_data[i]=data[i-1]; bsplines_data[i]=data[i-1];
} }
/* Evaluate the actual bspline moduli for X/Y/Z */ /* Evaluate the actual bspline moduli for X/Y/Z */
for(d=0;d<3;d++) for (d=0;d<3;d++)
{ {
ndata = pme->ngrid[d]; ndata = pme->ngrid[d];
for(i=0;i<ndata;i++) for (i=0;i<ndata;i++)
{ {
sc=ss=0; sc=ss=0;
for(j=0;j<ndata;j++) for (j=0;j<ndata;j++)
{ {
arg=(RealOpenMM) ((2.0*M_PI*i*j)/ndata); arg=(RealOpenMM) ((2.0*M_PI*i*j)/ndata);
sc+=bsplines_data[j]*cos(arg); sc+=bsplines_data[j]*cos(arg);
...@@ -176,9 +176,9 @@ pme_calculate_bsplines_moduli(pme_t pme) ...@@ -176,9 +176,9 @@ pme_calculate_bsplines_moduli(pme_t pme)
} }
pme->bsplines_moduli[d][i]=sc*sc+ss*ss; pme->bsplines_moduli[d][i]=sc*sc+ss*ss;
} }
for(i=0;i<ndata;i++) for (i=0;i<ndata;i++)
{ {
if(pme->bsplines_moduli[d][i]<1.0e-7) if (pme->bsplines_moduli[d][i]<1.0e-7)
{ {
pme->bsplines_moduli[d][i]=(pme->bsplines_moduli[d][i-1]+pme->bsplines_moduli[d][i+1])/2; pme->bsplines_moduli[d][i]=(pme->bsplines_moduli[d][i-1]+pme->bsplines_moduli[d][i+1])/2;
} }
...@@ -213,7 +213,7 @@ pme_update_grid_index_and_fraction(pme_t pme, ...@@ -213,7 +213,7 @@ pme_update_grid_index_and_fraction(pme_t pme,
RealOpenMM t; RealOpenMM t;
int ti; int ti;
for(i=0;i<pme->natoms;i++) for (i=0;i<pme->natoms;i++)
{ {
/* Index calculation (Look mom, no conditionals!): /* Index calculation (Look mom, no conditionals!):
* *
...@@ -252,7 +252,7 @@ pme_update_grid_index_and_fraction(pme_t pme, ...@@ -252,7 +252,7 @@ pme_update_grid_index_and_fraction(pme_t pme,
* (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!) * (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!)
*/ */
RealVec coord = atomCoordinates[i]; RealVec coord = atomCoordinates[i];
for(d=0;d<3;d++) for (d=0;d<3;d++)
{ {
t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d]; t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d];
t = (t-floor(t))*pme->ngrid[d]; t = (t-floor(t))*pme->ngrid[d];
...@@ -281,9 +281,9 @@ pme_update_bsplines(pme_t pme) ...@@ -281,9 +281,9 @@ pme_update_bsplines(pme_t pme)
order = pme->order; order = pme->order;
for(i=0; (i<pme->natoms); i++) for (i=0; (i<pme->natoms); i++)
{ {
for(j=0; j<3; j++) for (j=0; j<3; j++)
{ {
/* dr is relative offset from lower cell limit */ /* dr is relative offset from lower cell limit */
dr = pme->particlefraction[i][j]; dr = pme->particlefraction[i][j];
...@@ -294,11 +294,11 @@ pme_update_bsplines(pme_t pme) ...@@ -294,11 +294,11 @@ pme_update_bsplines(pme_t pme)
data[1] = dr; data[1] = dr;
data[0] = 1-dr; data[0] = 1-dr;
for(k=3; k<order; k++) for (k=3; k<order; k++)
{ {
div = (RealOpenMM) (1.0/(k-1.0)); div = (RealOpenMM) (1.0/(k-1.0));
data[k-1] = div*dr*data[k-2]; data[k-1] = div*dr*data[k-2];
for(l=1; l<(k-1); l++) for (l=1; l<(k-1); l++)
{ {
data[k-l-1] = div*((dr+l)*data[k-l-2]+(k-l-dr)*data[k-l-1]); data[k-l-1] = div*((dr+l)*data[k-l-2]+(k-l-dr)*data[k-l-1]);
} }
...@@ -308,7 +308,7 @@ pme_update_bsplines(pme_t pme) ...@@ -308,7 +308,7 @@ pme_update_bsplines(pme_t pme)
/* differentiate */ /* differentiate */
ddata[0] = -data[0]; ddata[0] = -data[0];
for(k=1; k<order; k++) for (k=1; k<order; k++)
{ {
ddata[k] = data[k-1]-data[k]; ddata[k] = data[k-1]-data[k];
} }
...@@ -316,7 +316,7 @@ pme_update_bsplines(pme_t pme) ...@@ -316,7 +316,7 @@ pme_update_bsplines(pme_t pme)
div = (RealOpenMM) (1.0/(order-1)); div = (RealOpenMM) (1.0/(order-1));
data[order-1] = div*dr*data[order-2]; data[order-1] = div*dr*data[order-2];
for(l=1; l<(order-1); l++) for (l=1; l<(order-1); l++)
{ {
data[order-l-1] = div*((dr+l)*data[order-l-2]+(order-l-dr)*data[order-l-1]); data[order-l-1] = div*((dr+l)*data[order-l-2]+(order-l-dr)*data[order-l-1]);
} }
...@@ -343,12 +343,12 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges) ...@@ -343,12 +343,12 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges)
order = pme->order; order = pme->order;
/* Reset the grid */ /* Reset the grid */
for(i=0;i<pme->ngrid[0]*pme->ngrid[1]*pme->ngrid[2];i++) for (i=0;i<pme->ngrid[0]*pme->ngrid[1]*pme->ngrid[2];i++)
{ {
pme->grid[i].re = pme->grid[i].im = 0; pme->grid[i].re = pme->grid[i].im = 0;
} }
for(i=0;i<pme->natoms;i++) for (i=0;i<pme->natoms;i++)
{ {
q = charges[i]; q = charges[i];
...@@ -380,16 +380,16 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges) ...@@ -380,16 +380,16 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges)
* 3) When we parallelize things, we only need to communicate in one direction instead of two! * 3) When we parallelize things, we only need to communicate in one direction instead of two!
*/ */
for(ix=0;ix<order;ix++) for (ix=0;ix<order;ix++)
{ {
/* Calculate index, apply PBC so we spread to index 0/1/2 when a particle is close to the upper limit of the grid */ /* Calculate index, apply PBC so we spread to index 0/1/2 when a particle is close to the upper limit of the grid */
xindex = (x0index + ix) % pme->ngrid[0]; xindex = (x0index + ix) % pme->ngrid[0];
for(iy=0;iy<order;iy++) for (iy=0;iy<order;iy++)
{ {
yindex = (y0index + iy) % pme->ngrid[1]; yindex = (y0index + iy) % pme->ngrid[1];
for(iz=0;iz<order;iz++) for (iz=0;iz<order;iz++)
{ {
/* Can be optimized, but we keep it simple here */ /* Can be optimized, but we keep it simple here */
zindex = (z0index + iz) % pme->ngrid[2]; zindex = (z0index + iz) % pme->ngrid[2];
...@@ -448,21 +448,21 @@ pme_reciprocal_convolution(pme_t pme, ...@@ -448,21 +448,21 @@ pme_reciprocal_convolution(pme_t pme,
maxky = (RealOpenMM) ((ny+1)/2); maxky = (RealOpenMM) ((ny+1)/2);
maxkz = (RealOpenMM) ((nz+1)/2); maxkz = (RealOpenMM) ((nz+1)/2);
for(kx=0;kx<nx;kx++) for (kx=0;kx<nx;kx++)
{ {
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */ /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mx = (RealOpenMM) ((kx<maxkx) ? kx : (kx-nx)); mx = (RealOpenMM) ((kx<maxkx) ? kx : (kx-nx));
mhx = mx*recipBoxVectors[0][0]; mhx = mx*recipBoxVectors[0][0];
bx = boxfactor*pme->bsplines_moduli[0][kx]; bx = boxfactor*pme->bsplines_moduli[0][kx];
for(ky=0;ky<ny;ky++) for (ky=0;ky<ny;ky++)
{ {
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */ /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
my = (RealOpenMM) ((ky<maxky) ? ky : (ky-ny)); my = (RealOpenMM) ((ky<maxky) ? ky : (ky-ny));
mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1]; mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
by = pme->bsplines_moduli[1][ky]; by = pme->bsplines_moduli[1][ky];
for(kz=0;kz<nz;kz++) for (kz=0;kz<nz;kz++)
{ {
/* If the net charge of the system is 0.0, there will not be any DC (direct current, zero frequency) component. However, /* If the net charge of the system is 0.0, there will not be any DC (direct current, zero frequency) component. However,
* we can still handle charged systems through a charge correction, in which case the DC * we can still handle charged systems through a charge correction, in which case the DC
...@@ -546,7 +546,7 @@ pme_grid_interpolate_force(pme_t pme, ...@@ -546,7 +546,7 @@ pme_grid_interpolate_force(pme_t pme,
/* This is almost identical to the charge spreading routine! */ /* This is almost identical to the charge spreading routine! */
for(i=0;i<pme->natoms;i++) for (i=0;i<pme->natoms;i++)
{ {
fx = fy = fz = 0; fx = fy = fz = 0;
...@@ -570,21 +570,21 @@ pme_grid_interpolate_force(pme_t pme, ...@@ -570,21 +570,21 @@ pme_grid_interpolate_force(pme_t pme,
/* Since we will add order^3 (typically 4*4*4=64) terms to the force on each particle, we use temporary fx/fy/fz /* Since we will add order^3 (typically 4*4*4=64) terms to the force on each particle, we use temporary fx/fy/fz
* variables, and only add it to memory forces[] at the end. * variables, and only add it to memory forces[] at the end.
*/ */
for(ix=0;ix<order;ix++) for (ix=0;ix<order;ix++)
{ {
xindex = (x0index + ix) % pme->ngrid[0]; xindex = (x0index + ix) % pme->ngrid[0];
/* Get both the bspline factor and its derivative with respect to the x coordinate! */ /* Get both the bspline factor and its derivative with respect to the x coordinate! */
tx = thetax[ix]; tx = thetax[ix];
dtx = dthetax[ix]; dtx = dthetax[ix];
for(iy=0;iy<order;iy++) for (iy=0;iy<order;iy++)
{ {
yindex = (y0index + iy) % pme->ngrid[1]; yindex = (y0index + iy) % pme->ngrid[1];
/* bspline + derivative wrt y */ /* bspline + derivative wrt y */
ty = thetay[iy]; ty = thetay[iy];
dty = dthetay[iy]; dty = dthetay[iy];
for(iz=0;iz<order;iz++) for (iz=0;iz<order;iz++)
{ {
/* Can be optimized, but we keep it simple here */ /* Can be optimized, but we keep it simple here */
zindex = (z0index + iz) % pme->ngrid[2]; zindex = (z0index + iz) % pme->ngrid[2];
...@@ -633,7 +633,7 @@ pme_init(pme_t * ppme, ...@@ -633,7 +633,7 @@ pme_init(pme_t * ppme,
pme->ewaldcoeff = ewaldcoeff; pme->ewaldcoeff = ewaldcoeff;
pme->natoms = natoms; pme->natoms = natoms;
for(d=0;d<3;d++) for (d=0;d<3;d++)
{ {
pme->ngrid[d] = ngrid[d]; pme->ngrid[d] = ngrid[d];
pme->bsplines_theta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms); pme->bsplines_theta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
...@@ -712,7 +712,7 @@ pme_destroy(pme_t pme) ...@@ -712,7 +712,7 @@ pme_destroy(pme_t pme)
free(pme->grid); free(pme->grid);
for(d=0;d<3;d++) for (d=0;d<3;d++)
{ {
free(pme->bsplines_moduli[d]); free(pme->bsplines_moduli[d]);
free(pme->bsplines_theta[d]); free(pme->bsplines_theta[d]);
......
platforms/reference/src/SimTKReference/ReferencePairIxn.cpp
View file @ 7fb10336
...@@ -38,7 +38,7 @@ using namespace OpenMM; ...@@ -38,7 +38,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferencePairIxn::ReferencePairIxn( ){ ReferencePairIxn::ReferencePairIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -54,7 +54,7 @@ ReferencePairIxn::ReferencePairIxn( ){ ...@@ -54,7 +54,7 @@ ReferencePairIxn::ReferencePairIxn( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferencePairIxn::~ReferencePairIxn( ){ ReferencePairIxn::~ReferencePairIxn() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceProperDihedralBond.cpp
View file @ 7fb10336
...@@ -39,7 +39,7 @@ using namespace OpenMM; ...@@ -39,7 +39,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceProperDihedralBond::ReferenceProperDihedralBond( ){ ReferenceProperDihedralBond::ReferenceProperDihedralBond() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -55,7 +55,7 @@ ReferenceProperDihedralBond::ReferenceProperDihedralBond( ){ ...@@ -55,7 +55,7 @@ ReferenceProperDihedralBond::ReferenceProperDihedralBond( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceProperDihedralBond::~ReferenceProperDihedralBond( ){ ReferenceProperDihedralBond::~ReferenceProperDihedralBond() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -79,11 +79,11 @@ ReferenceProperDihedralBond::~ReferenceProperDihedralBond( ){ ...@@ -79,11 +79,11 @@ ReferenceProperDihedralBond::~ReferenceProperDihedralBond( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices, void ReferenceProperDihedralBond::calculateBondIxn(int* atomIndices,
vector<RealVec>& atomCoordinates, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, RealOpenMM* parameters,
vector<RealVec>& forces, vector<RealVec>& forces,
RealOpenMM* totalEnergy ) const { RealOpenMM* totalEnergy) const {
static const std::string methodName = "\nReferenceProperDihedralBond::calculateBondIxn"; static const std::string methodName = "\nReferenceProperDihedralBond::calculateBondIxn";
...@@ -111,9 +111,9 @@ void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -111,9 +111,9 @@ void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices,
int atomBIndex = atomIndices[1]; int atomBIndex = atomIndices[1];
int atomCIndex = atomIndices[2]; int atomCIndex = atomIndices[2];
int atomDIndex = atomIndices[3]; int atomDIndex = atomIndices[3];
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] ); ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1] ); ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);
ReferenceForce::getDeltaR( atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2] ); ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
RealOpenMM dotDihedral; RealOpenMM dotDihedral;
RealOpenMM signOfAngle; RealOpenMM signOfAngle;
...@@ -127,35 +127,35 @@ void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -127,35 +127,35 @@ void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices,
// get dihedral angle // get dihedral angle
RealOpenMM dihedralAngle = getDihedralAngleBetweenThreeVectors( deltaR[0], deltaR[1], deltaR[2], RealOpenMM dihedralAngle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &dotDihedral, deltaR[0], crossProduct, &dotDihedral, deltaR[0],
&signOfAngle, hasREntry ); &signOfAngle, hasREntry);
// evaluate delta angle, dE/d(angle) // evaluate delta angle, dE/d(angle)
RealOpenMM deltaAngle = parameters[2]*dihedralAngle - parameters[1]; RealOpenMM deltaAngle = parameters[2]*dihedralAngle - parameters[1];
RealOpenMM sinDeltaAngle = SIN( deltaAngle ); RealOpenMM sinDeltaAngle = SIN(deltaAngle);
RealOpenMM dEdAngle = -parameters[0]*parameters[2]*sinDeltaAngle; RealOpenMM dEdAngle = -parameters[0]*parameters[2]*sinDeltaAngle;
RealOpenMM energy = parameters[0]*(one + COS( deltaAngle ) ); RealOpenMM energy = parameters[0]*(one + COS(deltaAngle));
// compute force // compute force
RealOpenMM internalF[4][3]; RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4]; RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3( crossProduct[0], crossProduct[0] ); RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex]; RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1; forceFactors[0] = (-dEdAngle*normBC)/normCross1;
RealOpenMM normCross2 = DOT3( crossProduct[1], crossProduct[1] ); RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2; forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3( deltaR[0], deltaR[1] ); forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index]; forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3( deltaR[2], deltaR[1] ); forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index]; forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii]; internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii]; internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
...@@ -168,7 +168,7 @@ void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -168,7 +168,7 @@ void ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices,
// accumulate forces // accumulate forces
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forces[atomAIndex][ii] += internalF[0][ii]; forces[atomAIndex][ii] += internalF[0][ii];
forces[atomBIndex][ii] -= internalF[1][ii]; forces[atomBIndex][ii] -= internalF[1][ii];
forces[atomCIndex][ii] -= internalF[2][ii]; forces[atomCIndex][ii] -= internalF[2][ii];
......
platforms/reference/src/SimTKReference/ReferenceRbDihedralBond.cpp
View file @ 7fb10336
...@@ -39,7 +39,7 @@ using namespace OpenMM; ...@@ -39,7 +39,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceRbDihedralBond::ReferenceRbDihedralBond( ){ ReferenceRbDihedralBond::ReferenceRbDihedralBond() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -55,7 +55,7 @@ ReferenceRbDihedralBond::ReferenceRbDihedralBond( ){ ...@@ -55,7 +55,7 @@ ReferenceRbDihedralBond::ReferenceRbDihedralBond( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceRbDihedralBond::~ReferenceRbDihedralBond( ){ ReferenceRbDihedralBond::~ReferenceRbDihedralBond() {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -77,11 +77,11 @@ ReferenceRbDihedralBond::~ReferenceRbDihedralBond( ){ ...@@ -77,11 +77,11 @@ ReferenceRbDihedralBond::~ReferenceRbDihedralBond( ){
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices, void ReferenceRbDihedralBond::calculateBondIxn(int* atomIndices,
vector<RealVec>& atomCoordinates, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, RealOpenMM* parameters,
vector<RealVec>& forces, vector<RealVec>& forces,
RealOpenMM* totalEnergy ) const { RealOpenMM* totalEnergy) const {
static const std::string methodName = "\nReferenceRbDihedralBond::calculateBondIxn"; static const std::string methodName = "\nReferenceRbDihedralBond::calculateBondIxn";
...@@ -113,9 +113,9 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -113,9 +113,9 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
int atomBIndex = atomIndices[1]; int atomBIndex = atomIndices[1];
int atomCIndex = atomIndices[2]; int atomCIndex = atomIndices[2];
int atomDIndex = atomIndices[3]; int atomDIndex = atomIndices[3];
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] ); ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1] ); ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);
ReferenceForce::getDeltaR( atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2] ); ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
RealOpenMM cosPhi; RealOpenMM cosPhi;
RealOpenMM signOfAngle; RealOpenMM signOfAngle;
...@@ -126,13 +126,13 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -126,13 +126,13 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
RealOpenMM* crossProduct[2]; RealOpenMM* crossProduct[2];
crossProduct[0] = crossProductMemory; crossProduct[0] = crossProductMemory;
crossProduct[1] = crossProductMemory + 3; crossProduct[1] = crossProductMemory + 3;
RealOpenMM dihederalAngle = getDihedralAngleBetweenThreeVectors( deltaR[0], deltaR[1], deltaR[2], RealOpenMM dihederalAngle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &cosPhi, deltaR[0], crossProduct, &cosPhi, deltaR[0],
&signOfAngle, hasREntry ); &signOfAngle, hasREntry);
// Gromacs: use polymer convention // Gromacs: use polymer convention
if( dihederalAngle < zero ){ if (dihederalAngle < zero) {
dihederalAngle += PI_M; dihederalAngle += PI_M;
} else { } else {
dihederalAngle -= PI_M; dihederalAngle -= PI_M;
...@@ -141,36 +141,36 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -141,36 +141,36 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
// Ryckaert-Bellemans: // Ryckaert-Bellemans:
// V = sum over i: { C_i*cos( psi )**i }, where psi = phi - PI, // V = sum over i: { C_i*cos(psi)**i }, where psi = phi - PI,
// C_i is ith RB coefficient // C_i is ith RB coefficient
RealOpenMM dEdAngle = zero; RealOpenMM dEdAngle = zero;
RealOpenMM energy = parameters[0]; RealOpenMM energy = parameters[0];
RealOpenMM cosFactor = one; RealOpenMM cosFactor = one;
for( int ii = 1; ii < numberOfParameters; ii++ ){ for (int ii = 1; ii < numberOfParameters; ii++) {
dEdAngle -= ((RealOpenMM) ii)*parameters[ii]*cosFactor; dEdAngle -= ((RealOpenMM) ii)*parameters[ii]*cosFactor;
cosFactor *= cosPhi; cosFactor *= cosPhi;
energy += cosFactor*parameters[ii]; energy += cosFactor*parameters[ii];
} }
dEdAngle *= SIN( dihederalAngle ); dEdAngle *= SIN(dihederalAngle);
RealOpenMM internalF[4][3]; RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4]; RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3( crossProduct[0], crossProduct[0] ); RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex]; RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1; forceFactors[0] = (-dEdAngle*normBC)/normCross1;
RealOpenMM normCross2 = DOT3( crossProduct[1], crossProduct[1] ); RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2; forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3( deltaR[0], deltaR[1] ); forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index]; forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3( deltaR[2], deltaR[1] ); forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index]; forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii]; internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii]; internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
...@@ -183,7 +183,7 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices, ...@@ -183,7 +183,7 @@ void ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
// accumulate forces // accumulate forces
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forces[atomAIndex][ii] += internalF[0][ii]; forces[atomAIndex][ii] += internalF[0][ii];
forces[atomBIndex][ii] -= internalF[1][ii]; forces[atomBIndex][ii] -= internalF[1][ii];
forces[atomCIndex][ii] -= internalF[2][ii]; forces[atomCIndex][ii] -= internalF[2][ii];
......
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