Library to check that forces agree across platforms

libraries/validate/include/ValidateOpenMM.h

+#ifndef VALIDATE_OPENMM_H_
+#define VALIDATE_OPENMM_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenMM.h"
+#include "../../../platforms/reference/include/ReferencePlatform.h"
+#include "openmm/Context.h"
+#include "openmm/System.h"
+
+// free-energy plugin includes
+
+//#define	INCLUDE_FREE_ENERGY_PLUGIN
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+#include "../../../plugins/freeEnergy/openmmapi/include/OpenMMFreeEnergy.h"
+#include "../../../plugins/freeEnergy/openmmapi/include/openmm/freeEnergyKernels.h"
+#include "../../../plugins/freeEnergy/platforms/reference/include/ReferenceFreeEnergyPlatform.h"
+#endif
+
+#include <sstream>
+#include <typeinfo>
+
+#ifdef _MSC_VER
+   #define isinf !_finite
+   #define isnan _isnan
+#endif
+
+namespace OpenMM {
+
+typedef std::map< std::string, int > StringIntMap;
+typedef StringIntMap::iterator StringIntMapI;
+typedef StringIntMap::const_iterator StringIntMapCI;
+
+typedef std::map< std::string, double > StringDoubleMap;
+typedef StringDoubleMap::iterator StringDoubleMapI;
+typedef StringDoubleMap::const_iterator StringDoubleMapCI;
+
+typedef std::map< std::string, unsigned int > StringUIntMap;
+typedef StringUIntMap::iterator StringUIntMapI;
+typedef StringUIntMap::const_iterator StringUIntMapCI;
+
+typedef std::vector< std::string > StringVector;
+typedef StringVector::iterator StringVectorI;
+typedef StringVector::const_iterator StringVectorCI;
+
+typedef std::vector< int > IntVector;
+typedef IntVector::iterator IntVectorI;
+typedef IntVector::const_iterator IntVectorCI;
+
+typedef std::map< std::string, std::vector<std::string> > StringStringVectorMap;
+typedef StringStringVectorMap::iterator StringStringVectorMapI;
+typedef StringStringVectorMap::const_iterator StringStringVectorMapCI;
+
+/**
+ * Base class w/ common functionality 
+ */
+class ValidateOpenMM {
+public:
+
+    ValidateOpenMM( void );
+    ~ValidateOpenMM();
+
+    // force names
+    
+    static const std::string HARMONIC_BOND_FORCE;
+    static const std::string HARMONIC_ANGLE_FORCE;
+    static const std::string PERIODIC_TORSION_FORCE;
+    static const std::string RB_TORSION_FORCE;
+    
+    static const std::string NB_FORCE;
+    static const std::string NB_SOFTCORE_FORCE;
+    
+    static const std::string NB_EXCEPTION_FORCE;
+    static const std::string NB_EXCEPTION_SOFTCORE_FORCE;
+    
+    static const std::string GBSA_OBC_FORCE;
+    static const std::string GBSA_OBC_SOFTCORE_FORCE;
+    
+    static const std::string GBVI_FORCE;
+    static const std::string GBVI_SOFTCORE_FORCE;
+
+    static const std::string CM_MOTION_REMOVER;
+    static const std::string ANDERSEN_THERMOSTAT;
+    static const std::string CUSTOM_BOND_FORCE;
+    static const std::string CUSTOM_EXTERNAL_FORCE;
+    static const std::string CUSTOM_NONBONDED_FORCE;
+
+    /**
+     * Return true if input number is nan
+     *
+     * @param number   number to test
+     *
+     * @return true if number is nan
+     */
+    
+    static int isNan( double number );
+    
+    /**
+     * Get force name
+     * 
+     * @param force      OpenMM system force 
+     *
+     * @return force name or "NA" if force is not recognized
+     */
+    std::string getForceName(const Force& force ) const;
+
+    /**
+     * Copy force 
+     * 
+     * @param force      OpenMM system force to copy 
+     *
+     * @return force or NULL if not recognized
+     */
+    Force* copyForce(const Force& force) const;
+
+    /**
+     * Get copy of input system, but omit forces
+     * 
+     * @param systemToCopy   system to copy
+     *
+     * @return copy of system but w/o forces
+     */
+    System* copySystemExcludingForces( const System& systemToCopy ) const;
+
+    /**
+     *
+     * Set the velocities/positions of context2 to those of context1
+     * 
+     * @param context1                 context1 
+     * @param context2                 context2 
+     *
+     * @return 0
+     */
+    void synchContexts( const Context& context1, Context& context2 ) const;
+    
+    /**
+     *
+     * Get log FILE* reference
+     * 
+     * @return log
+     *
+     */
+    FILE* getLog( ) const;
+    
+    /**
+     *
+     * Set log FILE* reference
+     * 
+     * @param log  log
+     *
+     */
+    void setLog( FILE* log );
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy harmonic bond force
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyHarmonicBondForce( const HarmonicBondForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy harmonic angle
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyHarmonicAngleForce( const HarmonicAngleForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy PeriodicTorsionForce  
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyPeriodicTorsionForce( const PeriodicTorsionForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy RBTorsionForce
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyRBTorsionForce( const RBTorsionForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy NonbondedException force
+    
+       @param forceToCopy          force to copy
+       @param nonbondedForce       force to copy execeptions to
+       @param log                  log file pointer -- may be NULL
+    
+       --------------------------------------------------------------------------------------- */
+    
+    void copyNonbondedExceptions( const NonbondedForce& forceToCopy,
+                                  NonbondedForce* nonbondedForce, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy NonbondedSoftcoreExceptions force (free energy plugin force)
+    
+       @param forceToCopy          force to copy
+       @param nonbondedForce       force to copy execeptions to
+       @param log                  log file pointer -- may be NULL
+    
+       --------------------------------------------------------------------------------------- */
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+    void copyNonbondedSoftcoreExceptions( const NonbondedSoftcoreForce& forceToCopy,
+                                          NonbondedSoftcoreForce* nonbondedForce, FILE* log = NULL ) const;
+#endif
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy NonbondedForce
+    
+       @param forceToCopy          force to copy
+       @param nonbondedForce       force to copy execeptions to
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyNonbondedForce( const NonbondedForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy NonbondedSoftcoreForce (free energy plugin force)
+    
+       @param forceToCopy          force to copy
+       @param nonbondedForce       force to copy execeptions to
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+    Force* copyNonbondedSoftcoreForce( const NonbondedSoftcoreForce& forceToCopy, FILE* log = NULL ) const;
+#endif
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy GBSAOBCForce 
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyGBSAOBCForce( const GBSAOBCForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy GBSAOBCSoftcoreForce (free energy plugin force)
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+    Force* copyGBSAOBCSoftcoreForce( const GBSAOBCSoftcoreForce& forceToCopy, FILE* log = NULL ) const;
+#endif
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy GBVIForce
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+    Force* copyGBVIForce( const GBVIForce& forceToCopy, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy GBVISoftcoreForce (free energy plugin force)
+    
+       @param forceToCopy          force to copy
+       @param log                  log file pointer -- may be NULL
+    
+       @return copy of force
+    
+       --------------------------------------------------------------------------------------- */
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+    Force* copyGBVISoftcoreForce( const GBVISoftcoreForce& forceToCopy, FILE* log = NULL ) const;
+#endif
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Copy constraints
+    
+       @param systemToCopy         system whose constraints are to be copied
+       @param system               system to add constraints to
+       @param log                  log file pointer -- may be NULL
+    
+       --------------------------------------------------------------------------------------- */
+    
+    void copyConstraints( const System& systemToCopy, System* system, FILE* log = NULL ) const;
+    
+    /**---------------------------------------------------------------------------------------
+    
+       Get force dependencies
+    
+       @param forceName            force to check if there exist any dependencies
+       @param returnVector         vector of forces the input force is dependent on (example: GBSAOBC force requires Nonbonded force since
+                                   on Cuda platofrm they are computed in same loop and hence are inseparable)
+    
+       --------------------------------------------------------------------------------------- */
+    
+    void getForceDependencies( std::string forceName, StringVector& returnVector ) const;
+    
+    /**
+     * Write masses to parameter file
+     *
+     * @param filePtr    file to write masses to
+     * @param system     write masses in system
+     */
+    
+    void writeMasses( FILE* filePtr, const System& system ) const;
+    
+    /**
+     * Write constraints to parameter file
+     *
+     * @param filePtr    file to write constraints to
+     * @param system     write constraints in system
+     *
+     */
+    
+    void writeConstraints( FILE* filePtr, const System& system ) const;
+    
+    /**
+     * Write harmonicBondForce parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param harmonicBondForce     write harmonicBondForce parameters
+     *
+     */
+    
+    void writeHarmonicBondForce( FILE* filePtr, const HarmonicBondForce& harmonicBondForce ) const;
+
+    /**
+     * Write harmonicAngleForce parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param harmonicAngleForce     write harmonicAngleForce parameters
+     *
+     */
+    
+    void writeHarmonicAngleForce( FILE* filePtr, const HarmonicAngleForce& harmonicAngleForce ) const;
+    
+    /**
+     * Write rbTorsionForce parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param rbTorsionForce         write rbTorsionForce parameters
+     *
+     */
+    
+    void writeRbTorsionForce( FILE* filePtr, const RBTorsionForce& rbTorsionForce ) const;
+    
+    /**
+     * Write periodicTorsionForce parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param periodicTorsionForce   write periodicTorsionForce parameters
+     *
+     */
+    
+    void writePeriodicTorsionForce( FILE* filePtr, const PeriodicTorsionForce& periodicTorsionForce ) const;
+    
+    /**
+     * Write nonbonded parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param nonbondedForce         write nonbondedForce parameters
+     *
+     */
+    
+    void writeNonbondedForce( FILE* filePtr, const NonbondedForce & nonbondedForce ) const;
+    
+    /**
+     * Write GBSAOBCForce parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param gbsaObcForce           write gbsaObcForce parameters
+     *
+     */
+    
+    void writeGbsaObcForce( FILE* filePtr, const GBSAOBCForce& gbsaObcForce ) const;
+    
+    /**
+     * Write GBSA GB/VI Force parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param gbviObcForce           write gbviObcForce parameters
+     *
+     */
+    
+    void writeGBVIForce( FILE* filePtr, const GBVIForce& gbviForce ) const;
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+    
+    /**
+     * Write nonbonded softcore parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param nonbondedForce         write nonbondedForce parameters
+     */
+    
+    void writeNonbondedSoftcoreForce( FILE* filePtr, const NonbondedSoftcoreForce & nonbondedSoftcoreForce ) const;
+
+    /**
+     * Write GBSAOBCSoftcoreForce parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param gbsaObcForce           write gbsaObcForce parameters
+     *
+     */
+    
+    void writeGbsaObcSoftcoreForce( FILE* filePtr, const GBSAOBCSoftcoreForce& gbsaObcForce ) const;
+    
+    /**
+     * Write GBSA GB/VI softcore force parameters to file
+     *
+     * @param filePtr                file to write forces to
+     * @param gbviObcForce           write gbviObcForce parameters
+     *
+     */
+    
+    void writeGBVISoftcoreForce( FILE* filePtr, const GBVISoftcoreForce& gbviSoftcoreForce ) const;
+    
+#endif
+    
+    /**
+     * Write coordinates, velocities, ... to file
+     *
+     * @param filePtr                file to write Vec3 entries to
+     * @param vect3Array             write array of Vec3
+     *
+     */
+    
+    void writeVec3( FILE* filePtr, const std::vector<Vec3>& vect3Array ) const;
+    
+    /**
+     * Write context info to file (positions, velocities, forces, energies)
+     *
+     * @param filePtr                file to write entries to
+     * @param context                write context positions, velocities, forces, energies to file
+     *
+     */
+    
+    void writeContext( FILE* filePtr, const Context& context ) const;
+    
+    /**
+     * Write integrator
+     *
+     * @param filePtr                file to write integrator info to
+     * @param integrator             write integrator info (time step, seed, ... as applicable)
+     *
+     */
+    
+    void writeIntegrator( FILE* filePtr, const Integrator& integrator ) const;
+    
+    /**
+     * Write parameter file
+     * @param context                context whose entries are to be written to file
+     * @param parameterFileName      file name
+     *
+     */
+    
+    void writeParameterFile( const Context& context, const std::string& parameterFileName ) const;
+    
+private:
+
+    FILE* _log;
+
+    // map of force dependencies (e.g., GBSAObc requires NB force on CudaPlatform)
+
+    StringStringVectorMap _forceDependencies;
+};
+
+
+} // namespace OpenMM
+
+#endif /*VALIDATE_OPENMM_H_*/

libraries/validate/include/ValidateOpenMMForces.h

+#ifndef VALIDATE_OPENMM_FORCES_H_
+#define VALIDATE_OPENMM_FORCES_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Mark Friedrichs                                                   *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "ValidateOpenMM.h"
+
+namespace OpenMM {
+
+typedef std::map< int, int > MapIntInt;
+typedef MapIntInt::iterator MapIntIntI;
+typedef MapIntInt::const_iterator MapIntIntCI;
+
+// Helper class for ValidateOpenMMForces class used to store force results and facitilitate comparisons of
+// resulting forces
+
+class ForceValidationResult {
+public:
+
+    ForceValidationResult( const Context& context1, const Context& context2, StringUIntMap& forceNamesMap );
+    ~ForceValidationResult();
+
+    /**
+     * Get potential energy at specified platform index (0 || 1)
+     * 
+     * @return  potential energy for spercifed platform
+     *
+     * @throws OpenMMException if energyIndex is not 0 or 1
+     */
+    double getPotentialEnergy( int energyIndex  ) const; 
+
+    /**
+     * Get array of forces at specified platform index (0 || 1)
+     * 
+     * @return array of force norms
+     *
+     * @throws OpenMMException if forceIndex is not 0 or 1
+     */
+    std::vector<double> getForceNorms( int forceIndex  ) const; 
+
+    /**
+     * Get array of forces at platform index (0 || 1)
+     * 
+     * @return force array
+     *
+     * @throws OpenMMException if forceIndex is not 0 or 1
+     */
+    std::vector<Vec3> getForces( int forceIndex  ) const; 
+
+    /**
+     * Get maximum delta in force norm
+     * 
+     * @param maxIndex  return atom index of entry with maximum delta norm (optional)
+     *
+     * @return max delta in norm of forces
+     */
+    double getMaxDeltaForceNorm( int* maxIndex = NULL  ) const;
+
+    /**
+     * Get maximum relative delta in force norm
+     * 
+     * @param maxIndex  return atom index of entry w/ maximum relative delta norm (optional)
+     *
+     * @return max relative delta in norm of forces
+     */
+    double getMaxRelativeDeltaForceNorm( int* maxIndex = NULL  ) const;
+
+    /**
+     * Get maximum dot product between forces
+     * 
+     * @param maxIndex  return atom index of entry w/ maximum dot product between forces (optional)
+     *
+     * @return max dot product between forces
+     */
+    double getMaxDotProduct( int* maxIndex = NULL  ) const;
+
+    /**
+     * Get name of force associated w/ computed results
+     * 
+     * @return force name(s); if more than one force active in computation, 
+     * then names are concatenated and separated by '::' (e.g., 'NB_FORCE::GBSA_OBC_FORCE')
+     */
+    std::string getForceName( void ) const;
+
+    /**
+     * Get platform name
+     * 
+     * @param index index of platform (0 or 1)
+     *
+     * @return platform name
+     *
+     * @throws OpenMMException if index is not 0 or 1
+     */
+    std::string getPlatformName( int index ) const;
+
+    /**
+     * Register index of two entries that differ by a specified tolerance
+     * 
+     * @param index inconsistent index
+     *
+     */
+    void registerInconsistentForceIndex( int index, int value = 1 );
+
+    /**
+     * Clear list of entries that differ by a specified tolerance
+     * 
+     */
+    void clearInconsistentForceIndexList( void );
+
+    /**
+     * Get list of entries that differ by a specified tolerance
+     * 
+     */
+    void getInconsistentForceIndexList( std::vector<int>& inconsistentIndices ) const;
+
+    /**
+     * Get number of entries in inconsistent index list
+     * 
+     */
+    int getNumberOfInconsistentForceEntries( void ) const;
+
+    /**
+     * Return true if nans were detected
+     * 
+     * @return true if nans were detected
+     */
+    int nansDetected( void ) const;
+
+    /**
+     * Determine if force norms are valid
+     * 
+     * @param tolerance              tolerance
+     */
+    void compareForceNorms( double tolerance );
+
+    /**
+     * Determine if forces are valid
+     * 
+     * @param tolerance   tolerance
+     */
+    void compareForces( double tolerance );
+
+private:
+
+    // computed potential energies and forces fror two platforms
+
+    double _potentialEnergies[2];
+    std::vector<Vec3> _forces[2];
+
+    // platform and force names
+
+    std::string _platforms[2];
+    std::vector<std::string> _forceNames;
+
+    // force norms and stat entries
+
+    std::vector<double> _norms[2];
+    std::vector<double> _normStatVectors[2];
+
+    // map of indicies w/ inconsistent force entries
+
+    std::map<int, int> _inconsistentForceIndicies;
+
+    // if set, then nans detected
+
+    int _nansDetected;
+
+    /**
+     * Calculate norms of vectors
+     * 
+     */
+    void _calculateNorms( void );
+
+    /**
+     * Calculate norms of specified vector
+     * 
+     */
+    void _calculateNormOfForceVector( int forceIndex );
+
+    // stat indices
+
+    static const int STAT_AVG = 0;
+    static const int STAT_STD = 1;
+    static const int STAT_MIN = 2;
+    static const int STAT_ID1 = 3;
+    static const int STAT_MAX = 4;
+    static const int STAT_ID2 = 5;
+    static const int STAT_CNT = 6;
+    static const int STAT_END = 7;
+
+    /**
+     * Find vector stats
+     * 
+     */
+    void _findStatsForDouble( const std::vector<double>& array, std::vector<double>& statVector ) const;
+};
+
+// Class used to compare forces/potential energies on two platforms
+
+class ValidateOpenMMForces : public ValidateOpenMM {
+public:
+
+    ValidateOpenMMForces( void );
+    ~ValidateOpenMMForces();
+
+    /**
+     * Validate force/energy by comparing the results between the forces/energies computed on user-provided context platform
+     * with Reference platform
+     * 
+     * @param context          context reference
+     * @param summaryString    output summary string of results of comparison (optional)
+     *
+     * @return number of inconsistent entries
+     */
+     int compareWithReferencePlatform(Context& context, std::string* summaryString = NULL );
+
+    /**
+     * Validate force/energy by comparing the results between the forces/energies computed on two different platforms
+     * 
+     * @param context          context reference
+     * @param compareForces    indices of force to be tested
+     * @param platform1        first platform to compute forces
+     * @param platform2        second platform to compute forces
+     *
+     * @return ForceValidationResult reference containing results of force/energy computations
+     *         on the two input platforms
+     */
+     ForceValidationResult* compareForce(Context& context, std::vector<int>& compareForces,
+                                          Platform& platform1, Platform& platform2 ) const;
+
+    /**
+     * Compare individual forces by comparing calculations across two platforms (platform associated w/ input context and
+     * comparisonPlatform)
+     * 
+     * @param context                 context reference
+     * @param platform                comparsion platform reference 
+     * @param forceValidationResults  output vector of ForceValidationResult ptrs (user is responsible for deleting
+     *                                individual ForceValidationResult objects)
+     */
+    void compareOpenMMForces(Context& context, Platform& comparisonPlatform, std::vector<ForceValidationResult*>& forceValidationResults ) const;
+
+    /**
+     * Determine if results are consistent
+     * 
+     * @param forceValidationResults  vector of ForceValidationResult ptrs to check if forces are consistent
+     */
+    void checkForInconsistentForceEntries( std::vector<ForceValidationResult*>& forceValidationResults ) const;
+
+    /**
+     * Get total number of force entries that are inconsistent
+     * 
+     * @param forceValidationResults  vector of ForceValidationResult ptrs to check if forces are consistent
+     */
+    int getTotalNumberOfInconsistentForceEntries( std::vector<ForceValidationResult*>& forceValidationResults ) const;
+
+    /**
+     * Get summary string of results
+     * 
+     * @param forceValidationResults  vector of ForceValidationResult ptrs
+     */
+    std::string getSummary( std::vector<ForceValidationResult*>& forceValidationResults ) const;
+
+    /**
+     * Set force tolerance
+     * 
+     * @param tolerance     force tolerance
+     */
+    void setForceTolerance( double tolerance );
+
+    /**
+     * Get force tolerance
+     * 
+     * @return force tolerance
+     */
+    double getForceTolerance( void ) const;
+
+    /* 
+     * Get force tolerance for specified force
+     *
+     * @param forceName   name of force
+     *
+     * @return force tolerance
+     *
+     * */
+      
+    double getForceTolerance( const std::string& forceName ) const;
+       
+    /* 
+     * Get max errors to print in summary string
+     *
+     * @return max errors to print
+     *
+     * */
+      
+    int getMaxErrorsToPrint( void ) const;
+       
+    /* 
+     * Set max errors to print in summary string
+     *
+     * @param maxErrorsToPrint max errors to print
+     *
+     * */
+      
+    void setMaxErrorsToPrint( int maxErrorsToPrint );
+       
+    /* 
+     * Return true if force is not to be validated (Andersen thermostat, CM motion remover, ...)
+     *
+     * @param forceName   force name
+     *
+     * @return true if force is not currently validated
+     **/
+      
+    int isExcludedForce( std::string forceName ) const;
+       
+private:
+
+     // initialize class entries
+
+     void _initialize( void );
+
+      /* 
+       * Format output line
+       *
+       * @param tab         tab
+       * @param description description
+       * @param value       value
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string _getLine( const std::string& tab,
+                            const std::string& description,
+                            const std::string& value ) const;
+      
+     std::vector<ForceValidationResult*> _forceValidationResults;
+
+     // max errors to print
+
+     int _maxErrorsToPrint;
+
+     // tolerence
+
+     double _forceTolerance;
+
+     // map of force tolerances to type (name)
+
+     StringDoubleMap _forceTolerances;
+     
+     // forces to be excluded from validation
+
+     StringIntMap _forcesToBeExcluded;
+};
+
+} // namespace OpenMM
+
+#endif /*VALIDATE_OPENMM_FORCES_H_*/

libraries/validate/src/ValidateOpenMM.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Authors: Mark Friedrichs                                                   *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "ValidateOpenMM.h" 
+
+using namespace OpenMM;
+
+// fixed force names
+    
+const std::string ValidateOpenMM::HARMONIC_BOND_FORCE             = "HarmonicBond";
+const std::string ValidateOpenMM::HARMONIC_ANGLE_FORCE            = "HarmonicAngle";
+const std::string ValidateOpenMM::PERIODIC_TORSION_FORCE          = "PeriodicTorsion";
+const std::string ValidateOpenMM::RB_TORSION_FORCE                = "RbTorsion";
+
+const std::string ValidateOpenMM::NB_FORCE                        = "Nb";
+const std::string ValidateOpenMM::NB_SOFTCORE_FORCE               = "NbSoftcore";
+
+const std::string ValidateOpenMM::NB_EXCEPTION_FORCE              = "NbException";
+const std::string ValidateOpenMM::NB_EXCEPTION_SOFTCORE_FORCE     = "NbSoftcoreException";
+
+const std::string ValidateOpenMM::GBSA_OBC_FORCE                  = "Obc";
+const std::string ValidateOpenMM::GBSA_OBC_SOFTCORE_FORCE         = "ObcSoftcore";
+
+const std::string ValidateOpenMM::GBVI_FORCE                      = "GBVI";
+const std::string ValidateOpenMM::GBVI_SOFTCORE_FORCE             = "GBVISoftcore";
+
+const std::string ValidateOpenMM::CM_MOTION_REMOVER               = "CMMotionRemover";
+const std::string ValidateOpenMM::ANDERSEN_THERMOSTAT             = "AndersenThermostat";
+const std::string ValidateOpenMM::CUSTOM_BOND_FORCE               = "CustomBond";
+const std::string ValidateOpenMM::CUSTOM_EXTERNAL_FORCE           = "CustomExternal";
+const std::string ValidateOpenMM::CUSTOM_NONBONDED_FORCE          = "CustomNonBonded";
+
+
+ValidateOpenMM::ValidateOpenMM( void ) {
+
+    _log                                          = NULL;
+
+    // force dependencies 
+    // these may need to specialized depending on platform since
+    // currently apply to Cuda platform, but may not apply to OpenCL platform
+
+    std::vector< std::string > nbVector;
+    nbVector.push_back( NB_FORCE );
+    _forceDependencies[GBSA_OBC_FORCE]            = nbVector;
+    _forceDependencies[GBVI_FORCE]                = nbVector;
+
+    std::vector< std::string > nbSoftcoreVector;
+    nbSoftcoreVector.push_back( NB_SOFTCORE_FORCE );
+    _forceDependencies[GBSA_OBC_SOFTCORE_FORCE]   = nbSoftcoreVector;
+    _forceDependencies[GBVI_SOFTCORE_FORCE]       = nbSoftcoreVector;
+    
+}
+
+ValidateOpenMM::~ValidateOpenMM() {
+}
+
+int ValidateOpenMM::isNan( double number ){
+    return isinf( number ) || isnan( number ) ? 1 : 0;
+}
+
+FILE* ValidateOpenMM::getLog( void ) const {
+    return _log;
+}
+
+void ValidateOpenMM::setLog(  FILE* log ){
+    _log = log;
+}
+
+std::string ValidateOpenMM::getForceName(const Force& force) const { 
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ValidateForce::getForceName";
+
+// ---------------------------------------------------------------------------------------
+
+    std::string forceName  = "NA";
+
+    // bond force
+
+    try {
+        const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
+        return HARMONIC_BOND_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // angle force
+
+    try {
+        const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
+        return HARMONIC_ANGLE_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // periodic torsion force
+
+    try {
+        const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
+       return PERIODIC_TORSION_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // RB torsion force
+
+    try {
+        const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
+       return RB_TORSION_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // nonbonded force
+
+    try {
+        const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
+        return  NB_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // GBSA OBC
+
+    try {
+       const GBSAOBCForce& obcForce       = dynamic_cast<const GBSAOBCForce&>(force);
+       return GBSA_OBC_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // GB/VI
+
+    try {
+        const GBVIForce& obcForce  = dynamic_cast<const GBVIForce&>(force);
+        return GBVI_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+
+    // free energy plugin forces
+
+    // nonbonded softcore
+
+    try {
+        const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
+        return NB_SOFTCORE_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // GBSA OBC softcore
+
+    try {
+        const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
+        return GBSA_OBC_SOFTCORE_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // GB/VI softcore
+
+    try {
+        const GBVISoftcoreForce& gbviForce = dynamic_cast<const GBVISoftcoreForce&>(force);
+        return  GBVI_SOFTCORE_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+#endif
+
+    // CMMotionRemover
+
+    try {
+        const CMMotionRemover& cMMotionRemover  = dynamic_cast<const CMMotionRemover&>(force);
+        return CM_MOTION_REMOVER;
+    } catch( std::bad_cast ){
+    }
+
+    // AndersenThermostat
+
+    try {
+        const AndersenThermostat & andersenThermostat = dynamic_cast<const AndersenThermostat&>(force);
+        return ANDERSEN_THERMOSTAT;
+    } catch( std::bad_cast ){
+    }
+
+    // CustomBondForce
+
+    try {
+        const CustomBondForce & customBondForce = dynamic_cast<const CustomBondForce&>(force);
+        return CUSTOM_BOND_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // CustomExternalForce
+
+    try {
+        const CustomExternalForce& customExternalForce = dynamic_cast<const CustomExternalForce&>(force);
+        return CUSTOM_EXTERNAL_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    // CustomNonbondedForce
+
+    try {
+        const CustomNonbondedForce& customNonbondedForce = dynamic_cast<const CustomNonbondedForce&>(force);
+        return CUSTOM_NONBONDED_FORCE;
+    } catch( std::bad_cast ){
+    }
+
+    return forceName;
+}
+
+Force* ValidateOpenMM::copyForce(const Force& force) const { 
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ValidateForce::copyForce";
+
+// ---------------------------------------------------------------------------------------
+
+    // bond force
+
+    try {
+        const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
+        return copyHarmonicBondForce( harmonicBondForce );
+    } catch( std::bad_cast ){
+    }
+
+    // angle force
+
+    try {
+        const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
+        return copyHarmonicAngleForce( harmonicAngleForce );
+    } catch( std::bad_cast ){
+    }
+
+    // periodic torsion force
+
+    try {
+        const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
+        return copyPeriodicTorsionForce( periodicTorsionForce );
+    } catch( std::bad_cast ){
+    }
+
+    // RB torsion force
+
+    try {
+        const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
+        return copyRBTorsionForce( rBTorsionForce );
+    } catch( std::bad_cast ){
+    }
+
+    // nonbonded force
+
+    try {
+        const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
+        return copyNonbondedForce( nbForce ); 
+    } catch( std::bad_cast ){
+    }
+
+    // GBSA OBC
+
+    try {
+        const GBSAOBCForce& obcForce       = dynamic_cast<const GBSAOBCForce&>(force);
+        return copyGBSAOBCForce( obcForce );
+    } catch( std::bad_cast ){
+    }
+
+    // GB/VI
+
+    try {
+        const GBVIForce& gbviForce  = dynamic_cast<const GBVIForce&>(force);
+        return copyGBVIForce( gbviForce );
+    } catch( std::bad_cast ){
+    }
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+
+    // free energy plugin forces
+
+    // nonbonded softcore
+    try {
+         const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
+         return copyNonbondedSoftcoreForce( nbForce );
+    } catch( std::bad_cast ){
+    }
+
+    // GBSA OBC softcore
+
+    try {
+        const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
+        return copyGBSAOBCSoftcoreForce( obcForce );
+    } catch( std::bad_cast ){
+    }
+
+    // GB/VI softcore
+
+    try {
+        const GBVISoftcoreForce& gbviForce = dynamic_cast<const GBVISoftcoreForce&>(force);
+        return copyGBVISoftcoreForce( gbviForce );
+    } catch( std::bad_cast ){
+    }
+#endif
+
+    // CMMotionRemover
+
+    try {
+        const CMMotionRemover& cMMotionRemover  = dynamic_cast<const CMMotionRemover&>(force);
+        return NULL;
+    } catch( std::bad_cast ){
+    }
+
+    // AndersenThermostat
+
+    try {
+        const AndersenThermostat & andersenThermostat = dynamic_cast<const AndersenThermostat&>(force);
+        return NULL;
+    } catch( std::bad_cast ){
+    }
+
+    // CustomBondForce
+
+    try {
+        const CustomBondForce & customBondForce = dynamic_cast<const CustomBondForce&>(force);
+        return NULL;
+    } catch( std::bad_cast ){
+    }
+
+    // CustomExternalForce
+
+    try {
+        const CustomExternalForce& customExternalForce = dynamic_cast<const CustomExternalForce&>(force);
+        return NULL;
+    } catch( std::bad_cast ){
+    }
+
+    // CustomNonbondedForce
+
+    try {
+        const CustomNonbondedForce& customNonbondedForce = dynamic_cast<const CustomNonbondedForce&>(force);
+        return NULL;
+    } catch( std::bad_cast ){
+    }
+
+    return NULL;
+}
+
+/** 
+ * Get copy of input system, but leave out forces
+ * 
+ * @param systemToCopy   system to copy
+ * @return system  w/o forces
+ */
+System* ValidateOpenMM::copySystemExcludingForces( const System& systemToCopy ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ValidateOpenMM::copySystemExcludingForces";
+
+// ---------------------------------------------------------------------------------------
+
+    // create new system and add masses and box dimensions
+
+    System* system = new System();
+    for(int ii = 0; ii < systemToCopy.getNumParticles(); ii++ ){
+        double mass = systemToCopy.getParticleMass( ii );
+        system->addParticle( mass );
+    }   
+
+    Vec3 a, b, c;
+    systemToCopy.getPeriodicBoxVectors( a, b, c ); 
+    system->setPeriodicBoxVectors( a, b, c );
+
+    copyConstraints( systemToCopy, system );
+
+    return system;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set the velocities/positions of context2 to those of context1
+
+   @param context1                 context1 
+   @param context2                 context2 
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+void ValidateOpenMM::synchContexts( const Context& context1, Context& context2 ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const char* methodName  = "ValidateOpenMM::synchContexts: ";
+
+// ---------------------------------------------------------------------------------------
+
+    const State state                       = context1.getState(State::Positions | State::Velocities);
+    const std::vector<Vec3>& positions      = state.getPositions();
+    const std::vector<Vec3>& velocities     = state.getVelocities();
+ 
+    context2.setPositions( positions );
+    context2.setVelocities( velocities );
+ 
+    return;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy harmonic bond force
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyHarmonicBondForce( const HarmonicBondForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "copyHarmonicBondForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    HarmonicBondForce* bondForce = new HarmonicBondForce();
+    for( int ii = 0; ii < forceToCopy.getNumBonds(); ii++ ){
+        int particle1, particle2;
+        double length, k;
+        forceToCopy.getBondParameters( ii, particle1, particle2, length, k ); 
+        bondForce->addBond( particle1, particle2, length, k );
+    }
+
+    return bondForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy harmonic angle
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyHarmonicAngleForce( const HarmonicAngleForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "copyHarmonicAngleForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    HarmonicAngleForce* bondForce = new HarmonicAngleForce();
+    for( int ii = 0; ii < forceToCopy.getNumAngles(); ii++ ){
+        int particle1, particle2, particle3;
+        double angle, k;
+        forceToCopy.getAngleParameters( ii, particle1, particle2, particle3, angle, k ); 
+        bondForce->addAngle( particle1, particle2, particle3, angle, k );
+    }
+
+    return bondForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy PeriodicTorsionForce  
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyPeriodicTorsionForce( const PeriodicTorsionForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "copyPeriodicTorsionForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    PeriodicTorsionForce* bondForce = new PeriodicTorsionForce();
+    for( int ii = 0; ii < forceToCopy.getNumTorsions(); ii++ ){
+        int particle1, particle2, particle3, particle4, periodicity;
+        double phase, k;
+        forceToCopy.getTorsionParameters( ii, particle1, particle2, particle3, particle4, periodicity, phase, k );
+        bondForce->addTorsion( particle1, particle2, particle3, particle4, periodicity, phase, k );
+    }
+
+    return bondForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy RBTorsionForce
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyRBTorsionForce( const RBTorsionForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyRBTorsionForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    RBTorsionForce* bondForce = new RBTorsionForce();
+    for( int ii = 0; ii < forceToCopy.getNumTorsions(); ii++ ){
+        int particle1, particle2, particle3, particle4;
+        double c0, c1, c2, c3, c4, c5;
+        forceToCopy.getTorsionParameters( ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
+        bondForce->addTorsion( particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
+    }
+
+    return bondForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy NonbondedException force
+
+   @param forceToCopy          force to copy
+   @param nonbondedForce       force to copy execeptions to
+   @param log                  log file pointer -- may be NULL
+
+   --------------------------------------------------------------------------------------- */
+
+void ValidateOpenMM::copyNonbondedExceptions( const NonbondedForce& forceToCopy, NonbondedForce* nonbondedForce, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "copyNonbondedExceptions";
+   
+// ---------------------------------------------------------------------------------------
+
+    for( int ii = 0; ii < forceToCopy.getNumExceptions(); ii++ ){
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        forceToCopy.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon );
+        nonbondedForce->addException( particle1, particle2, chargeProd, sigma, epsilon );
+    }
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy NonbondedSoftcoreExceptions force
+
+   @param forceToCopy          force to copy
+   @param nonbondedForce       force to copy execeptions to
+   @param log                  log file pointer -- may be NULL
+
+   --------------------------------------------------------------------------------------- */
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+void ValidateOpenMM::copyNonbondedSoftcoreExceptions( const NonbondedSoftcoreForce& forceToCopy,
+                                                      NonbondedSoftcoreForce* nonbondedForce, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyNonbondedSoftcoreExceptions";
+   
+// ---------------------------------------------------------------------------------------
+
+    for( int ii = 0; ii < forceToCopy.getNumExceptions(); ii++ ){
+         int particle1, particle2;
+         double chargeProd, sigma, epsilon, softcoreLJLambda;
+         forceToCopy.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon, softcoreLJLambda );
+         nonbondedForce->addException( particle1, particle2, chargeProd, sigma, epsilon, softcoreLJLambda );
+    }
+
+    return;
+}
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Copy NonbondedForce
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyNonbondedForce( const NonbondedForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyNonbondedForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    NonbondedForce* nonbondedForce         = new NonbondedForce();
+
+    for( int ii = 0; ii < forceToCopy.getNumParticles(); ii++ ){
+        double charge, sigma, epsilon;
+        forceToCopy.getParticleParameters( ii, charge, sigma, epsilon );
+        nonbondedForce->addParticle( charge, sigma, epsilon );
+    }
+
+    nonbondedForce->setNonbondedMethod( forceToCopy.getNonbondedMethod() );
+    nonbondedForce->setCutoffDistance( forceToCopy.getCutoffDistance() );
+    nonbondedForce->setReactionFieldDielectric( forceToCopy.getReactionFieldDielectric() );
+    nonbondedForce->setEwaldErrorTolerance( forceToCopy.getEwaldErrorTolerance() );
+
+    copyNonbondedExceptions( forceToCopy, nonbondedForce, log );
+
+    return nonbondedForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy NonbondedSoftcoreForce
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+Force* ValidateOpenMM::copyNonbondedSoftcoreForce( const NonbondedSoftcoreForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyNonbondedSoftcoreForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    NonbondedSoftcoreForce* nonbondedSoftcore  = new NonbondedSoftcoreForce();
+
+    for( int ii = 0; ii < forceToCopy.getNumParticles(); ii++ ){
+         double charge, sigma, epsilon, softcoreLJLambda;
+         forceToCopy.getParticleParameters( ii, charge, sigma, epsilon, softcoreLJLambda );
+         nonbondedSoftcore->addParticle( charge, sigma, epsilon, softcoreLJLambda );
+    }
+
+    nonbondedSoftcore->setNonbondedMethod( forceToCopy.getNonbondedMethod() );
+    nonbondedSoftcore->setCutoffDistance( forceToCopy.getCutoffDistance() );
+    nonbondedSoftcore->setReactionFieldDielectric( forceToCopy.getReactionFieldDielectric() );
+    nonbondedSoftcore->setEwaldErrorTolerance( forceToCopy.getEwaldErrorTolerance() );
+
+    copyNonbondedSoftcoreExceptions( forceToCopy, nonbondedSoftcore, log );
+
+    return nonbondedSoftcore;
+}
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Copy GBSAOBCForce 
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyGBSAOBCForce( const GBSAOBCForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyGBSAOBCForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    GBSAOBCForce* gbsaObcForce = new GBSAOBCForce();
+    gbsaObcForce->setSoluteDielectric( forceToCopy.getSoluteDielectric() );
+    gbsaObcForce->setSolventDielectric( forceToCopy.getSolventDielectric() );
+
+    for( int ii = 0; ii < forceToCopy.getNumParticles(); ii++ ){
+         double charge, radius, scalingFactor;
+         forceToCopy.getParticleParameters( ii, charge, radius, scalingFactor );
+         gbsaObcForce->addParticle( charge, radius, scalingFactor );
+    }
+
+    return gbsaObcForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy GBSAOBCSoftcoreForce
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+Force* ValidateOpenMM::copyGBSAOBCSoftcoreForce( const GBSAOBCSoftcoreForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyGBSAOBCSoftcoreForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    GBSAOBCSoftcoreForce* gbsaObcSoftcoreForce   = new GBSAOBCSoftcoreForce();
+
+    gbsaObcSoftcoreForce->setSoluteDielectric( forceToCopy.getSoluteDielectric() );
+    gbsaObcSoftcoreForce->setSolventDielectric( forceToCopy.getSolventDielectric() );
+
+    for( int ii = 0; ii < forceToCopy.getNumParticles(); ii++ ){
+        double charge, radius, scalingFactor, nonpolarScaleFactor;
+        forceToCopy.getParticleParameters( ii, charge, radius, scalingFactor, nonpolarScaleFactor );
+        gbsaObcSoftcoreForce->addParticle( charge, radius, scalingFactor, nonpolarScaleFactor );
+    }
+
+    return gbsaObcSoftcoreForce;
+}
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Copy GBVIForce
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+Force* ValidateOpenMM::copyGBVIForce( const GBVIForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyGBVIForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    GBVIForce* gbviForce          = new GBVIForce();
+    gbviForce->setSoluteDielectric( forceToCopy.getSoluteDielectric() );
+    gbviForce->setSolventDielectric( forceToCopy.getSolventDielectric() );
+
+    for( int ii = 0; ii < forceToCopy.getNumParticles(); ii++ ){
+        double charge, radius, gamma;
+        forceToCopy.getParticleParameters( ii, charge, radius, gamma );
+        gbviForce->addParticle( charge, radius, gamma );
+    }
+    for( int ii = 0; ii < forceToCopy.getNumBonds(); ii++ ){
+        int atomI, atomJ;
+        double bondLength;
+        forceToCopy.getBondParameters( ii, atomI, atomJ, bondLength );
+        gbviForce->addBond( atomI, atomJ, bondLength );
+    }
+
+    return gbviForce;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Copy GBVISoftcoreForce
+
+   @param forceToCopy          force to copy
+   @param log                  log file pointer -- may be NULL
+
+   @return copy of force
+
+   --------------------------------------------------------------------------------------- */
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+Force* ValidateOpenMM::copyGBVISoftcoreForce( const GBVISoftcoreForce& forceToCopy, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyGBVISoftcoreForce";
+   
+// ---------------------------------------------------------------------------------------
+
+    GBVISoftcoreForce* gbviForce          = new GBVISoftcoreForce();
+    gbviForce->setSoluteDielectric( forceToCopy.getSoluteDielectric() );
+    gbviForce->setSolventDielectric( forceToCopy.getSolventDielectric() );
+    gbviForce->setBornRadiusScalingMethod( forceToCopy.getBornRadiusScalingMethod() );
+    gbviForce->setQuinticLowerLimitFactor( forceToCopy.getQuinticLowerLimitFactor() );
+    gbviForce->setQuinticUpperBornRadiusLimit( forceToCopy.getQuinticUpperBornRadiusLimit() );
+
+    for( int ii = 0; ii < forceToCopy.getNumParticles(); ii++ ){
+        double charge, radius, gamma, bornRadiusScaleFactor;
+        forceToCopy.getParticleParameters( ii, charge, radius, gamma, bornRadiusScaleFactor );
+        gbviForce->addParticle( charge, radius, gamma, bornRadiusScaleFactor );
+    }
+    for( int ii = 0; ii < forceToCopy.getNumBonds(); ii++ ){
+         int atomI, atomJ;
+         double bondLength;
+         forceToCopy.getBondParameters( ii, atomI, atomJ, bondLength );
+         gbviForce->addBond( atomI, atomJ, bondLength );
+    }
+
+    return gbviForce;
+}
+#endif
+
+/**---------------------------------------------------------------------------------------
+
+   Copy constraints
+
+   @param systemToCopy         system whose constraints are to be copied
+   @param system               system to add constraints to
+   @param log                  log file pointer -- may be NULL
+
+   --------------------------------------------------------------------------------------- */
+
+void ValidateOpenMM::copyConstraints( const System& systemToCopy, System* system, FILE* log ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "copyConstraints";
+   
+// ---------------------------------------------------------------------------------------
+
+    for( unsigned int ii = 0; ii < systemToCopy.getNumConstraints(); ii++ ){
+        int particle1, particle2;
+        double distance;
+        systemToCopy.getConstraintParameters( ii, particle1, particle2, distance ); 
+        system->addConstraint( particle1, particle2, distance );
+    }
+
+}
+
+/**---------------------------------------------------------------------------------------
+    
+   Get force dependencies
+    
+   @param forceName            force to check if there exist any dependencies
+   @param returnVector         vector of forces the input force is dependent on (example: GBSAOBC force requires Nonbonded force)
+    
+   --------------------------------------------------------------------------------------- */
+    
+void ValidateOpenMM::getForceDependencies( std::string forceName, StringVector& returnVector ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "getForceDependencies";
+   
+// ---------------------------------------------------------------------------------------
+
+    StringStringVectorMapCI dependencyVector = _forceDependencies.find( forceName );
+    if( dependencyVector != _forceDependencies.end() ){
+        StringVector dependices = (*dependencyVector).second;
+        for( unsigned int ii = 0; ii < dependices.size(); ii++ ){
+            returnVector.push_back( dependices[ii] );
+        }
+    }
+}
+
+/**
+ * Write masses and box dimensions to parameter file
+ */
+
+void ValidateOpenMM::writeMasses( FILE* filePtr, const System& system ) const {
+
+    (void) fprintf( filePtr, "Masses %d\n", system.getNumParticles() );
+    for(int ii = 0; ii < system.getNumParticles(); ii++ ){
+       double mass = system.getParticleMass( ii );
+       (void) fprintf( filePtr, "%8d %14.7e\n", ii, mass );
+    }
+
+    Vec3 a, b, c;
+    system.getPeriodicBoxVectors( a, b, c);
+    (void) fprintf( filePtr, "Box %14.6f %14.6f %14.6f  %14.6f %14.6f %14.6f   %14.6f %14.6f %14.6f\n",
+                    a[0], a[1], a[2], b[0], b[1], b[2], c[0], c[1], c[2] );
+}
+
+/**
+ * Write constraints to parameter file
+ */
+
+void ValidateOpenMM::writeConstraints( FILE* filePtr, const System& system ) const {
+
+    (void) fprintf( filePtr, "Constraints %d\n", system.getNumConstraints() );
+    for(int ii = 0; ii < system.getNumConstraints(); ii++ ){
+       int particle1, particle2;
+       double distance;
+       system.getConstraintParameters( ii, particle1, particle2, distance );
+       (void) fprintf( filePtr, "%8d %8d %8d %14.7e\n", ii, particle1, particle2, distance );
+    }
+}
+
+/**
+ * Write harmonicBondForce parameters to file
+ */
+
+void ValidateOpenMM::writeHarmonicBondForce( FILE* filePtr, const HarmonicBondForce& harmonicBondForce ) const {
+
+    (void) fprintf( filePtr, "HarmonicBondForce %d\n", harmonicBondForce.getNumBonds() );
+    for(int ii = 0; ii < harmonicBondForce.getNumBonds(); ii++ ){
+       int particle1, particle2;
+       double length, k;
+       harmonicBondForce.getBondParameters( ii, particle1, particle2, length, k );
+       (void) fprintf( filePtr, "%8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, length, k);
+    }
+}
+
+/**
+ * Write harmonicAngleForce parameters to file
+ */
+
+void ValidateOpenMM::writeHarmonicAngleForce( FILE* filePtr, const HarmonicAngleForce& harmonicAngleForce ) const {
+
+    (void) fprintf( filePtr, "HarmonicAngleForce %d\n", harmonicAngleForce.getNumAngles() );
+    for(int ii = 0; ii < harmonicAngleForce.getNumAngles(); ii++ ){
+       int particle1, particle2, particle3;
+       double angle, k;
+       harmonicAngleForce.getAngleParameters( ii, particle1, particle2, particle3, angle, k );
+       (void) fprintf( filePtr, "%8d %8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, particle3, angle, k);
+    }
+}
+
+/**
+ * Write rbTorsionForce parameters to file
+ */
+
+void ValidateOpenMM::writeRbTorsionForce( FILE* filePtr, const RBTorsionForce& rbTorsionForce ) const {
+
+    (void) fprintf( filePtr, "RBTorsionForce %d\n", rbTorsionForce.getNumTorsions() );
+    for(int ii = 0; ii < rbTorsionForce.getNumTorsions(); ii++ ){
+       int particle1, particle2, particle3, particle4;
+       double c0, c1, c2, c3, c4, c5;
+       rbTorsionForce.getTorsionParameters( ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
+       (void) fprintf( filePtr, "%8d %8d %8d %8d %8d %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
+    }
+}
+
+/**
+ * Write periodicTorsionForce parameters to file
+ */
+
+void ValidateOpenMM::writePeriodicTorsionForce( FILE* filePtr, const PeriodicTorsionForce& periodicTorsionForce ) const {
+
+    (void) fprintf( filePtr, "PeriodicTorsionForce %d\n", periodicTorsionForce.getNumTorsions() );
+    for(int ii = 0; ii < periodicTorsionForce.getNumTorsions(); ii++ ){
+       int particle1, particle2, particle3, particle4, periodicity;
+       double phase, k;
+       periodicTorsionForce.getTorsionParameters( ii, particle1, particle2, particle3, particle4, periodicity, phase, k );
+       (void) fprintf( filePtr, "%8d %8d %8d %8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, particle3, particle4, periodicity, phase, k );
+    }
+}
+
+/**
+ * Write GBSAOBCForce parameters to file
+ */
+
+void ValidateOpenMM::writeGbsaObcForce( FILE* filePtr, const GBSAOBCForce& gbsaObcForce ) const {
+
+    (void) fprintf( filePtr, "GBSAOBCForce %d\n", gbsaObcForce.getNumParticles() );
+    for(int ii = 0; ii <  gbsaObcForce.getNumParticles(); ii++ ){
+       double charge, radius, scalingFactor;
+       gbsaObcForce.getParticleParameters( ii, charge, radius, scalingFactor );
+       (void) fprintf( filePtr, "%8d  %14.7e %14.7e %14.7e\n", ii, charge, radius, scalingFactor );
+    }
+    (void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbsaObcForce.getSoluteDielectric() );
+    (void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbsaObcForce.getSolventDielectric() );
+}
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+
+/**
+ * Write GBSAOBCSoftcoreForce parameters to file
+ */
+
+void ValidateOpenMM::writeGbsaObcSoftcoreForce( FILE* filePtr, const GBSAOBCSoftcoreForce& gbsaObcForce ) const {
+
+    (void) fprintf( filePtr, "GBSAOBCSoftcoreForce %d\n", gbsaObcForce.getNumParticles() );
+    for(int ii = 0; ii <  gbsaObcForce.getNumParticles(); ii++ ){
+       double charge, radius, scalingFactor;
+       double nonPolarScalingFactor;
+       gbsaObcForce.getParticleParameters( ii, charge, radius, scalingFactor, nonPolarScalingFactor );
+       (void) fprintf( filePtr, "%8d  %14.7e %14.7e %14.7e %14.7e\n", ii, charge, radius, scalingFactor, nonPolarScalingFactor );
+    }
+    (void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbsaObcForce.getSoluteDielectric() );
+    (void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbsaObcForce.getSolventDielectric() );
+}
+
+#endif
+
+/**
+ * Write GBSA GB/VI Force parameters to file
+ */
+
+void ValidateOpenMM::writeGBVIForce( FILE* filePtr, const GBVIForce& gbviForce ) const {
+
+    (void) fprintf( filePtr, "GBVIForce %d\n", gbviForce.getNumParticles() );
+    for(int ii = 0; ii <  gbviForce.getNumParticles(); ii++ ){
+       double charge, radius, gamma;
+       gbviForce.getParticleParameters( ii, charge, radius, gamma );
+       (void) fprintf( filePtr, "%8d  %14.7e %14.7e %14.7e\n", ii, charge, radius, gamma );
+    }
+
+    (void) fprintf( filePtr, "GBVIBonds  %d\n", gbviForce.getNumBonds() );
+    for(int ii = 0; ii <  gbviForce.getNumBonds(); ii++ ){
+       int atomI, atomJ;
+       double bondLength;
+       gbviForce.getBondParameters( ii, atomI, atomJ, bondLength );
+       (void) fprintf( filePtr, "%8d %8d %8d %14.7e\n", ii, atomI, atomJ, bondLength );
+    }
+    (void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbviForce.getSoluteDielectric() );
+    (void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbviForce.getSolventDielectric() );
+
+}
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+
+/**
+ * Write GBSA GB/VI softcore force parameters to file
+ */
+
+void ValidateOpenMM::writeGBVISoftcoreForce( FILE* filePtr, const GBVISoftcoreForce& gbviSoftcoreForce ) const {
+
+    (void) fprintf( filePtr, "GBVISoftcoreForce %d\n", gbviSoftcoreForce.getNumParticles() );
+    for(int ii = 0; ii <  gbviSoftcoreForce.getNumParticles(); ii++ ){
+       double charge, radius, gamma;
+       double nonPolarScalingFactor;
+       gbviSoftcoreForce.getParticleParameters( ii, charge, radius, gamma, nonPolarScalingFactor );
+       (void) fprintf( filePtr, "%8d  %14.7e %14.7e %14.7e %14.7e\n", ii, charge, radius, gamma, nonPolarScalingFactor );
+    }
+
+    (void) fprintf( filePtr, "GBVISoftcoreBonds  %d\n", gbviSoftcoreForce.getNumBonds() );
+    for(int ii = 0; ii <  gbviSoftcoreForce.getNumBonds(); ii++ ){
+       int atomI, atomJ;
+       double bondLength;
+       gbviSoftcoreForce.getBondParameters( ii, atomI, atomJ, bondLength );
+       (void) fprintf( filePtr, "%8d %8d %8d %14.7e\n", ii, atomI, atomJ, bondLength );
+    }
+    (void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbviSoftcoreForce.getSoluteDielectric() );
+    (void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbviSoftcoreForce.getSolventDielectric() );
+
+    (void) fprintf( filePtr, "BornRadiusScalingMethod %d\n", gbviSoftcoreForce.getBornRadiusScalingMethod() );
+    (void) fprintf( filePtr, "QuinticLowerLimitFactor %14.7e\n", gbviSoftcoreForce.getQuinticLowerLimitFactor() );
+    (void) fprintf( filePtr, "QuinticUpperBornRadiusLimit %14.7e\n", gbviSoftcoreForce.getQuinticUpperBornRadiusLimit() );
+}
+
+#endif
+
+/**
+ * Write coordinates to file
+ */
+
+void ValidateOpenMM::writeVec3( FILE* filePtr, const std::vector<Vec3>& vect3Array ) const {
+
+    for( unsigned int ii = 0; ii < vect3Array.size(); ii++ ){
+       (void) fprintf( filePtr, "%8d  %14.7e %14.7e %14.7e\n", ii, 
+                       vect3Array[ii][0], vect3Array[ii][1], vect3Array[ii][2] );
+    }
+}
+
+/**
+ * Write context info to file
+ */
+
+void ValidateOpenMM::writeContext( FILE* filePtr, const Context& context ) const {
+
+    State state                  = context.getState( State::Positions | State::Velocities | State::Forces | State::Energy );
+    std::vector<Vec3> positions  = state.getPositions();
+    std::vector<Vec3> velocities = state.getVelocities();
+    std::vector<Vec3> forces     = state.getForces();
+
+    (void) fprintf( filePtr, "Positions %u\n", positions.size() );
+    writeVec3( filePtr, positions );
+
+    (void) fprintf( filePtr, "Velocities %u\n", velocities.size() );
+    writeVec3( filePtr, velocities );
+
+    (void) fprintf( filePtr, "Forces %u\n", forces.size() );
+    writeVec3( filePtr, forces );
+
+    (void) fprintf( filePtr, "KineticEnergy %14.7e\n", state.getKineticEnergy() );
+    (void) fprintf( filePtr, "PotentialEnergy %14.7e\n", state.getPotentialEnergy() );
+}
+
+/**
+ * Write nonbonded parameters to file
+ */
+
+void ValidateOpenMM::writeNonbondedForce( FILE* filePtr, const NonbondedForce & nonbondedForce ) const {
+
+    // charge and vdw parameters
+
+    (void) fprintf( filePtr, "NonbondedForce %d\n", nonbondedForce.getNumParticles() );
+    for(int ii = 0; ii < nonbondedForce.getNumParticles(); ii++ ){
+       double charge, sigma, epsilon;
+       nonbondedForce.getParticleParameters( ii, charge, sigma, epsilon );
+       (void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e\n", ii, charge, sigma, epsilon );
+    }
+
+    // cutoff, dielectric, Ewald tolerance
+
+    (void) fprintf( filePtr, "CutoffDistance %14.7e\n", nonbondedForce.getCutoffDistance() );
+    (void) fprintf( filePtr, "RFDielectric %14.7e\n", nonbondedForce.getReactionFieldDielectric() );
+    (void) fprintf( filePtr, "EwaldRTolerance %14.7e\n", nonbondedForce.getEwaldErrorTolerance() );
+
+    // cutoff mode
+
+    std::string nonbondedForceMethod;
+    switch( nonbondedForce.getNonbondedMethod() ){
+        case NonbondedForce::NoCutoff:
+            nonbondedForceMethod = "NoCutoff";
+            break;
+        case NonbondedForce::CutoffNonPeriodic:
+            nonbondedForceMethod = "CutoffNonPeriodic";
+            break;
+        case NonbondedForce::CutoffPeriodic:
+            nonbondedForceMethod = "CutoffPeriodic";
+            break;
+        case NonbondedForce::Ewald:
+            nonbondedForceMethod = "Ewald";
+            break;
+        case NonbondedForce::PME:
+            nonbondedForceMethod = "PME";
+            break;
+        default:
+            nonbondedForceMethod = "Unknown";
+    }
+    (void) fprintf( filePtr, "NonbondedForceMethod %s\n", nonbondedForceMethod.c_str() );
+
+    (void) fprintf( filePtr, "NonbondedForceExceptions %d\n", nonbondedForce.getNumExceptions() );
+    for(int ii = 0; ii < nonbondedForce.getNumExceptions(); ii++ ){
+       int particle1, particle2;
+       double chargeProd, sigma, epsilon;
+       nonbondedForce.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon );
+       (void) fprintf( filePtr, "%8d %8d %8d %14.7e %14.7e %14.7e\n", ii, particle1, particle2, chargeProd, sigma, epsilon );
+    }
+
+}
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+
+/**
+ * Write nonbonded softcore parameters to file
+ */
+
+void ValidateOpenMM::writeNonbondedSoftcoreForce( FILE* filePtr, const NonbondedSoftcoreForce & nonbondedSoftcoreForce ) const {
+
+    (void) fprintf( filePtr, "NonbondedSoftcoreForce %d\n", nonbondedSoftcoreForce.getNumParticles() );
+    for(int ii = 0; ii < nonbondedSoftcoreForce.getNumParticles(); ii++ ){
+       double charge, sigma, epsilon, softCoreLJ;
+       nonbondedSoftcoreForce.getParticleParameters( ii, charge, sigma, epsilon, softCoreLJ );
+       (void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e %14.7e\n", ii, charge, sigma, epsilon, softCoreLJ );
+    }
+
+    (void) fprintf( filePtr, "CutoffDistance %14.7e\n", nonbondedSoftcoreForce.getCutoffDistance() );
+    (void) fprintf( filePtr, "RFDielectric %14.7e\n", nonbondedSoftcoreForce.getReactionFieldDielectric() );
+    (void) fprintf( filePtr, "EwaldRTolerance %14.7e\n", nonbondedSoftcoreForce.getEwaldErrorTolerance() );
+
+    std::string nonbondedSoftcoreForceMethod;
+    switch( nonbondedSoftcoreForce.getNonbondedMethod() ){
+        case NonbondedSoftcoreForce::NoCutoff:
+            nonbondedSoftcoreForceMethod = "NoCutoff";
+            break;
+/*
+        case NonbondedSoftcoreForce::CutoffNonPeriodic:
+            nonbondedSoftcoreForceMethod = "CutoffNonPeriodic";
+            break;
+        case NonbondedSoftcoreForce::CutoffPeriodic:
+            nonbondedSoftcoreForceMethod = "CutoffPeriodic";
+            break;
+        case NonbondedSoftcoreForce::Ewald:
+            nonbondedSoftcoreForceMethod = "Ewald";
+            break;
+        case NonbondedSoftcoreForce::PME:
+            nonbondedSoftcoreForceMethod = "PME";
+            break;
+*/
+        default:
+            nonbondedSoftcoreForceMethod = "Unknown";
+    }
+    (void) fprintf( filePtr, "NonbondedSoftcoreForceMethod %s\n", nonbondedSoftcoreForceMethod.c_str() );
+
+    (void) fprintf( filePtr, "NonbondedSoftcoreForceExceptions %d\n", nonbondedSoftcoreForce.getNumExceptions() );
+    for(int ii = 0; ii < nonbondedSoftcoreForce.getNumExceptions(); ii++ ){
+       int particle1, particle2;
+       double chargeProd, sigma, epsilon, softCoreLJ;
+       nonbondedSoftcoreForce.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon, softCoreLJ );
+       (void) fprintf( filePtr, "%8d %8d %8d %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, chargeProd, sigma, epsilon, softCoreLJ );
+    }
+
+    return;
+}
+#endif
+
+void ValidateOpenMM::writeIntegrator( FILE* filePtr, const Integrator& integrator ) const {
+
+    // LangevinIntegrator
+ 
+    try {
+        const LangevinIntegrator langevinIntegrator = dynamic_cast<const LangevinIntegrator&>(integrator);
+        (void) fprintf( filePtr, "Integrator LangevinIntegrator\n" );
+        (void) fprintf( filePtr, "StepSize %14.7e\n", langevinIntegrator.getStepSize() );
+        (void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", langevinIntegrator.getConstraintTolerance() );
+        (void) fprintf( filePtr, "Temperature %14.7e\n", langevinIntegrator.getTemperature() );
+        (void) fprintf( filePtr, "Friction %14.7e\n", langevinIntegrator.getFriction() );
+        (void) fprintf( filePtr, "RandomNumberSeed %d\n", langevinIntegrator.getRandomNumberSeed() );
+        return;
+    } catch( std::bad_cast ){
+    }
+ 
+    // VariableLangevinIntegrator
+ 
+    try {
+        const VariableLangevinIntegrator& langevinIntegrator = dynamic_cast<const VariableLangevinIntegrator&>(integrator);
+        (void) fprintf( filePtr, "Integrator VariableLangevinIntegrator\n" );
+        (void) fprintf( filePtr, "StepSize %14.7e\n", langevinIntegrator.getStepSize() );
+        (void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", langevinIntegrator.getConstraintTolerance() );
+        (void) fprintf( filePtr, "Temperature %14.7e\n", langevinIntegrator.getTemperature() );
+        (void) fprintf( filePtr, "Friction %14.7e\n", langevinIntegrator.getFriction() );
+        (void) fprintf( filePtr, "RandomNumberSeed %d\n", langevinIntegrator.getRandomNumberSeed() );
+        (void) fprintf( filePtr, "ErrorTolerance %14.7e\n", langevinIntegrator.getErrorTolerance() );
+        return;
+    } catch( std::bad_cast ){
+    }
+ 
+    // VerletIntegrator
+ 
+    try {
+        const VerletIntegrator& verletIntegrator = dynamic_cast<const VerletIntegrator&>(integrator);
+        (void) fprintf( filePtr, "Integrator VerletIntegrator\n" );
+        (void) fprintf( filePtr, "StepSize %14.7e\n", verletIntegrator.getStepSize() );
+        (void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", verletIntegrator.getConstraintTolerance() );
+        return;
+    } catch( std::bad_cast ){
+    }
+     
+    // VariableVerletIntegrator
+ 
+    try {
+        const VariableVerletIntegrator & variableVerletIntegrator = dynamic_cast<const VariableVerletIntegrator&>(integrator);
+        (void) fprintf( filePtr, "Integrator VariableVerletIntegrator\n" );
+        (void) fprintf( filePtr, "StepSize %14.7e\n", variableVerletIntegrator.getStepSize() );
+        (void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", variableVerletIntegrator.getConstraintTolerance() );
+        (void) fprintf( filePtr, "ErrorTolerance %14.7e\n", variableVerletIntegrator.getErrorTolerance() );
+        return;
+    } catch( std::bad_cast ){
+    }
+     
+    // BrownianIntegrator
+ 
+    try {
+        const BrownianIntegrator& brownianIntegrator = dynamic_cast<const BrownianIntegrator&>(integrator);
+        (void) fprintf( filePtr, "Integrator BrownianIntegrator\n" );
+        (void) fprintf( filePtr, "StepSize %14.7e\n", brownianIntegrator.getStepSize() );
+        (void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", brownianIntegrator.getConstraintTolerance() );
+        (void) fprintf( filePtr, "Temperature %14.7e\n", brownianIntegrator.getTemperature() );
+        (void) fprintf( filePtr, "Friction %14.7e\n", brownianIntegrator.getFriction() );
+        (void) fprintf( filePtr, "RandomNumberSeed %d\n", brownianIntegrator.getRandomNumberSeed() );
+        return;
+    } catch( std::bad_cast ){
+    }
+     
+    if( getLog() ){
+       (void) fprintf( getLog(), "Integrator not recognized." );
+       (void) fflush( getLog() );
+    }
+
+    return;
+}
+
+/**
+ * Write parameter file
+ */
+
+void ValidateOpenMM::writeParameterFile( const Context& context, const std::string& parameterFileName ) const {
+
+    // open file
+ 
+    FILE* filePtr = fopen( parameterFileName.c_str(), "w" );
+
+    const System& system          = context.getSystem();
+    const Integrator& integrator  = context.getIntegrator();
+
+    // (void) fprintf( filePtr, "Version %s\n", versionString.c_str() );
+    (void) fprintf( filePtr, "Particles %8d\n", system.getNumParticles() );
+    writeMasses( filePtr, system );
+    (void) fprintf( filePtr, "NumberOfForces    %8d\n", system.getNumForces() );
+
+    // print active forces and relevant parameters
+
+    for(int i = 0; i < system.getNumForces(); ++i) {
+
+        int hit                       = 0;
+        const Force& force            = system.getForce(i);
+
+        // bond
+
+        if( !hit ){
+
+            try {
+               const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
+               writeHarmonicBondForce( filePtr, harmonicBondForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+    
+        // angle
+
+        if( !hit ){
+    
+            try {
+               const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
+               writeHarmonicAngleForce( filePtr, harmonicAngleForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+
+        // PeriodicTorsionForce
+    
+        if( !hit ){
+    
+            try {
+               const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
+               writePeriodicTorsionForce( filePtr, periodicTorsionForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+    
+        // RBTorsionForce
+    
+        if( !hit ){
+            try {
+               const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
+               writeRbTorsionForce(  filePtr, rBTorsionForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+    
+        // nonbonded
+    
+        if( !hit ){
+            try {
+               const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
+               writeNonbondedForce( filePtr, nbForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        } 
+
+        // nonbonded softcore
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+        if( !hit ){
+            try {
+               const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
+               writeNonbondedSoftcoreForce( filePtr, nbForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        } 
+#endif
+
+        // GBSA OBC
+    
+        if( !hit ){
+            try {
+               const GBSAOBCForce& obcForce = dynamic_cast<const GBSAOBCForce&>(force);
+               writeGbsaObcForce(  filePtr, obcForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+    
+
+        // GBSA OBC softcore
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+        if( !hit ){
+            try {
+               const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
+               writeGbsaObcSoftcoreForce(  filePtr, obcForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+#endif
+    
+        // GB/VI
+    
+        if( !hit ){
+            try {
+               const GBVIForce& obcForce = dynamic_cast<const GBVIForce&>(force);
+               writeGBVIForce(  filePtr, obcForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+    
+        // GB/VI softcore
+    
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+        if( !hit ){
+            try {
+               const GBVISoftcoreForce& obcForce = dynamic_cast<const GBVISoftcoreForce&>(force);
+               writeGBVISoftcoreForce(  filePtr, obcForce );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+#endif
+    
+        // COM
+
+        if( !hit ){
+    
+            try {
+               const CMMotionRemover& cMMotionRemover = dynamic_cast<const CMMotionRemover&>(force);
+               (void) fprintf( filePtr, "CMMotionRemover %d\n", cMMotionRemover.getFrequency() );
+               hit++;
+            } catch( std::bad_cast ){
+            }
+        }
+
+        if( !hit ){
+           (void) fprintf( stderr, "   %2d force not recognized.\n", i );
+           exit(-1);
+        }
+
+    }
+
+    // constraints
+
+    writeConstraints( filePtr, system );
+
+    // context
+
+    writeContext( filePtr, context );
+
+    // integrator 
+
+    writeIntegrator( filePtr, integrator );
+
+    // close file
+
+    (void) fclose( filePtr );
+}

libraries/validate/src/ValidateOpenMMForces.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "ValidateOpenMMForces.h"
+
+using namespace OpenMM;
+
+ForceValidationResult::ForceValidationResult( const Context& context1, const Context& context2, StringUIntMap& forceNamesMap  ){
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::ForceValidationResult";
+
+// ---------------------------------------------------------------------------------------
+
+    _nansDetected                             = 0;
+
+    // calculate forces/energies
+
+    int types                                 = State::Forces | State::Energy;
+    State state1                              = context1.getState( types );
+    State state2                              = context2.getState( types );
+
+    _potentialEnergies[0]                     = state1.getPotentialEnergy();
+    _potentialEnergies[1]                     = state2.getPotentialEnergy();
+    if( ValidateOpenMM::isNan( _potentialEnergies[0] ) || ValidateOpenMM::isNan( _potentialEnergies[1] ) ){
+       _nansDetected++;
+    }
+
+    _forces[0]                                = state1.getForces();
+    _forces[1]                                = state2.getForces();
+
+    _platforms[0]                             = context1.getPlatform().getName();
+    _platforms[1]                             = context2.getPlatform().getName();
+
+    for( StringUIntMapI ii = forceNamesMap.begin(); ii != forceNamesMap.end(); ii++ ){
+       _forceNames.push_back( (*ii).first );
+    }
+    _calculateNorms();
+}
+
+ForceValidationResult::~ForceValidationResult(  ){
+}
+
+void ForceValidationResult::_calculateNorms( void ){
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::_calculateNorms";
+
+// ---------------------------------------------------------------------------------------
+
+    for( int jj = 0; jj < 2; jj++ ){
+        if( _norms[jj].size() < 1 ){
+            _calculateNormOfForceVector( jj );
+            _findStatsForDouble( _norms[jj], _normStatVectors[jj] );
+        }
+    }
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Find stats for double array
+
+   @param array                   array 
+   @param statVector              vector of stats 
+
+   @return 0
+
+   --------------------------------------------------------------------------------------- */
+
+void ForceValidationResult::_findStatsForDouble( const std::vector<double>& array, std::vector<double>& statVector ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const char* methodName  = "ForceValidationResult::_findStatsForDouble";
+
+// ---------------------------------------------------------------------------------------
+
+    statVector.resize( STAT_END );
+ 
+    double avgValue   =  0.0;
+    double stdValue   =  0.0;
+    double minValue   =  1.0e+30;
+    double maxValue   = -1.0e+30;
+    int minValueIndex = 0;
+    int maxValueIndex = 0;
+ 
+    for( unsigned int ii = 0; ii < array.size(); ii++ ){
+ 
+       double norm  =  array[ii];
+ 
+       avgValue    += norm;
+       stdValue    += norm*norm;
+ 
+       if( norm > maxValue ){
+          maxValue       = norm;
+          maxValueIndex  = ii;
+       }
+       if( norm < minValue ){
+          minValue       = norm;
+          minValueIndex  = ii;
+       }
+    }
+ 
+    double count  = static_cast<double>(array.size());
+    double iCount = count > 0.0 ? 1.0/count : 0.0;
+ 
+    statVector[STAT_AVG] = avgValue*iCount;
+    statVector[STAT_STD] = stdValue - avgValue*avgValue*count;
+    if( count > 1.0 ){
+       statVector[STAT_STD] = std::sqrt( stdValue/( count - 1.0 ) );
+    }
+    statVector[STAT_MIN] = minValue;
+    statVector[STAT_ID1] = static_cast<double>(minValueIndex);
+    statVector[STAT_MAX] = maxValue;
+    statVector[STAT_ID2] = static_cast<double>(maxValueIndex);
+    statVector[STAT_CNT] = count;
+ 
+    return;
+}
+
+void ForceValidationResult::_calculateNormOfForceVector( int forceIndex ){
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::_calculateNormOfForceVector";
+
+// ---------------------------------------------------------------------------------------
+
+    for( unsigned int ii = 0; ii < _forces[forceIndex].size(); ii++ ){
+        Vec3 f1 = _forces[forceIndex][ii];
+        _norms[forceIndex].push_back( std::sqrt( (f1[0]*f1[0]) + (f1[1]*f1[1]) + (f1[2]*f1[2]) ) );
+        if( ValidateOpenMM::isNan( f1[0] ) || ValidateOpenMM::isNan( f1[1] ) || ValidateOpenMM::isNan( f1[2] ) ){
+            _nansDetected++;
+        }
+    }
+
+    return;
+}
+
+double ForceValidationResult::getPotentialEnergy( int energyIndex ) const {
+    if( energyIndex == 0 || energyIndex == 1 ){
+        return _potentialEnergies[energyIndex];
+    } else {
+       char buffer[1024];
+       (void) sprintf( buffer, "getPotentialEnergy: energyIndex=%d is inconsistent.", energyIndex );
+       throw OpenMMException( buffer );
+   }
+}
+
+std::vector<Vec3> ForceValidationResult::getForces( int forceIndex  ) const {
+    if( forceIndex == 0 || forceIndex == 1 ){
+        return _forces[forceIndex];
+    } else {
+       char buffer[1024];
+       (void) sprintf( buffer, "getForces: forceIndex=%d is inconsistent.", forceIndex );
+       throw OpenMMException( buffer );
+   }
+}
+
+std::vector<double> ForceValidationResult::getForceNorms( int forceIndex  ) const {
+    if( forceIndex == 0 || forceIndex == 1 ){
+        return _norms[forceIndex];
+    } else {
+       char buffer[1024];
+       (void) sprintf( buffer, "getForceNorms: forceIndex=%d is inconsistent.", forceIndex );
+       throw OpenMMException( buffer );
+   }
+}
+
+int ForceValidationResult::nansDetected( void ) const {
+    return _nansDetected;
+}
+
+void ForceValidationResult::registerInconsistentForceIndex( int index, int value ){
+    _inconsistentForceIndicies[index] = value;
+}
+
+void ForceValidationResult::clearInconsistentForceIndexList( void ){
+    _inconsistentForceIndicies.clear();
+}
+
+void ForceValidationResult::getInconsistentForceIndexList( std::vector<int>& inconsistentIndices ) const {
+    for( MapIntIntCI ii = _inconsistentForceIndicies.begin(); ii != _inconsistentForceIndicies.end(); ii++ ){
+        inconsistentIndices.push_back( (*ii).first );
+    }
+}
+
+int ForceValidationResult::getNumberOfInconsistentForceEntries( void ) const {
+    return  static_cast<int>(_inconsistentForceIndicies.size() );
+}
+
+double ForceValidationResult::getMaxDeltaForceNorm( int* maxIndex  ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::calculateNorm";
+
+// ---------------------------------------------------------------------------------------
+
+    double maxDelta                           = -1.0e+90;
+    int maxDeltaIndex                         = -1;
+
+    for( unsigned int ii = 0; ii < _forces[0].size(); ii++ ){
+ 
+        Vec3 f1                = _forces[0][ii];
+        double normF1          = _norms[0][ii];
+     
+        Vec3 f2                = _forces[1][ii];
+        double normF2          = _norms[1][ii];
+  
+        double delta           = std::sqrt( (f1[0]-f2[0])*(f1[0]-f2[0]) + (f1[1]-f2[1])*(f1[1]-f2[1]) + (f1[2]-f2[2])*(f1[2]-f2[2]) );
+        if( maxDelta < delta ){
+           maxDelta      = delta;
+           maxDeltaIndex = static_cast<int>(ii);
+        }
+    }
+
+    if( maxIndex ){
+        *maxIndex = maxDeltaIndex;
+    }
+         
+    return maxDelta;
+}
+
+double ForceValidationResult::getMaxRelativeDeltaForceNorm( int* maxIndex  ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::getMaxRelativeDeltaForceNorm";
+
+// ---------------------------------------------------------------------------------------
+
+    double maxRelativeDelta    = -1.0e+90;
+    int maxRelativeDeltaIndex  = -1;
+
+    for( unsigned int ii = 0; ii < _forces[0].size(); ii++ ){
+ 
+        Vec3 f1                = _forces[0][ii];
+        double normF1          = _norms[0][ii];
+
+        Vec3 f2                = _forces[1][ii];
+        double normF2          = _norms[1][ii];
+
+        double delta           = std::sqrt( (f1[0]-f2[0])*(f1[0]-f2[0]) + (f1[1]-f2[1])*(f1[1]-f2[1]) + (f1[2]-f2[2])*(f1[2]-f2[2]) );
+        double forceSum        = 0.5*(normF1 + normF2);
+        if( forceSum > 0.0 ){
+            delta /= forceSum;
+        }
+  
+        if( maxRelativeDelta < delta ){
+           maxRelativeDelta      = delta;
+           maxRelativeDeltaIndex = static_cast<int>(ii);
+        }
+    }
+
+    if( maxIndex ){
+        *maxIndex = maxRelativeDeltaIndex;
+    }
+         
+    return maxRelativeDelta;
+}
+
+double ForceValidationResult::getMaxDotProduct( int* maxIndex  ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::getMaxDotProduct";
+
+// ---------------------------------------------------------------------------------------
+
+    double maxDotProduct       = -1.0e+90;
+    int maxDotProductIndex     = -1;
+
+    for( unsigned int ii = 0; ii < _forces[0].size(); ii++ ){
+ 
+        Vec3 f1                = _forces[0][ii];
+        double normF1          = _norms[0][ii];
+
+        Vec3 f2                = _forces[1][ii];
+        double normF2          = _norms[1][ii];
+
+        double dotProduct      = f1[0]*f2[0] + f1[1]*f2[1] + f1[2]*f2[2];
+
+        if(  (normF1*normF2) > 0.0 ){
+            dotProduct     /= (normF1*normF2);
+            dotProduct      = 1.0 - dotProduct;
+        } else {
+            dotProduct      = 0.0;
+        }
+ 
+        if( maxDotProduct < dotProduct ){
+           maxDotProduct      = dotProduct;
+           maxDotProductIndex = static_cast<int>(ii);
+        }
+    }
+
+    if( maxIndex ){
+        *maxIndex = maxDotProductIndex;
+    }
+         
+    return maxDotProduct;
+}
+
+std::string ForceValidationResult::getForceName( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::getForceName";
+
+// ---------------------------------------------------------------------------------------
+
+    std::stringstream forceName;
+    for( unsigned int ii = 0; ii < _forceNames.size(); ii++ ){
+       forceName << _forceNames[ii];
+       if( ii < (_forceNames.size()-1) )forceName << "::";
+    }
+    return forceName.str();
+}
+
+std::string ForceValidationResult::getPlatformName( int index ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ForceValidationResult::getPlatformName";
+
+// ---------------------------------------------------------------------------------------
+
+    if( index == 0 || index == 1 ){
+        return _platforms[index];
+    } else {
+       char buffer[1024];
+       (void) sprintf( buffer, "getPlatformName: index=%d is inconsistent.", index );
+       throw OpenMMException( buffer );
+   }
+}
+
+void ForceValidationResult::compareForces( double tolerance ){
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = " ForceValidationResult::compareForces";
+
+// ---------------------------------------------------------------------------------------
+
+    std::vector<Vec3> forces1          = getForces( 0 );
+    std::vector<double> forceNorms1    = getForceNorms( 0 );
+
+    std::vector<Vec3> forces2          = getForces( 1 );
+    std::vector<double> forceNorms2    = getForceNorms( 1 );
+
+    clearInconsistentForceIndexList( );
+    for( unsigned int jj = 0; jj < forces1.size(); jj++ ){
+        if( ValidateOpenMM::isNan( forceNorms1[jj] ) ||  ValidateOpenMM::isNan( forceNorms2[jj] ) ){
+            registerInconsistentForceIndex( jj );
+        } else {
+            double delta   = std::sqrt( (forces1[jj][0]-forces2[jj][0])*(forces1[jj][0]-forces2[jj][0]) +
+                                        (forces1[jj][1]-forces2[jj][1])*(forces1[jj][1]-forces2[jj][1]) +
+                                        (forces1[jj][2]-forces2[jj][2])*(forces1[jj][2]-forces2[jj][2]) );
+
+            double sum     = 0.5*( fabs( forceNorms1[jj] ) + fabs( forceNorms2[jj] ) );
+            double diff    = delta > 0.0 ? delta/sum : 0.0;
+            if( diff > tolerance && sum > tolerance ){
+                registerInconsistentForceIndex( jj );
+            }
+        }
+    }
+
+}
+
+void ForceValidationResult::compareForceNorms( double tolerance ){
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "ForceValidationResult::compareForceNorms";
+
+// ---------------------------------------------------------------------------------------
+
+    std::vector<double> forceNorms0    = getForceNorms( 0 );
+    std::vector<double> forceNorms1    = getForceNorms( 1 );
+
+    clearInconsistentForceIndexList( );
+    for( unsigned int jj = 0; jj < forceNorms0.size(); jj++ ){
+        if( ValidateOpenMM::isNan( forceNorms0[jj] ) || ValidateOpenMM::isNan( forceNorms1[jj] ) ){
+            registerInconsistentForceIndex( jj );
+        } else {
+            double sum  = 0.5*( fabs( forceNorms0[jj] ) + fabs( forceNorms1[jj] ) );
+            double diff = sum > 0.0 ? fabs( forceNorms0[jj] - forceNorms1[jj] )/sum : 0.0;
+            if( diff > tolerance && sum > tolerance ){
+                registerInconsistentForceIndex( jj );
+            }
+        }
+    }
+}
+
+ValidateOpenMMForces::ValidateOpenMMForces( void ) {
+    _initialize();
+}
+
+void ValidateOpenMMForces::_initialize( void ){
+
+    _forceTolerance                               = 1.0e-02;
+    _maxErrorsToPrint                             = 25;
+
+    // force tolerances by name
+
+    _forceTolerances[HARMONIC_BOND_FORCE]         = _forceTolerance;
+    _forceTolerances[HARMONIC_ANGLE_FORCE]        = _forceTolerance;
+    //_forceTolerances[PERIODIC_TORSION_FORCE]      = _forceTolerance;
+    _forceTolerances[PERIODIC_TORSION_FORCE]      = 0.3;
+    //_forceTolerances[RB_TORSION_FORCE]            = _forceTolerance;
+    _forceTolerances[RB_TORSION_FORCE]            = 0.3;
+    _forceTolerances[NB_FORCE]                    = _forceTolerance;
+    _forceTolerances[NB_SOFTCORE_FORCE]           = _forceTolerance;
+    _forceTolerances[GBSA_OBC_FORCE]              = _forceTolerance;
+    _forceTolerances[GBSA_OBC_SOFTCORE_FORCE]     = _forceTolerance;
+    _forceTolerances[GBVI_FORCE]                  = _forceTolerance;
+    _forceTolerances[GBVI_SOFTCORE_FORCE]         = _forceTolerance;
+
+    // forces to be excluded from validation
+
+    _forcesToBeExcluded[CM_MOTION_REMOVER]        = 1;
+    _forcesToBeExcluded[ANDERSEN_THERMOSTAT]      = 1;
+    _forcesToBeExcluded[CUSTOM_BOND_FORCE]        = 1;
+    _forcesToBeExcluded[CUSTOM_EXTERNAL_FORCE]    = 1;
+    _forcesToBeExcluded[CUSTOM_NONBONDED_FORCE]   = 1;
+}
+
+ValidateOpenMMForces::~ValidateOpenMMForces( ){
+
+    for( unsigned int ii = 0; ii < _forceValidationResults.size(); ii++ ){
+         delete _forceValidationResults[ii];
+    }
+    _forceValidationResults.resize( 0 );
+    
+}
+
+double ValidateOpenMMForces::getForceTolerance( void ) const {
+    return _forceTolerance;
+}
+
+void ValidateOpenMMForces::setForceTolerance( double forceTolerance ){
+    _forceTolerance = forceTolerance;
+}
+
+int ValidateOpenMMForces::compareWithReferencePlatform( Context& context, std::string* summaryString ){
+
+// ---------------------------------------------------------------------------------------
+
+    static const std::string methodName      = "ValidateOpenMMForces::compareWithReferencePlatform";
+
+// ---------------------------------------------------------------------------------------
+
+#ifdef INCLUDE_FREE_ENERGY_PLUGIN
+    ReferenceFreeEnergyPlatform referencePlatform;
+#else
+    ReferencePlatform referencePlatform;
+#endif
+
+    compareOpenMMForces( context, referencePlatform, _forceValidationResults );
+
+    checkForInconsistentForceEntries( _forceValidationResults );
+
+    if( summaryString ){
+        *summaryString = getSummary( _forceValidationResults );
+    }
+
+    return getTotalNumberOfInconsistentForceEntries( _forceValidationResults );
+}
+
+void ValidateOpenMMForces::compareOpenMMForces( Context& context, Platform& platform,
+                                                std::vector<ForceValidationResult*>& forceValidationResults ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "ValidateOpenMMForces::compareOpenMMForces";
+
+// ---------------------------------------------------------------------------------------
+
+    // loop over forces
+ 
+    const System& system         = context.getSystem();
+    std::vector<int> compareAllForces;
+    for( unsigned int ii = 0; ii < system.getNumForces(); ii++ ){
+        std::vector<int> compareForces;
+        std::string forceName  = getForceName( system.getForce( ii ) );
+        if( isExcludedForce( forceName ) == 0 ){
+            compareForces.push_back( ii );
+            compareAllForces.push_back( ii );
+            ForceValidationResult* forceValidationResult = compareForce(context, compareForces, platform, context.getPlatform() );
+            forceValidationResults.push_back( forceValidationResult );
+        } else if( getLog() ){
+           (void) fprintf( getLog(), "Force %s is not being validated.\n", forceName.c_str() );
+        }
+    }
+    ForceValidationResult* forceValidationResult = compareForce(context, compareAllForces, platform, context.getPlatform() );
+    forceValidationResults.push_back( forceValidationResult );
+}
+
+ForceValidationResult* ValidateOpenMMForces::compareForce(Context& context, std::vector<int>& compareForces,
+                                                          Platform& platform1, Platform& platform2 ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    //static const std::string methodName      = "ValidateOpenMMForces::compareForce";
+
+// ---------------------------------------------------------------------------------------
+
+    // note if platforms are identical
+
+    if( getLog() && platform1.getName().compare( platform2.getName() ) == 0 ){
+        (void) fprintf( getLog(), "Note: Platforms to compares %s are identical.\n", platform1.getName().c_str() );
+        (void) fflush( getLog() );
+    }
+
+    const System& system         = context.getSystem();
+
+    // collect systemForceNameMap[forceName] = index in system
+    //         systemForceNameIndex[index]   = force name
+
+    StringIntMap systemForceNameMap;
+    StringVector systemForceNameIndex;
+    systemForceNameIndex.resize( system.getNumForces() );
+    for( int ii = 0; ii < system.getNumForces(); ii++ ){
+        std::string forceName         = getForceName( system.getForce( ii ) );
+        if( forceName.compare( "NA" ) == 0 ){
+            (void) fprintf( stderr, "Force at index=%d not found -- aborting!\n", ii );
+            exit(-1);
+        }
+        systemForceNameMap[forceName] = ii;
+        systemForceNameIndex[ii]      = forceName;
+    }
+
+    // diagnostics
+
+    if( 0 && getLog() ){
+        for( StringIntMapI ii = systemForceNameMap.begin(); ii != systemForceNameMap.end(); ii++ ){
+            int index = (*ii).second;
+            (void) fprintf( getLog(), "  System force map %s index=%d reverse map=%s\n", (*ii).first.c_str(), index, systemForceNameIndex[index].c_str() );
+        }
+        for( unsigned int ii = 0; ii < compareForces.size(); ii++ ){
+           (void) fprintf( getLog(), "   ValidateOpenMMForces %u %s\n", ii, systemForceNameIndex[compareForces[ii]].c_str() );
+        }
+        (void) fflush( getLog() );
+    }
+
+    // get system copy and add forces to system
+
+    System* validationSystem     = copySystemExcludingForces( system );
+    StringUIntMap forceNamesMap;
+    for( unsigned int ii = 0; ii < compareForces.size(); ii++ ){
+        const Force& forceToCopy = system.getForce( compareForces[ii] );
+        Force* force             = copyForce( forceToCopy );
+        validationSystem->addForce( force );
+        forceNamesMap[systemForceNameIndex[compareForces[ii]]] = ii;
+    }
+
+    // include any missing dependencies (e.g, OBC force requires NB force for Cuda platform)
+
+    for( StringUIntMapI ii = forceNamesMap.begin(); ii != forceNamesMap.end(); ii++ ){
+       std::string forceName = (*ii).first;
+       StringVector dependencyVector;
+       getForceDependencies( forceName, dependencyVector ); 
+       for( unsigned int jj = 0; jj < dependencyVector.size(); jj++ ){
+           std::string dependentForceName = dependencyVector[jj];
+           StringUIntMapCI dependent      = forceNamesMap.find( dependentForceName );
+           if( dependent == forceNamesMap.end() ){
+              forceNamesMap[dependentForceName] = 1;
+              int forceIndex                    = systemForceNameMap[dependentForceName];
+              const Force& forceToCopy          = system.getForce( forceIndex );
+              validationSystem->addForce( copyForce( forceToCopy ) ); 
+           }
+       }
+    }
+
+    // create contexts
+
+    VerletIntegrator verletIntegrator( 0.001 );
+    Context* validationContext1  = new Context( *validationSystem, verletIntegrator, platform1);
+    Context* validationContext2  = new Context( *validationSystem, verletIntegrator, platform2);
+
+    // set positions
+
+    synchContexts( context, *validationContext1 );
+    synchContexts( context, *validationContext2 );
+
+    // diagnostics
+
+    if( 0 && getLog() ){
+        std::stringstream forceNames;
+        (void) fprintf( getLog(), "    Validating system forces=%d\n", validationSystem->getNumForces() );
+        for( int ii = 0; ii < validationSystem->getNumForces(); ii++ ){
+            std::string forceName         = getForceName( validationSystem->getForce( ii ) );
+            forceNames << forceName;
+            if( ii < (validationSystem->getNumForces()-1) ){
+                forceNames << "_";
+            } else {
+                forceNames << "Parameters.txt";
+            }
+            (void) fprintf( getLog(), "       force %d %s\n", ii, forceName.c_str() );
+        }
+        writeParameterFile( *validationContext1, forceNames.str() ); 
+        (void) fflush( getLog() );
+    }
+
+    // calculate forces & build return result
+
+    ForceValidationResult* forceValidationResult = new ForceValidationResult( *validationContext1, *validationContext2, forceNamesMap );
+
+    delete validationContext1;
+    delete validationContext2;
+    delete validationSystem;
+        
+    return forceValidationResult;
+}
+
+
+void ValidateOpenMMForces::checkForInconsistentForceEntries( std::vector<ForceValidationResult*>& forceValidationResults  ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "ValidateOpenMMForces::checkForInconsistentForceEntries";
+
+// ---------------------------------------------------------------------------------------
+
+    for( unsigned int ii = 0; ii < forceValidationResults.size(); ii++ ){
+        std::string forceName              = forceValidationResults[ii]->getForceName();
+        double tolerance                   = getForceTolerance( forceName );
+        forceValidationResults[ii]->compareForces( tolerance ); 
+    }
+
+}
+
+int ValidateOpenMMForces::getTotalNumberOfInconsistentForceEntries( std::vector<ForceValidationResult*>& forceValidationResults  ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "ValidateOpenMMForces::getTotalNumberOfInconsistentForceEntries";
+
+// ---------------------------------------------------------------------------------------
+
+    int inconsistentEntries = 0;
+    for( unsigned int ii = 0; ii < forceValidationResults.size(); ii++ ){
+        inconsistentEntries  += forceValidationResults[ii]->getNumberOfInconsistentForceEntries();
+    }
+
+    return inconsistentEntries;
+}
+
+int ValidateOpenMMForces::isExcludedForce( std::string forceName ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "ValidateOpenMMForces::isExcludedForce";
+
+// ---------------------------------------------------------------------------------------
+
+    StringIntMapCI isPresent = _forcesToBeExcluded.find( forceName );
+    return isPresent != _forcesToBeExcluded.end() ? 1 : 0;
+}
+
+
+/* 
+ * Get force tolerance for specified force
+ *
+ * @param forceName   name of force
+ *
+ * @return force tolerance
+ *
+ * */
+
+double ValidateOpenMMForces::getForceTolerance( const std::string& forceName ) const {
+    
+    StringDoubleMapCI forceIsPresent = _forceTolerances.find( forceName );
+    if( forceIsPresent != _forceTolerances.end() ){
+        return (*forceIsPresent).second;
+    }
+    return _forceTolerance;
+}
+
+int ValidateOpenMMForces::getMaxErrorsToPrint( void ) const {
+    return _maxErrorsToPrint;
+}
+
+void ValidateOpenMMForces::setMaxErrorsToPrint( int maxErrorsToPrint ){
+    _maxErrorsToPrint = maxErrorsToPrint;
+}
+
+/* 
+ * Format line
+ *
+ * @param tab         tab
+ * @param description description
+ * @param value       value
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string ValidateOpenMMForces::_getLine( const std::string& tab, 
+                                            const std::string& description,
+                                            const std::string& value ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    static const std::string methodName      = "GromacsOpenMMTest::_getLine";
+    static const unsigned int MAX_LINE_CHARS = 256; 
+    char line[MAX_LINE_CHARS];
+ 
+ // ---------------------------------------------------------------------------------------
+ 
+    std::stringstream message;
+    memset( line, ' ', MAX_LINE_CHARS );   
+#ifdef WIN32
+    (void) sprintf_s( line, MAX_LINE_CHARS, "%s %-40s %s", tab.c_str(), description.c_str(), value.c_str() );
+#else
+    (void) sprintf( line, "%s %-40s %s", tab.c_str(), description.c_str(), value.c_str() );
+#endif
+    message << std::string( line ) << std::endl;
+ 
+    return message.str();
+
+}
+
+std::string ValidateOpenMMForces::getSummary( std::vector<ForceValidationResult*>& forceValidationResults ) const {
+
+// ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName      = "ValidateOpenMMForces::getSummary";
+   static const unsigned int MAX_LINE_CHARS = 256; 
+   char value[MAX_LINE_CHARS];
+
+// ---------------------------------------------------------------------------------------
+
+#ifdef WIN32
+#define LOCAL_SPRINTF0(a,b) sprintf_s( (a), MAX_LINE_CHARS, (b) );   
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#define LOCAL_SPRINTF4(a,b,c,d) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d) );   
+#define LOCAL_SPRINTF5(a,b,c,d,e) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d), (e) );   
+#define LOCAL_SPRINTF6(a,b,c,d,e,f) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d), (e), (f) );   
+#define LOCAL_SPRINTF12(a,b,c,d,e,f,g,h,i,j,k,l) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d), (e), (f), (g), (h), (i), (j), (k), (l) );   
+#else
+#define LOCAL_SPRINTF0(a,b) sprintf( (a), (b) );   
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#define LOCAL_SPRINTF4(a,b,c,d) sprintf( (a), (b), (c), (d) );   
+#define LOCAL_SPRINTF5(a,b,c,d,e) sprintf( (a), (b), (c), (d), (e) );   
+#define LOCAL_SPRINTF6(a,b,c,d,e,f) sprintf( (a), (b), (c), (d), (e), (f) );   
+#define LOCAL_SPRINTF12(a,b,c,d,e,f,g,h,i,j,k,l) sprintf( (a), (b), (c), (d), (e), (f), (g), (h), (i), (j), (k), (l) );   
+#endif
+
+    int maxMissesToPrint = getMaxErrorsToPrint();
+    std::stringstream summary;
+    std::string tab      = "   ";
+    int useForces        = 1;
+    for( unsigned int ii = 0; ii < forceValidationResults.size(); ii++ ){
+
+        // platforms
+
+        if( ii == 0 ){
+            std::string platform1 = forceValidationResults[ii]->getPlatformName( 0 );
+            std::string platform2 = forceValidationResults[ii]->getPlatformName( 1 );
+            (void) LOCAL_SPRINTF4( value, "%10s  %10s\n", platform1.c_str(), platform2.c_str());
+            summary << _getLine( tab, "Platforms", value );
+        }
+
+        // force name
+        std::string forceName              = forceValidationResults[ii]->getForceName();
+        double tolerance                   = getForceTolerance( forceName );
+        summary << _getLine( tab, "Force", forceName );
+        (void) LOCAL_SPRINTF( value, "%.3e ", tolerance );
+        summary << _getLine( tab, "Tolerance", value );
+        int index;
+
+        // delta
+
+        double delta = forceValidationResults[ii]->getMaxDeltaForceNorm( &index );
+        (void) LOCAL_SPRINTF4( value, "%.3e at index %6d", delta, index );
+        summary << _getLine( tab, "Max Delta", value );
+
+        // relative delta
+
+        delta        = forceValidationResults[ii]->getMaxRelativeDeltaForceNorm( &index );
+        (void) LOCAL_SPRINTF4( value, "%.3e at index %6d", delta, index );
+        summary << _getLine( tab, "Max Relative Delta", value );
+
+        // PE delta
+
+        double pe1   = forceValidationResults[ii]->getPotentialEnergy( 0 );
+        double pe2   = forceValidationResults[ii]->getPotentialEnergy( 1 );
+        double sum   = 0.5*(fabs( pe1 ) + fabs( pe2 ) );
+        delta        = sum > 0.0 ? fabs( pe1 - pe2 )/sum : 0.0;
+        (void) LOCAL_SPRINTF5( value, "%10.4e PE[%14.6e %14.6e]", delta, pe1, pe2 );
+        summary << _getLine( tab, "Potential energies relative delta", value );
+
+        // misses
+
+        std::vector<int> inconsistentIndices;
+        forceValidationResults[ii]->getInconsistentForceIndexList( inconsistentIndices );
+        if( inconsistentIndices.size() > 0 ){ 
+    
+            std::vector<Vec3> forces1          = forceValidationResults[ii]->getForces( 0 );
+            std::vector<Vec3> forces2          = forceValidationResults[ii]->getForces( 1 );
+    
+            std::vector<double> forceNorms1    = forceValidationResults[ii]->getForceNorms( 0 );
+            std::vector<double> forceNorms2    = forceValidationResults[ii]->getForceNorms( 1 );
+            for( unsigned int kk = 0; kk < inconsistentIndices.size() && kk < maxMissesToPrint; kk++ ){
+                int jj = inconsistentIndices[kk];
+                if( isinf( forceNorms1[jj] ) || isnan( forceNorms1[jj] ) ||
+                    isinf( forceNorms2[jj] ) || isnan( forceNorms2[jj] ) ){
+                     (void) LOCAL_SPRINTF5( value, "         nan at index %6d  norms: [%12.5e  %12.5e]",  jj, forceNorms1[jj], forceNorms2[jj] );
+                     summary << _getLine( tab, "Error", value );
+                } else {
+    
+                    double diff, sum;
+                    if( useForces ){
+                        double delta   = std::sqrt( (forces1[jj][0]-forces2[jj][0])*(forces1[jj][0]-forces2[jj][0]) +
+                                                    (forces1[jj][1]-forces2[jj][1])*(forces1[jj][1]-forces2[jj][1]) +
+                                                    (forces1[jj][2]-forces2[jj][2])*(forces1[jj][2]-forces2[jj][2]) );
+    
+                        sum            = 0.5*( fabs( forceNorms1[jj] ) + fabs( forceNorms2[jj] ) );
+                        diff           = delta > 0.0 ? delta/sum : 0.0;
+                    } else {
+                        sum            = 0.5*( fabs( forceNorms1[jj] ) + fabs( forceNorms2[jj] ) );
+                        diff           = sum > 0.0 ? fabs( forceNorms1[jj] - forceNorms2[jj] )/sum : 0.0;
+                    }
+        
+                    (void) LOCAL_SPRINTF12( value, "%11.5e  at index %6d  norms: [%11.5e  %11.5e] forces: [%12.5e  %12.5e %12.5e] [%12.5e  %12.5e %12.5e]",
+                                            diff, jj, forceNorms1[jj], forceNorms2[jj],
+                                            forces1[jj][0], forces1[jj][1], forces1[jj][2],
+                                            forces2[jj][0], forces2[jj][1], forces2[jj][2] );
+                    summary << _getLine( tab, "Error", value );
+                    if( kk == (maxMissesToPrint-1) ){
+                        value[0] = '\0';
+                        summary << _getLine( tab, "No more errors will be printed: call ValidateOpenMMForces::setMaxErrorsToPrint() to modifiy this limit.", value );
+                    }
+                }
+            }
+        }
+
+        (void) LOCAL_SPRINTF0( value, "\n" );
+        summary << _getLine( tab, " ", value );
+        if( inconsistentIndices.size() > 0 ){
+            (void) LOCAL_SPRINTF( value, "%u\n", inconsistentIndices.size() );
+            summary << _getLine( tab, "Total errors", value );
+        }
+        if( forceValidationResults[ii]->nansDetected() ){
+            value[0] = '\0';
+            summary << _getLine( tab, "Nans detected.", value );
+        }
+    }
+
+    return summary.str();
+}