Moving the reference code into the platforms folder

platforms/reference/src/ReferenceProperDihedralBond.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceProperDihedralBond.h"
+#include "ReferenceForce.h"
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceProperDihedralBond constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceProperDihedralBond::ReferenceProperDihedralBond( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceProperDihedralBond::ReferenceProperDihedralBond";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceProperDihedralBond destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceProperDihedralBond::~ReferenceProperDihedralBond( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceProperDihedralBond::~ReferenceProperDihedralBond";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate proper dihedral bond ixn
+
+   @param atomIndices      atom indices of 4 atoms in bond
+   @param atomCoordinates  atom coordinates
+   @param parameters       3 parameters: parameters[0] = k
+                                         parameters[1] = ideal bond angle in degrees
+                                         parameters[2] = multiplicity
+   @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
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceProperDihedralBond::calculateBondIxn( int* atomIndices,
+                                                   RealOpenMM** atomCoordinates,
+                                                   RealOpenMM* parameters,
+                                                   RealOpenMM** forces,
+                                                   RealOpenMM* energiesByBond,
+                                                   RealOpenMM* energiesByAtom ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceProperDihedralBond::calculateBondIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "\nReferenceProperDihedralBond::calculateBondIxn";
+
+   // constants -- reduce Visual Studio warnings regarding conversions between float & double
+
+   static const RealOpenMM zero        =  0.0;
+   static const RealOpenMM one         =  1.0;
+   static const RealOpenMM two         =  2.0;
+   static const RealOpenMM three       =  3.0;
+   static const RealOpenMM oneM        = -1.0;
+
+   static const int threeI             = 3;
+
+   // debug flag
+
+   static const int debug              = 0;
+
+   static const int LastAtomIndex      = 4;
+
+   RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
+
+   RealOpenMM crossProductMemory[6];
+
+   // ---------------------------------------------------------------------------------------
+
+   // get deltaR, R2, and R between three pairs of atoms: [j,i], [j,k], [l,k]
+
+   int atomAIndex = atomIndices[0];
+   int atomBIndex = atomIndices[1];
+   int atomCIndex = atomIndices[2];
+   int atomDIndex = atomIndices[3];
+   ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] );  
+   ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1] );  
+   ReferenceForce::getDeltaR( atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2] );  
+
+   RealOpenMM dotDihedral;
+   RealOpenMM signOfAngle;
+   int hasREntry             = 1;
+
+   // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
+
+   RealOpenMM* crossProduct[2];
+   crossProduct[0]           = crossProductMemory;
+   crossProduct[1]           = crossProductMemory + 3;
+
+   // get dihedral angle
+
+   RealOpenMM dihedralAngle  =  getDihedralAngleBetweenThreeVectors( deltaR[0], deltaR[1], deltaR[2],
+                                                                     crossProduct, &dotDihedral, deltaR[0], 
+                                                                     &signOfAngle, hasREntry );
+
+   // evaluate delta angle, dE/d(angle) 
+
+   RealOpenMM deltaAngle     = parameters[2]*dihedralAngle - (parameters[1]*DEGREE_TO_RADIAN); 
+   RealOpenMM sinDeltaAngle  = SIN( deltaAngle );
+   RealOpenMM dEdAngle       = -parameters[0]*parameters[2]*sinDeltaAngle;
+   RealOpenMM energy         =  parameters[0]*(one + COS( deltaAngle ) );
+   
+   // compute force
+
+   RealOpenMM internalF[4][3];
+   RealOpenMM forceFactors[4];
+   RealOpenMM normCross1         = DOT3( crossProduct[0], crossProduct[0] );
+   RealOpenMM normBC             = deltaR[1][ReferenceForce::RIndex];
+              forceFactors[0]    = (-dEdAngle*normBC)/normCross1;
+
+   RealOpenMM normCross2         = DOT3( crossProduct[1], crossProduct[1] );
+              forceFactors[3]    = (dEdAngle*normBC)/normCross2;
+  
+              forceFactors[1]    = DOT3( deltaR[0], deltaR[1] );
+              forceFactors[1]   /= deltaR[1][ReferenceForce::R2Index];
+
+              forceFactors[2]    = DOT3( deltaR[2], deltaR[1] );
+              forceFactors[2]   /= deltaR[1][ReferenceForce::R2Index];
+
+   for( int ii = 0; ii < 3; ii++ ){
+
+      internalF[0][ii]  = forceFactors[0]*crossProduct[0][ii];
+      internalF[3][ii]  = forceFactors[3]*crossProduct[1][ii];
+
+      RealOpenMM s      = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii];
+
+      internalF[1][ii]  = internalF[0][ii] - s;
+      internalF[2][ii]  = internalF[3][ii] + s;
+   }
+  
+   // accumulate forces
+
+   for( int ii = 0; ii < 3; ii++ ){
+      forces[atomAIndex][ii] += internalF[0][ii];
+      forces[atomBIndex][ii] -= internalF[1][ii];
+      forces[atomCIndex][ii] -= internalF[2][ii];
+      forces[atomDIndex][ii] += internalF[3][ii];
+   }
+
+   // accumulate energies
+
+   updateEnergy( energy, energiesByBond, LastAtomIndex, atomIndices, energiesByAtom );
+
+   // debug 
+
+   if( debug ){
+      static bool printHeader = false;
+      std::stringstream message;
+      message << methodName;
+      message << std::endl;
+      if( !printHeader  ){  
+         printHeader = true;
+         message << std::endl;
+         message << methodName.c_str() << " a0 k [c q p s] r1 r2  angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
+      }   
+
+      message << std::endl;
+      for( int ii = 0; ii < 4; ii++ ){
+         message << " Atm " << atomIndices[ii] << " [" << atomCoordinates[atomIndices[ii]][0] << " " << atomCoordinates[atomIndices[ii]][1] << " " << atomCoordinates[atomIndices[ii]][2] << "] ";
+      }
+      message << std::endl << " Delta:";
+      for( int ii = 0; ii < (LastAtomIndex - 1); ii++ ){
+         message << " [";
+         for( int jj = 0; jj < ReferenceForce::LastDeltaRIndex; jj++ ){
+            message << deltaR[ii][jj] << " ";
+         }
+         message << "]";
+      }
+      message << std::endl;
+
+      message << std::endl << " Cross:";
+      for( int ii = 0; ii < 2; ii++ ){
+         message << " [";
+         for( int jj = 0; jj < 3; jj++ ){
+            message << crossProduct[ii][jj] << " ";
+         }
+         message << "]";
+      }
+      message << std::endl;
+
+      message << " k="     << parameters[0];
+      message << " a="     << parameters[1];
+      message << " m="     << parameters[2];
+      message << " ang="   << dihedralAngle;
+      message << " dotD="  << dotDihedral;
+      message << " sign="  << signOfAngle;
+      message << std::endl << "  ";
+
+      message << " deltaAngle=" << deltaAngle;
+      message << " dEdAngle=" << dEdAngle;
+      message << " E=" << energy << " force factors: [";
+      for( int ii = 0; ii < 4; ii++ ){
+         message << forceFactors[ii] << " ";
+      }
+      message << "] F=compute force; f=cumulative force";
+
+      message << std::endl << "  ";
+      for( int ii = 0; ii < LastAtomIndex; ii++ ){
+         message << " F" << (ii+1) << "[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, internalF[ii], threeI );
+         message << "]";
+      }   
+      message << std::endl << "  ";
+
+      for( int ii = 0; ii < LastAtomIndex; ii++ ){
+         message << " f" << (ii+1) << "[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, forces[atomIndices[ii]], threeI );
+         message << "]";
+      }
+
+      SimTKOpenMMLog::printMessage( message );
+   }   
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName << " DONE";
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   return ReferenceForce::DefaultReturn;
+}

platforms/reference/src/ReferenceProperDihedralBond.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceProperDihedralBond_H__
+#define __ReferenceProperDihedralBond_H__
+
+#include "ReferenceBondIxn.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceProperDihedralBond : public ReferenceBondIxn {
+
+   private:
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceProperDihedralBond( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceProperDihedralBond( );
+
+      /**---------------------------------------------------------------------------------------
+      
+        Calculate proper dihedral bond ixn
+      
+         @param atomIndices      atom indices of 4 atoms in bond
+         @param atomCoordinates  atom coordinates
+         @param parameters       3 parameters: parameters[0] = k
+                                               parameters[1] = ideal bond angle in degrees
+                                               parameters[2] = multiplicity
+         @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
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+                            RealOpenMM* parameters, RealOpenMM** forces,
+                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+      
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceProperDihedralBond_H__

platforms/reference/src/ReferenceRbDihedralBond.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceRbDihedralBond.h"
+#include "ReferenceForce.h"
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceRbDihedralBond constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceRbDihedralBond::ReferenceRbDihedralBond( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceRbDihedralBond::ReferenceRbDihedralBond";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceRbDihedralBond destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceRbDihedralBond::~ReferenceRbDihedralBond( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceRbDihedralBond::~ReferenceRbDihedralBond";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Ryckaert-Bellemans bond ixn
+
+   @param atomIndices      atom indices of 4 atoms in bond
+   @param atomCoordinates  atom coordinates
+   @param parameters       six RB parameters
+   @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
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceRbDihedralBond::calculateBondIxn( int* atomIndices,
+                                               RealOpenMM** atomCoordinates,
+                                               RealOpenMM* parameters,
+                                               RealOpenMM** forces,
+                                               RealOpenMM* energiesByBond,
+                                               RealOpenMM* energiesByAtom ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceRbDihedralBond::calculateBondIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "\nReferenceRbDihedralBond::calculateBondIxn";
+
+   // constants -- reduce Visual Studio warnings regarding conversions between float & double
+
+   static const RealOpenMM zero        =  0.0;
+   static const RealOpenMM one         =  1.0;
+   static const RealOpenMM two         =  2.0;
+   static const RealOpenMM three       =  3.0;
+   static const RealOpenMM oneM        = -1.0;
+
+   static const int threeI             = 3;
+
+   // number of parameters
+
+   static const int numberOfParameters = 6;
+
+   // debug flag
+
+   static const int debug              = 0;
+
+   static const int LastAtomIndex      = 4;
+
+   RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
+
+   RealOpenMM crossProductMemory[6];
+
+   // ---------------------------------------------------------------------------------------
+
+   // get deltaR, R2, and R between 2 atoms
+
+   int atomAIndex = atomIndices[0];
+   int atomBIndex = atomIndices[1];
+   int atomCIndex = atomIndices[2];
+   int atomDIndex = atomIndices[3];
+   ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] );  
+   ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1] );  
+   ReferenceForce::getDeltaR( atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2] );  
+
+   RealOpenMM cosPhi;
+   RealOpenMM signOfAngle;
+   int hasREntry             = 1;
+
+   // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
+
+   RealOpenMM* crossProduct[2];
+   crossProduct[0]           = crossProductMemory;
+   crossProduct[1]           = crossProductMemory + 3;
+   RealOpenMM dihederalAngle = getDihedralAngleBetweenThreeVectors( deltaR[0], deltaR[1], deltaR[2],
+                                                                    crossProduct, &cosPhi, deltaR[0], 
+                                                                    &signOfAngle, hasREntry );
+
+   // Gromacs: use polymer convention
+
+   if( dihederalAngle < zero ){
+      dihederalAngle += PI_M;
+   } else {
+      dihederalAngle -= PI_M;
+   }
+   cosPhi *= -one;
+
+   // Ryckaert-Bellemans:
+
+   // V = sum over i: { C_i*cos( psi )**i }, where psi = phi - PI, 
+   //                                              C_i is ith RB coefficient
+
+   RealOpenMM dEdAngle       = zero;
+   RealOpenMM energy         = parameters[0];
+   RealOpenMM cosFactor      = one;
+   for( int ii = 1; ii < numberOfParameters; ii++ ){
+      dEdAngle  -= ((RealOpenMM) ii)*parameters[ii]*cosFactor;
+      cosFactor *= cosPhi;
+      energy    += cosFactor*parameters[ii];
+   }
+
+   dEdAngle *= SIN( dihederalAngle );
+
+   RealOpenMM internalF[4][3];
+   RealOpenMM forceFactors[4];
+   RealOpenMM normCross1         = DOT3( crossProduct[0], crossProduct[0] );
+   RealOpenMM normBC             = deltaR[1][ReferenceForce::RIndex];
+              forceFactors[0]    = (-dEdAngle*normBC)/normCross1;
+
+   RealOpenMM normCross2         = DOT3( crossProduct[1], crossProduct[1] );
+              forceFactors[3]    = (dEdAngle*normBC)/normCross2;
+  
+              forceFactors[1]    = DOT3( deltaR[0], deltaR[1] );
+              forceFactors[1]   /= deltaR[1][ReferenceForce::R2Index];
+
+              forceFactors[2]    = DOT3( deltaR[2], deltaR[1] );
+              forceFactors[2]   /= deltaR[1][ReferenceForce::R2Index];
+
+   for( int ii = 0; ii < 3; ii++ ){
+
+      internalF[0][ii]  = forceFactors[0]*crossProduct[0][ii];
+      internalF[3][ii]  = forceFactors[3]*crossProduct[1][ii];
+
+      RealOpenMM s      = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii]; 
+
+      internalF[1][ii]  = internalF[0][ii] - s;
+      internalF[2][ii]  = internalF[3][ii] + s;
+   }
+
+   // accumulate forces
+
+   for( int ii = 0; ii < 3; ii++ ){
+      forces[atomAIndex][ii] += internalF[0][ii];
+      forces[atomBIndex][ii] -= internalF[1][ii];
+      forces[atomCIndex][ii] -= internalF[2][ii];
+      forces[atomDIndex][ii] += internalF[3][ii];
+   }
+
+   // accumulate energies
+
+   updateEnergy( energy, energiesByBond, LastAtomIndex, atomIndices, energiesByAtom );
+
+   // debug 
+
+   if( debug ){
+      static bool printHeader = false;
+      std::stringstream message;
+      message << methodName;
+      message << std::endl;
+      if( !printHeader  ){  
+         printHeader = true;
+         message << std::endl;
+         message << methodName.c_str() << " a0 k [c q p s] r1 r2  angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
+      }   
+
+      message << std::endl;
+      for( int ii = 0; ii < 4; ii++ ){
+         message << " Atm " << atomIndices[ii] << " [" << atomCoordinates[atomIndices[ii]][0] << " ";
+         message <<            atomCoordinates[atomIndices[ii]][1] << "  " << atomCoordinates[atomIndices[ii]][2] << "] ";
+      }
+      message << std::endl << " Delta:";
+      for( int ii = 0; ii < (LastAtomIndex - 1); ii++ ){
+         message << " [";
+         for( int jj = 0; jj < ReferenceForce::LastDeltaRIndex; jj++ ){
+            message << deltaR[ii][jj] << " ";
+         }
+         message << "]";
+      }
+      message << std::endl;
+
+      message << std::endl << " Cross:";
+      for( int ii = 0; ii < 2; ii++ ){
+         message << " [";
+         for( int jj = 0; jj < 3; jj++ ){
+            message << crossProduct[ii][jj] << " ";
+         }
+         message << "]";
+      }
+      message << std::endl;
+
+      message << " k="     << parameters[0];
+      message << " a="     << parameters[1];
+      message << " m="     << parameters[2];
+      message << " ang="   << dihederalAngle;
+      message << " dotD="  << cosPhi;
+      message << " sign="  << signOfAngle;
+      message << std::endl << "  ";
+
+      message << " dEdAngle=" << dEdAngle;
+      message << " E=" << energy << " force factors: [";
+      for( int ii = 0; ii < 4; ii++ ){
+         message << forceFactors[ii] << " ";
+      }
+      message << "] F=compute force; f=cumulative force";
+
+      message << std::endl << "  ";
+      for( int ii = 0; ii < LastAtomIndex; ii++ ){
+         message << " F" << (ii+1) << "[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, internalF[ii], threeI );
+         message << "]";
+      }   
+      message << std::endl << "  ";
+
+      for( int ii = 0; ii < LastAtomIndex; ii++ ){
+         message << " f" << (ii+1) << "[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, forces[atomIndices[ii]], threeI );
+         message << "]";
+      }
+
+      SimTKOpenMMLog::printMessage( message );
+   }   
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName << " DONE";
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   return ReferenceForce::DefaultReturn;
+}

platforms/reference/src/ReferenceRbDihedralBond.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceRbDihedralBond_H__
+#define __ReferenceRbDihedralBond_H__
+
+#include "ReferenceBondIxn.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceRbDihedralBond : public ReferenceBondIxn {
+
+   private:
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceRbDihedralBond( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceRbDihedralBond( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Calculate Ryckaert-Bellemans bond ixn
+      
+         @param atomIndices      atom indices of 4 atoms in bond
+         @param atomCoordinates  atom coordinates
+         @param parameters       six RB parameters
+         @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
+            
+         --------------------------------------------------------------------------------------- */
+      
+      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+                            RealOpenMM* parameters, RealOpenMM** forces,
+                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+      
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceRbDihedralBond_H__

platforms/reference/src/ReferenceShakeAlgorithm.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMUtilities.h"
+#include "SimTKOpenMMLog.h"
+#include "ReferenceShakeAlgorithm.h"
+#include "ReferenceDynamics.h"
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceShakeAlgorithm constructor
+
+   @param numberOfConstraints      number of constraints
+   @param atomIndices              block of atom indices
+   @param shakeParameters          Shake parameters
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceShakeAlgorithm::ReferenceShakeAlgorithm( int numberOfConstraints,
+                                                  int** atomIndices,
+                                                  RealOpenMM** shakeParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::ReferenceShakeAlgorithm";
+
+   static const RealOpenMM zero        =  0.0;
+   static const RealOpenMM one         =  1.0;
+   static const RealOpenMM two         =  2.0;
+   static const RealOpenMM three       =  3.0;
+   static const RealOpenMM oneM        = -1.0;
+
+   static const int threeI             =  3;
+
+   // ---------------------------------------------------------------------------------------
+
+   _numberOfConstraints        = numberOfConstraints;
+   _atomIndices                = atomIndices;
+   _shakeParameters            = shakeParameters;
+
+   _maximumNumberOfIterations  = 15;
+   _tolerance                  = (RealOpenMM) 1.0e-04;
+
+   // work arrays
+
+   _r_ij                       = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfConstraints, threeI, NULL, 
+                                                                                    1, zero, "r_ij" );
+
+   _d_ij2                      = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "dij_2" );
+   _distanceTolerance          = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "distanceTolerance" );
+   _reducedMasses              = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "reducedMasses" );
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceShakeAlgorithm destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceShakeAlgorithm::~ReferenceShakeAlgorithm( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::~ReferenceShakeAlgorithm";
+
+   // ---------------------------------------------------------------------------------------
+
+   SimTKOpenMMUtilities::freeTwoDRealOpenMMArray( _r_ij,  "r_ij" );
+
+   SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _d_ij2, "d_ij2" );
+   SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _distanceTolerance, "distanceTolerance" );
+   SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _reducedMasses, "reducedMasses" );
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get number of constraints
+
+   @return number of constraints
+
+   --------------------------------------------------------------------------------------- */
+   
+int ReferenceShakeAlgorithm::getNumberOfConstraints( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::getNumberOfConstraints";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _numberOfConstraints;
+}
+   
+/**---------------------------------------------------------------------------------------
+
+   Get maximum number of iterations
+
+   @return maximum number of iterations
+
+   --------------------------------------------------------------------------------------- */
+   
+int ReferenceShakeAlgorithm::getMaximumNumberOfIterations( void ) const {
+   
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::getMaximumNumberOfIterations";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _maximumNumberOfIterations;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set maximum number of iterations
+
+   @param maximumNumberOfIterations   new maximum number of iterations
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+   
+int ReferenceShakeAlgorithm::setMaximumNumberOfIterations( int maximumNumberOfIterations ){
+   
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::setMaximumNumberOfIterations";
+
+   // ---------------------------------------------------------------------------------------
+
+   _maximumNumberOfIterations = maximumNumberOfIterations;
+
+   return ReferenceDynamics::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get tolerance
+
+   @return tolerance
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM ReferenceShakeAlgorithm::getTolerance( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::getTolerance";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _tolerance;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set tolerance
+
+   @param tolerance new tolerance
+
+   @return tolerance
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceShakeAlgorithm::setTolerance( RealOpenMM tolerance ){
+
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceShakeAlgorithm::setTolerance";
+
+   // ---------------------------------------------------------------------------------------
+
+   _tolerance = tolerance;;
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Apply Shake algorithm
+
+   @param numberOfAtoms    number of atoms
+   @param atomCoordinates  atom coordinates
+   @param atomCoordinatesP atom coordinates prime
+   @param inverseMasses    1/mass
+
+   @return ReferenceDynamics::DefaultReturn if converge; else 
+    return ReferenceDynamics::ErrorReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceShakeAlgorithm::applyShake( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                         RealOpenMM** atomCoordinatesP,
+                                         RealOpenMM* inverseMasses ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\nReferenceShakeAlgorithm::applyShake";
+
+   static const RealOpenMM zero        =  0.0;
+   static const RealOpenMM one         =  1.0;
+   static const RealOpenMM two         =  2.0;
+   static const RealOpenMM three       =  3.0;
+   static const RealOpenMM oneM        = -1.0;
+   static const RealOpenMM half        =  0.5;
+
+   static const RealOpenMM epsilon6    = (RealOpenMM) 1.0e-06;
+
+   static int firstPass                = 1;
+   static int debug                    = 0;
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfConstraints = getNumberOfConstraints();
+
+   // temp arrays
+
+   RealOpenMM** r_ij                   = _r_ij;
+   RealOpenMM* d_ij2                   = _d_ij2;
+   RealOpenMM* distanceTolerance       = _distanceTolerance;
+   RealOpenMM* reducedMasses           = _reducedMasses;
+
+   // calculate reduced masses on 1st pass
+
+   if( firstPass ){
+      firstPass = 0;
+      for( int ii = 0; ii < _numberOfConstraints; ii++ ){
+         int atomI          = _atomIndices[ii][0];
+         int atomJ          = _atomIndices[ii][1];
+         reducedMasses[ii]  = half/( inverseMasses[atomI] + inverseMasses[atomJ] );
+      }
+   }
+
+   // setup: r_ij for each (i,j) constraint
+
+   RealOpenMM tolerance     = getTolerance();
+              tolerance    *= two;
+   for( int ii = 0; ii < _numberOfConstraints; ii++ ){
+
+      int atomI   = _atomIndices[ii][0];
+      int atomJ   = _atomIndices[ii][1];
+      for( int jj = 0; jj < 3; jj++ ){
+         r_ij[ii][jj]   = atomCoordinates[atomI][jj]  - atomCoordinates[atomJ][jj];
+      }
+      d_ij2[ii]              = DOT3( r_ij[ii], r_ij[ii] );
+      distanceTolerance[ii]  = d_ij2[ii]*tolerance;
+      if( distanceTolerance[ii] > zero ){
+         distanceTolerance[ii] = one/distanceTolerance[ii];
+      }
+   }
+
+   // main loop
+
+   int done                 = 0;
+   int iterations           = 0;
+   int numberConverged      = 0;
+   while( !done && iterations++ < getMaximumNumberOfIterations() ){
+      numberConverged  = 0; 
+      for( int ii = 0; ii < _numberOfConstraints; ii++ ){
+
+         int atomI   = _atomIndices[ii][0];
+         int atomJ   = _atomIndices[ii][1];
+
+         RealOpenMM rp_ij[3];
+         for( int jj = 0; jj < 3; jj++ ){
+            rp_ij[jj] = atomCoordinatesP[atomI][jj] - atomCoordinatesP[atomJ][jj]; 
+         }
+
+         RealOpenMM rp2  = DOT3( rp_ij, rp_ij );
+         RealOpenMM diff = d_ij2[ii] - rp2;
+         int iconv       = (int) ( FABS( diff )*distanceTolerance[ii] );
+         if( iconv ){
+            RealOpenMM rrpr  = DOT3(  rp_ij, r_ij[ii] );
+            RealOpenMM acor;
+            if( rrpr <  d_ij2[ii]*epsilon6 ){
+               std::stringstream message;
+               message << iterations << " Error: sign of rrpr < 0?\n";
+               SimTKOpenMMLog::printMessage( message );
+            } else {
+               acor  = reducedMasses[ii]*diff/rrpr;
+               for( int jj = 0; jj < 3; jj++ ){
+                  RealOpenMM dr                = acor*r_ij[ii][jj];
+                  atomCoordinatesP[atomI][jj] += inverseMasses[atomI]*dr;
+                  atomCoordinatesP[atomJ][jj] -= inverseMasses[atomJ]*dr;
+               }
+            }
+
+/*
+if( ii < -3 ){
+   std::stringstream message;
+   message << iterations << " C0 it=" << ii << " [" << atomI << " " << atomJ << "]";
+   message <<  " rp2=" << rp2 << " tol=" << d_ij2[ii];
+   message <<  " diff=" << diff << " acor=" << acor;
+   message <<  " m2=" << reducedMasses[ii];
+   message <<  " rm[" << inverseMasses[atomI] << " " << inverseMasses[atomJ];
+   message << " [" << atomCoordinatesP[atomI][0] << " " << atomCoordinatesP[atomI][1] << " " << atomCoordinatesP[atomI][2] << "] ";
+   message << " [" << atomCoordinatesP[atomJ][0] << " " << atomCoordinatesP[atomJ][1] << " " << atomCoordinatesP[atomJ][2] << "] ";
+   message <<  " rrpr=" << rrpr << " rijx=" << r_ij[ii][0];
+   message <<  " \n";
+   SimTKOpenMMLog::printMessage( message );
+} */
+
+         } else {
+            numberConverged++;
+         }
+      }
+      if( numberConverged == _numberOfConstraints ){
+         done = true;
+      }
+   }
+          
+   // diagnostics
+
+   if( debug || !done ){
+      std::stringstream message;
+      message << methodName;
+      message << " iterations=" << iterations << " no. converged=" << numberConverged << " out of " << _numberOfConstraints; 
+      if( done ){
+         message << " SUCCESS";
+      } else {
+         message << " FAILED";
+      }
+      message << "\n";
+      int errors = reportShake( numberOfAtoms, atomCoordinatesP, message );
+      if( !errors ){
+         message << "*** no errors recorded in explicit check ***";
+      }
+      message << "\n";
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   return (done ? ReferenceDynamics::DefaultReturn : ReferenceDynamics::ErrorReturn);
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Report any violated constriants
+
+   @param numberOfAtoms    number of atoms
+   @param atomCoordinates  atom coordinates
+   @param message          report
+
+   @return number of violated constraints
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceShakeAlgorithm::reportShake( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                          std::stringstream& message ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\nReferenceShakeAlgorithm::reportShake";
+
+   static const RealOpenMM zero        =  0.0;
+   static const RealOpenMM one         =  1.0;
+   static const RealOpenMM two         =  2.0;
+   static const RealOpenMM three       =  3.0;
+   static const RealOpenMM oneM        = -1.0;
+   static const RealOpenMM half        =  0.5;
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfConstraints = getNumberOfConstraints();
+
+   // loop over constraints calculating distance and comparing to
+   // expected distance -- report any contraints that are violated
+
+   int numberConverged  = 0; 
+   RealOpenMM tolerance = getTolerance();
+   for( int ii = 0; ii < _numberOfConstraints; ii++ ){
+
+      int atomI   = _atomIndices[ii][0];
+      int atomJ   = _atomIndices[ii][1];
+
+      RealOpenMM rp2  = zero;
+      for( int jj = 0; jj < 3; jj++ ){
+         rp2 += (atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj])*(atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj]); 
+      }
+      RealOpenMM diff = FABS( rp2 - (_shakeParameters[ii][0]*_shakeParameters[ii][0]) );
+      if( diff > tolerance ){
+         message << ii << " constraint violated: " << atomI << " " << atomJ << "] d=" << SQRT( rp2 ) << " " << rp2 << " d0=" << _shakeParameters[ii][0];
+         message << " diff=" << diff;
+         message << " [" << atomCoordinates[atomI][0] << " " << atomCoordinates[atomI][1] << " " << atomCoordinates[atomI][2] << "] ";
+         message << " [" << atomCoordinates[atomJ][0] << " " << atomCoordinates[atomJ][1] << " " << atomCoordinates[atomJ][2] << "] ";
+         message << "\n";
+      } else {
+         numberConverged++;
+      }
+   }
+          
+   return (numberOfConstraints-numberConverged);
+
+}

platforms/reference/src/ReferenceShakeAlgorithm.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceShakeAlgorithm_H__
+#define __ReferenceShakeAlgorithm_H__
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceShakeAlgorithm {
+
+   protected:
+
+      int _maximumNumberOfIterations;
+      RealOpenMM _tolerance;
+
+      int _numberOfConstraints;
+      int** _atomIndices;
+      RealOpenMM** _shakeParameters;
+
+      RealOpenMM** _r_ij;
+      RealOpenMM* _d_ij2;
+      RealOpenMM* _distanceTolerance;
+      RealOpenMM* _reducedMasses;
+      
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         ReferenceShakeAlgorithm constructor
+      
+         @param numberOfConstraints      number of constraints 
+         @param atomIndices              atom indices for contraints
+         @param shakeParameters          Shake parameters
+      
+         --------------------------------------------------------------------------------------- */
+      
+      ReferenceShakeAlgorithm( int numberOfConstraints, int** atomIndices, RealOpenMM** shakeParameters );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceShakeAlgorithm( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get number of constraints
+      
+         @return number of constraints
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int getNumberOfConstraints( void ) const;
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get maximum number of iterations
+      
+         @return maximum number of iterations
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int getMaximumNumberOfIterations( void ) const;
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Set maximum number of iterations
+      
+         @param maximumNumberOfIterations   new maximum number of iterations
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setMaximumNumberOfIterations( int maximumNumberOfIterations );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get tolerance
+      
+         @return tolerance
+      
+         --------------------------------------------------------------------------------------- */
+      
+      RealOpenMM getTolerance( void ) const;
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Set tolerance
+      
+         @param tolerance new tolerance
+
+         @return tolerance
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setTolerance( RealOpenMM tolerance );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Print parameters
+      
+         @param message message
+
+         @return ReferenceShakeAlgorithm::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+          
+      int printParameters( std::stringstream& message ) const;
+
+      /**---------------------------------------------------------------------------------------
+      
+         Apply Shake algorithm
+      
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param atomCoordinatesP atom coordinates prime
+         @param inverseMasses    1/mass
+      
+         @return ReferenceDynamics::DefaultReturn if converge; else 
+          return ReferenceDynamics::ErrorReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int applyShake( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                       RealOpenMM** atomCoordinatesP, RealOpenMM* inverseMasses );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Report any violated constraints
+      
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param message          report
+      
+         @return number of violated constraints
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int reportShake( int numberOfAtoms, RealOpenMM** atomCoordinates, std::stringstream& message );
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceShakeAlgorithm_H__

platforms/reference/src/ReferenceStochasticDynamics.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceStochasticDynamics.h"
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceStochasticDynamics constructor
+
+   @param numberOfAtoms  number of atoms
+   @param deltaT         delta t for dynamics
+   @param tau            viscosity(?)
+   @param temperature    temperature
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceStochasticDynamics::ReferenceStochasticDynamics( int numberOfAtoms,
+                                                          RealOpenMM deltaT, RealOpenMM tau,
+                                                          RealOpenMM temperature ) : 
+           ReferenceDynamics( numberOfAtoms, deltaT, temperature ), _tau( tau ) {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName      = "\nReferenceStochasticDynamics::ReferenceStochasticDynamics";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // insure tau is not zero -- if it is print warning message
+
+   if( _tau == zero ){
+
+      std::stringstream message;
+      message << methodName;
+      message << " input tau value=" << tau << " is invalid -- setting to 1.";
+      SimTKOpenMMLog::printError( message );
+
+      _tau = one;
+     
+   }
+   _setFixedParameters( );
+
+   allocate2DArrays( numberOfAtoms, 3, Max2DArrays );
+   allocate1DArrays( numberOfAtoms, Max1DArrays );
+   
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceStochasticDynamics destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceStochasticDynamics::~ReferenceStochasticDynamics( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceStochasticDynamics::~ReferenceStochasticDynamics";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set fixed parameters
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceStochasticDynamics::_setFixedParameters( void ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceStochasticDynamics::_setFixedParameters";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+   static const RealOpenMM two        =  2.0;
+   static const RealOpenMM three      =  3.0;
+   static const RealOpenMM four       =  4.0;
+   static const RealOpenMM half       =  0.5;
+
+   // ---------------------------------------------------------------------------------------
+
+   _fixedParameters[GDT]      = getDeltaT()/getTau();
+   _fixedParameters[EPH]      = EXP(  half*_fixedParameters[GDT] );
+   _fixedParameters[EMH]      = EXP( -half*_fixedParameters[GDT] );
+   _fixedParameters[EM]       = EXP(      -_fixedParameters[GDT] );
+   _fixedParameters[EP]       = EXP(       _fixedParameters[GDT] );
+
+   if( _fixedParameters[GDT] >= (RealOpenMM) 0.1 ){
+
+      RealOpenMM term1        = _fixedParameters[EPH] - one;
+                 term1       *= term1;
+      _fixedParameters[B]     = _fixedParameters[GDT]*(_fixedParameters[EP] - one) - four*term1;
+
+      _fixedParameters[C]     = _fixedParameters[GDT] - three + four*_fixedParameters[EMH] - _fixedParameters[EM];
+      _fixedParameters[D]     = two - _fixedParameters[EPH] - _fixedParameters[EMH];
+
+    } else {
+
+      // this has not been debugged
+
+      RealOpenMM term1        = half*_fixedParameters[GDT];
+      RealOpenMM term2        = term1*term1;
+      RealOpenMM term4        = term2*term2;
+
+      RealOpenMM third        = (RealOpenMM) ( 1.0/3.0 );
+      RealOpenMM o7_9         = (RealOpenMM) ( 7.0/9.0 );
+      RealOpenMM o1_12        = (RealOpenMM) ( 1.0/12.0 );
+      RealOpenMM o17_90       = (RealOpenMM) ( 17.0/90.0 );
+      RealOpenMM o7_30        = (RealOpenMM) ( 7.0/30.0 );
+      RealOpenMM o31_1260     = (RealOpenMM) ( 31.0/1260.0 );
+      RealOpenMM o_360        = (RealOpenMM) ( 1.0/360.0 );
+
+      _fixedParameters[B]     = term4*( third  + term1*( third + term1*( o17_90 + term1*o7_9 )));
+      _fixedParameters[C]     = term2*term1*( two*third + term1*( -half + term1*( o7_30 + term1*(-o1_12 + term1*o31_1260 ))));
+      _fixedParameters[D]     = term2*( -one + term2*(-o1_12 - term2*o_360));
+   }    
+
+   RealOpenMM kT        = ((RealOpenMM) BOLTZ)*getTemperature();
+
+   _fixedParameters[V]  = SQRT( kT*( one - _fixedParameters[EM]) );
+   _fixedParameters[X]  = getTau()*SQRT( kT*_fixedParameters[C] );
+   _fixedParameters[Yv] = SQRT( kT*_fixedParameters[B]/_fixedParameters[C] );
+   _fixedParameters[Yx] = getTau()*SQRT( kT*_fixedParameters[B]/(one - _fixedParameters[EM]) );
+
+   return ReferenceDynamics::DefaultReturn;
+
+};
+
+/**---------------------------------------------------------------------------------------
+
+   Get tau
+
+   @return tau
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM ReferenceStochasticDynamics::getTau( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nReferenceStochasticDynamics::getTau";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _tau;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get array of fixed parameters indexed by 'FixedParameters' enums
+
+   @return array
+
+   --------------------------------------------------------------------------------------- */
+   
+const RealOpenMM* ReferenceStochasticDynamics::getFixedParameters( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nReferenceStochasticDynamics::getFixedParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _fixedParameters;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Print parameters
+
+   @param message             message
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceStochasticDynamics::printParameters( std::stringstream& message ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nReferenceStochasticDynamics::printParameters";
+   static const char* parameterNames[MaxFixedParameters] = { "gdt", "ep", "eph", "emh", "em", "B", "C", "D",
+                                                             "V", "X", "Yv", "Yx" };
+
+   // ---------------------------------------------------------------------------------------
+
+   // print parameters
+
+   ReferenceDynamics::printParameters( message );
+   message << " tau=" << getTau();
+   message << " T=" << getTemperature();
+   int cut = 3;
+   for( int ii = 0; ii < MaxFixedParameters; ii++ ){
+      message << " " << parameterNames[ii] << "=" << _fixedParameters[ii];
+      if( cut++ > 5 ){
+         cut = 0;
+         message << std::endl;
+      }
+   }
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   First SD update; based on code in update.c do_update_sd() Gromacs 3.1.4
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param inverseMasses       inverse atom masses
+   @param xPrime              xPrime
+   @param oldVelocities       previous velocities
+   @param xVector             xVector
+   @param vVector             vVector
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                              RealOpenMM** velocities,
+                                              RealOpenMM** forces, RealOpenMM* inverseMasses,
+                                              RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                                              RealOpenMM** xVector, RealOpenMM** vVector ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nReferenceStochasticDynamics::updatePart1";
+
+   static const RealOpenMM one    =  1.0;
+   static int debug               = 0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // perform first update
+
+   const RealOpenMM*  fixedParameters = getFixedParameters();
+   RealOpenMM   tau                   = getTau();
+   RealOpenMM   fix1                  = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
+
+      for( int jj = 0; jj < 3; jj++ ){
+
+         oldVelocities[ii][jj] = velocities[ii][jj];
+
+         RealOpenMM Vmh        = xVector[ii][jj]*fixedParameters[D]/(tau*fixedParameters[C]) +
+                                 sqrtInvMass*fixedParameters[Yv]*getNormallyDistributedRandomNumber();
+
+           vVector[ii][jj]     = sqrtInvMass*fixedParameters[V]*getNormallyDistributedRandomNumber();
+          
+           velocities[ii][jj]  = oldVelocities[ii][jj]*fixedParameters[EM] +
+                                 inverseMasses[ii]*forces[ii][jj]*tau*( one - fixedParameters[EM]) +
+                                 vVector[ii][jj] - fixedParameters[EM]*Vmh;
+            
+           xPrime[ii][jj]      = atomCoordinates[ii][jj] + velocities[ii][jj]*fix1;
+      }
+   }
+
+   // diagnostics
+
+   if( debug ){
+      int maxPrint = 5;
+      std::stringstream message;
+      message << methodName << " Post SD1 atoms=" << numberOfAtoms << "\n";
+      RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+      for( int ii = 0; ii < maxPrint; ii++ ){
+         message << " mI=" << inverseMasses[ii];
+        	RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
+         message << " sdpc[" <<	sqrtInvMass*fixedParameters[Yv] << " " << sqrtInvMass*fixedParameters[V];
+         message << " x[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
+         message << "] xp[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
+         message << "] v[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
+         message << "] vV[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, vVector[ii], 3, one );
+         message << "] ov[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
+         message << "] xV[";
+         message << "] f[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, forces[ii], 3, one );
+         message << "] xV[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, xVector[ii], 3, one );
+         message << "]\n";
+      }
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Second update; based on code in update.c do_update_sd() w/ bFirstHalf = false in Gromacs 3.1.4
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                              RealOpenMM** velocities,
+                                              RealOpenMM** forces, RealOpenMM* inverseMasses,
+                                              RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                                              RealOpenMM** xVector, RealOpenMM** vVector ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nReferenceStochasticDynamics::updatePart2";
+
+   static const RealOpenMM one    =  1.0;
+   static int debug               = 0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // perform second update
+
+   const RealOpenMM*  fixedParameters = getFixedParameters();
+   RealOpenMM   tau                   = getTau();
+   RealOpenMM   fix1                  = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
+                fix1                  = one/fix1;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
+
+      for( int jj = 0; jj < 3; jj++ ){
+
+           velocities[ii][jj]  = (xPrime[ii][jj] - atomCoordinates[ii][jj])*fix1;
+         RealOpenMM Xmh        = vVector[ii][jj]*tau*fixedParameters[D]/(fixedParameters[EM] - one) +
+                                 sqrtInvMass*fixedParameters[Yx]*getNormallyDistributedRandomNumber();
+
+           xVector[ii][jj]     = sqrtInvMass*fixedParameters[X]*getNormallyDistributedRandomNumber();
+            
+           xPrime[ii][jj]     += xVector[ii][jj] - Xmh;
+      }
+   }
+
+   // diagnostics
+
+   if( debug ){
+      int maxPrint = 5;
+      std::stringstream message;
+      message << methodName << " Post SD2 atoms=" << numberOfAtoms << "\n";
+      RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+      for( int ii = 0; ii < maxPrint; ii++ ){
+         message << " mI=" << inverseMasses[ii];
+        	RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
+         message << " sdpc[" <<	sqrtInvMass*fixedParameters[Yx] << " " << sqrtInvMass*fixedParameters[X];
+         message << " x[";
+         message << " x[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
+         message << "] xp[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
+         message << "] v[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
+         message << "] vV[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, vVector[ii], 3, one );
+         message << "] oV[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
+         message << "] xV[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, xVector[ii], 3, one );
+         message << "]\n";
+      }
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Update -- driver routine for performing stochastic dynamics update of coordinates
+   and velocities
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                          RealOpenMM** velocities,
+                                          RealOpenMM** forces, RealOpenMM* masses ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName      = "\nReferenceStochasticDynamics::update";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+
+   static int debug                   =  0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get work arrays
+
+   RealOpenMM** xPrime          = get2DArrayAtIndex( xPrime2D );
+   RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+   RealOpenMM** xVector         = get2DArrayAtIndex( X2D );
+   RealOpenMM** vVector         = get2DArrayAtIndex( V2D );
+
+   RealOpenMM* inverseMasses    = get1DArrayAtIndex( InverseMasses );
+
+   // first-time-through initialization
+
+   if( getTimeStep() == 0 ){
+
+      std::stringstream message;
+      message << methodName;
+      int errors = 0;
+
+      // invert masses
+
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         if( masses[ii] <= zero ){
+            message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
+            errors++;
+         } else {
+            inverseMasses[ii] = one/masses[ii];
+         }
+      }
+
+      // set xVector 
+
+      const RealOpenMM*  fixedParameters = getFixedParameters();
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
+         for( int jj = 0; jj < 3; jj++ ){
+            xVector[ii][jj] = sqrtInverseMass*getNormallyDistributedRandomNumber();
+         }
+      }
+
+      // exit if errors
+
+      if( errors ){
+         SimTKOpenMMLog::printError( message );
+      }
+   }
+
+   // 1st update
+
+   updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses,
+                xPrime, oldVelocities, xVector, vVector );
+
+   //writeState( numberOfAtoms, atomCoordinates, velocities, forces, masses, -1 , "Sd1" );
+
+   ReferenceShakeAlgorithm* referenceShakeAlgorithm = getReferenceShakeAlgorithm();
+   if( referenceShakeAlgorithm ){
+
+/*
+      std::stringstream message;
+      message << methodName;
+      message << " calling shake1\n";
+      SimTKOpenMMLog::printMessage( message );
+*/
+
+      referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime,
+                                           inverseMasses );
+
+   }
+
+   // 2nd update
+
+   if( debug ){
+      int maxPrint = 5;
+      std::stringstream message;
+      message << methodName << " Pre SD2 atoms=" << numberOfAtoms << std::endl;
+      RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+      for( int ii = 0; ii < maxPrint; ii++ ){
+         message << " x[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
+         message << "] xp[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
+         message << "] v[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
+         message << "] ov[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
+         message << "]\n";
+      }
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses,
+                xPrime, oldVelocities, xVector, vVector );
+
+   if( debug ){
+      int maxPrint = 5;
+      std::stringstream message;
+      message << methodName << " Post SD2 atoms=" << numberOfAtoms << "\n";
+      RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+      for( int ii = 0; ii < maxPrint; ii++ ){
+         message << " x[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
+         message << "] xp[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
+         message << "] v[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
+         message << "] ov[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
+         message << "]\n";
+      }
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   if( referenceShakeAlgorithm ){
+
+/*
+      std::stringstream message;
+      message << methodName;
+      message << " calling shake2\n";
+      SimTKOpenMMLog::printMessage( message );
+*/
+
+      referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime,
+                                           inverseMasses );
+
+      // copy xPrime -> atomCoordinates
+
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         atomCoordinates[ii][0] = xPrime[ii][0];
+         atomCoordinates[ii][1] = xPrime[ii][1];
+         atomCoordinates[ii][2] = xPrime[ii][2];
+      }
+   }
+
+   incrementTimeStep();
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Write state
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+   @param state               0 if initial state; otherwise nonzero
+   @param baseFileName        base file name
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceStochasticDynamics::writeState( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                             RealOpenMM** velocities,
+                                             RealOpenMM** forces, RealOpenMM* masses,
+                                             int state, const std::string& baseFileName ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName      = "\nReferenceStochasticDynamics::writeState";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+   static const int threeI            =  3;
+
+   // ---------------------------------------------------------------------------------------
+
+   std::stringstream stateFileName;
+
+   stateFileName << baseFileName;
+   stateFileName << "_Step" << getTimeStep();
+   // stateFileName << "_State" << state;
+   stateFileName << ".txt";
+
+   // ---------------------------------------------------------------------------------------
+
+   // open file -- return if unsuccessful
+
+   FILE* stateFile = NULL;
+#ifdef WIN32
+   fopen_s( &stateFile, stateFileName.str().c_str(), "w" );
+#else
+   stateFile = fopen( stateFileName.str().c_str(), "w" );
+#endif
+
+   // ---------------------------------------------------------------------------------------
+
+   // diagnostics
+
+   if( stateFile != NULL ){
+      std::stringstream message;
+      message << methodName;
+      message << " Opened file=<" << stateFileName.str() << ">.\n";
+      SimTKOpenMMLog::printMessage( message );
+   } else {
+      std::stringstream message;
+      message << methodName;
+      message << " could not open file=<" << stateFileName.str() << "> -- abort output.\n";
+      SimTKOpenMMLog::printMessage( message );
+      return ReferenceDynamics::ErrorReturn;
+   }   
+
+   // ---------------------------------------------------------------------------------------
+
+   StringVector scalarNameI;
+   IntVector scalarI;
+
+   StringVector scalarNameR;
+   RealOpenMMVector scalarR;
+
+   StringVector scalarNameR1;
+   RealOpenMMPtrVector scalarR1;
+
+   StringVector scalarNameR2;
+   RealOpenMMPtrPtrVector scalarR2;
+
+   scalarI.push_back( getNumberOfAtoms() );
+   scalarNameI.push_back( "Atoms" );
+
+   scalarI.push_back( getTimeStep() );
+   scalarNameI.push_back( "Timestep" );
+
+   if( state == 0 || state == -1 ){
+
+      scalarR.push_back( getDeltaT() );
+      scalarNameR.push_back( "delta_t" );
+
+      scalarR.push_back( getTemperature() );
+      scalarNameR.push_back( "T" );
+
+      scalarR.push_back( getTau() );
+      scalarNameR.push_back( "Tau" );
+
+      scalarR1.push_back( masses );
+      scalarNameR1.push_back( "mass" );
+
+      scalarR2.push_back( atomCoordinates );
+      scalarNameR2.push_back( "coord" );
+
+      scalarR2.push_back( velocities );
+      scalarNameR2.push_back( "velocities" );
+
+      scalarR2.push_back( forces );
+      scalarNameR2.push_back( "forces" );
+
+      if( state == -1 ){
+
+         RealOpenMM** xPrime          = get2DArrayAtIndex( xPrime2D );
+         RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+         RealOpenMM** xVector         = get2DArrayAtIndex( X2D );
+         RealOpenMM** vVector         = get2DArrayAtIndex( V2D );
+
+         scalarR2.push_back( xPrime );
+         scalarNameR2.push_back( "xPrime" );
+
+         scalarR2.push_back( oldVelocities);
+         scalarNameR2.push_back( "vold" );
+
+         scalarR2.push_back( xVector );
+         scalarNameR2.push_back( "xVector" );
+
+         scalarR2.push_back( vVector );
+         scalarNameR2.push_back( "vVector" );
+      }
+      
+   } else {
+
+      scalarR2.push_back( atomCoordinates );
+      scalarNameR2.push_back( "coord" );
+
+      scalarR2.push_back( velocities );
+      scalarNameR2.push_back( "velocities" );
+
+   }
+
+   writeStateToFile( stateFile, scalarNameI, scalarI, scalarNameR, scalarR, getNumberOfAtoms(), scalarNameR1, scalarR1, threeI, scalarNameR2, scalarR2 ); 
+
+   (void) fclose( stateFile );
+
+   return ReferenceDynamics::DefaultReturn;
+
+}

platforms/reference/src/ReferenceStochasticDynamics.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceStochasticDynamics_H__
+#define __ReferenceStochasticDynamics_H__
+
+#include "ReferenceDynamics.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceStochasticDynamics : public ReferenceDynamics {
+
+   private:
+
+      enum FixedParameters   { GDT, EPH, EMH, EP, EM, B, C, D, V, X, Yv, Yx, MaxFixedParameters };
+      enum TwoDArrayIndicies { X2D, V2D, OldV, xPrime2D, vPrime2D, Max2DArrays };
+      enum OneDArrayIndicies { InverseMasses, Max1DArrays };
+
+      RealOpenMM _tau;
+      RealOpenMM _fixedParameters[MaxFixedParameters];
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set fixed values
+      
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int _setFixedParameters( void );
+      
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+
+         @param numberOfAtoms  number of atoms
+         @param deltaT         delta t for dynamics
+         @param tau            viscosity
+         @param temperature    temperature
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceStochasticDynamics( int numberOfAtoms, RealOpenMM deltaT, RealOpenMM tau, RealOpenMM temperature );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceStochasticDynamics( );
+
+      /**---------------------------------------------------------------------------------------
+      
+         Get tau
+      
+         @return tau
+      
+         --------------------------------------------------------------------------------------- */
+      
+      RealOpenMM getTau( void ) const;
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Get array of fixed parameters indexed by 'FixedParameters' enums
+      
+         @return array
+      
+         --------------------------------------------------------------------------------------- */
+      
+      const RealOpenMM* getFixedParameters( void ) const;
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Print parameters
+      
+         @param message message
+
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int printParameters( std::stringstream& message ) const;
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Update
+      
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+      
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+     
+      int update( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                  RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses );
+     
+      /**---------------------------------------------------------------------------------------
+      
+         First update; based on code in update.c do_update_sd() Gromacs 3.1.4
+      
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param inverseMasses       inverse atom masses
+         @param xPrime              xPrime
+         @param oldVelocities       previous velocities
+         @param xVector             xVector
+         @param vVector             vVector
+      
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
+                       RealOpenMM** forces, RealOpenMM* inverseMasses,
+                       RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                       RealOpenMM** xVector, RealOpenMM** vVector );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Second update
+      
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+      
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
+                       RealOpenMM** forces, RealOpenMM* inverseMasses,
+                       RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                       RealOpenMM** xVector, RealOpenMM** vVector );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Write state
+      
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+         @param state               0 if initial state; otherwise nonzero
+         @param baseFileName        base file name
+      
+         @return ReferenceDynamics::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int writeState( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                      RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses,
+                      int state, const std::string& baseFileName ) const;
+      
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceStochasticDynamics_H__

platforms/reference/src/SimTKReference-old-2005.vcproj

+<?xml version="1.0" encoding="UTF-8"?>
+<VisualStudioProject
+	ProjectType="Visual C++"
+	Version="8.00"
+	Name="SimTKReference"
+	ProjectGUID="{B8F36313-772E-46F3-8C07-C59D6B1772A0}"
+	Keyword="Win32Proj"
+	>
+	<Platforms>
+		<Platform
+			Name="Win32"
+		/>
+	</Platforms>
+	<ToolFiles>
+	</ToolFiles>
+	<Configurations>
+		<Configuration
+			Name="Debug|Win32"
+			OutputDirectory="Debug"
+			IntermediateDirectory="Debug"
+			ConfigurationType="4"
+			>
+			<Tool
+				Name="VCPreBuildEventTool"
+			/>
+			<Tool
+				Name="VCCustomBuildTool"
+			/>
+			<Tool
+				Name="VCXMLDataGeneratorTool"
+			/>
+			<Tool
+				Name="VCWebServiceProxyGeneratorTool"
+			/>
+			<Tool
+				Name="VCMIDLTool"
+			/>
+			<Tool
+				Name="VCCLCompilerTool"
+				Optimization="0"
+				AdditionalIncludeDirectories="../SimTKUtilities;../gromacs/include"
+				PreprocessorDefinitions="WIN32;_DEBUG;_LIB;"
+				MinimalRebuild="true"
+				BasicRuntimeChecks="3"
+				RuntimeLibrary="1"
+				UsePrecompiledHeader="0"
+				WarningLevel="3"
+				DebugInformationFormat="4"
+			/>
+			<Tool
+				Name="VCManagedResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCPreLinkEventTool"
+			/>
+			<Tool
+				Name="VCLibrarianTool"
+			/>
+			<Tool
+				Name="VCALinkTool"
+			/>
+			<Tool
+				Name="VCXDCMakeTool"
+			/>
+			<Tool
+				Name="VCBscMakeTool"
+			/>
+			<Tool
+				Name="VCFxCopTool"
+			/>
+			<Tool
+				Name="VCPostBuildEventTool"
+			/>
+		</Configuration>
+		<Configuration
+			Name="Release|Win32"
+			OutputDirectory="Release"
+			IntermediateDirectory="Release"
+			ConfigurationType="4"
+			>
+			<Tool
+				Name="VCPreBuildEventTool"
+			/>
+			<Tool
+				Name="VCCustomBuildTool"
+			/>
+			<Tool
+				Name="VCXMLDataGeneratorTool"
+			/>
+			<Tool
+				Name="VCWebServiceProxyGeneratorTool"
+			/>
+			<Tool
+				Name="VCMIDLTool"
+			/>
+			<Tool
+				Name="VCCLCompilerTool"
+				AdditionalIncludeDirectories="../SimTKUtilities; ../gromacs/include"
+				PreprocessorDefinitions="WIN32;NDEBUG;_LIB;"
+				RuntimeLibrary="0"
+				UsePrecompiledHeader="0"
+				WarningLevel="3"
+				DebugInformationFormat="3"
+			/>
+			<Tool
+				Name="VCManagedResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCPreLinkEventTool"
+			/>
+			<Tool
+				Name="VCLibrarianTool"
+			/>
+			<Tool
+				Name="VCALinkTool"
+			/>
+			<Tool
+				Name="VCXDCMakeTool"
+			/>
+			<Tool
+				Name="VCBscMakeTool"
+			/>
+			<Tool
+				Name="VCFxCopTool"
+			/>
+			<Tool
+				Name="VCPostBuildEventTool"
+			/>
+		</Configuration>
+		<Configuration
+			Name="DebugPurify|Win32"
+			OutputDirectory="$(ConfigurationName)"
+			IntermediateDirectory="$(ConfigurationName)"
+			ConfigurationType="4"
+			>
+			<Tool
+				Name="VCPreBuildEventTool"
+			/>
+			<Tool
+				Name="VCCustomBuildTool"
+			/>
+			<Tool
+				Name="VCXMLDataGeneratorTool"
+			/>
+			<Tool
+				Name="VCWebServiceProxyGeneratorTool"
+			/>
+			<Tool
+				Name="VCMIDLTool"
+			/>
+			<Tool
+				Name="VCCLCompilerTool"
+				Optimization="0"
+				AdditionalIncludeDirectories="../SimTKUtilities;../gromacs/include"
+				PreprocessorDefinitions="WIN32;_DEBUG;_LIB;"
+				MinimalRebuild="true"
+				BasicRuntimeChecks="3"
+				RuntimeLibrary="1"
+				UsePrecompiledHeader="0"
+				WarningLevel="3"
+				DebugInformationFormat="4"
+			/>
+			<Tool
+				Name="VCManagedResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCPreLinkEventTool"
+			/>
+			<Tool
+				Name="VCLibrarianTool"
+			/>
+			<Tool
+				Name="VCALinkTool"
+			/>
+			<Tool
+				Name="VCXDCMakeTool"
+			/>
+			<Tool
+				Name="VCBscMakeTool"
+			/>
+			<Tool
+				Name="VCFxCopTool"
+			/>
+			<Tool
+				Name="VCPostBuildEventTool"
+			/>
+		</Configuration>
+	</Configurations>
+	<References>
+	</References>
+	<Files>
+		<Filter
+			Name="Header Files"
+			Filter="h;hpp;hxx;hm;inl;inc;xsd"
+			UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"
+			>
+			<File
+				RelativePath=".\gromacsReferenceInterface.h"
+				>
+			</File>
+			<File
+				RelativePath=".\gromacsReferenceInterfaceCpp.h"
+				>
+			</File>
+			<File
+				RelativePath=".\gromacsReferenceParameterPrint.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceAngleBondIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondForce.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceConstraint.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceDynamics.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceForce.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceHarmonicBondIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJ14.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJCoulombIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferencePairIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceProperDihedralBond.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceRbDihedralBond.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceShakeAlgorithm.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceStochasticDynamics.h"
+				>
+			</File>
+		</Filter>
+		<Filter
+			Name="Resource Files"
+			Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx"
+			UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"
+			>
+		</Filter>
+		<Filter
+			Name="Source Files"
+			Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"
+			UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"
+			>
+			<File
+				RelativePath=".\gromacsReferenceInterface.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\gromacsReferenceParameterPrint.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceAngleBondIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondForce.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceConstraint.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceDynamics.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceForce.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceHarmonicBondIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJ14.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJCoulombIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferencePairIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceProperDihedralBond.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceRbDihedralBond.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceShakeAlgorithm.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceStochasticDynamics.cpp"
+				>
+			</File>
+		</Filter>
+		<File
+			RelativePath=".\debug\BuildLog.htm"
+			DeploymentContent="true"
+			>
+		</File>
+	</Files>
+	<Globals>
+	</Globals>
+</VisualStudioProject>

platforms/reference/src/SimTKReference.sln

+
+Microsoft Visual Studio Solution File, Format Version 9.00
+# Visual Studio 2005
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "SimTKReference", "SimTKReference.vcproj", "{B8F36313-772E-46F3-8C07-C59D6B1772A0}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|Win32 = Debug|Win32
+		Release|Win32 = Release|Win32
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{B8F36313-772E-46F3-8C07-C59D6B1772A0}.Debug|Win32.ActiveCfg = Debug|Win32
+		{B8F36313-772E-46F3-8C07-C59D6B1772A0}.Debug|Win32.Build.0 = Debug|Win32
+		{B8F36313-772E-46F3-8C07-C59D6B1772A0}.Release|Win32.ActiveCfg = Release|Win32
+		{B8F36313-772E-46F3-8C07-C59D6B1772A0}.Release|Win32.Build.0 = Release|Win32
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+EndGlobal

platforms/reference/src/SimTKReference.vcproj

+<?xml version="1.0" encoding="UTF-8"?>
+<VisualStudioProject
+	ProjectType="Visual C++"
+	Version="8.00"
+	Name="SimTKReference"
+	ProjectGUID="{B8F36313-772E-46F3-8C07-C59D6B1772A0}"
+	Keyword="Win32Proj"
+	>
+	<Platforms>
+		<Platform
+			Name="Win32"
+		/>
+	</Platforms>
+	<ToolFiles>
+	</ToolFiles>
+	<Configurations>
+		<Configuration
+			Name="Debug|Win32"
+			OutputDirectory="Debug"
+			IntermediateDirectory="Debug"
+			ConfigurationType="4"
+			>
+			<Tool
+				Name="VCPreBuildEventTool"
+			/>
+			<Tool
+				Name="VCCustomBuildTool"
+			/>
+			<Tool
+				Name="VCXMLDataGeneratorTool"
+			/>
+			<Tool
+				Name="VCWebServiceProxyGeneratorTool"
+			/>
+			<Tool
+				Name="VCMIDLTool"
+			/>
+			<Tool
+				Name="VCCLCompilerTool"
+				Optimization="0"
+				AdditionalIncludeDirectories="../SimTKUtilities;../gromacs/include"
+				PreprocessorDefinitions="WIN32;_DEBUG;_LIB;"
+				MinimalRebuild="true"
+				BasicRuntimeChecks="3"
+				RuntimeLibrary="1"
+				UsePrecompiledHeader="0"
+				WarningLevel="3"
+				DebugInformationFormat="4"
+			/>
+			<Tool
+				Name="VCManagedResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCPreLinkEventTool"
+			/>
+			<Tool
+				Name="VCLibrarianTool"
+			/>
+			<Tool
+				Name="VCALinkTool"
+			/>
+			<Tool
+				Name="VCXDCMakeTool"
+			/>
+			<Tool
+				Name="VCBscMakeTool"
+			/>
+			<Tool
+				Name="VCFxCopTool"
+			/>
+			<Tool
+				Name="VCPostBuildEventTool"
+			/>
+		</Configuration>
+		<Configuration
+			Name="Release|Win32"
+			OutputDirectory="Release"
+			IntermediateDirectory="Release"
+			ConfigurationType="4"
+			>
+			<Tool
+				Name="VCPreBuildEventTool"
+			/>
+			<Tool
+				Name="VCCustomBuildTool"
+			/>
+			<Tool
+				Name="VCXMLDataGeneratorTool"
+			/>
+			<Tool
+				Name="VCWebServiceProxyGeneratorTool"
+			/>
+			<Tool
+				Name="VCMIDLTool"
+			/>
+			<Tool
+				Name="VCCLCompilerTool"
+				AdditionalIncludeDirectories="../SimTKUtilities; ../gromacs/include"
+				PreprocessorDefinitions="WIN32;NDEBUG;_LIB;"
+				RuntimeLibrary="0"
+				UsePrecompiledHeader="0"
+				WarningLevel="3"
+				DebugInformationFormat="3"
+			/>
+			<Tool
+				Name="VCManagedResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCPreLinkEventTool"
+			/>
+			<Tool
+				Name="VCLibrarianTool"
+			/>
+			<Tool
+				Name="VCALinkTool"
+			/>
+			<Tool
+				Name="VCXDCMakeTool"
+			/>
+			<Tool
+				Name="VCBscMakeTool"
+			/>
+			<Tool
+				Name="VCFxCopTool"
+			/>
+			<Tool
+				Name="VCPostBuildEventTool"
+			/>
+		</Configuration>
+		<Configuration
+			Name="DebugPurify|Win32"
+			OutputDirectory="$(ConfigurationName)"
+			IntermediateDirectory="$(ConfigurationName)"
+			ConfigurationType="4"
+			>
+			<Tool
+				Name="VCPreBuildEventTool"
+			/>
+			<Tool
+				Name="VCCustomBuildTool"
+			/>
+			<Tool
+				Name="VCXMLDataGeneratorTool"
+			/>
+			<Tool
+				Name="VCWebServiceProxyGeneratorTool"
+			/>
+			<Tool
+				Name="VCMIDLTool"
+			/>
+			<Tool
+				Name="VCCLCompilerTool"
+				Optimization="0"
+				AdditionalIncludeDirectories="../SimTKUtilities;../gromacs/include"
+				PreprocessorDefinitions="WIN32;_DEBUG;_LIB;"
+				MinimalRebuild="true"
+				BasicRuntimeChecks="3"
+				RuntimeLibrary="1"
+				UsePrecompiledHeader="0"
+				WarningLevel="3"
+				DebugInformationFormat="4"
+			/>
+			<Tool
+				Name="VCManagedResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCResourceCompilerTool"
+			/>
+			<Tool
+				Name="VCPreLinkEventTool"
+			/>
+			<Tool
+				Name="VCLibrarianTool"
+			/>
+			<Tool
+				Name="VCALinkTool"
+			/>
+			<Tool
+				Name="VCXDCMakeTool"
+			/>
+			<Tool
+				Name="VCBscMakeTool"
+			/>
+			<Tool
+				Name="VCFxCopTool"
+			/>
+			<Tool
+				Name="VCPostBuildEventTool"
+			/>
+		</Configuration>
+	</Configurations>
+	<References>
+	</References>
+	<Files>
+		<Filter
+			Name="Header Files"
+			Filter="h;hpp;hxx;hm;inl;inc;xsd"
+			UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"
+			>
+			<File
+				RelativePath=".\gromacsReferenceInterface.h"
+				>
+			</File>
+			<File
+				RelativePath=".\gromacsReferenceInterfaceCpp.h"
+				>
+			</File>
+			<File
+				RelativePath=".\gromacsReferenceParameterPrint.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceAngleBondIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondForce.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceConstraint.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceDynamics.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceForce.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceHarmonicBondIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJ14.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJCoulombIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferencePairIxn.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceProperDihedralBond.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceRbDihedralBond.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceShakeAlgorithm.h"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceStochasticDynamics.h"
+				>
+			</File>
+		</Filter>
+		<Filter
+			Name="Resource Files"
+			Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx"
+			UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"
+			>
+		</Filter>
+		<Filter
+			Name="Source Files"
+			Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"
+			UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"
+			>
+			<File
+				RelativePath=".\gromacsReferenceInterface.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\gromacsReferenceParameterPrint.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceAngleBondIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondForce.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceBondIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceConstraint.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceDynamics.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceForce.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceHarmonicBondIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJ14.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceLJCoulombIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferencePairIxn.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceProperDihedralBond.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceRbDihedralBond.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceShakeAlgorithm.cpp"
+				>
+			</File>
+			<File
+				RelativePath=".\ReferenceStochasticDynamics.cpp"
+				>
+			</File>
+		</Filter>
+		<File
+			RelativePath=".\debug\BuildLog.htm"
+			>
+		</File>
+	</Files>
+	<Globals>
+	</Globals>
+</VisualStudioProject>

platforms/reference/src/gromacsReferenceInterface.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+// #include "SimTKOpenMMGromacsUtilities.h"
+
+#include "SimTKOpenMMGromacsUtilities.h"
+#include "gromacsReferenceInterface.h"
+#include "gromacsReferenceInterfaceCpp.h"
+
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+
+#include "ReferenceBondForce.h"
+#include "ReferenceHarmonicBondIxn.h"
+#include "ReferenceAngleBondIxn.h"
+#include "ReferenceProperDihedralBond.h"
+#include "ReferenceRbDihedralBond.h"
+#include "ReferenceLJ14.h"
+#include "ReferenceLJCoulombIxn.h"
+
+#include "ReferenceStochasticDynamics.h"
+#include "ReferenceConstraint.h"
+
+/**---------------------------------------------------------------------------------------
+
+   Allocate memory for atom indices participating in bonds and the bond parameters
+
+   @param numberOfBonds            number of bonds
+   @param numberOfAtomIndices      (number of atoms)/bond
+   @param outputBondIndices        allocated memory for atom indices outputBondIndices[bondIndex][atomIndex]
+   @param numberParameters         number of parameters/bond
+   @param outputBondParameters     allocated memory for parameters outputBondParameters[bondIndex][parameterIndex]
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsAllocateBondIxnArrays( int numberOfBonds, int numberOfAtomIndices,
+                                  int*** outputBondIndices, int numberParameters,
+                                      RealOpenMM*** outputBondParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsAllocateBondIxnArrays: ";
+
+   // ---------------------------------------------------------------------------------------
+
+   int** atomIndices                   = (int**) malloc( sizeof( int* )*numberOfBonds );
+   int* atomIndicesBlock               = (int*)  malloc( sizeof( int  )*numberOfBonds*numberOfAtomIndices );
+   memset( atomIndicesBlock, 0, sizeof( int  )*numberOfBonds*numberOfAtomIndices );
+
+   RealOpenMM** bondParameters         = (RealOpenMM**) malloc( sizeof( RealOpenMM* )*numberOfBonds );
+   RealOpenMM* bondParametersBlock     = (RealOpenMM*)  malloc( sizeof( RealOpenMM )*numberOfBonds*numberParameters );
+   memset( bondParametersBlock, 0, sizeof( RealOpenMM )*numberOfBonds*numberParameters );
+
+   for( int ii = 0; ii < numberOfBonds; ii++ ){ 
+
+      atomIndices[ii]       = atomIndicesBlock;
+      atomIndicesBlock     += numberOfAtomIndices;
+
+      bondParameters[ii]    = bondParametersBlock;
+      bondParametersBlock  += numberParameters;
+
+   }
+
+   *outputBondIndices    = atomIndices;
+   *outputBondParameters = bondParameters;
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Free memory for arrays
+
+   @param indexOneDArraysToFree    number of 1D RealOpenMM arrays
+   @param indexTwoDArraysToFree    number of 2D RealOpenMM arrays
+   @param indexTwoDIntArraysToFree number of 2D int arrays
+   @param MaxFreeArrayIndex        max index on arrays
+   @param oneDArraysToFree         1D RealOpenMM arrays to free
+   @param twoDArraysToFree         2D RealOpenMM arrays to free
+   @param twoDIntArraysToFree      2D int arrays to free
+   @param whoIsCalling             callee
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsFreeSundryArrays( int indexOneDArraysToFree, int indexTwoDArraysToFree,
+                             int indexTwoDIntArraysToFree, int MaxFreeArrayIndex,
+                             RealOpenMM** oneDArraysToFree,
+                             RealOpenMM*** twoDArraysToFree, int*** twoDIntArraysToFree,
+                             const char* whoIsCalling ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsFreeSundryArrays ";
+
+   // ---------------------------------------------------------------------------------------
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   if( indexTwoDArraysToFree    >= MaxFreeArrayIndex ||
+       indexTwoDIntArraysToFree >= MaxFreeArrayIndex ||
+       indexOneDArraysToFree    >= MaxFreeArrayIndex ){
+
+      std::stringstream message;
+      message << methodName << " " << whoIsCalling;
+      message << " Free array index too large: [" << indexTwoDArraysToFree << " " << indexTwoDIntArraysToFree << " " << indexOneDArraysToFree << "]";
+      message << " should be less than " << MaxFreeArrayIndex;
+
+      SimTKOpenMMLog::printMessage( message );
+
+      indexTwoDArraysToFree     = indexTwoDArraysToFree    >= MaxFreeArrayIndex ? MaxFreeArrayIndex : indexTwoDArraysToFree;
+      indexTwoDIntArraysToFree  = indexTwoDIntArraysToFree >= MaxFreeArrayIndex ? MaxFreeArrayIndex : indexTwoDIntArraysToFree;
+      indexOneDArraysToFree     = indexOneDArraysToFree    >= MaxFreeArrayIndex ? MaxFreeArrayIndex : indexOneDArraysToFree;
+   }
+
+   for( int ii = 0; ii < indexTwoDArraysToFree; ii++ ){
+      free( twoDArraysToFree[ii][0] );
+      free(  twoDArraysToFree[ii] );
+   }
+
+   for( int ii = 0; ii < indexTwoDIntArraysToFree; ii++ ){
+      free( twoDIntArraysToFree[ii][0] );
+      free(  twoDIntArraysToFree[ii] );
+   }
+
+   for( int ii = 0; ii < indexOneDArraysToFree; ii++ ){
+      free( oneDArraysToFree[ii] );
+   }
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for harmonic bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetHarmonicBondParameters( const t_topology* top, int*** atomIndices,
+                                      RealOpenMM*** bondParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetHarmonicBondParameters: ";
+   static const int twoI           = 2;
+   static const RealOpenMM two     = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef   = &top->idef;
+   int numberOfBonds    = idef->il[ F_BONDS ].nr/3; 
+   int* atoms           = (int*) idef->il[ F_BONDS ].iatoms;
+
+   gromacsAllocateBondIxnArrays( numberOfBonds, twoI, atomIndices, twoI, bondParameters );
+
+   int offsetIndex = 0;
+   for( int ii = 0; ii < numberOfBonds; ii++ ) { 
+
+      int type                 = atoms[ offsetIndex++ ];  
+      (*atomIndices)[ii][0]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][1]    = atoms[ offsetIndex++ ];
+
+      (*bondParameters)[ii][0] = idef->iparams[ type ].harmonic.rA;
+      (*bondParameters)[ii][1] = idef->iparams[ type ].harmonic.krA;
+   }
+
+   return numberOfBonds;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate harmonic bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param harmonicBondForces       output array of forces: harmonicBondForces[atomIndex][3]
+   @param harmonicBondEnergies     output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsHarmonicBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                       RealOpenMM** harmonicBondForces,
+                                       RealOpenMM* harmonicBondAtomEnergies,
+                                       RealOpenMM* totalEnergy ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName         = "\ngromacsHarmonicReferenceForce: ";
+   static int debug                      = true;
+   static int startAtom                  = 0;
+   static int stopAtom                   = 3;
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+
+   static const int  MaxFreeArrayIndex   = 20;
+
+   int indexOneDArraysToFree             = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree             = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree          = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get bond atom indices & parameters
+
+   int** harmonicBondAtomIndices;
+   RealOpenMM** harmonicBondParameters;
+   int numberOfBonds                                 = gromacsGetHarmonicBondParameters( top, &harmonicBondAtomIndices, &harmonicBondParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = harmonicBondAtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = harmonicBondParameters;
+
+   // bond energy array
+
+   RealOpenMM* harmonicBondBondEnergies              = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfBonds );
+   memset( harmonicBondBondEnergies, 0, sizeof( RealOpenMM )*numberOfBonds );
+   oneDArraysToFree[indexOneDArraysToFree++]         = harmonicBondBondEnergies;
+
+   ReferenceHarmonicBondIxn referenceHarmonicBondIxn;
+   ReferenceBondForce referenceBondForce;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " harmonic bonds=" << numberOfBonds;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // harmonic bond forces
+
+   referenceBondForce.calculateForce( numberOfBonds, harmonicBondAtomIndices, atomCoordinates,
+                                      harmonicBondParameters, harmonicBondForces,
+                                      harmonicBondBondEnergies, harmonicBondAtomEnergies, 
+                                      totalEnergy, referenceHarmonicBondIxn );
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " Harmonic Total energy=" << totalEnergy << std::endl;
+      for( int ii = startAtom; ii < stopAtom; ii++ ){
+         message << ii << " F[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, harmonicBondForces[ii], 3, one );
+         message << "] E=" << harmonicBondAtomEnergies[ii] << std::endl;
+      } 
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for angle bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetAngleBondParameters( const t_topology* top, int*** atomIndices,
+                                   RealOpenMM*** bondParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetAngleBondParameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const RealOpenMM two     = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef   = &top->idef;
+   int numberOfBonds    = idef->il[ F_ANGLES ].nr/4; 
+   int* atoms           = (int*) idef->il[ F_ANGLES ].iatoms;
+
+   gromacsAllocateBondIxnArrays( numberOfBonds, threeI, atomIndices, twoI, bondParameters );
+
+   int offsetIndex = 0;
+   for( int ii = 0; ii < numberOfBonds; ii++ ) { 
+
+      int type                 = atoms[ offsetIndex++ ];  
+      (*atomIndices)[ii][0]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][1]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][2]    = atoms[ offsetIndex++ ];
+
+      (*bondParameters)[ii][0] = idef->iparams[ type ].harmonic.rA;
+      (*bondParameters)[ii][1] = idef->iparams[ type ].harmonic.krA;
+   }
+
+   return numberOfBonds;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate angle bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param angleBondForces       output array of forces: angleBondForces[atomIndex][3]
+   @param angleBondEnergies     output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsAngleBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                    RealOpenMM** angleBondForces,
+                                    RealOpenMM* angleBondAtomEnergies,
+                                    RealOpenMM* totalEnergy ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName         = "\ngromacsAngleReferenceForce: ";
+   static int debug                      = true;
+   static int startAtom                  = 0;
+   static int stopAtom                   = 3;
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+
+   static const int  MaxFreeArrayIndex   = 20;
+
+   int indexOneDArraysToFree             = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree             = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree          = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get bond atom indices & parameters
+
+   int** angleBondAtomIndices;
+   RealOpenMM** angleBondParameters;
+   int numberOfBonds                                 = gromacsGetAngleBondParameters( top, &angleBondAtomIndices, &angleBondParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = angleBondAtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = angleBondParameters;
+
+   // bond energy array
+
+   RealOpenMM* angleBondBondEnergies                 = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfBonds );
+   memset( angleBondBondEnergies, 0, sizeof( RealOpenMM )*numberOfBonds );
+   oneDArraysToFree[indexOneDArraysToFree++]         = angleBondBondEnergies;
+
+   ReferenceAngleBondIxn referenceAngleBondIxn;
+   ReferenceBondForce referenceBondForce;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " angle bonds=" << numberOfBonds;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // angle bond forces
+
+   referenceBondForce.calculateForce( numberOfBonds, angleBondAtomIndices, atomCoordinates,
+                                      angleBondParameters, angleBondForces, 
+                                      angleBondBondEnergies, angleBondAtomEnergies, 
+                                      totalEnergy, referenceAngleBondIxn );
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << "Angle  Total energy=" << *totalEnergy << std::endl;
+      for( int ii = startAtom; ii < stopAtom; ii++ ){
+         message << ii << " F[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, angleBondForces[ii], 3, one );
+         message << "] E=" << angleBondAtomEnergies[ii] << std::endl;
+      } 
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   // ---------------------------------------------------------------------------------------
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for properDihedral bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetProperDihedralBondParameters( const t_topology* top, int*** atomIndices,
+                                            RealOpenMM*** bondParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetProperDihedralBondParameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const int fourI          = 4;
+   static const RealOpenMM two     = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef   = &top->idef;
+   int numberOfBonds    = idef->il[ F_PDIHS ].nr/5; 
+   int* atoms           = (int*) idef->il[ F_PDIHS ].iatoms;
+
+   gromacsAllocateBondIxnArrays( numberOfBonds, fourI, atomIndices, threeI, bondParameters );
+
+   int offsetIndex = 0;
+   for( int ii = 0; ii < numberOfBonds; ii++ ) { 
+
+      int type                 = atoms[ offsetIndex++ ];  
+      (*atomIndices)[ii][0]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][1]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][2]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][3]    = atoms[ offsetIndex++ ];
+
+      (*bondParameters)[ii][0] = idef->iparams[ type ].pdihs.cpA;
+      (*bondParameters)[ii][1] = idef->iparams[ type ].pdihs.phiA;
+      (*bondParameters)[ii][2] = (RealOpenMM) idef->iparams[ type ].pdihs.mult;
+   }
+
+   return numberOfBonds;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate properDihedral bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param properDihedralBondForces       output array of forces: angleBondForces[atomIndex][3]
+   @param properDihedralBondEnergies     output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsProperDihedralBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                             RealOpenMM** properDihedralBondForces,
+                                             RealOpenMM* properDihedralBondAtomEnergies,
+                                             RealOpenMM* totalEnergy ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName         = "\ngromacsProperDihedralReferenceForce: ";
+   static int debug                      = true;
+   static int startAtom                  = 0;
+   static int stopAtom                   = 3;
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+
+   static const int  MaxFreeArrayIndex   = 20;
+
+   int indexOneDArraysToFree             = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree             = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree          = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get bond atom indices & parameters
+
+   int** properDihedralBondAtomIndices;
+   RealOpenMM** properDihedralBondParameters;
+   int numberOfBonds                                 = gromacsGetProperDihedralBondParameters( top, &properDihedralBondAtomIndices, &properDihedralBondParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = properDihedralBondAtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = properDihedralBondParameters;
+
+   // bond energy array
+
+   RealOpenMM* properDihedralBondBondEnergies        = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfBonds );
+   memset( properDihedralBondBondEnergies, 0, sizeof( RealOpenMM )*numberOfBonds );
+   oneDArraysToFree[indexOneDArraysToFree++]         = properDihedralBondBondEnergies;
+
+   ReferenceProperDihedralBond referenceProperDihedralBond;
+   ReferenceBondForce referenceBondForce;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " properDihedral bonds=" << numberOfBonds << std::endl;
+/*
+      for( int ii  = 0; ii < numberOfBonds; ii++ ){
+         message << "   bnd=" << ii << " a[";
+         for( int jj  = 0; jj < 4; jj++ ){
+            message << properDihedralBondAtomIndices[ii][jj] << " ";
+         }
+         message << "] p[";
+         for( int jj  = 0; jj < 3; jj++ ){
+            message << properDihedralBondParameters[ii][jj] << " ";
+         }
+         message << "]" << std::endl; 
+      }
+*/
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // properDihedral bond forces
+
+   referenceBondForce.calculateForce( numberOfBonds, properDihedralBondAtomIndices, atomCoordinates,
+                                      properDihedralBondParameters, properDihedralBondForces, 
+                                      properDihedralBondBondEnergies, properDihedralBondAtomEnergies, 
+                                      totalEnergy, referenceProperDihedralBond);
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << "ProperDihedral  Total energy=" << *totalEnergy << std::endl;
+      for( int ii = startAtom; ii < stopAtom; ii++ ){
+         message << ii << " F[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, properDihedralBondForces[ii], 3, one );
+         message << "] E=" << properDihedralBondAtomEnergies[ii] << std::endl;
+      } 
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   // ---------------------------------------------------------------------------------------
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for RB Dihedral bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetRbDihedralBondParameters( const t_topology* top, int*** atomIndices,
+                                            RealOpenMM*** bondParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetRbDihedralBondParameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const int fourI          = 4;
+   static const int sixI           = 6;
+
+   static const RealOpenMM two     = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef   = &top->idef;
+   int numberOfBonds    = idef->il[ F_RBDIHS ].nr/5; 
+   int* atoms           = (int*) idef->il[ F_RBDIHS ].iatoms;
+
+   gromacsAllocateBondIxnArrays( numberOfBonds, fourI, atomIndices, sixI, bondParameters );
+
+   int offsetIndex = 0;
+   for( int ii = 0; ii < numberOfBonds; ii++ ) { 
+
+      int type                 = atoms[ offsetIndex++ ];  
+      (*atomIndices)[ii][0]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][1]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][2]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][3]    = atoms[ offsetIndex++ ];
+
+      for( int jj = 0; jj < sixI; jj++ ) { 
+         (*bondParameters)[ii][jj] = idef->iparams[ type ].rbdihs.rbc[jj];
+      }
+   }
+
+   return numberOfBonds;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate RB dihedral bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param rbDihedralBondForces     output array of forces: angleBondForces[atomIndex][3]
+   @param rbDihedralBondEnergies   output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsRbDihedralBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                         RealOpenMM** rbDihedralBondForces,
+                                         RealOpenMM* rbDihedralBondAtomEnergies,
+                                         RealOpenMM* totalEnergy ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName         = "\ngromacsRbDihedralReferenceForce: ";
+   static int debug                      = true;
+   static int startAtom                  = 0;
+   static int stopAtom                   = 3;
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+
+   static const int  MaxFreeArrayIndex   = 20;
+
+   int indexOneDArraysToFree             = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree             = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree          = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get bond atom indices & parameters
+
+   int** rbDihedralBondAtomIndices;
+   RealOpenMM** rbDihedralBondParameters;
+   int numberOfBonds                                 = gromacsGetRbDihedralBondParameters( top, &rbDihedralBondAtomIndices, &rbDihedralBondParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = rbDihedralBondAtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = rbDihedralBondParameters;
+
+   // bond energy array
+
+   RealOpenMM* rbDihedralBondBondEnergies            = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfBonds );
+   memset( rbDihedralBondBondEnergies, 0, sizeof( RealOpenMM )*numberOfBonds );
+   oneDArraysToFree[indexOneDArraysToFree++]         = rbDihedralBondBondEnergies;
+
+   ReferenceRbDihedralBond referenceRbDihedralBond;
+   ReferenceBondForce referenceBondForce;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " rbDihedral bonds=" << numberOfBonds << std::endl;
+/*
+      for( int ii  = 0; ii < numberOfBonds; ii++ ){
+         message << "   bnd=" << ii << " a[";
+         for( int jj  = 0; jj < 4; jj++ ){
+            message << rbDihedralBondAtomIndices[ii][jj] << " ";
+         }
+         message << "] p[";
+         for( int jj  = 0; jj < 3; jj++ ){
+            message << rbDihedralBondParameters[ii][jj] << " ";
+         }
+         message << "]" << std::endl; 
+      }
+*/
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // rbDihedral bond forces
+
+   referenceBondForce.calculateForce( numberOfBonds, rbDihedralBondAtomIndices, atomCoordinates,
+                                      rbDihedralBondParameters, rbDihedralBondForces, 
+                                      rbDihedralBondBondEnergies, rbDihedralBondAtomEnergies, 
+                                      totalEnergy, referenceRbDihedralBond);
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << "RbDihedral  Total energy=" << *totalEnergy << std::endl;
+      for( int ii = startAtom; ii < stopAtom; ii++ ){
+         message << ii << " F[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, rbDihedralBondForces[ii], 3, one );
+         message << "] E=" << rbDihedralBondAtomEnergies[ii] << std::endl;
+      } 
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   // ---------------------------------------------------------------------------------------
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for LJ 14 bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetLJ14Parameters( const t_topology* top, int*** atomIndices,
+                              RealOpenMM*** bondParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetLJ14Parameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const int fourI          = 4;
+   static const int sixI           = 6;
+
+   static const RealOpenMM two     = 2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef   = &top->idef;
+   int numberOfBonds    = idef->il[ F_LJ14 ].nr/3; 
+   int* atoms           = (int*) idef->il[ F_LJ14 ].iatoms;
+
+   gromacsAllocateBondIxnArrays( numberOfBonds, twoI, atomIndices, threeI, bondParameters );
+
+   int offsetIndex = 0;
+   for( int ii = 0; ii < numberOfBonds; ii++ ) { 
+
+      int type                 = atoms[ offsetIndex++ ];  
+      (*atomIndices)[ii][0]    = atoms[ offsetIndex++ ];
+      (*atomIndices)[ii][1]    = atoms[ offsetIndex++ ];
+
+      (*bondParameters)[ii][0] = idef->iparams[ type ].lj14.c6A;
+      (*bondParameters)[ii][1] = idef->iparams[ type ].lj14.c12A;
+   }
+
+   return numberOfBonds;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate LJ 14 forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param fr                       Gromacs t_forcerec struct
+   @param lj14Forces               output array of forces: lj14Forces[atomIndex][3]
+   @param lj14AtomEnergies         output array of LJ 14 energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsLJ14ReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                               t_forcerec *fr, const RealOpenMM* charges, RealOpenMM** lj14Forces,
+                               RealOpenMM* lj14AtomEnergies,
+                               RealOpenMM* totalEnergy ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName         = "\ngromacsLJ14ReferenceForce: ";
+   static int debug                      = true;
+   static int startAtom                  = 0;
+   static int stopAtom                   = 3;
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+
+   static const int  MaxFreeArrayIndex   = 20;
+
+   int indexOneDArraysToFree             = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree             = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree          = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get bond atom indices & parameters
+
+   int** lj14AtomIndices;
+   RealOpenMM** lj14Parameters;
+   int numberOfBonds                                 = gromacsGetLJ14Parameters( top, &lj14AtomIndices, &lj14Parameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = lj14AtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = lj14Parameters;
+
+   // bond energy array
+
+   RealOpenMM* lj14BondEnergies                      = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfBonds );
+   memset( lj14BondEnergies, 0, sizeof( RealOpenMM )*numberOfBonds );
+   oneDArraysToFree[indexOneDArraysToFree++]         = lj14BondEnergies;
+
+   ReferenceLJ14 referenceLJ14;
+   ReferenceBondForce referenceBondForce;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " lj 14 bonds=" << numberOfBonds << std::endl;
+/*
+      for( int ii  = 0; ii < numberOfBonds; ii++ ){
+         message << "   bnd=" << ii << " a[";
+         for( int jj  = 0; jj < 4; jj++ ){
+            message << lj14AtomIndices[ii][jj] << " ";
+         }
+         message << "] p[";
+         for( int jj  = 0; jj < 3; jj++ ){
+            message << lj14Parameters[ii][jj] << " ";
+         }
+         message << "]" << std::endl; 
+      }
+*/
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // calculate derived parameters
+
+   for( int ii  = 0; ii < numberOfBonds; ii++ ){
+      referenceLJ14.getDerivedParameters( lj14Parameters[ii][0], lj14Parameters[ii][1],
+                                          charges[lj14AtomIndices[ii][0]],
+                                          charges[lj14AtomIndices[ii][1]], fr->epsfac*fr->fudgeQQ, lj14Parameters[ii] );
+   }
+
+   // lj14 bond forces
+
+   referenceBondForce.calculateForce( numberOfBonds, lj14AtomIndices, atomCoordinates,
+                                      lj14Parameters, lj14Forces, 
+                                      lj14BondEnergies, lj14AtomEnergies, 
+                                      totalEnergy, referenceLJ14);
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << "RbDihedral  Total energy=" << *totalEnergy << std::endl;
+      for( int ii = startAtom; ii < stopAtom; ii++ ){
+         message << ii << " F[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, lj14Forces[ii], 3, one );
+         message << "] E=" << lj14AtomEnergies[ii] << std::endl;
+      } 
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   // ---------------------------------------------------------------------------------------
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for LJ Coulomb bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetLJCoulombParameters( const t_topology* top, t_forcerec *fr, const t_mdatoms *md,
+                                   int*** ljCoulombExclusions, RealOpenMM*** outputAtomParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetLJCoulombParameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const int fourI          = 4;
+   static const int sixI           = 6;
+
+   static const RealOpenMM two     = 2.0;
+
+   int numberOfAtoms, exclusionBlockSize, numberOfParameters;
+   int** exclusions;
+   int* exclusionBlock;
+
+   RealOpenMM** atomParameters;
+   RealOpenMM* atomParametersBlock;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_block *gmx_excl             = &top->atoms.excl;
+
+   int ntypes                          = fr->ntype;
+   int* types                          = md->typeA;
+   real* nbfp                          = fr->nbfp;
+   
+   // ---------------------------------------------------------------------------------------
+
+   numberOfAtoms                       = top->atoms.nr;
+   numberOfParameters                  = 4;
+
+   // allocate & initialize memory for exclusions and parameters
+
+   exclusionBlockSize                  = gmx_excl->index[numberOfAtoms] - gmx_excl->index[0] + numberOfAtoms;
+   exclusions                          = (int**) malloc( sizeof( int* )*numberOfAtoms );
+   exclusionBlock                      = (int*)  malloc( sizeof( int  )*exclusionBlockSize );
+   memset( exclusionBlock, 0, sizeof( int  )*exclusionBlockSize );
+   *ljCoulombExclusions                = exclusions;
+
+   atomParameters                      = (RealOpenMM**) malloc( sizeof( RealOpenMM* )*numberOfAtoms );
+   atomParametersBlock                 = (RealOpenMM*)  malloc( sizeof( RealOpenMM )*numberOfAtoms*numberOfParameters );
+   memset( atomParametersBlock, 0, sizeof( RealOpenMM )*numberOfAtoms*numberOfParameters );
+   *outputAtomParameters               = atomParameters;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){ 
+
+      exclusions[ii]        = exclusionBlock;
+      exclusionBlock       += (gmx_excl->index[ii+1] - gmx_excl->index[ii] + 1);
+
+      atomParameters[ii]    = atomParametersBlock;
+      atomParametersBlock  += numberOfParameters;
+
+   }
+   
+   // ---------------------------------------------------------------------------------------
+
+   // load parameters
+
+   RealOpenMM* charges = md->chargeA;
+   for( int ii = 0; ii < numberOfAtoms; ii++ ) { 
+
+      exclusions[ii][0]        = gmx_excl->index[ii+1] - gmx_excl->index[ii];
+      int offsetIndex          = 1;
+      for ( int jj = gmx_excl->index[ii]; jj < gmx_excl->index[ii+1]; jj++ ) {
+         exclusions[ii][offsetIndex++] = gmx_excl->a[jj];
+      }    
+ 
+      atomParameters[ii][0]     = nbfp[ types[ii]*2*ntypes + types[ii]*2 ];
+      atomParameters[ii][1]     = nbfp[ types[ii]*2*ntypes + types[ii]*2 + 1 ]; 
+      atomParameters[ii][2]     = charges[ii]; 
+
+   }
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate LJ Coulomb forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param fr                       Gromacs t_forcerec struct
+   @param md                       Gromacs t_mdatoms struct
+   @param ljCoulombForces          output array of forces: ljCoulombForces[atomIndex][3]
+   @param ljCoulombAtomEnergies    output array of LJ Coulomb energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsLJCoulombReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                    t_forcerec *fr, const t_mdatoms *md, RealOpenMM** ljCoulombForces,
+                                    RealOpenMM* ljCoulombAtomEnergies,
+                                    RealOpenMM* totalEnergy ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName         = "\ngromacsLJCoulombReferenceForce: ";
+   static int debug                      = true;
+   static int startAtom                  = 0;
+   static int stopAtom                   = 3;
+
+   static const RealOpenMM zero          = 0.0;
+   static const RealOpenMM one           = 1.0;
+
+   static const int  MaxFreeArrayIndex   = 20;
+
+   int indexOneDArraysToFree             = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree             = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   // currently not used!
+
+   RealOpenMM fixedParameters[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree          = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get exclusions & parameters
+
+   int** ljCoulombExclusions;
+   RealOpenMM** ljCoulombParameters;
+   int numberOfBonds                                 = gromacsGetLJCoulombParameters( top, fr, md, &ljCoulombExclusions, &ljCoulombParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = ljCoulombExclusions;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = ljCoulombParameters;
+
+   ReferenceLJCoulombIxn referenceLJCoulomb;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " lj Coulomb atoms=" << numberOfAtoms << std::endl;
+/*
+      for( int ii  = 0; ii < numberOfAtoms; ii++ ){
+
+         message << "   " << ii << " P[";
+         for( int jj  = 0; jj < 3; jj++ ){
+            message << ljCoulombParameters[ii][jj] << " ";
+         }
+         message << "]";
+
+         message << "  X " << ljCoulombExclusions[ii][0] << " [";
+         for( int jj  = 1; jj <= ljCoulombExclusions[ii][0]; jj++ ){
+            message << ljCoulombExclusions[ii][jj] << " ";
+         }
+         message << "]" << std::endl; 
+      }
+*/
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // calculate derived parameters
+
+   FILE* paramFile = fopen( "LJCoulombDerived.txt", "w" );
+   RealOpenMM epsFacSqrt = SQRT( fr->epsfac );
+   for( int ii  = 0; ii < numberOfAtoms; ii++ ){
+      (void) fprintf( paramFile, "%d %12.5e %12.5e %12.5e", ii, ljCoulombParameters[ii][0], ljCoulombParameters[ii][1], ljCoulombParameters[ii][2] );
+      referenceLJCoulomb.getDerivedParameters( ljCoulombParameters[ii][0], ljCoulombParameters[ii][1],
+                                               ljCoulombParameters[ii][2],
+                                               epsFacSqrt, ljCoulombParameters[ii] );
+      (void) fprintf( paramFile, "   %12.5e %12.5e %12.5e\n", ljCoulombParameters[ii][0], ljCoulombParameters[ii][1], ljCoulombParameters[ii][2] );
+   }
+   (void) fclose( paramFile );
+
+   // ljCoulomb bond forces
+
+   referenceLJCoulomb.calculatePairIxn( numberOfAtoms, atomCoordinates,
+                                        ljCoulombParameters, ljCoulombExclusions,
+                                        fixedParameters, ljCoulombForces, 
+                                        ljCoulombAtomEnergies, totalEnergy );
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << "LJCoulomb Total energy=" << *totalEnergy << std::endl;
+      for( int ii = startAtom; ii < stopAtom; ii++ ){
+         message << ii << " F[";
+         SimTKOpenMMUtilities::formatRealStringStream( message, ljCoulombForces[ii], 3, one );
+         message << "] E=" << ljCoulombAtomEnergies[ii] << std::endl;
+      } 
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   // ---------------------------------------------------------------------------------------
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate forces for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param gromacAtomCoordinates    atom configuration
+   @param md                       Gromacs t_mdatoms data struct
+   @param partialChargesIn         array of partial charges
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsReferenceForce( const t_topology* top, const rvec* gromacAtomCoordinates,
+                           const t_mdatoms *md, t_forcerec* fr,
+                           char* baseFileName, FILE* log ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName               = "\ngromacsReferenceForce: ";
+   static int debug                            = true;
+   static int startAtom                        = 0;
+   static int stopAtom                         = 3;
+
+   static const RealOpenMM zero                = 0.0;
+   static const RealOpenMM one                 = 1.0;
+
+   static const int  zeroI                     = 0;
+   static const int  twoI                      = 2;
+   static const int  threeI                    = 3;
+   static const int  fourI                     = 4;
+
+   static const int  MaxFreeArrayIndex         = 20;
+
+   static const int  includeHarmonic           = 1;
+   static const int  includeAngle              = 1;
+   static const int  includeProperDihedral     = 1;
+   static const int  includeRbDihedral         = 1;
+   static const int  includeLJ14               = 1;
+   static const int  includeLJCoulomb          = 1;
+
+   static const int  outputParameterFile       = 0;
+   static const int  outputForceFile           = 1;
+
+   RealOpenMM totalEnergy                      = zero;
+
+   int indexOneDArraysToFree                   = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree                   = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree                = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   // set log file, if available
+
+   if( log ){
+      SimTKOpenMMLog::setSimTKOpenMMLog( log );
+   }
+
+   int numberOfAtoms = top->atoms.nr;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " atoms=" << numberOfAtoms;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   RealOpenMM** atomCoordinates              = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = atomCoordinates;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      atomCoordinates[ii][0] = gromacAtomCoordinates[ii][0];
+      atomCoordinates[ii][1] = gromacAtomCoordinates[ii][1];
+      atomCoordinates[ii][2] = gromacAtomCoordinates[ii][2];
+   }
+
+   // accumulate bonded force contributions
+
+   RealOpenMM** bondedForces                 = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = bondedForces;
+
+   // ---------------------------------------------------------------------------------------
+
+   // harmonic bond forces
+
+   if( includeHarmonic ){
+      RealOpenMM** harmonicBondForces                   = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]         = harmonicBondForces;
+   
+      RealOpenMM* harmonicBondEnergies                  = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      memset( harmonicBondEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]         = harmonicBondEnergies;
+   
+      gromacsHarmonicBondReferenceForce( top, atomCoordinates, harmonicBondForces, harmonicBondEnergies, &totalEnergy ); 
+
+      if( debug ){
+         std::stringstream message;
+         message << methodName;
+         message << " Bond Total energy=" << totalEnergy << std::endl;
+         for( int ii = startAtom; ii < stopAtom; ii++ ){
+            message << ii << " F[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, harmonicBondForces[ii], 3, one );
+            message << "] E=" << harmonicBondEnergies[ii] << std::endl;
+         } 
+         SimTKOpenMMLog::printMessage( message );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream bondForceFileName;
+         bondForceFileName << baseFileName;
+         bondForceFileName << "HarmonicBond.txt";
+         ReferenceForce::writeForces( numberOfAtoms, twoI, atomCoordinates, harmonicBondForces,
+                                      harmonicBondEnergies, bondForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( bondedForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, harmonicBondForces, bondedForces );
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // angle bond forces
+
+   if( includeAngle ){
+
+      RealOpenMM** angleBondForces                      = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]         = angleBondForces;
+   
+      RealOpenMM* angleBondAtomEnergies                 = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      memset( angleBondAtomEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]         = angleBondAtomEnergies;
+   
+      gromacsAngleBondReferenceForce( top, atomCoordinates, angleBondForces, angleBondAtomEnergies, &totalEnergy ); 
+   
+      if( debug ){
+         std::stringstream message;
+         message << methodName;
+         message << " Angle Total energy=" << totalEnergy << std::endl;
+         for( int ii = startAtom; ii < stopAtom; ii++ ){
+            message << ii << " F[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, angleBondForces[ii], 3, one );
+            message << "] E=" << angleBondAtomEnergies[ii] << std::endl;
+         } 
+         SimTKOpenMMLog::printMessage( message );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream angleForceFileName;
+         angleForceFileName << baseFileName;
+         angleForceFileName << "AngleBond.txt";
+         ReferenceForce::writeForces( numberOfAtoms, threeI, atomCoordinates, angleBondForces,
+                                      angleBondAtomEnergies, angleForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( bondedForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, angleBondForces, bondedForces );
+      }
+
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // properDihedral bond forces
+
+   if( includeProperDihedral ){
+
+      RealOpenMM** properDihedralBondForces             = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]         = properDihedralBondForces;
+   
+      RealOpenMM* properDihedralBondAtomEnergies        = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      memset( properDihedralBondAtomEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]         = properDihedralBondAtomEnergies;
+   
+      gromacsProperDihedralBondReferenceForce( top, atomCoordinates, properDihedralBondForces, properDihedralBondAtomEnergies, &totalEnergy ); 
+   
+      if( debug ){
+         std::stringstream message;
+         message << methodName;
+         message << " ProperDihedral Total energy=" << totalEnergy << std::endl;
+         for( int ii = startAtom; ii < stopAtom; ii++ ){
+            message << ii << " F[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, properDihedralBondForces[ii], 3, one );
+            message << "] E=" << properDihedralBondAtomEnergies[ii] << std::endl;
+         } 
+         SimTKOpenMMLog::printMessage( message );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream properDihedralForceFileName;
+         properDihedralForceFileName << baseFileName;
+         properDihedralForceFileName << "ProperDihedral.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, properDihedralBondForces,
+                                      properDihedralBondAtomEnergies, properDihedralForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( bondedForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, properDihedralBondForces, bondedForces );
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // rbDihedral bond forces
+
+   if( includeRbDihedral ){
+
+      RealOpenMM** rbDihedralBondForces             = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]     = rbDihedralBondForces;
+   
+      RealOpenMM* rbDihedralBondAtomEnergies        = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      memset( rbDihedralBondAtomEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]     = rbDihedralBondAtomEnergies;
+   
+      gromacsRbDihedralBondReferenceForce( top, atomCoordinates, rbDihedralBondForces, rbDihedralBondAtomEnergies, &totalEnergy ); 
+  
+      if( debug ){
+         std::stringstream message;
+         message << methodName;
+         message << " RbDihedral Total energy=" << totalEnergy << std::endl;
+         for( int ii = startAtom; ii < stopAtom; ii++ ){
+            message << ii << " F[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, rbDihedralBondForces[ii], 3, one );
+            message << "] E=" << rbDihedralBondAtomEnergies[ii] << std::endl;
+         } 
+         SimTKOpenMMLog::printMessage( message );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream rbDihedralForceFileName;
+         rbDihedralForceFileName << baseFileName;
+         rbDihedralForceFileName << "RbDihedral.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, rbDihedralBondForces,
+                                      rbDihedralBondAtomEnergies, rbDihedralForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( bondedForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, rbDihedralBondForces, bondedForces );
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // lj14 bond forces
+
+   if( includeLJ14 ){
+
+      RealOpenMM** lj14BondForces                   = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]     = lj14BondForces;
+   
+      RealOpenMM* lj14BondAtomEnergies              = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      memset( lj14BondAtomEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]     = lj14BondAtomEnergies;
+   
+      gromacsLJ14ReferenceForce( top, atomCoordinates, fr, md->chargeA, lj14BondForces, lj14BondAtomEnergies, &totalEnergy ); 
+  
+      if( debug ){
+         std::stringstream message;
+         message << methodName;
+         message << " LJ14 Total energy=" << totalEnergy << std::endl;
+         for( int ii = startAtom; ii < stopAtom; ii++ ){
+            message << ii << " F[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, lj14BondForces[ii], 3, one );
+            message << "] E=" << lj14BondAtomEnergies[ii] << std::endl;
+         } 
+         SimTKOpenMMLog::printMessage( message );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream lj14ForceFileName;
+         lj14ForceFileName << baseFileName;
+         lj14ForceFileName << "LJ14.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, lj14BondForces,
+                                      lj14BondAtomEnergies, lj14ForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( bondedForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, lj14BondForces, bondedForces );
+      }
+   }
+
+   // output file
+
+   if( outputForceFile && bondedForces ){
+      std::stringstream bondedForceFileName;
+      bondedForceFileName << baseFileName;
+      bondedForceFileName << "Bonded.txt";
+      ReferenceForce::writeForces( numberOfAtoms, zeroI, atomCoordinates, bondedForces,
+                                   NULL, bondedForceFileName.str(), one, one, one );
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // ljCoulomb bond forces
+
+   if( includeLJCoulomb ){
+
+      RealOpenMM** ljCoulombBondForces              = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]     = ljCoulombBondForces;
+   
+      RealOpenMM* ljCoulombBondAtomEnergies         = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      memset( ljCoulombBondAtomEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]     = ljCoulombBondAtomEnergies;
+   
+      gromacsLJCoulombReferenceForce( top, atomCoordinates, fr, md, ljCoulombBondForces, ljCoulombBondAtomEnergies, &totalEnergy ); 
+  
+      if( debug ){
+         std::stringstream message;
+         message << methodName;
+         message << " LJCoulomb Total energy=" << totalEnergy << std::endl;
+         for( int ii = startAtom; ii < stopAtom; ii++ ){
+            message << ii << " F[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, ljCoulombBondForces[ii], 3, one );
+            message << "] E=" << ljCoulombBondAtomEnergies[ii] << std::endl;
+         } 
+         SimTKOpenMMLog::printMessage( message );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream ljCoulombForceFileName;
+         ljCoulombForceFileName << baseFileName;
+         ljCoulombForceFileName << "LJCoulomb.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, ljCoulombBondForces,
+                                      ljCoulombBondAtomEnergies, ljCoulombForceFileName.str(), one, one, one );
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get Shake parameters and atom indices (simple version)
+
+   @param top                      Gromacs t_topology data struct
+   @param md                       Gromacs t_mdatoms data struct
+   @param outputAtomIndices        array of atomIndices -- see below
+   @param outputShakeParameters    array of Shake parameters -- see below
+
+   The arrays atomIndices & outputShakeParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of blocks)*(second dimension)
+
+   The memory should be freed by the callee via calls:
+		free( outputAtomIndices[0] ); free( outputAtomIndices );
+		free( outputShakeParameters[0] ); free( outputShakeParameters );
+
+
+   outputAtomIndices has the following structure:
+
+   outputAtomIndices[blockIndex][i], where
+
+			i = 0:   index of atom1
+			i = 1:   index of atom2
+
+   outputShakeParameters has the following structure:
+
+   outputShakeParameters[blockIndex][i], where
+
+			i = 0:  constraint distance
+
+   @return number of constraints
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetShakeParameters( const t_topology* top, const t_mdatoms *md,
+                               int*** outputAtomIndices, RealOpenMM*** outputShakeParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetShakeParameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const int fourI          = 4;
+   static const int fiveI          = 5;
+   static const int sixI           = 6;
+
+   static const int debug          = 1;
+
+   static const RealOpenMM zero    = 0.0;
+   static const RealOpenMM one     = 1.0;
+   static const RealOpenMM two     = 2.0;
+   static const RealOpenMM epsilon = (RealOpenMM) 1.0e-05;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef           = &top->idef;
+   const t_iatom* atoms         = idef->il[ F_SHAKE ].iatoms;
+   int numberOfConstraints      = idef->il[ F_SHAKE ].nr/3;
+
+   gromacsAllocateBondIxnArrays( numberOfConstraints, twoI, outputAtomIndices, twoI, outputShakeParameters );
+   int** atomIndices            = *outputAtomIndices;
+   RealOpenMM** shakeParameters = *outputShakeParameters;
+
+   // loop over constraints 
+
+   int offset     = 0;
+   for( int ii = 0; ii < numberOfConstraints; ii++ ) { 
+
+      int type                    = atoms[ offset     ];  
+      int atomI                   = atoms[ offset + 1 ];  
+      int atomJ                   = atoms[ offset + 2 ];  
+
+      offset                     += 3;
+
+      atomIndices[ii][0]          = atomI;
+      atomIndices[ii][1]          = atomJ;
+
+      shakeParameters[ii][0]      = idef->iparams[type].shake.dA;
+      shakeParameters[ii][1]      = zero; // unused
+   }
+
+   return numberOfConstraints;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get Shake parameters and atom indices in blocks
+
+   @param top                      Gromacs t_topology data struct
+   @param md                       Gromacs t_mdatoms data struct
+   @param outputAtomIndices        array of atomIndices -- see below
+   @param outputShakeParameters    array of Shake parameters -- see below
+
+   The arrays atomIndices & outputShakeParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of blocks)*(second dimension [2 or 5] )
+
+   The memory should be freed by the callee via calls:
+		free( outputAtomIndices[0] ); free( outputAtomIndices );
+		free( outputShakeParameters[0] ); free( outputShakeParameters );
+
+
+
+   outputAtomIndices has the following structure:
+
+   outputAtomIndices[blockIndex][i], where
+
+			i = 0:   number of atoms in block (max is 4)
+			i = 1:   index of heavy atom in constraint block
+			i = 2-4: index of light atoms in constraint block
+
+   If max > 4, then program aborts
+
+
+
+   outputShakeParameters has the following structure:
+
+   outputShakeParameters[blockIndex][i], where
+
+			i = 0:  inverse mass of heavy atom
+			i = 2:  1.0/2*[ inverse mass of heavy atom + inverse mass of light atom ]
+			i = 3:  distance constraint**2
+
+   The following assumptions are made:
+
+			all light atoms have same mass
+			constraint distance is the same for all constraints in a block
+
+   If these are violated, the program should abort
+
+
+   @return number of constraint blocks
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsGetBlockShakeParameters( const t_topology* top, const t_mdatoms *md,
+                                    int*** outputAtomIndices, RealOpenMM*** outputShakeParameters ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName   = "\ngromacsGetBlockShakeParameters: ";
+   static const int twoI           = 2;
+   static const int threeI         = 3;
+   static const int fourI          = 4;
+   static const int fiveI          = 5;
+   static const int sixI           = 6;
+
+   static const int debug          = 1;
+
+   static const RealOpenMM zero    = 0.0;
+   static const RealOpenMM one     = 1.0;
+   static const RealOpenMM two     = 2.0;
+   static const RealOpenMM epsilon = (RealOpenMM) 1.0e-05;
+
+   // ---------------------------------------------------------------------------------------
+
+   const t_idef* idef          = &top->idef;
+   const t_iatom* atoms        = idef->il[ F_SHAKE ].iatoms;
+   int numberOfConstraints     = idef->il[ F_SHAKE ].nr/3;
+
+   // loop over constraints 
+   // create a ReferenceShakeConstraint object w/ info for that constraint
+   // add the constraint to a map of constraints based on the index of the heavy atom
+   //      shakeMap[heavyAtomIndex] = Vector{ ReferenceShakeConstraints } containing that atom
+
+   int offset = 0;
+   IntShakeMap shakeMap;
+   for( int ii = 0; ii < numberOfConstraints; ii++ ) { 
+
+      int type       = atoms[ offset     ];  
+      int atomI      = atoms[ offset + 1 ];  
+      int atomJ      = atoms[ offset + 2 ];  
+      offset        += 3;
+      ReferenceShakeConstraint* referenceShakeConstraint = new ReferenceShakeConstraint( atomI, atomJ, idef->iparams[type].shake.dA,
+                                                                                         md->invmass[atomI], md->invmass[atomJ] );
+      int heavyAtomIndex = referenceShakeConstraint->getHeavyAtomIndex();
+
+      IntShakeMapI shakeMapI= shakeMap.find( heavyAtomIndex );
+      ShakeVector* shakeVector;
+      if( shakeMapI == shakeMap.end() ){
+         shakeVector              = new ShakeVector();
+         shakeMap[heavyAtomIndex] = shakeVector;
+      } else {
+         shakeVector = (*shakeMapI).second;
+      }
+      shakeVector->push_back( referenceShakeConstraint );
+   }
+
+   // validate blocks
+   //    (1) at most 3 elements/block
+   //    (2) distances the same
+   //    (3) masses the same for light atoms
+
+   int numberOfConstraintBlocks = (int) shakeMap.size();
+   int blockCount               = 1;
+   int errors                   = 0;
+   std::stringstream message;
+   message << "Number of Shake constraint blocks=" << numberOfConstraintBlocks;
+   for( IntShakeMapI ii = shakeMap.begin(); ii != shakeMap.end(); ii++ ){ 
+      ShakeVector* shakeVector = (*ii).second;
+      message << "\nBlock " << blockCount++;
+      if( shakeVector->size() > 3 ){
+         message << "\n      ERROR: Block has too many elements(=" << shakeVector->size() << std::endl;
+         errors++;
+      }
+      RealOpenMM blockDistance     = -one;
+      RealOpenMM blockInverseMass  = -one;
+      for( ShakeVectorI jj = shakeVector->begin(); jj != shakeVector->end(); jj++ ) { 
+
+         (*jj)->printState( message );
+         message << "   ";
+
+         // distance check
+
+         if( blockDistance < zero ){
+            blockDistance = (*jj)->getConstraintDistance();
+         } else {
+            RealOpenMM constraintDistance = (*jj)->getConstraintDistance();
+            RealOpenMM diff               = FABS( blockDistance - constraintDistance );
+            if( diff > epsilon ){
+               errors++;
+               message << "\n      ERROR: distances not equal for block[" << blockDistance << " " << constraintDistance << "] diff=" << diff << std::endl;
+            }
+         }
+
+         // mass check
+
+         if( blockInverseMass < zero ){
+            blockInverseMass = (*jj)->getLightAtomInverseMass();
+         } else {
+            RealOpenMM lightMass          = (*jj)->getLightAtomInverseMass();
+            RealOpenMM diff               = FABS( blockInverseMass - lightMass );
+            if( diff > epsilon ){
+               errors++;
+               message << "\n      ERROR: light atom masses not equal for block[" << blockInverseMass << " " << lightMass << "] diff=" << diff << std::endl;
+            }
+         }
+      }
+   }
+
+   // abort if errors
+
+   if( errors ){
+      message << "\n\nErrors=" << errors << " -- aborting.";
+      SimTKOpenMMLog::printError( message );
+   } else {
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // allocate & initialize memory for Shake atom indices & parameters
+   // and load data into arrays
+
+   gromacsAllocateBondIxnArrays( numberOfConstraintBlocks, fiveI, outputAtomIndices, threeI, outputShakeParameters );
+   int** atomIndices            = *outputAtomIndices;
+   RealOpenMM** shakeParameters = *outputShakeParameters;
+
+   blockCount = 0;
+   for( IntShakeMapI ii = shakeMap.begin(); ii != shakeMap.end(); ii++, blockCount++ ){ 
+      ShakeVector* shakeVector   = (*ii).second;
+      atomIndices[blockCount][0] = (int) (shakeVector->size() + 1);
+      int lightAtomCount         = 2;
+      for( ShakeVectorI jj = shakeVector->begin(); jj != shakeVector->end(); jj++ ) { 
+
+         if( lightAtomCount == 2 ){
+            atomIndices[blockCount][1]                = (*jj)->getHeavyAtomIndex();
+
+            shakeParameters[blockCount][0]            = (*jj)->getHeavyAtomInverseMass();
+   
+            shakeParameters[blockCount][1]            = two*( (*jj)->getHeavyAtomInverseMass() +  (*jj)->getLightAtomInverseMass() );
+            shakeParameters[blockCount][1]            = one/shakeParameters[blockCount][1];
+   
+            shakeParameters[blockCount][2]            = (*jj)->getConstraintDistance()* (*jj)->getConstraintDistance();
+         }
+         atomIndices[blockCount][lightAtomCount++]    = (*jj)->getLightAtomIndex();
+      }
+   }
+         
+   // free allocated memory
+
+   for( IntShakeMapI ii = shakeMap.begin(); ii != shakeMap.end(); ii++ ) { 
+      ShakeVector* shakeVector = (*ii).second;
+      for( ShakeVectorI jj = shakeVector->begin(); jj != shakeVector->end(); jj++ ) { 
+         delete *jj;
+      }
+     delete shakeVector;
+   }
+
+   return numberOfConstraintBlocks;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Update coordinates using stochastic dynamics integrator
+
+   @param top                      Gromacs t_topology data struct
+   @param gromacAtomCoordinates    atom configuration
+   @param md                       Gromacs t_mdatoms data struct
+   @param fr                       Gromacs t_forcerec data struct
+   @param gromacsVelocities        velocities
+   @param gromacsForces            forces
+   @param masses                   masses
+   @param deltaT                   delta t
+   @param tau                      viscosity(?)
+   @param temperature              temperature
+   @param baseFileName             base file name
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsReferenceSdUpdate( int numberOfAtoms, const t_topology* top, const rvec* gromacAtomCoordinates,
+                              const t_mdatoms *md, t_forcerec* fr,
+                              const rvec* gromacsVelocities, const rvec* gromacsForces,
+                              real* masses, real deltaT, real tau, real temperature,
+                              char* baseFileName, FILE* log ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName               = "\ngromacsReferenceSdUpdate: ";
+   static int debug                            = true;
+
+   static const RealOpenMM zero                = 0.0;
+   static const RealOpenMM one                 = 1.0;
+
+   static const int  zeroI                     = 0;
+   static const int  twoI                      = 2;
+   static const int  threeI                    = 3;
+   static const int  fourI                     = 4;
+
+   static const int  MaxFreeArrayIndex         = 20;
+
+   static const int  outputParameterFile       = 0;
+   static const int  outputForceFile           = 1;
+
+   RealOpenMM totalEnergy                      = zero;
+
+   int indexOneDArraysToFree                   = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree                   = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree                = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   // set log file, if available
+
+   if( log ){
+      SimTKOpenMMLog::setSimTKOpenMMLog( log );
+   }
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " atoms=" << numberOfAtoms;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // load coordinates, velocities, and forces
+
+   RealOpenMM** atomCoordinates              = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = atomCoordinates;
+
+   RealOpenMM** velocities                   = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = velocities;
+
+   RealOpenMM** forces                       = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = forces;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      atomCoordinates[ii][0] = gromacAtomCoordinates[ii][0];
+      atomCoordinates[ii][1] = gromacAtomCoordinates[ii][1];
+      atomCoordinates[ii][2] = gromacAtomCoordinates[ii][2];
+
+      velocities[ii][0]      = gromacsVelocities[ii][0];
+      velocities[ii][1]      = gromacsVelocities[ii][1];
+      velocities[ii][2]      = gromacsVelocities[ii][2];
+
+      forces[ii][0]          = gromacsForces[ii][0];
+      forces[ii][1]          = gromacsForces[ii][1];
+      forces[ii][2]          = gromacsForces[ii][2];
+
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // perform update
+
+   // create referenceStochasticDynamics object
+   // get Shake parameters
+
+   ReferenceStochasticDynamics referenceStochasticDynamics( numberOfAtoms, deltaT, tau, temperature ); 
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName << " Parameters:\n   ";
+      referenceStochasticDynamics.printParameters( message );
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   int** shakeAtomIndices;
+   RealOpenMM** shakeParameters;
+   //int numberOfConstraintBlocks                      = gromacsGetBlockShakeParameters( top, md, &shakeAtomIndices, &shakeParameters );
+   int numberOfConstraints                           = gromacsGetShakeParameters( top, md, &shakeAtomIndices, &shakeParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = shakeAtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = shakeParameters;
+
+   if( debug ){
+/*
+      // block constraints
+
+      std::stringstream message;
+      message << methodName << " Shake Parameters: " << numberOfConstraintBlocks << "\n";
+      for( int ii = 0; ii < numberOfConstraintBlocks; ii++ ){
+         message << " " << (ii+1) << " " << shakeParameters[ii][0] << " " << shakeParameters[ii][1] << " " << shakeParameters[ii][2] << " ";
+         message << shakeAtomIndices[ii][0] << " [";
+         for( int jj = 0; jj < shakeAtomIndices[ii][0]; jj++ ){
+            message << shakeAtomIndices[ii][jj+1] << " ";
+         }
+         message << "]\n";
+      }
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+*/
+
+      
+      std::stringstream message;
+      message << methodName << " Shake Parameters: " << numberOfConstraints << "\n";
+      for( int ii = 0; ii < numberOfConstraints && ii < 10; ii++ ){
+         message << " " << (ii+1) << " " << shakeParameters[ii][0] << " [" << shakeAtomIndices[ii][0] << " " << shakeAtomIndices[ii][1] << "]\n";
+      }
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // add Shake to dynamics
+
+   ReferenceShakeAlgorithm referenceShakeAlgorithm( numberOfConstraints, shakeAtomIndices, shakeParameters );
+   referenceStochasticDynamics.setReferenceShakeAlgorithm( &referenceShakeAlgorithm );
+
+	referenceStochasticDynamics.writeState( numberOfAtoms, atomCoordinates,
+                                           velocities, forces, masses,
+                                           zeroI, baseFileName );
+
+   referenceStochasticDynamics.update( numberOfAtoms, atomCoordinates, velocities, forces, masses );
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName << " Called update: " << referenceStochasticDynamics.getTimeStep() << "\n";
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+	referenceStochasticDynamics.writeState( numberOfAtoms, atomCoordinates,
+                                           velocities, forces, masses,
+                                           zeroI, baseFileName );
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Update coordinates using stochastic dynamics integrator
+
+   @param top                      Gromacs t_topology data struct
+   @param gromacAtomCoordinates    atom configuration
+   @param md                       Gromacs t_mdatoms data struct
+   @param fr                       Gromacs t_forcerec data struct
+   @param gromacsVelocities        velocities
+   @param gromacsForces            forces
+   @param masses                   masses
+   @param timesteps                number of timesteps
+   @param deltaT                   delta t
+   @param tau                      viscosity(?)
+   @param temperature              temperature
+   @param baseFileName             base file name
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsReferenceSimulation( int numberOfAtoms, const t_topology* top, const rvec* gromacAtomCoordinates,
+                                const t_mdatoms *md, t_forcerec* fr,
+                                const rvec* gromacsVelocities, const rvec* gromacsForces,
+                                real* masses, int timesteps, real deltaT, real tau, real temperature,
+                                char* baseFileName, FILE* log ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName               = "\ngromacsReferenceSimulation: ";
+   static int debug                            = false;
+
+   static const RealOpenMM zero                = 0.0;
+   static const RealOpenMM one                 = 1.0;
+
+   static const int  zeroI                     = 0;
+   static const int  twoI                      = 2;
+   static const int  threeI                    = 3;
+   static const int  fourI                     = 4;
+
+   static const int  MaxFreeArrayIndex         = 100;
+
+   static const int  outputParameterFile       = 0;
+   static const int  outputForceFile           = 1;
+
+   RealOpenMM fixedParameters[2]; // unused!
+
+   RealOpenMM totalEnergy                      = zero;
+
+   int indexOneDArraysToFree                   = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree                   = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree                = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   // set log file, if available
+
+   if( log ){
+      SimTKOpenMMLog::setSimTKOpenMMLog( log );
+   }
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " atoms=" << numberOfAtoms;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // load coordinates, velocities, and forces
+
+   RealOpenMM** atomCoordinates              = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = atomCoordinates;
+
+   RealOpenMM** velocities                   = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = velocities;
+
+   RealOpenMM** forces                       = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = forces;
+
+   RealOpenMM* energies                      = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfAtoms, NULL, true, (RealOpenMM) 0.0 );
+   oneDArraysToFree[indexOneDArraysToFree++] = energies;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      atomCoordinates[ii][0] = gromacAtomCoordinates[ii][0];
+      atomCoordinates[ii][1] = gromacAtomCoordinates[ii][1];
+      atomCoordinates[ii][2] = gromacAtomCoordinates[ii][2];
+
+      velocities[ii][0]      = gromacsVelocities[ii][0];
+      velocities[ii][1]      = gromacsVelocities[ii][1];
+      velocities[ii][2]      = gromacsVelocities[ii][2];
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // force setup
+
+   enum ForceTypes { HarmonicBond, AngleBond, ProperDihedral, RbDihedral, LJ14, LJCoulomb, MaxForceTypes };
+ 
+   int** bondAtomIndices[MaxForceTypes];
+   RealOpenMM** bondParameters[MaxForceTypes];
+   int numberOfBonds[MaxForceTypes];
+   ReferenceBondIxn* referenceBondIxn[MaxForceTypes];
+
+   ReferenceBondForce referenceBondForce;
+
+   for( int ii = 0; ii < MaxForceTypes; ii++ ){
+      bondAtomIndices[ii]     = NULL;
+      bondParameters[ii]      = NULL;
+      numberOfBonds[ii]       = 0;
+      referenceBondIxn[ii]    = NULL;
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // harmonic bonds
+
+   numberOfBonds[HarmonicBond]                       = gromacsGetHarmonicBondParameters( top, &bondAtomIndices[HarmonicBond], &bondParameters[HarmonicBond] );
+   referenceBondIxn[HarmonicBond]                    = new ReferenceHarmonicBondIxn();
+
+   // angle bonds
+
+   numberOfBonds[AngleBond]                          = gromacsGetAngleBondParameters( top, &bondAtomIndices[AngleBond], &bondParameters[AngleBond] );
+   referenceBondIxn[AngleBond]                       = new ReferenceAngleBondIxn();
+   
+   // proper dihedral
+
+   numberOfBonds[ProperDihedral]                     = gromacsGetProperDihedralBondParameters( top, &bondAtomIndices[ProperDihedral], &bondParameters[ProperDihedral] );
+   referenceBondIxn[ProperDihedral]                  = new ReferenceProperDihedralBond();
+
+   // rb dihedral
+
+   numberOfBonds[RbDihedral]                         = gromacsGetRbDihedralBondParameters( top, &bondAtomIndices[RbDihedral], &bondParameters[RbDihedral] );
+   referenceBondIxn[RbDihedral]                      = new ReferenceRbDihedralBond();
+
+   // 1-4 ixns
+
+   numberOfBonds[LJ14]                               = gromacsGetLJ14Parameters( top, &bondAtomIndices[LJ14], &bondParameters[LJ14] );
+   ReferenceLJ14* referenceLJ14                      = new ReferenceLJ14();
+   referenceBondIxn[LJ14]                            = referenceLJ14;
+
+   // calculate derived parameters for 1-4 LJ
+
+   RealOpenMM* charges = md->chargeA;
+   for( int ii  = 0; ii < numberOfBonds[LJ14]; ii++ ){
+      referenceLJ14->getDerivedParameters( bondParameters[LJ14][ii][0],
+                                           bondParameters[LJ14][ii][1],
+                                           charges[bondAtomIndices[LJ14][ii][0]],
+                                           charges[bondAtomIndices[LJ14][ii][1]], fr->epsfac*fr->fudgeQQ, bondParameters[LJ14][ii] );
+   }
+
+   // LJ Coulomb
+
+   gromacsGetLJCoulombParameters( top, fr, md, &bondAtomIndices[LJCoulomb], &bondParameters[LJCoulomb] );
+   ReferenceLJCoulombIxn referenceLJCoulomb;
+
+   // calculate derived parameters
+
+   RealOpenMM epsFacSqrt = SQRT( fr->epsfac );
+   for( int ii  = 0; ii < numberOfAtoms; ii++ ){
+      referenceLJCoulomb.getDerivedParameters( bondParameters[LJCoulomb][ii][0], bondParameters[LJCoulomb][ii][1],
+                                               bondParameters[LJCoulomb][ii][2],
+                                               epsFacSqrt, bondParameters[LJCoulomb][ii] );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // record arrays to be freed
+
+   for( int ii = 0; ii < MaxForceTypes; ii++ ){
+      if( bondAtomIndices[ii] != NULL ){
+         twoDIntArraysToFree[indexTwoDIntArraysToFree++] = bondAtomIndices[ii];
+      }
+      if( bondParameters[ii] != NULL ){
+         twoDArraysToFree[indexTwoDArraysToFree++]       = bondParameters[ii];
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // update() setup
+
+   // create referenceStochasticDynamics object
+   // get Shake parameters
+
+   ReferenceStochasticDynamics referenceStochasticDynamics( numberOfAtoms, deltaT, tau, temperature ); 
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName << " Parameters:\n   ";
+      referenceStochasticDynamics.printParameters( message );
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   int** shakeAtomIndices;
+   RealOpenMM** shakeParameters;
+   int numberOfConstraints                           = gromacsGetShakeParameters( top, md, &shakeAtomIndices, &shakeParameters );
+   twoDIntArraysToFree[indexTwoDIntArraysToFree++]   = shakeAtomIndices;
+   twoDArraysToFree[indexTwoDArraysToFree++]         = shakeParameters;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName << " Shake Parameters: " << numberOfConstraints << "\n";
+      for( int ii = 0; ii < numberOfConstraints && ii < 10; ii++ ){
+         message << " " << (ii+1) << " " << shakeParameters[ii][0] << " [" << shakeAtomIndices[ii][0] << " " << shakeAtomIndices[ii][1] << "]\n";
+      }
+      message << std::endl;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // add Shake to dynamics
+
+   ReferenceShakeAlgorithm referenceShakeAlgorithm( numberOfConstraints, shakeAtomIndices, shakeParameters );
+   referenceStochasticDynamics.setReferenceShakeAlgorithm( &referenceShakeAlgorithm );
+
+	referenceStochasticDynamics.writeState( numberOfAtoms, atomCoordinates,
+                                           velocities, forces, masses,
+                                           zeroI, baseFileName );
+
+   // ---------------------------------------------------------------------------------------
+
+   // time steps
+
+   for( int step = 0; step < timesteps; step++ ){
+
+      // zero forces
+
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         forces[ii][0]   = zero;
+         forces[ii][1]   = zero;
+         forces[ii][2]   = zero;
+         energies[ii]    = zero;
+      }
+
+      // calculate bond forces
+
+      for( int ii = 0; ii < MaxForceTypes; ii++ ){
+         if( numberOfBonds[ii] > 0 ){
+            referenceBondForce.calculateForce( numberOfBonds[ii], bondAtomIndices[ii], atomCoordinates,
+                                               bondParameters[ii], forces,
+                                               NULL, energies, &totalEnergy, *(referenceBondIxn[ii]) );
+         }
+      }
+
+      if( debug ){
+         std::stringstream message;
+         int maxAtom = 5;
+         message << methodName << " Bonded forces:\n";
+         message << " Bond counts: [";
+         for( int ii = 0; ii < MaxForceTypes; ii++ ){
+            message <<  numberOfBonds[ii] << " ";
+         }
+         message << "]";
+         message << std::endl;
+
+         for( int ii = 0; ii < maxAtom; ii++ ){
+            message << ii;
+            message << " x[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
+            message << "] f[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, forces[ii], 3, one );
+            message << "]";
+            message << std::endl;
+         }
+         SimTKOpenMMLog::printMessage( message );
+      }
+   
+      // LJ Coulomb forces
+      // bondAtomIndices contains exclusions
+
+      referenceLJCoulomb.calculatePairIxn( numberOfAtoms, atomCoordinates,
+                                           bondParameters[LJCoulomb], bondAtomIndices[LJCoulomb],
+                                           fixedParameters, forces, NULL, &totalEnergy );
+
+      if( debug ){
+         std::stringstream message;
+         int maxAtom = 5;
+         message << methodName << " Bonded + LJ/Coulomb forces:\n   ";
+         for( int ii = 0; ii < maxAtom; ii++ ){
+            message << ii;
+            message << " x[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
+            message << "] f[";
+            SimTKOpenMMUtilities::formatRealStringStream( message, forces[ii], 3, one );
+            message << "]";
+            message << std::endl;
+         }
+         SimTKOpenMMLog::printMessage( message );
+      }
+   
+      // do update
+
+      referenceStochasticDynamics.update( numberOfAtoms, atomCoordinates, velocities, forces, masses );
+
+      if( debug ){
+         std::stringstream message;
+         message << methodName << " Called update: " << referenceStochasticDynamics.getTimeStep() << "\n";
+         SimTKOpenMMLog::printMessage( message );
+      }
+   }
+
+	referenceStochasticDynamics.writeState( numberOfAtoms, atomCoordinates,
+                                           velocities, forces, masses,
+                                           zeroI, baseFileName );
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   for( int ii = 0; ii < MaxForceTypes; ii++ ){
+      if( referenceBondIxn[ii] ){
+         delete referenceBondIxn[ii];
+      }
+   }
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Remove linear momentum from velocities
+
+   @param numberOfAtoms            number of atoms
+   @param gromacsMasses            masses
+   @param gromacsVelocities        velocities
+   @param baseFileName             base file name
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsReferenceRemoveLinearMomentum( int numberOfAtoms, const real* gromacsMasses,
+                                          rvec* gromacsVelocities,
+                                          char* baseFileName, FILE* log ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName               = "\ngromacsReferenceRemoveLinearMomentum: ";
+   static int debug                            = false;
+
+   static const RealOpenMM zero                = 0.0;
+   static const RealOpenMM one                 = 1.0;
+
+   static const int  zeroI                     = 0;
+   static const int  twoI                      = 2;
+   static const int  threeI                    = 3;
+   static const int  fourI                     = 4;
+
+   static const int  MaxFreeArrayIndex         = 100;
+
+   int indexOneDArraysToFree                   = 0;
+   RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDArraysToFree                   = 0;
+   RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   int indexTwoDIntArraysToFree                = 0;
+   int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   // set log file, if available
+
+   if( log ){
+      SimTKOpenMMLog::setSimTKOpenMMLog( log );
+   }
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " atoms=" << numberOfAtoms;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // load masses and velocities
+
+   RealOpenMM*  masses                       = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfAtoms, NULL, true, (RealOpenMM) 0.0 );
+   oneDArraysToFree[indexOneDArraysToFree++] = masses;
+
+   RealOpenMM** velocities                   = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+   twoDArraysToFree[indexTwoDArraysToFree++] = velocities;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      masses[ii]             = gromacsMasses[ii];
+
+      velocities[ii][0]      = gromacsVelocities[ii][0];
+      velocities[ii][1]      = gromacsVelocities[ii][1];
+      velocities[ii][2]      = gromacsVelocities[ii][2];
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // remove linear momentum, writing out velocities before and after
+
+   ReferenceDynamics referenceDynamics( numberOfAtoms, zero, zero );
+	referenceDynamics.writeState( numberOfAtoms, NULL, velocities, NULL, masses, zeroI, baseFileName );
+   referenceDynamics.removeTotalLinearMomentum( numberOfAtoms, masses, velocities );
+   referenceDynamics.incrementTimeStep();
+	referenceDynamics.writeState( numberOfAtoms, NULL, velocities, NULL, masses, zeroI, baseFileName );
+
+   // ---------------------------------------------------------------------------------------
+
+   // free memory
+
+   gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                            MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                            methodName );
+
+   // ---------------------------------------------------------------------------------------
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate forces for given configuration and topology
+   Used in self-test; parameter arrays, ... only allocated once (static); freed
+   if debug is -1
+
+   @param top                      Gromacs t_topology data struct
+   @param gromacAtomCoordinates    atom configuration
+   @param md                       Gromacs t_mdatoms data struct
+   @param fr                       Gromacs t_forcerec data struct
+   @param forces                   output forces
+   @param includeNonBonded         include nonbonded ixn
+   @param baseFileName             Base file name
+   @param debug                    debug flag (== -1, free arrays and return)
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+extern "C"
+int gromacsReferenceForceSelfTest( 
+                           const t_topology* top, const rvec* gromacAtomCoordinates,
+                           const t_mdatoms *md, t_forcerec* fr, rvec* forces,
+                           int includeNonBonded, char* baseFileName, int debug, FILE* log ){ 
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName                     = "\ngromacsReferenceForceSelfTest: ";
+
+   // ---------------------------------------------------------------------------------------
+
+   static const RealOpenMM zero                      = 0.0;
+   static const RealOpenMM one                       = 1.0;
+
+   static const int  zeroI                           = 0;
+   static const int  twoI                            = 2;
+   static const int  threeI                          = 3;
+   static const int  fourI                           = 4;
+
+   // ---------------------------------------------------------------------------------------
+
+   static const int  includeHarmonic                 = 1;
+   static const int  includeAngle                    = 1;
+   static const int  includeProperDihedral           = 1;
+   static const int  includeRbDihedral               = 1;
+   static const int  includeLJ14                     = 1;
+   static const int  includeLJCoulomb                = 1;
+
+   // ---------------------------------------------------------------------------------------
+
+   static const int  HarmonicIndex                   = 0;
+   static const int  AngleIndex                      = 1;
+   static const int  ProperDihedralIndex             = 2;
+   static const int  RbDihedralIndex                 = 3;
+   static const int  LJ14Index                       = 4;
+   static const int  LJCoulombIndex                  = 5;
+   static const int  MaxForceIndex                   = 6;
+   static int bondCount[MaxForceIndex];
+
+   static const int  outputParameterFile             = 0;
+   static int        outputForceFile                 = 1;
+
+   RealOpenMM totalEnergy                            = zero;
+
+   // ---------------------------------------------------------------------------------------
+
+   static const int  MaxFreeArrayIndex               = 20;
+   static int indexOneDArraysToFree                  = 0;
+   static RealOpenMM* oneDArraysToFree[MaxFreeArrayIndex];
+
+   static int indexTwoDArraysToFree                  = 0;
+   static RealOpenMM** twoDArraysToFree[MaxFreeArrayIndex];
+
+   static int indexTwoDIntArraysToFree               = 0;
+   static int** twoDIntArraysToFree[MaxFreeArrayIndex];
+
+   // ---------------------------------------------------------------------------------------
+
+   static RealOpenMM** atomCoordinates               = NULL;
+   static RealOpenMM** atomForces                    = NULL;
+   static RealOpenMM** tempForces                    = NULL;
+   static RealOpenMM*  tempEnergies                  = NULL;
+
+   static int** harmonicBondAtomIndices              = NULL;
+   static RealOpenMM** harmonicBondParameters        = NULL;
+
+   static int** angleBondAtomIndices                 = NULL;
+   static RealOpenMM** angleBondParameters           = NULL;
+
+   static int** properDihedralBondAtomIndices        = NULL;
+   static RealOpenMM** properDihedralBondParameters  = NULL;
+
+   static int** rbDihedralBondAtomIndices            = NULL;
+   static RealOpenMM** rbDihedralBondParameters      = NULL;
+
+   static int** lj14AtomIndices                      = NULL;
+   static RealOpenMM** lj14Parameters                = NULL;
+
+   static int** ljCoulombExclusions                  = NULL;
+   static RealOpenMM** ljCoulombParameters           = NULL;
+
+   // ---------------------------------------------------------------------------------------
+
+   // set log file, if available
+
+   if( log ){
+      SimTKOpenMMLog::setSimTKOpenMMLog( log );
+   }
+
+   if( !baseFileName ){
+      outputForceFile = 0;
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // free arrays and return
+
+   if( debug == -1 ){
+
+      gromacsFreeSundryArrays( indexOneDArraysToFree, indexTwoDArraysToFree, indexTwoDIntArraysToFree,
+                               MaxFreeArrayIndex, oneDArraysToFree, twoDArraysToFree, twoDIntArraysToFree,
+                               methodName );
+
+      atomCoordinates               = NULL;
+      atomForces                    = NULL;
+      tempForces                    = NULL;
+      tempEnergies                  = NULL;
+   
+      harmonicBondAtomIndices       = NULL;
+      harmonicBondParameters        = NULL;
+   
+      angleBondAtomIndices          = NULL;
+      angleBondParameters           = NULL;
+   
+      properDihedralBondAtomIndices = NULL;
+      properDihedralBondParameters  = NULL;
+   
+      rbDihedralBondAtomIndices     = NULL;
+      rbDihedralBondParameters      = NULL;
+   
+      lj14AtomIndices               = NULL;
+      lj14Parameters                = NULL;
+   
+      ljCoulombExclusions           = NULL;
+      ljCoulombParameters           = NULL;
+   
+      return 0;
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms = top->atoms.nr;
+
+   if( debug ){
+      std::stringstream message;
+      message << methodName;
+      message << " atoms=" << numberOfAtoms;
+      SimTKOpenMMLog::printMessage( message );
+   }
+
+   // setup
+
+   if( atomCoordinates == NULL ){
+
+      atomCoordinates                                     = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]           = atomCoordinates;
+
+      atomForces                                          = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]           = atomForces;
+
+
+      tempForces                                          = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfAtoms, 3, NULL, true, (RealOpenMM) 0.0 );
+      twoDArraysToFree[indexTwoDArraysToFree++]           = tempForces;
+   
+      tempEnergies                                        = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
+      oneDArraysToFree[indexOneDArraysToFree++]           = tempEnergies;
+   
+      // harmonic 
+
+      bondCount[HarmonicIndex]                            = gromacsGetHarmonicBondParameters( top, &harmonicBondAtomIndices, &harmonicBondParameters );
+      twoDIntArraysToFree[indexTwoDIntArraysToFree++]     = harmonicBondAtomIndices;
+      twoDArraysToFree[indexTwoDArraysToFree++]           = harmonicBondParameters;
+
+      // angle
+
+      bondCount[AngleIndex]                               = gromacsGetAngleBondParameters( top, &angleBondAtomIndices, &angleBondParameters );
+      twoDIntArraysToFree[indexTwoDIntArraysToFree++]     = angleBondAtomIndices;
+      twoDArraysToFree[indexTwoDArraysToFree++]           = angleBondParameters;
+
+      // proper dihedral
+
+      bondCount[ProperDihedralIndex]                      = gromacsGetProperDihedralBondParameters( top, &properDihedralBondAtomIndices, &properDihedralBondParameters );
+      twoDIntArraysToFree[indexTwoDIntArraysToFree++]     = properDihedralBondAtomIndices;
+      twoDArraysToFree[indexTwoDArraysToFree++]           = properDihedralBondParameters;
+
+      // rb dihedral
+
+      bondCount[RbDihedralIndex]                          = gromacsGetRbDihedralBondParameters( top, &rbDihedralBondAtomIndices, &rbDihedralBondParameters );
+      twoDIntArraysToFree[indexTwoDIntArraysToFree++]     = rbDihedralBondAtomIndices;
+      twoDArraysToFree[indexTwoDArraysToFree++]           = rbDihedralBondParameters;
+   
+      // LJ14
+
+      bondCount[LJ14Index]                                = gromacsGetLJ14Parameters( top, &lj14AtomIndices, &lj14Parameters );
+      twoDIntArraysToFree[indexTwoDIntArraysToFree++]     = lj14AtomIndices;
+      twoDArraysToFree[indexTwoDArraysToFree++]           = lj14Parameters;
+
+      // calculate derived parameters
+
+      ReferenceLJ14 referenceLJ14;
+      RealOpenMM* charges = md->chargeA;
+      for( int ii  = 0; ii < bondCount[LJ14Index]; ii++ ){
+         referenceLJ14.getDerivedParameters( lj14Parameters[ii][0], lj14Parameters[ii][1],
+                                             charges[lj14AtomIndices[ii][0]],
+                                             charges[lj14AtomIndices[ii][1]], fr->epsfac*fr->fudgeQQ, lj14Parameters[ii] );
+      }
+
+      // LJ Coulomb
+
+      if( includeNonBonded ){
+         bondCount[LJCoulombIndex]                           = gromacsGetLJCoulombParameters( top, fr, md, &ljCoulombExclusions, &ljCoulombParameters );
+         twoDIntArraysToFree[indexTwoDIntArraysToFree++]     = ljCoulombExclusions;
+         twoDArraysToFree[indexTwoDArraysToFree++]           = ljCoulombParameters;
+      
+         ReferenceLJCoulombIxn referenceLJCoulomb;
+      
+         // calculate derived parameters
+      
+         RealOpenMM epsFacSqrt = SQRT( fr->epsfac );
+         for( int ii  = 0; ii < numberOfAtoms; ii++ ){
+            referenceLJCoulomb.getDerivedParameters( ljCoulombParameters[ii][0], ljCoulombParameters[ii][1],
+                                                     ljCoulombParameters[ii][2],
+                                                     epsFacSqrt, ljCoulombParameters[ii] );
+         }
+      }
+   
+   }
+   memset( atomForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+   memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      atomCoordinates[ii][0] = gromacAtomCoordinates[ii][0];
+      atomCoordinates[ii][1] = gromacAtomCoordinates[ii][1];
+      atomCoordinates[ii][2] = gromacAtomCoordinates[ii][2];
+   }
+
+   // accumulate bonded force contributions
+
+   // ---------------------------------------------------------------------------------------
+
+   ReferenceBondForce referenceBondForce;
+
+   // harmonic bond forces
+
+   if( includeHarmonic ){
+
+      ReferenceHarmonicBondIxn referenceHarmonicBondIxn;
+
+      // harmonic bond forces
+
+      memset( tempForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+      memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+      referenceBondForce.calculateForce( bondCount[HarmonicIndex], harmonicBondAtomIndices, atomCoordinates,
+                                         harmonicBondParameters, tempForces,
+                                         NULL, tempEnergies, 
+                                         &totalEnergy, referenceHarmonicBondIxn );
+
+      if( atomForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, tempForces, atomForces );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream harmonicBondForceFileName;
+         harmonicBondForceFileName << baseFileName;
+         harmonicBondForceFileName << "SelfHarmonicBond.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, tempForces,
+                                      tempEnergies, harmonicBondForceFileName.str(), one, one, one );
+      }
+
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // angle bond forces
+
+   if( includeAngle ){
+
+      ReferenceAngleBondIxn referenceAngleBondIxn;
+
+      // angle bond forces
+
+      memset( tempForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+      memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+      referenceBondForce.calculateForce( bondCount[AngleIndex], angleBondAtomIndices, atomCoordinates,
+                                         angleBondParameters, tempForces, 
+                                         NULL, tempEnergies, 
+                                         &totalEnergy, referenceAngleBondIxn );
+
+      if( atomForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, tempForces, atomForces );
+      }
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream angleForceFileName;
+         angleForceFileName << baseFileName;
+         angleForceFileName << "SelfAngleBond.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, tempForces,
+                                      tempEnergies, angleForceFileName.str(), one, one, one );
+      }
+
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // properDihedral bond forces
+
+   if( includeProperDihedral ){
+
+      ReferenceProperDihedralBond referenceProperDihedralBond;
+
+      memset( tempForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+      memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+      referenceBondForce.calculateForce( bondCount[ProperDihedralIndex], properDihedralBondAtomIndices, atomCoordinates,
+                                         properDihedralBondParameters, tempForces, 
+                                         NULL, tempEnergies, 
+                                         &totalEnergy, referenceProperDihedralBond );
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream properDihedralForceFileName;
+         properDihedralForceFileName << baseFileName;
+         properDihedralForceFileName << "SelfProperDihedral.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, tempForces,
+                                      tempEnergies, properDihedralForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( atomForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, tempForces, atomForces );
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // rbDihedral bond forces
+
+   if( includeRbDihedral ){
+
+      ReferenceRbDihedralBond referenceRbDihedralBond;
+
+      memset( tempForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+      memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+      referenceBondForce.calculateForce( bondCount[RbDihedralIndex], rbDihedralBondAtomIndices, atomCoordinates,
+                                         rbDihedralBondParameters, tempForces, 
+                                         NULL, tempEnergies, 
+                                         &totalEnergy, referenceRbDihedralBond );
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream rbDihedralForceFileName;
+         rbDihedralForceFileName << baseFileName;
+         rbDihedralForceFileName << "SelfRbDihedral.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, tempForces,
+                                      tempEnergies, rbDihedralForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( atomForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, tempForces, atomForces );
+      }
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // lj14 bond forces
+
+   if( includeLJ14 ){
+
+      ReferenceLJ14 referenceLJ14;
+
+      memset( tempForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+      memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+      referenceBondForce.calculateForce( bondCount[LJ14Index], lj14AtomIndices, atomCoordinates,
+                                         lj14Parameters, tempForces, 
+                                         NULL, tempEnergies, 
+                                         &totalEnergy, referenceLJ14 );
+
+      // output file
+   
+      if( outputForceFile ){
+         std::stringstream lj14ForceFileName;
+         lj14ForceFileName << baseFileName;
+         lj14ForceFileName << "SelfLJ14.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, tempForces,
+                                      tempEnergies, lj14ForceFileName.str(), one, one, one );
+      }
+
+      // accumulate bonded forces
+
+      if( atomForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, tempForces, atomForces );
+      }
+   }
+
+   // output file
+
+   if( outputForceFile && atomForces ){
+      std::stringstream bondedForceFileName;
+      bondedForceFileName << baseFileName;
+      bondedForceFileName << "SelfBonded.txt";
+      ReferenceForce::writeForces( numberOfAtoms, zeroI, atomCoordinates, atomForces,
+                                   NULL, bondedForceFileName.str(), one, one, one );
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // ljCoulomb forces
+
+   if( includeNonBonded && includeLJCoulomb ){
+
+      ReferenceLJCoulombIxn referenceLJCoulomb;
+
+      // currently not used!
+
+      RealOpenMM fixedParameters[MaxFreeArrayIndex];
+
+      memset( tempForces[0], 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
+      memset( tempEnergies, 0, sizeof( RealOpenMM )*numberOfAtoms );
+
+      // ljCoulomb bond forces
+   
+      referenceLJCoulomb.calculatePairIxn( numberOfAtoms, atomCoordinates,
+                                           ljCoulombParameters, ljCoulombExclusions,
+                                           fixedParameters, tempForces, 
+                                           tempEnergies, &totalEnergy );
+   
+      if( outputForceFile ){
+         std::stringstream lJCoulombForceFileName;
+         lJCoulombForceFileName << baseFileName;
+         lJCoulombForceFileName << "SelfLJCoulomb.txt";
+         ReferenceForce::writeForces( numberOfAtoms, fourI, atomCoordinates, tempForces,
+                                      tempEnergies, lJCoulombForceFileName.str(), one, one, one );
+      }
+
+      // accumulate forces
+
+      if( atomForces ){
+         SimTKOpenMMUtilities::addTwoDimArray( numberOfAtoms, 3, tempForces, atomForces );
+      }
+   }
+
+   // output file
+
+   if( outputForceFile && atomForces ){
+      std::stringstream forceFileName;
+      forceFileName << baseFileName;
+      forceFileName << "SelfTotalF.txt";
+      ReferenceForce::writeForces( numberOfAtoms, zeroI, atomCoordinates, atomForces,
+                                   NULL, forceFileName.str(), one, one, one );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // accumulate forces
+
+   if( forces ){
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         forces[ii][0] = atomForces[ii][0];
+         forces[ii][1] = atomForces[ii][1];
+         forces[ii][2] = atomForces[ii][2];
+      }
+   }
+
+   return 0;
+}

platforms/reference/src/gromacsReferenceInterface.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __GromacsReferenceInterface_H__
+#define __GromacsReferenceInterface_H__
+
+#ifdef __cplusplus
+#define externC extern "C"
+#else
+#define externC extern
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate forces for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param gromacAtomCoordinates    atom configuration
+   @param partialChargesIn         array of partial charges
+   @param fr                       Gromac's force_record?????
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsReferenceForce( const t_topology* top, const rvec* gromacAtomCoordinates,
+                           const t_mdatoms *md, t_forcerec* fr,
+                           char* baseFileName, FILE* log );
+
+externC
+int gromacsReferenceSdUpdate( int numberOfAtoms, const t_topology* top, const rvec* gromacAtomCoordinates,
+                              const t_mdatoms *md, t_forcerec* fr,
+                              const rvec* gromacsVelocities, const rvec* gromacsForces,
+                              real* masses, real deltaT, real tau, real temperature,
+                              char* baseFileName, FILE* log );
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate forces for given configuration and topology
+   Used in self-test; parameter arrays, ... only allocated once (static); freed
+   if debug is -1
+
+   @param top                      Gromacs t_topology data struct
+   @param gromacAtomCoordinates    atom configuration
+   @param md                       Gromacs t_mdatoms data struct
+   @param fr                       Gromacs t_forcerec data struct
+   @param includeNonBonded         include nonbonded ixn
+   @param baseFileName             Base file name
+   @param debug                    debug flag (== -1, free arrays and return)
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsReferenceForceSelfTest( 
+                                   const t_topology* top, const rvec* gromacAtomCoordinates,
+                                   const t_mdatoms *md, t_forcerec* fr,
+                                   rvec* forces, int includeNonBonded, char* baseFileName,
+                                   int debug, FILE* log );
+
+/**---------------------------------------------------------------------------------------
+
+   Remove linear momentum from velocities
+
+   @param numberOfAtoms            number of atoms
+   @param gromacsMasses            masses
+   @param gromacsVelocities        velocities
+   @param baseFileName             base file name
+   @param log                      log reference (stdlog in md.c)
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsReferenceRemoveLinearMomentum( int numberOfAtoms, const real* gromacsMasses,
+                                          rvec* gromacsVelocities, char* baseFileName, FILE* log );
+
+#endif

platforms/reference/src/gromacsReferenceInterfaceCpp.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __GromacsReferenceInterfaceCpp_H__
+#define __GromacsReferenceInterfaceCpp_H__
+
+#ifdef __cplusplus
+#define externC extern "C"
+#else
+#define externC extern
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Allocate memory for atom indices participating in bonds and the bond parameters
+
+   @param numberOfBonds            number of bonds
+   @param numberOfAtomIndices      (number of atoms)/bond
+   @param outputBondIndices        allocated memory for atom indices outputBondIndices[bondIndex][atomIndex]
+   @param numberParameters         number of parameters/bond
+   @param outputBondParameters     allocated memory for parameters outputBondParameters[bondIndex][parameterIndex]
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsAllocateBondIxnArrays( int numberOfBonds, int numberOfAtomIndices,
+                                  int*** outputBondIndices, int numberParameters,
+                                   RealOpenMM*** outputBondParameters );
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for harmonic bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsGetHarmonicBondParameters( const t_topology* top, int*** atomIndices,
+                                      RealOpenMM*** bondParameters );
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate harmonic bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param harmonicBondForces       output array of forces: harmonicBondForces[atomIndex][3]
+   @param harmonicBondEnergies     output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsHarmonicBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                       RealOpenMM** harmonicBondForces,
+                                       RealOpenMM* harmonicBondAtomEnergies,
+                                       RealOpenMM* totalEnergy );
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for angle bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsGetAngleBondParameters( const t_topology* top, int*** atomIndices,
+                                   RealOpenMM*** bondParameters );
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate angle bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param angleBondForces       output array of forces: angleBondForces[atomIndex][3]
+   @param angleBondEnergies     output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsAngleBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                    RealOpenMM** angleBondForces,
+                                    RealOpenMM* angleBondAtomEnergies,
+                                    RealOpenMM* totalEnergy ); 
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for properDihedral bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsGetProperDihedralBondParameters( const t_topology* top, int*** atomIndices,
+                                            RealOpenMM*** bondParameters );
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate properDihedral bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param properDihedralBondForces       output array of forces: angleBondForces[atomIndex][3]
+   @param properDihedralBondEnergies     output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsProperDihedralBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                             RealOpenMM** properDihedralBondForces,
+                                             RealOpenMM* properDihedralBondAtomEnergies,
+                                             RealOpenMM* totalEnergy );
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for RB Dihedral bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+		free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsGetRbDihedralBondParameters( const t_topology* top, int*** atomIndices,
+                                            RealOpenMM*** bondParameters );
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate RB dihedral bond forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param rbDihedralBondForces     output array of forces: angleBondForces[atomIndex][3]
+   @param rbDihedralBondEnergies   output array of bond energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsRbDihedralBondReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                                         RealOpenMM** rbDihedralBondForces,
+                                         RealOpenMM* rbDihedralBondAtomEnergies,
+                                         RealOpenMM* totalEnergy );
+
+/**---------------------------------------------------------------------------------------
+
+   Get bond atom indices and parameters for LJ 14 bonds
+
+   @param top                      Gromacs t_topology data struct
+   @param atomIndices              array of atomIndices[bondIndex][atomIndex]
+   @param bondParameters           array of bond parameters[bondIndex][parameterIndex]
+
+   The arrays atomIndices & bondParameters are allocated off the heap;
+   The memory for the second index is one block of size (number of bonds)*(number of atoms in bond)
+
+   The memory should be freed by the callee via calls:
+      free( atomIndices[0] ); free( atomIndices );
+
+   @return number of bonds
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsGetLJ14Parameters( const t_topology* top, int*** atomIndices,
+                              RealOpenMM*** bondParameters );
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate LJ 14 forces/energies for given configuration and topology
+
+   @param top                      Gromacs t_topology data struct
+   @param atomCoordinates          atom configuration
+   @param fr                       Gromacs t_forcerec struct
+   @param lj14Forces               output array of forces: lj14Forces[atomIndex][3]
+   @param lj14AtomEnergies         output array of LJ 14 energies[atomIndex]
+   @param totalEnergy              output total energy
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+externC
+int gromacsLJ14ReferenceForce( const t_topology* top, RealOpenMM** atomCoordinates,
+                               t_forcerec *fr, const RealOpenMM* charges, RealOpenMM** lj14Forces,
+                               RealOpenMM* lj14AtomEnergies,
+                               RealOpenMM* totalEnergy );
+
+#endif

platforms/reference/src/gromacsReferenceParameterPrint.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include "SimTKOpenMMGromacsUtilities.h"
+#include "gromacsReferenceParameterPrint.h"
+
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+
+/**---------------------------------------------------------------------------------------
+      
+   Print bond parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printBondParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* bondFile ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintBondParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   t_iatom *atoms;
+   int  ii, numberOfBonds, type, i, j, bondIndex;
+  
+   t_idef *idef   = &top->idef;
+   atoms          = idef->il[ F_BONDS ].iatoms;
+   numberOfBonds  = idef->il[ F_BONDS ].nr/3; 
+   
+  (void) fprintf( bondFile, "%d Bond parameters: [atom indices] idealBond k\n", numberOfBonds );
+
+   bondIndex      = 1;
+   for( ii = 0; ii < idef->il[ F_BONDS ].nr; ii += 3 ) {
+
+      type = atoms[ ii     ]; 
+      i    = atoms[ ii + 1 ]; 
+      j    = atoms[ ii + 2 ]; 
+
+      (void) fprintf( bondFile, "%6d %6d %6d %13.6e %13.6e\n",
+                                 bondIndex++, i, j, idef->iparams[ type ].harmonic.rA, idef->iparams[ type ].harmonic.krA );
+
+   }
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print angle parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printAngleParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* angleFile ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintAngleParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   t_iatom *atoms;
+   int  ii, numberOfAngles, type, i, j, k, angleIndex;
+  
+   t_idef *idef   = &top->idef;
+   atoms          = idef->il[ F_ANGLES ].iatoms;
+   numberOfAngles = idef->il[ F_ANGLES ].nr/4; 
+   
+  (void) fprintf( angleFile, "%d Angle parameters: [atom indices] idealAngle k\n", numberOfAngles );
+
+   angleIndex     = 1;
+   for( ii = 0; ii < idef->il[ F_ANGLES ].nr; ii += 4 ) {
+
+      type = atoms[ ii     ]; 
+      i    = atoms[ ii + 1 ]; 
+      j    = atoms[ ii + 2 ]; 
+      k    = atoms[ ii + 3 ]; 
+
+      (void) fprintf( angleFile, "%6d %6d %6d %6d %13.6e %13.6e\n",
+                                 angleIndex++, i, j, k, idef->iparams[ type ].harmonic.rA, idef->iparams[ type ].harmonic.krA );
+
+   }
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print rbDihedral parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printRbDihedralParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* rbDihedralFile ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintRbDihedralParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   t_iatom *atoms;
+   int  ii, numberOfRbDihedrals, type, i, j, k, l, rbDihedralIndex;
+  
+   t_idef *idef        = &top->idef;
+   atoms               = idef->il[ F_RBDIHS ].iatoms;
+   numberOfRbDihedrals = idef->il[ F_RBDIHS ].nr/5; 
+   
+  (void) fprintf( rbDihedralFile, "%d RbDihedral parameters: [atom indices] RbDihedral parameters\n", numberOfRbDihedrals );
+
+   rbDihedralIndex     = 1;
+   for( ii = 0; ii < idef->il[ F_RBDIHS ].nr; ii += 5 ) {
+
+      type = atoms[ ii     ]; 
+      i    = atoms[ ii + 1 ]; 
+      j    = atoms[ ii + 2 ]; 
+      k    = atoms[ ii + 3 ]; 
+      l    = atoms[ ii + 4 ]; 
+
+      (void) fprintf( rbDihedralFile, "%6d %6d %6d %6d %6d %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e\n",
+                      rbDihedralIndex++, i, j, k, l, 
+                      idef->iparams[ type ].rbdihs.rbc[0],
+                      idef->iparams[ type ].rbdihs.rbc[1],
+                      idef->iparams[ type ].rbdihs.rbc[2],
+                      idef->iparams[ type ].rbdihs.rbc[3],
+                      idef->iparams[ type ].rbdihs.rbc[4],
+                      idef->iparams[ type ].rbdihs.rbc[5] );
+
+   }
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print properDihedral parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printProperDihedralParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* properDihedralFile ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintProperDihedralParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   t_iatom *atoms;
+   int  ii, numberOfProperDihedrals, type, i, j, k, l, properDihedralIndex;
+  
+   t_idef *idef            = &top->idef;
+   atoms                   = idef->il[ F_PDIHS ].iatoms;
+   numberOfProperDihedrals = idef->il[ F_PDIHS ].nr/5; 
+   
+  (void) fprintf( properDihedralFile, "%d ProperDihedral parameters: [atom indices] ProperDihedral parameters\n", numberOfProperDihedrals );
+
+   properDihedralIndex     = 1;
+   for( ii = 0; ii < idef->il[ F_PDIHS ].nr; ii += 5 ) {
+
+      type = atoms[ ii     ]; 
+      i    = atoms[ ii + 1 ]; 
+      j    = atoms[ ii + 2 ]; 
+      k    = atoms[ ii + 3 ]; 
+      l    = atoms[ ii + 4 ]; 
+
+      (void) fprintf( properDihedralFile, "%6d %6d %6d %6d %6d %13.6e %13.6e %3d\n",
+                      properDihedralIndex++, i, j, k, l, 
+                      idef->iparams[ type ].pdihs.cpA,
+                      idef->iparams[ type ].pdihs.phiA,
+                      idef->iparams[ type ].pdihs.mult );
+
+   }
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print LJ 14 parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printLJ14Parameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* lj14File ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintLj14Parameters";
+
+   real c6, c12;
+   real oneSixth = ((real) 1.0)/ ( (real) 6.0);
+
+   real zero     = 0.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   t_iatom *atoms;
+   int  ii, numberOfLJ14, type, i, j, lj14Index;
+  
+   t_idef *idef            = &top->idef;
+   real* charges           = md->chargeA;
+
+   atoms                   = idef->il[ F_LJ14 ].iatoms;
+   numberOfLJ14            = idef->il[ F_LJ14 ].nr/3; 
+   
+  (void) fprintf( lj14File, "%d LJ 14 parameters: [atom indices] c6 c12 q1 q2 where epsfac=%.5e fudge=%.5e\n", numberOfLJ14,
+                  fr->epsfac, fr->fudgeQQ );
+
+   lj14Index     = 1;
+   for( ii = 0; ii < idef->il[ F_LJ14 ].nr; ii += 3 ) {
+
+      type = atoms[ ii     ]; 
+      i    = atoms[ ii + 1 ]; 
+      j    = atoms[ ii + 2 ]; 
+
+      c6   = idef->iparams[ type ].lj14.c6A;
+      c12  = idef->iparams[ type ].lj14.c12A;
+
+
+/*
+      if( c12 > zero && c6 > zero ){
+         sig  =  c6*c6/c12;
+         eps  =  (real) pow( (c12/c6), oneSixth );
+      } else {
+         sig  = zero;
+         eps  = zero;
+      }
+
+      coulomb = fr->epsfac*fr->fudgeQQ*charges[i]*charges[j];
+*/
+
+      (void) fprintf( lj14File, "%6d %6d %6d %14.6e %14.6e %9.5f %9.5f\n",
+                      lj14Index++, i, j, c6, c12, charges[i], charges[j] );
+
+   }
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print LJ Coulomb parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printLJCoulombParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* ljCoulombFile ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintLjCoulombParameters";
+
+   real c6, c12;
+
+   real zero      = 0.0;
+   real half      = ((real) 1.0)/ ( (real) 2.0);
+   real oneSixth  = ((real) 1.0)/ ( (real) 6.0);
+   real onetwelth = ((real) 1.0)/ ( (real) 12.0);
+
+   real mass;
+   real scaleFactor;
+
+   int errors    = 0;
+
+   // ---------------------------------------------------------------------------------------
+
+   int  ii, jj;
+  
+   t_idef *idef            = &top->idef;
+   real* charges           = md->chargeA;
+   t_block *gmx_excl       = &top->atoms.excl;
+   int ntypes              = fr->ntype;
+   int* types              = md->typeA;
+   real* nbfp              = fr->nbfp;
+   
+  (void) fprintf( ljCoulombFile, "%d epsfac=%.5e fudge=%.5e  LJ Coulomb parameters: { atom c6 c12 q numberOfExclusions Exclsions[] }\n", md->nr,
+                  fr->epsfac, fr->fudgeQQ );
+
+   // loop over atoms
+
+   for( ii = 0; ii < md->nr; ii++ ) {
+
+      c6              = nbfp[ types[ii]*2*ntypes + types[ii]*2 ];
+      c12             = nbfp[ types[ii]*2*ntypes + types[ii]*2 + 1 ];
+
+      char* atomTypeC = *(top->atoms.atomtype[ii]);
+
+/*
+      if( c12 > zero && c6 > zero ){
+         sig     =  half*c6*c6/c12;
+         eps     =  (real) pow( (c12/c6), onetwelth );
+      } else {
+         sig     = zero;
+         eps     = zero;
+      }
+*/
+
+      mass = md->massA[ii];
+      if ( mass < 1.2f && mass >= 1.0f ){        // hydrogen
+         scaleFactor  = 0.85f; 
+      } else if( mass > 11.8f && mass < 12.2f ){ // carbon
+         scaleFactor  = 0.72f; 
+      } else if( mass > 14.0f && mass < 15.0f ){ // nitrogen
+         scaleFactor  = 0.79f;
+      } else if( mass > 15.5f && mass < 16.5f ){ // oxygen
+         scaleFactor  = 0.85f; 
+      } else if( mass > 31.5f && mass < 32.5f ){ // sulphur
+         scaleFactor  = 0.96f;
+      } else if( mass > 29.5f && mass < 30.5f ){ // phosphorus
+         scaleFactor  = 0.86f;
+      } else {
+         errors++;
+         (void) fprintf( ljCoulombFile, "Warning: mass for atom %d mass=%.4f not recognized.",
+                         ii, mass );
+      }   
+
+      (void) fprintf( ljCoulombFile, "%6d %14.6e %14.6e %9.5f %12s %9.4f %d",
+                      ii, c6, c12, charges[ii], atomTypeC, scaleFactor, (gmx_excl->index[ii+1] - gmx_excl->index[ii]) );
+
+      for ( jj = gmx_excl->index[ii]; jj < gmx_excl->index[ii+1]; jj++ ) {
+         (void) fprintf( ljCoulombFile, " %d", gmx_excl->a[jj] );
+      }    
+      (void) fprintf( ljCoulombFile, "\n" );
+
+      // symmetry check
+
+      for ( jj = gmx_excl->index[ii]; jj < gmx_excl->index[ii+1]; jj++ ){
+         int checkI = gmx_excl->a[jj];
+         int hit    = 0;
+         for( int kk = gmx_excl->index[checkI]; kk < gmx_excl->index[checkI+1] && !hit; kk++ ){
+            if( gmx_excl->a[kk] == ii ){
+               hit = 1;
+            }
+         }
+         if( !hit ){
+            (void) fprintf( ljCoulombFile, " Missing %d in list for %d\n", ii, checkI );
+         }
+      }    
+
+   }
+
+   if( errors ){
+      (void) fprintf( stdout, "printLJ14Parameters errors -- aborting!\n" );
+      (void) fprintf( stderr, "printLJ14Parameters errors -- aborting!\n" );
+   }
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print Shake parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" int printShakeParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* shakeFile ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nprintShakeParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   t_idef* idef                = &top->idef;
+   t_iatom* atoms              = idef->il[ F_SHAKE ].iatoms;
+   int numberOfConstraints     = idef->il[ F_SHAKE ].nr/3;
+
+  (void) fprintf( shakeFile, "# %d atomI atomJ d0 invMassI invMassJ\n", numberOfConstraints );
+
+   // loop over constraints 
+
+   int offset = 0;
+   for( int ii = 0; ii < numberOfConstraints; ii++ ) {
+
+      int type       = atoms[ offset     ]; 
+      int atomI      = atoms[ offset + 1 ]; 
+      int atomJ      = atoms[ offset + 2 ]; 
+      offset        += 3;
+
+      // atom indicies, distance constraint, inverse mass
+
+      (void) fprintf( shakeFile, "%6d %6d %6d %14.6e %14.6e %14.6e\n", ii, atomI, atomJ, idef->iparams[type].shake.dA, md->invmass[atomI], md->invmass[atomJ] );
+   }
+
+   return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print Gromacs parameters 
+   
+   @param baseFileName        base file name
+   @param top                 Gromacs parameter data struct
+   @param fr                  Gromacs parameter data struct
+   @param md                  Gromacs parameter data struct
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" void printReferenceParameters( char* baseFileName, t_topology *top, t_forcerec *fr,
+                                          t_mdatoms *md, FILE *log ){
+
+   // ---------------------------------------------------------------------------------------
+
+   typedef int printParameterFunction( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* file );
+
+   typedef struct {
+      char* parameterFileName;
+      printParameterFunction* printP;
+   } printParameters;
+
+   static const int ParameterCount = 7;
+   static printParameters printParameterList[7] = {
+      { "HarmonicBondParameter",    printBondParameters             },
+      { "AngleBondParameter",       printAngleParameters            },
+      { "RbDihedralParameter",      printRbDihedralParameters       },
+      { "ProperDihedralParameter",  printProperDihedralParameters   },
+      { "LJ14Parameter",            printLJ14Parameters             },
+      { "ShakeParameters",          printShakeParameters            },
+      { "LJCoulombParameter",       printLJCoulombParameters        }
+                                                  };
+
+   // ---------------------------------------------------------------------------------------
+
+   int ii;
+   static const int bufferSz = 1024;
+   char  parameterFileName[1024];
+   FILE* parameterFile;
+
+   // ---------------------------------------------------------------------------------------
+
+   // loop over parameter types 
+
+   for( ii = 0; ii < ParameterCount; ii++ ){
+
+      // open file
+
+#ifdef   WIN32
+      (void) sprintf_s( parameterFileName, bufferSz, "%s%s.txt", baseFileName, printParameterList[ii].parameterFileName );
+      fopen_s( &parameterFile, parameterFileName, "w" );
+#else
+      (void) sprintf( parameterFileName, "%s%s.txt", baseFileName, printParameterList[ii].parameterFileName );
+      parameterFile = fopen( parameterFileName, "w" );
+#endif
+
+      if( parameterFile == NULL ){
+         (void) fprintf( log, "Parameter file=<%s> could not be opened.\n", parameterFileName );
+      } else {
+
+         // do it
+
+         (*printParameterList[ii].printP)( top, fr, md, parameterFile );
+         (void) fclose( parameterFile );
+      }
+   }
+
+   return;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print forces/energies
+   
+   @param numberOfAtoms       number of atoms
+   @param x                   atom coordinates
+   @param f                   atom forces
+   @param energies            energies by atom
+   @param totalEnergy         total energy
+   @param fileName            file name
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" void writeGromacsForceFile(  int numberOfAtoms, rvec x[], rvec f[], real* energies,
+                                         real totalEnergy, char* fileName ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\ngromacsReferenceParameterPrint::writeGromacsForceFile";
+
+   // ---------------------------------------------------------------------------------------
+
+   // open file
+
+   FILE* forceFile;
+#ifdef   WIN32
+   fopen_s( &forceFile, fileName, "w" );
+#else
+   forceFile = fopen( fileName, "w" );
+#endif
+
+   if( forceFile == NULL ){
+      std::stringstream message;
+      message << methodName << " force file=<" << fileName << " could not be opened.\n";
+      SimTKOpenMMLog::printMessage( message );
+   } else {
+
+      // do it
+
+      (void) fprintf( forceFile, "# %6d E=%14.6e x[i][3] f[i][3] E[i]\n", numberOfAtoms, totalEnergy );
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         (void) fprintf( forceFile, "%5d %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", ii,
+                         x[ii][0], x[ii][1], x[ii][2],
+                         f[ii][0], f[ii][1], f[ii][2], (energies ? energies[ii] : 0.0) );
+      }
+      (void) fclose( forceFile );
+   }
+
+   return;
+}
+
+/**---------------------------------------------------------------------------------------
+      
+   Print 'state' mass/coordinates/velocities/forces
+   
+   @param numberOfAtoms       number of atoms
+   @param x                   atom coordinates
+   @param v                   atom velocities
+   @param f                   atom forces
+   @param timestep            time step
+   @param deltaT              delta t
+   @param temperature         temperature
+   @param tau                 tau
+   @param fileName            file name
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+extern "C" void writeGromacsStateFile( int numberOfAtoms, real mass[], rvec x[], rvec v[],
+                                       rvec f[], int timestep, real deltaT,
+                                       real temperature, real tau, char* fileName ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\ngromacsReferenceParameterPrint::writeGromacsStateFile";
+
+   // ---------------------------------------------------------------------------------------
+
+   // open file
+
+   FILE* stateFile;
+#ifdef   WIN32
+   fopen_s( &stateFile, fileName, "w" );
+#else
+   stateFile = fopen( fileName, "w" );
+#endif
+
+   if( stateFile == NULL ){
+
+      std::stringstream message;
+      message << methodName << " state file=<" << fileName << " could not be opened.\n";
+      SimTKOpenMMLog::printMessage( message );
+
+   } else {
+
+      // do it
+
+      (void) fprintf( stateFile, "# mass coordinates velocities forces\n", numberOfAtoms );
+      (void) fprintf( stateFile, "%6d         # Atoms\n", numberOfAtoms );
+      (void) fprintf( stateFile, "%6d         # timestep\n", timestep);
+      (void) fprintf( stateFile, "%14.6e # delta_t\n", deltaT );
+      (void) fprintf( stateFile, "%14.6e # temperature\n", temperature );
+      (void) fprintf( stateFile, "%14.6e # tau\n", tau );
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+         (void) fprintf( stateFile, "%6d ", ii );
+
+         if( mass ){
+            (void) fprintf( stateFile, "%14.6e   ", mass[ii] );
+         }
+         if( x ){
+            (void) fprintf( stateFile, "%14.6e %14.6e %14.6e    ", x[ii][0], x[ii][1], x[ii][2] );
+         }
+         if( v ){
+            (void) fprintf( stateFile, "%14.6e %14.6e %14.6e    ", v[ii][0], v[ii][1], v[ii][2] );
+         }
+         if( f ){
+            (void) fprintf( stateFile, "%14.6e %14.6e %14.6e    ", f[ii][0], f[ii][1], f[ii][2] );
+         }
+         (void) fprintf( stateFile, "\n" );
+      }
+      (void) fclose( stateFile );
+
+   }
+
+   return;
+}
+

platforms/reference/src/gromacsReferenceParameterPrint.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __GromacsReferenceParameterPrint_H__
+#define __GromacsReferenceParameterPrint_H__
+
+#ifdef __cplusplus
+#define externC extern "C"
+#else
+#define externC extern
+#endif
+
+
+/**---------------------------------------------------------------------------------------
+      
+   Print bond parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC int printBondParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* bondFile );
+ 
+/**---------------------------------------------------------------------------------------
+      
+   Print angle parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC int printAngleParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* angleFile );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print rbDihedral parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+externC int printRbDihedralParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* rbDihedralFile );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print properDihedral parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC int printProperDihedralParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* properDihedralFile );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print LJ 14 parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC int printLJ14Parameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* lj14File );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print LJ Coulomb parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC int printLJCoulombParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* ljCoulombFile );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print Shake parameters 
+   
+   @param idef                Gromacs parameter data struct
+   @param fr                  force record?
+   @param charges             atom charges
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC int printShakeParameters( t_topology *top, t_forcerec *fr, t_mdatoms *md, FILE* shakeFile );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print Gromacs parameters 
+   
+   @param baseFileName        base file name
+   @param top                 Gromacs parameter data struct
+   @param fr                  Gromacs parameter data struct
+   @param md                  Gromacs parameter data struct
+   @param file                FILE* reference
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC void printReferenceParameters( char* baseFileName, t_topology *top, t_forcerec *fr,
+                                           t_mdatoms *md, FILE *log );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print forces/energies
+   
+   @param numberOfAtoms       number of atoms
+   @param x                   atom coordinates
+   @param f                   atom forces
+   @param energies            energies by atom
+   @param totalEnergy         total energy
+   @param fileName            file name
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC void writeGromacsForceFile(  int numberOfAtoms, rvec x[], rvec f[], real* energies,
+                                         real totalEnergy, char* fileName );
+
+/**---------------------------------------------------------------------------------------
+      
+   Print 'state' mass/coordinates/velocities/forces
+   
+   @param numberOfAtoms       number of atoms
+   @param x                   atom coordinates
+   @param v                   atom velocities
+   @param f                   atom forces
+   @param timestep            time step
+   @param deltaT              delta t
+   @param temperature         temperature
+   @param tau                 tau
+   @param fileName            file name
+
+   @return 0
+      
+   --------------------------------------------------------------------------------------- */
+
+externC void writeGromacsStateFile( int numberOfAtoms, real mass[], rvec x[], rvec v[],
+                                       rvec f[], int timestep, real deltaT,
+                                       real temperature, real tau, char* fileName );
+
+
+#endif