Simplified reference code for accumulating energies

platforms/reference/src/ReferenceKernels.cpp

@@ -320,21 +320,17 @@ void ReferenceCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
     RealOpenMM** forceData = extractForces(context);
     ReferenceBondForce refBondForce;
     ReferenceHarmonicBondIxn harmonicBond;
-    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, harmonicBond);
 }
 
 double ReferenceCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numBonds];
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceHarmonicBondIxn harmonicBond;
-    for (int i = 0; i < numBonds; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, &energy, harmonicBond);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -380,24 +376,20 @@ void ReferenceCalcCustomBondForceKernel::executeForces(ContextImpl& context) {
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
-    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, harmonicBond);
 }
 
 double ReferenceCalcCustomBondForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numBonds];
     RealOpenMM energy = 0;
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
-    for (int i = 0; i < numBonds; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, &energy, harmonicBond);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -427,21 +419,17 @@ void ReferenceCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context)
     RealOpenMM** forceData = extractForces(context);
     ReferenceBondForce refBondForce;
     ReferenceAngleBondIxn angleBond;
-    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, angleBond);
 }
 
 double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numAngles];
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceAngleBondIxn angleBond;
-    for (int i = 0; i < numAngles; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, angleBond);
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, &energy, angleBond);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -488,24 +476,20 @@ void ReferenceCalcCustomAngleForceKernel::executeForces(ContextImpl& context) {
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
-    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, customAngle);
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, customAngle);
 }
 
 double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numAngles];
     RealOpenMM energy = 0;
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
-    for (int i = 0; i < numAngles; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, customAngle);
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, &energy, customAngle);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -537,21 +521,17 @@ void ReferenceCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context
     RealOpenMM** forceData = extractForces(context);
     ReferenceBondForce refBondForce;
     ReferenceProperDihedralBond periodicTorsionBond;
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
+    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, periodicTorsionBond);
 }
 
 double ReferenceCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numTorsions];
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceProperDihedralBond periodicTorsionBond;
-    for (int i = 0; i < numTorsions; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, periodicTorsionBond);
+    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, &energy, periodicTorsionBond);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -586,21 +566,17 @@ void ReferenceCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
     RealOpenMM** forceData = extractForces(context);
     ReferenceBondForce refBondForce;
     ReferenceRbDihedralBond rbTorsionBond;
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
+    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, rbTorsionBond);
 }
 
 double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numTorsions];
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceRbDihedralBond rbTorsionBond;
-    for (int i = 0; i < numTorsions; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
+    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, &energy, rbTorsionBond);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -690,24 +666,20 @@ void ReferenceCalcCustomTorsionForceKernel::executeForces(ContextImpl& context)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, customTorsion);
+    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, customTorsion);
 }
 
 double ReferenceCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
     RealOpenMM** posData = extractPositions(context);
     RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
-    RealOpenMM* energyArray = new RealOpenMM[numTorsions];
     RealOpenMM energy = 0;
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
-    for (int i = 0; i < numTorsions; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, customTorsion);
+    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, &energy, customTorsion);
     disposeRealArray(forceData, context.getSystem().getNumParticles());
-    delete[] energyArray;
     return energy;
 }
 
@@ -812,7 +784,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
     clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
-    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
+    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, nonbonded14);
 }
 
 double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
@@ -836,12 +808,8 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
-    RealOpenMM* energyArray = new RealOpenMM[num14];
-    for (int i = 0; i < num14; ++i)
-        energyArray[i] = 0;
-    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
+    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, &energy, nonbonded14);
     disposeRealArray(forceData, numParticles);
-    delete[] energyArray;
     if (periodic || ewald || pme) {
         RealOpenMM* boxSize = extractBoxSize(context);
         energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);

platforms/reference/src/SimTKReference/ReferenceAngleBondIxn.cpp

@@ -117,27 +117,15 @@ int ReferenceAngleBondIxn::getPrefactorsGivenAngleCosine( RealOpenMM cosine, Rea
    @param parameters       parameters: parameters[0] = ideal bond length
                                        parameters[1] = bond k (includes factor of 2)
    @param forces           force array (forces added)
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-
-   @return ReferenceForce::DefaultReturn
+   @param totalEnergy      if not null, the energy will be added to this
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceAngleBondIxn::calculateBondIxn( int* atomIndices,
+void ReferenceAngleBondIxn::calculateBondIxn( int* atomIndices,
                                              RealOpenMM** atomCoordinates,
                                              RealOpenMM* parameters,
                                              RealOpenMM** forces,
-                                             RealOpenMM* energiesByBond,
-                                             RealOpenMM* energiesByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceAngleBondIxn::calculateBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
-   static const std::string methodName = "\nReferenceAngleBondIxn::calculateBondIxn";
+                                             RealOpenMM* totalEnergy ) const {
 
    // constants -- reduce Visual Studio warnings regarding conversions between float & double
 
@@ -200,7 +188,6 @@ int ReferenceAngleBondIxn::calculateBondIxn( int* atomIndices,
 
    // accumulate energies
 
-   updateEnergy( energy, energiesByBond, LastAtomIndex, atomIndices, energiesByAtom );
-
-   return ReferenceForce::DefaultReturn;
+   if (totalEnergy != NULL)
+       *totalEnergy += energy;
 }

platforms/reference/src/SimTKReference/ReferenceAngleBondIxn.h

@@ -77,16 +77,13 @@ class ReferenceAngleBondIxn : public ReferenceBondIxn {
          @param parameters       parameters: parameters[0] = ideal bond length
                                              parameters[1] = bond k (includes factor of 2)
          @param forces           force array (forces added)
-         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-      
-         @return ReferenceForce::DefaultReturn
+         @param totalEnergy      if not null, the energy will be added to this
       
          --------------------------------------------------------------------------------------- */
       
-      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+      void calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                             RealOpenMM* parameters, RealOpenMM** forces,
-                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+                            RealOpenMM* totalEnergy ) const;
       
 
 };

platforms/reference/src/SimTKReference/ReferenceBondForce.cpp

@@ -72,21 +72,15 @@ ReferenceBondForce::~ReferenceBondForce( ){
    @param parameters       parameters: parameters[bondIndex][*]; contents of array 
                            depend on ixn
    @param forces           force array (forces added to current values): forces[atomIndex][3]
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-   @param totalEnergy      totalEnergy: sum over { energies[atomIndex] }
+   @param totalEnergy      if not null, the energy will be added to this
    @param ReferenceBondIxn ixn to be calculated
 
-   @return  ReferenceForce::DefaultReturn
-
    --------------------------------------------------------------------------------------- */
 
-int ReferenceBondForce::calculateForce( int numberOfBonds, int** atomIndices,
+void ReferenceBondForce::calculateForce( int numberOfBonds, int** atomIndices,
                                         RealOpenMM** atomCoordinates,
                                         RealOpenMM** parameters,
                                         RealOpenMM** forces, 
-                                        RealOpenMM* energiesByBond, 
-                                        RealOpenMM* energiesByAtom, 
                                         RealOpenMM *totalEnergy, 
                                         ReferenceBondIxn& referenceBondIxn ){
 
@@ -102,13 +96,7 @@ int ReferenceBondForce::calculateForce( int numberOfBonds, int** atomIndices,
       // calculate bond ixn
 
       referenceBondIxn.calculateBondIxn( atomIndices[ii], atomCoordinates, parameters[ii], 
-                                         forces, (energiesByBond == NULL ? NULL : energiesByBond + ii), energiesByAtom ); 
-
-      if( energiesByBond != NULL ){
-         *totalEnergy += energiesByBond[ii];
-      }
+                                         forces, totalEnergy );
    }
-
-   return ReferenceForce::DefaultReturn;
 }
 

platforms/reference/src/SimTKReference/ReferenceBondForce.h

@@ -62,19 +62,14 @@ class OPENMM_EXPORT ReferenceBondForce : public ReferenceForce {
          @param parameters       parameters: parameters[bondIndex][*]; contents of array 
                                  depend on ixn
          @param forces           force array (forces added to current values): forces[atomIndex][3]
-         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-         @param totalEnergy      totalEnergy: sum over { energies[atomIndex] }
+         @param totalEnergy      if not null, the energy will be added to this
          @param ReferenceBondIxn ixn to be calculated
       
-         @return  ReferenceForce::DefaultReturn
-      
          --------------------------------------------------------------------------------------- */
       
-       int calculateForce( int numberOfBonds, int** atomIndices, 
+       void calculateForce( int numberOfBonds, int** atomIndices,
                            RealOpenMM** atomCoordinates,
                            RealOpenMM** parameters, RealOpenMM** forces, 
-                           RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom, 
                            RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn );
 
 };

platforms/reference/src/SimTKReference/ReferenceBondIxn.cpp

@@ -63,42 +63,6 @@ ReferenceBondIxn::~ReferenceBondIxn( ){
 
 }
 
-/**---------------------------------------------------------------------------------------
-
-   Update energy
-
-   @param  energy               energy value to update 
-   @param  energyByBond         ptr to energyByBond accumulator (may be null)
-   @param  numberOfAtomIndices  number of atoms in bond
-   @param  atomIndices          array of atom indices of size 'numberOfAtomIndices'
-   @param  energyByAtom         array of energies by atom (may be null)
-
-   @return ReferenceForce::DefaultReturn 
-
-   --------------------------------------------------------------------------------------- */
-
-int ReferenceBondIxn::updateEnergy( RealOpenMM energy, RealOpenMM* energyByBond,
-                                    int numberOfAtomIndices, int* atomIndices, RealOpenMM* energyByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName = "\nReferenceBondIxn::updateEnergy";
-
-   // ---------------------------------------------------------------------------------------
-
-   if( energyByBond ){
-      *energyByBond += energy;
-   }
-   if( energyByAtom ){
-      for( int ii = 0; ii < numberOfAtomIndices; ii++ ){
-         energyByAtom[atomIndices[ii]] += energy;
-      }
-   }
-
-   return ReferenceForce::DefaultReturn;
-
-}
-
 /**---------------------------------------------------------------------------------------
       
    Calculate Bond Ixn -- virtual method -- does nothing
@@ -107,22 +71,19 @@ int ReferenceBondIxn::updateEnergy( RealOpenMM energy, RealOpenMM* energyByBond,
    @param atomCoordinates  atom coordinates
    @param parameters       parameters
    @param forces           force array (forces added)
-   @param energyByBond     bond energy 
-   @param energy           atom energy
+   @param totalEnergy      if not null, the energy will be added to this
       
    --------------------------------------------------------------------------------------- */
      
-   int ReferenceBondIxn::calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+   void ReferenceBondIxn::calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                                            RealOpenMM* parameters, RealOpenMM** forces,
-                                           RealOpenMM* energyByBond, RealOpenMM* energyByAtom ) const {
+                                           RealOpenMM* totalEnergy ) const {
    // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName = "\nReferenceBondIxn::calculateBondIxn";
 
    // ---------------------------------------------------------------------------------------
 
-   return ReferenceForce::DefaultReturn;
-
 }
  
 /**---------------------------------------------------------------------------------------

platforms/reference/src/SimTKReference/ReferenceBondIxn.h

@@ -59,32 +59,13 @@ class OPENMM_EXPORT ReferenceBondIxn {
          @param atomCoordinates  atom coordinates
          @param parameters       parameters
          @param forces           force array (forces added)
-         @param energyByBond     bond energy 
-         @param energy           atom energy
+         @param totalEnergy      if not null, the energy will be added to this
       
          --------------------------------------------------------------------------------------- */
       
-      virtual int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+      virtual void calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                                     RealOpenMM* parameters, RealOpenMM** forces,
-                                    RealOpenMM* energyByBond, RealOpenMM* energyByAtom ) const;
-
-      /**---------------------------------------------------------------------------------------
-      
-         Update energy
-      
-         @param  energy               energy value to update 
-         @param  energyByBond         ptr to energyByBond accumulator (may be null)
-         @param  numberOfAtomIndices  number of atoms in bond
-         @param  atomIndices          array of atom indices of size 'numberOfAtomIndices'
-         @param  energyByAtom         array of energies by atom (may be null)
-      
-         @return ReferenceForce::DefaultReturn 
-      
-         --------------------------------------------------------------------------------------- */
-      
-      int updateEnergy( RealOpenMM energy, RealOpenMM* energyByBond,
-                        int numberOfAtomIndices, int* atomIndices, RealOpenMM* energyByAtom ) const;
-      
+                                    RealOpenMM* totalEnergy ) const;
       
       /**---------------------------------------------------------------------------------------
       

platforms/reference/src/SimTKReference/ReferenceCMAPTorsionIxn.cpp

@@ -188,7 +188,6 @@ void ReferenceCMAPTorsionIxn::calculateOneIxn(int index, RealOpenMM** atomCoordi
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceCMAPTorsionIxn::calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
-        RealOpenMM* parameters, RealOpenMM** forces, RealOpenMM* energyByBond, RealOpenMM* energyByAtom) const {
-    return ReferenceForce::DefaultReturn;
+void ReferenceCMAPTorsionIxn::calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
+        RealOpenMM* parameters, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
 }

platforms/reference/src/SimTKReference/ReferenceCMAPTorsionIxn.h

@@ -84,9 +84,9 @@ public:
 
        --------------------------------------------------------------------------------------- */
 
-    int calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
+    void calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
                          RealOpenMM* parameters, RealOpenMM** forces,
-                         RealOpenMM* energyByBond, RealOpenMM* energyByAtom) const;
+                         RealOpenMM* totalEnergy) const;
 
 // ---------------------------------------------------------------------------------------
 

platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp

@@ -74,25 +74,15 @@ ReferenceCustomAngleIxn::~ReferenceCustomAngleIxn( ){
    @param atomCoordinates  atom coordinates
    @param parameters       parameters values
    @param forces           force array (forces added to input values)
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-
-   @return ReferenceForce::DefaultReturn;
+   @param totalEnergy      if not null, the energy will be added to this
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices,
+void ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices,
                                                 RealOpenMM** atomCoordinates,
                                                 RealOpenMM* parameters,
                                                 RealOpenMM** forces,
-                                                RealOpenMM* energiesByBond,
-                                                RealOpenMM* energiesByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomAngleIxn::calculateAngleIxn";
-
-   // ---------------------------------------------------------------------------------------
+                                                RealOpenMM* totalEnergy ) const {
 
    static const std::string methodName = "\nReferenceCustomAngleIxn::calculateAngleIxn";
 
@@ -156,8 +146,7 @@ int ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices,
 
    // accumulate energies
 
-   updateEnergy(energy, energiesByBond, 3, atomIndices, energiesByAtom);
-
-   return ReferenceForce::DefaultReturn;
+   if (totalEnergy != NULL)
+       *totalEnergy += energy;
 }
 

platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.h

@@ -64,14 +64,13 @@ class ReferenceCustomAngleIxn : public ReferenceBondIxn {
          @param atomCoordinates  atom coordinates
          @param parameters       parameter values
          @param forces           force array (forces added)
-         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+         @param totalEnergy      if not null, the energy will be added to this
 
          --------------------------------------------------------------------------------------- */
 
-      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+      void calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                             RealOpenMM* parameters, RealOpenMM** forces,
-                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+                            RealOpenMM* totalEnergy) const;
 
 
 };

platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp

@@ -75,25 +75,15 @@ ReferenceCustomBondIxn::~ReferenceCustomBondIxn( ){
    @param atomCoordinates  atom coordinates
    @param parameters       parameters values
    @param forces           force array (forces added to input values)
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-
-   @return ReferenceForce::DefaultReturn;
+   @param totalEnergy      if not null, the energy will be added to this
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
+void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
                                                 RealOpenMM** atomCoordinates,
                                                 RealOpenMM* parameters,
                                                 RealOpenMM** forces,
-                                                RealOpenMM* energiesByBond,
-                                                RealOpenMM* energiesByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomBondIxn::calculateBondIxn";
-
-   // ---------------------------------------------------------------------------------------
+                                                RealOpenMM* totalEnergy ) const {
 
    static const std::string methodName = "\nReferenceCustomBondIxn::calculateBondIxn";
 
@@ -127,8 +117,6 @@ int ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
    forces[atomBIndex][1]     -= dEdR*deltaR[ReferenceForce::YIndex];
    forces[atomBIndex][2]     -= dEdR*deltaR[ReferenceForce::ZIndex];
 
-   RealOpenMM energy          = (RealOpenMM) energyExpression.evaluate(variables);
-   updateEnergy( energy, energiesByBond, twoI, atomIndices, energiesByAtom );
-
-   return ReferenceForce::DefaultReturn;
+   if (totalEnergy != NULL)
+       *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
 }

platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.h

@@ -65,14 +65,13 @@ class ReferenceCustomBondIxn : public ReferenceBondIxn {
          @param atomCoordinates  atom coordinates
          @param parameters       parameter values
          @param forces           force array (forces added)
-         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+         @param totalEnergy      if not null, the energy will be added to this
 
          --------------------------------------------------------------------------------------- */
 
-      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+      void calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                             RealOpenMM* parameters, RealOpenMM** forces,
-                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+                            RealOpenMM* totalEnergy ) const;
 
 
 };

platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp

@@ -79,22 +79,14 @@ ReferenceCustomExternalIxn::~ReferenceCustomExternalIxn( ){
    @param forces           force array (forces added to input values)
    @param energy           energy is added to this
 
-   @return ReferenceForce::DefaultReturn;
-
    --------------------------------------------------------------------------------------- */
 
-int ReferenceCustomExternalIxn::calculateForce( int atomIndex,
+void ReferenceCustomExternalIxn::calculateForce( int atomIndex,
                                                 RealOpenMM** atomCoordinates,
                                                 RealOpenMM* parameters,
                                                 RealOpenMM** forces,
                                                 RealOpenMM* energy ) const {
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
    static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
 
    map<string, double> variables = globalParameters;
@@ -111,6 +103,4 @@ int ReferenceCustomExternalIxn::calculateForce( int atomIndex,
    forces[atomIndex][2] -= (RealOpenMM) forceExpressionZ.evaluate(variables);
    if (energy != NULL)
        *energy += (RealOpenMM) energyExpression.evaluate(variables);
-
-   return ReferenceForce::DefaultReturn;
 }

platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.h

@@ -72,7 +72,7 @@ class ReferenceCustomExternalIxn {
 
          --------------------------------------------------------------------------------------- */
 
-      int calculateForce( int atomIndex, RealOpenMM** atomCoordinates,
+      void calculateForce( int atomIndex, RealOpenMM** atomCoordinates,
                             RealOpenMM* parameters, RealOpenMM** forces, RealOpenMM* energy ) const;
 
 

platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp

@@ -74,25 +74,15 @@ ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn( ){
    @param atomCoordinates  atom coordinates
    @param parameters       parameters values
    @param forces           force array (forces added to input values)
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-
-   @return ReferenceForce::DefaultReturn;
+   @param totalEnergy      if not null, the energy will be added to this
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
+void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
                                                 RealOpenMM** atomCoordinates,
                                                 RealOpenMM* parameters,
                                                 RealOpenMM** forces,
-                                                RealOpenMM* energiesByBond,
-                                                RealOpenMM* energiesByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomTorsionIxn::calculateTorsionIxn";
-
-   // ---------------------------------------------------------------------------------------
+                                                RealOpenMM* totalEnergy ) const {
 
    static const std::string methodName = "\nReferenceCustomTorsionIxn::calculateTorsionIxn";
 
@@ -133,7 +123,6 @@ int ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
 
    // evaluate delta angle, dE/d(angle)
 
-   RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
    RealOpenMM dEdAngle = (RealOpenMM) forceExpression.evaluate(variables);
 
    // compute force
@@ -175,8 +164,7 @@ int ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
 
    // accumulate energies
 
-   updateEnergy( energy, energiesByBond, 4, atomIndices, energiesByAtom );
-
-   return ReferenceForce::DefaultReturn;
+   if (totalEnergy != NULL)
+       *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
 }
 

platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.h

@@ -64,14 +64,13 @@ class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
          @param atomCoordinates  atom coordinates
          @param parameters       parameter values
          @param forces           force array (forces added)
-         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+         @param totalEnergy      if not null, the energy will be added to this
 
          --------------------------------------------------------------------------------------- */
 
-      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+      void calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                             RealOpenMM* parameters, RealOpenMM** forces,
-                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+                            RealOpenMM* totalEnergy ) const;
 
 
 };

platforms/reference/src/SimTKReference/ReferenceHarmonicBondIxn.cpp

@@ -72,25 +72,15 @@ ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn( ){
    @param parameters       parameters: parameters[0] = ideal bond length
                                        parameters[1] = bond k
    @param forces           force array (forces added to input values)
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-
-   @return ReferenceForce::DefaultReturn;
+   @param totalEnergy      if not null, the energy will be added to this
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceHarmonicBondIxn::calculateBondIxn( int* atomIndices,
+void ReferenceHarmonicBondIxn::calculateBondIxn( int* atomIndices,
                                                 RealOpenMM** atomCoordinates,
                                                 RealOpenMM* parameters,
                                                 RealOpenMM** forces,
-                                                RealOpenMM* energiesByBond, 
-                                                RealOpenMM* energiesByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceHarmonicBondIxn::calculateBondIxn";
-
-   // ---------------------------------------------------------------------------------------
+                                                RealOpenMM* totalEnergy ) const {
 
    static const std::string methodName = "\nReferenceHarmonicBondIxn::calculateBondIxn";
 
@@ -129,8 +119,6 @@ int ReferenceHarmonicBondIxn::calculateBondIxn( int* atomIndices,
    forces[atomBIndex][1]     -= dEdR*deltaR[ReferenceForce::YIndex];
    forces[atomBIndex][2]     -= dEdR*deltaR[ReferenceForce::ZIndex];
 
-   RealOpenMM energy          = half*parameters[1]*deltaIdeal2;
-   updateEnergy( energy, energiesByBond, twoI, atomIndices, energiesByAtom );
-
-   return ReferenceForce::DefaultReturn;
+   if (totalEnergy != NULL)
+       *totalEnergy += half*parameters[1]*deltaIdeal2;
 }

platforms/reference/src/SimTKReference/ReferenceHarmonicBondIxn.h

@@ -60,15 +60,13 @@ class ReferenceHarmonicBondIxn : public ReferenceBondIxn {
          @param parameters       parameters: parameters[0] = ideal bond length
                                              parameters[1] = bond k
          @param forces           force array (forces added)
-         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+         @param totalEnergy      if not null, the energy will be added to this
 
          --------------------------------------------------------------------------------------- */
       
-      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+      void calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                             RealOpenMM* parameters, RealOpenMM** forces,
-                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
-      
+                            RealOpenMM* totalEnergy ) const;
 
 };
 

platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp

@@ -74,23 +74,13 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){
 										parameters[1]= c6*c6/c12      (4*epsilon)
 										parameters[2]= epsfac*q1*q2
    @param forces           force array (forces added to current values)
-   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
-   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
-
-   @return ReferenceForce::DefaultReturn
+   @param totalEnergy      if not null, the energy will be added to this
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+void ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
                                      RealOpenMM* parameters, RealOpenMM** forces,
-                                     RealOpenMM* energiesByBond,
-                                     RealOpenMM* energiesByAtom ) const {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceLJCoulomb14::calculateBondIxn";
-
-   // ---------------------------------------------------------------------------------------
+                                     RealOpenMM* totalEnergy ) const {
 
    static const std::string methodName = "\nReferenceLJCoulomb14::calculateBondIxn";
 
@@ -140,12 +130,8 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
       forces[atomBIndex][ii] -= force;
    }
 
-   RealOpenMM energy = parameters[1]*( sig6 - one )*sig6;
-   energy += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR);
-
    // accumulate energies
 
-   updateEnergy( energy, energiesByBond, LastAtomIndex, atomIndices, energiesByAtom );
-
-   return ReferenceForce::DefaultReturn;
+   if (totalEnergy != NULL)
+       *totalEnergy += parameters[1]*( sig6 - one )*sig6 + (ONE_4PI_EPS0*parameters[2]*inverseR);
 }