Skip to content

GitLab

"platforms/reference/vscode:/vscode.git/clone" did not exist on "21847f4a12a3a2e38d03933f737e68a67559ad2c"
Commit 4fdf49d5 authored by peastman's avatar peastman
Browse files

Merge pull request #347 from peastman/master 

Removed some obsolete test code
parents b915aca5 b16bec99
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 3606 deletions
platforms/cuda/tests/CMakeLists.txt
View file @ 4fdf49d5
...@@ -31,35 +31,6 @@ FOREACH(TEST_PROG ${TEST_PROGS}) ...@@ -31,35 +31,6 @@ FOREACH(TEST_PROG ${TEST_PROGS})
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS}") SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS}")
ENDIF (APPLE) ENDIF (APPLE)
IF( ${TEST_ROOT} STREQUAL "TestCUDAGBSAOBCForce2" )
SET(DEFINE_STRING "-DTEST_PLATFORM=1 ")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE ")
ENDIF( INCLUDE_SERIALIZATION )
# nonbond
SET(NONBOND_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=0")
SET(NONBOND_TEST "TestCUDANonbondedForce2")
ADD_EXECUTABLE(${NONBOND_TEST} ${TEST_PROG})
SET_TARGET_PROPERTIES(${NONBOND_TEST} PROPERTIES COMPILE_FLAGS ${NONBOND_DEFINE_STRING} )
ADD_TEST(${NONBOND_TEST} ${EXECUTABLE_OUTPUT_PATH}/${NONBOND_TEST})
# OBC
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES COMPILE_FLAGS ${DEFINE_STRING} )
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
ENDIF( ${TEST_ROOT} STREQUAL "TestCUDAGBSAOBCForce2" )
ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single) ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
IF (OPENMM_BUILD_CUDA_DOUBLE_PRECISION_TESTS) IF (OPENMM_BUILD_CUDA_DOUBLE_PRECISION_TESTS)
ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed) ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed)
......
platforms/opencl/tests/CMakeLists.txt
View file @ 4fdf49d5
...@@ -27,35 +27,6 @@ FOREACH(TEST_PROG ${TEST_PROGS}) ...@@ -27,35 +27,6 @@ FOREACH(TEST_PROG ${TEST_PROGS})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET}) TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS}") SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS}")
IF( ${TEST_ROOT} STREQUAL "TestOpenCLGBSAOBCForce2" )
SET(DEFINE_STRING "-DTEST_PLATFORM=1 ")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE ")
ENDIF( INCLUDE_SERIALIZATION )
# nonbond
SET(NONBOND_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=0")
SET(NONBOND_TEST "TestOpenCLNonbondedForce2")
ADD_EXECUTABLE(${NONBOND_TEST} ${TEST_PROG})
SET_TARGET_PROPERTIES(${NONBOND_TEST} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} ${NONBOND_DEFINE_STRING}" )
ADD_TEST(${NONBOND_TEST} ${EXECUTABLE_OUTPUT_PATH}/${NONBOND_TEST})
# OBC
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} ${DEFINE_STRING}" )
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
ENDIF( ${TEST_ROOT} STREQUAL "TestOpenCLGBSAOBCForce2" )
ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single) ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
IF (OPENMM_BUILD_OPENCL_DOUBLE_PRECISION_TESTS) IF (OPENMM_BUILD_OPENCL_DOUBLE_PRECISION_TESTS)
ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed) ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed)
......
platforms/opencl/tests/TestOpenCLGBSAOBCForce2.cpp deleted 100644 → 0
View file @ b915aca5
/* -------------------------------------------------------------------------- *
* 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: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementations of GBVIForce, GBSAOBCForce and the softcore versions of
* these forces.
*/
#define TEST_NONBONDED 0
#define TEST_OBC 1
#define TEST_GBVI 2
#define TEST_CUDA_PLATFORM 0
#define TEST_OPENCL_PLATFORM 1
#include "openmm/GBVIForce.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/NonbondedForce.h"
#ifdef USE_SOFTCORE
#include "openmm/GBVISoftcoreForce.h"
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#endif
/**
* Utility methods shared across unit tests
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/System.h"
#include "SimTKOpenMMRealType.h"
#include "OpenMM.h"
#if TEST_PLATFORM == TEST_OPENCL_PLATFORM
#include "ReferencePlatform.h"
#include "OpenCLPlatform.h"
#endif
#if TEST_PLATFORM == TEST_CUDA_PLATFORM
#include "ReferencePlatform.h"
#include "CudaPlatform.h"
#endif
#ifdef USE_SOFTCORE
#include "OpenMMFreeEnergy.h"
#include "openmm/freeEnergyKernels.h"
#endif
#include "sfmt/SFMT.h"
#include "openmm/VerletIntegrator.h"
#ifdef OPENMM_SERIALIZE
#include "openmm/serialization/SerializationProxy.h"
#include "openmm/serialization/SerializationNode.h"
#include "openmm/serialization/XmlSerializer.h"
#endif
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <typeinfo>
#include <iomanip>
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
static const int NoCutoff_OpenMMTest = 0;
static const int CutoffNonPeriodic_OpenMMTest = 1;
static const int CutoffPeriodic_OpenMMTest = 2;
static const int Ewald_OpenMMTest = 3;
static const int PME_OpenMMTest = 4;
static const int ChargeIndex_OpenMMTest = 0;
static const int SigmaIndex_OpenMMTest = 1;
static const int EpsIndex_OpenMMTest = 2;
static const int GammaIndex_OpenMMTest = 3;
static const int LambdaIndex_OpenMMTest = 4;
static const int Reference_OpenMMTest = 0;
static const int Cuda_OpenMMTest = 1;
static const int OpenCL_OpenMMTest = 2;
class BondInfo_OpenMMTest {
public:
BondInfo_OpenMMTest( int particle1, int particle2, double distance );
int _particle1;
int _particle2;
double _distance;
};
BondInfo_OpenMMTest::BondInfo_OpenMMTest( int particle1, int particle2, double distance ){
_particle1 = particle1;
_particle2 = particle2;
_distance = distance;
}
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
typedef std::vector<StringVector> StringVectorVector;
typedef std::vector<int> IntVector;
typedef IntVector::iterator IntVectorI;
typedef IntVector::const_iterator IntVectorCI;
typedef std::vector<IntVector> IntVectorVector;
typedef std::vector<double> DoubleVector;
typedef DoubleVector::iterator DoubleVectorI;
typedef DoubleVector::const_iterator DoubleVectorCI;
typedef std::vector<DoubleVector> DoubleVectorVector;
// the following are used in parsing parameter file
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
typedef std::vector<StringVector> VectorStringVector;
typedef VectorStringVector::iterator VectorStringVectorI;
typedef VectorStringVector::const_iterator VectorStringVectorCI;
typedef std::vector<std::vector<double> > VectorOfDoubleVectors;
typedef VectorOfDoubleVectors::iterator VectorOfDoubleVectorsI;
typedef VectorOfDoubleVectors::const_iterator VectorOfDoubleVectorsCI;
typedef std::map< int, int> MapIntInt;
typedef MapIntInt::iterator MapIntIntI;
typedef MapIntInt::const_iterator MapIntIntCI;
typedef std::map< double, int> MapDoubleToInt;
typedef MapDoubleToInt::iterator MapDoubleToIntI;
typedef MapDoubleToInt::const_iterator MapDoubleToIntCI;
typedef std::map< std::string, VectorOfDoubleVectors > MapStringVectorOfDoubleVectors;
typedef MapStringVectorOfDoubleVectors::iterator MapStringVectorOfDoubleVectorsI;
typedef MapStringVectorOfDoubleVectors::const_iterator MapStringVectorOfDoubleVectorsCI;
typedef std::map< std::string, StringVector > MapStringStringVector;
typedef MapStringStringVector::iterator MapStringStringVectorI;
typedef MapStringStringVector::const_iterator MapStringStringVectorCI;
typedef std::map< std::string, std::string > MapStringString;
typedef MapStringString::iterator MapStringStringI;
typedef MapStringString::const_iterator MapStringStringCI;
typedef std::map< std::string, std::string > MapStringToInt;
typedef MapStringToInt::iterator MapStringToIntI;
typedef MapStringToInt::const_iterator MapStringToIntCI;
typedef std::vector< std::map< std::string, std::string > > VectorMapStringString;
typedef VectorMapStringString::iterator VectorMapStringStringI;
typedef VectorMapStringString::const_iterator VectorMapStringStringCI;
typedef std::map< std::string, int > MapStringInt;
typedef MapStringInt::iterator MapStringIntI;
typedef MapStringInt::const_iterator MapStringIntCI;
typedef std::map< std::string, std::vector<OpenMM::Vec3> > MapStringVec3;
typedef MapStringVec3::iterator MapStringVec3I;
typedef MapStringVec3::const_iterator MapStringVec3CI;
typedef std::map< std::string, double > MapStringToDouble;
typedef MapStringToDouble::iterator MapStringToDoubleI;
typedef MapStringToDouble::const_iterator MapStringToDoubleCI;
typedef std::vector< MapStringToDouble > VectorOfMapStringToDouble;
typedef std::map< std::string, DoubleVector> MapStringToDoubleVector;
typedef MapStringToDoubleVector::iterator MapStringToDoubleVectorI;
typedef MapStringToDoubleVector::const_iterator MapStringToDoubleVectorCI;
typedef std::map< std::string, DoubleVector > MapStringToDoubleVector;
typedef std::map< std::string, Force*> MapStringForce;
typedef MapStringForce::iterator MapStringForceI;
typedef MapStringForce::const_iterator MapStringForceCI;
typedef std::map< int, IntVector> MapIntIntVector;
typedef MapIntIntVector::const_iterator MapIntIntVectorCI;
typedef std::pair<int, int> IntIntPair;
typedef std::vector<IntIntPair> IntIntPairVector;
typedef IntIntPairVector::iterator IntIntPairVectorI;
typedef IntIntPairVector::const_iterator IntIntPairVectorCI;
typedef std::pair<int, double> IntDoublePair;
typedef std::vector<IntDoublePair> IntDoublePairVector;
typedef IntDoublePairVector::iterator IntDoublePairVectorI;
typedef IntDoublePairVector::const_iterator IntDoublePairVectorCI;
/**
* Predicate for sorting <int,double> pair
*
* @param d1 first IntDoublePair to compare
* @param d2 second IntDoublePair to compare
*
*/
bool TestIntDoublePair( const IntDoublePair& d1, const IntDoublePair& d2 ){
return d1.second < d2.second;
}
class PositionGenerator {
public:
enum GenerationMethod {
/**
* Random positions
*/
Random = 0,
/**
* On grid
*/
SimpleGrid = 1,
};
PositionGenerator( int numMolecules, int numParticlesPerMolecule, double boxSize );
~PositionGenerator( );
/**
* Set logging file reference
*
* @param log log
*
*/
void setLog( FILE* log );
/**
* Set positions
*
* @param method method placement
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int setPositions( GenerationMethod method, std::vector<Vec3>& positions ) const;
/**
* Set positions
*
* @param method method placement
* @param sfmt input random number generator
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int setPositions( GenerationMethod method, OpenMM_SFMT::SFMT& sfmt, std::vector<Vec3>& positions ) const;
/**
* Place particles on a grid
*
* @param origin origin
* @param boxDimensions box dimensions
* @param spacing spacing
* @param sfmt input random number generator
* @param bondDistance input bond distance
* @param array output vector of grid values
*
* @return -1 if particles will not fit on grid; 0 if they do
*/
int setParticlesOnGrid( const Vec3& origin, const Vec3& boxDimensions, const Vec3& spacing,
OpenMM_SFMT::SFMT& sfmt, double bondDistance, std::vector<Vec3>& array ) const;
/**
* Set bond distance
*
* @param bondDistance bond distance
*/
void setBondDistance( double bondDistance );
/**
* Get bond distance
*
* @return bond distance
*/
double getBondDistance( void ) const;
/**
* Get distance
*
* @param index1 index of first particle
* @param index2 index of second particle
* @param positions particle positions
*
* @return distance
*/
double getDistance( int index1, int index2, const std::vector<Vec3>& positions ) const;
/**
* Get distance assumming periodic boundary conditions
*
* @param index1 index of first particle
* @param index2 index of second particle
* @param positions particle positions
*
* @return distance
*/
double getPeriodicDistance( int index1, int index2, const std::vector<Vec3>& positions ) const;
/**
* Get settings
*
* @return info string
*/
std::string getSettings( void ) const;
/**
* Get enclosing box
*
* @param positions input vector of positions
* @param enclosingBox output vector of enclosing box dimensions
*
*/
void getEnclosingBox( const std::vector<Vec3>& positions, Vec3 enclosingBox[2] ) const;
/**
* Get sorted distances from particular position
*
* @param periodicBoundaryConditions if set, apply PBC
* @param positionIndex input position index
* @param positions input vector of positions
* @param sortVector on output sorted IntDoublePairVector
*
*/
void getSortedDistances( int periodicBoundaryConditions, int positionIndex, const std::vector<Vec3>& positions, IntDoublePairVector& sortVector ) const;
/**
* Show min/max distances between all positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
*
*/
void showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex );
/**
* Show min/max distances between a single position and all remaining positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
* @param positionIndexVector list of entries to show min/max distances from
*
*/
void showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex,
const IntVector& positionIndexVector );
/**
* Show distances between positions
*
* @param pairs particle indcies for which distance is to be reported
* @param positions input vector of positions
*
*/
void showDistances( const IntIntPairVector& pairs, const std::vector<Vec3>& positions ) const;
/**
* Show particles within a specified distance of a given particle
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param particleIndex particle to check
* @param distanceToCheckFor distance to check for
* @param tolerance distance tolerance
*
*/
void showParticlesWithinDistance( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, unsigned int particleIndex,
double distanceToCheckFor, double tolerance);
private:
int _numMolecules;
int _numParticlesPerMolecule;
int _numParticles;
int _seed;
double _boxSize;
double _bondDistance;
Vec3 _origin;
Vec3 _boxDimensions;
Vec3 _spacing;
FILE* _log;
};
PositionGenerator::PositionGenerator( int numMolecules, int numParticlesPerMolecule, double boxSize ) :
_numMolecules(numMolecules),
_seed(0),
_log(NULL),
_bondDistance(0.1),
_numParticlesPerMolecule(numParticlesPerMolecule),
_numParticles(numMolecules*numParticlesPerMolecule),
_boxSize(boxSize),
_boxDimensions(Vec3(boxSize,boxSize,boxSize)),
_origin(Vec3(0.0,0.0,0.0)) {
double particlesPerDimension = pow( static_cast<double>(_numParticles), (1.0/3.0) );
int particlesPerDimensionI = static_cast<int>(particlesPerDimension+0.999999);
double spacingPerDimension = _boxSize/particlesPerDimension;
_spacing[0] = spacingPerDimension;
_spacing[1] = spacingPerDimension;
_spacing[2] = spacingPerDimension;
}
PositionGenerator::~PositionGenerator( ){};
void PositionGenerator::setBondDistance( double bondDistance ){
_bondDistance = bondDistance;
}
void PositionGenerator::setLog( FILE* log ){
_log = log;
}
double PositionGenerator::getBondDistance( void ) const {
return _bondDistance;
}
double PositionGenerator::getDistance( int index1, int index2, const std::vector<Vec3>& positions ) const {
Vec3 delta = positions[index2] - positions[index1];
return sqrt( delta.dot( delta ) );
}
double PositionGenerator::getPeriodicDistance( int index1, int index2, const std::vector<Vec3>& positions ) const {
Vec3 delta = positions[index2] - positions[index1];
if( _boxSize > 0.0 ){
delta[0] -= floor(delta[0]/_boxSize+0.5f)*_boxSize;
delta[1] -= floor(delta[1]/_boxSize+0.5f)*_boxSize;
delta[2] -= floor(delta[2]/_boxSize+0.5f)*_boxSize;
}
return sqrt( delta.dot( delta ) );
}
/**
* Get positions
*
* @param method method placement
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int PositionGenerator::setPositions( GenerationMethod method, std::vector<Vec3>& positions ) const {
OpenMM_SFMT::SFMT sfmt;
init_gen_rand( _seed, sfmt);
return setPositions( method, sfmt, positions );
}
/**
* Get settings
*
* @return info string
*/
std::string PositionGenerator::getSettings( void ) const {
std::stringstream msg;
msg << "numMolecules " << _numMolecules << std::endl;
msg << "numParticlesPerMolecule " << _numParticlesPerMolecule << std::endl;
msg << "boxSize " << _boxSize << std::endl;
msg << "spacing " << _spacing[0] << std::endl;
msg << "seed " << _seed << std::endl;
return msg.str();
}
/**
* Get positions
*
* @param method method placement
* @param sfmt input random number generator
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int PositionGenerator::setPositions( GenerationMethod method, OpenMM_SFMT::SFMT& sfmt, std::vector<Vec3>& positions ) const {
int errorFlag = 0;
positions.resize( _numParticles );
if( method == Random ){
for( unsigned int ii = 0; ii < static_cast<unsigned int>(_numParticles); ii += _numParticlesPerMolecule ){
positions[ii] = Vec3(_boxSize*genrand_real2(sfmt), _boxSize*genrand_real2(sfmt), _boxSize*genrand_real2(sfmt));
for( unsigned int jj = 1; jj < static_cast<unsigned int>(_numParticlesPerMolecule); jj++) {
positions[ii+jj] = positions[ii] + Vec3(_bondDistance*genrand_real2(sfmt), _bondDistance*genrand_real2(sfmt), _bondDistance*genrand_real2(sfmt));
}
}
} else if( method == SimpleGrid ){
Vec3 origin, boxDimensions, spacing;
std::stringstream msg;
if( _numParticlesPerMolecule > 1 && _bondDistance > 0.0 ){
origin = Vec3(_bondDistance, _bondDistance, _bondDistance );
double particlesPerDimension = pow( static_cast<double>(_numMolecules), (1.0/3.0) );
int particlesPerDimensionI = static_cast<int>(particlesPerDimension+0.999999);
double boxSize = _boxSize;
double spacingPerDimension = (boxSize-_bondDistance)/(particlesPerDimension+1.0);
spacing = Vec3(spacingPerDimension, spacingPerDimension, spacingPerDimension );
boxDimensions = Vec3(boxSize, boxSize, boxSize );
msg << "Bond distance " << _bondDistance << std::endl;
msg << "particlesPerDimension " << particlesPerDimension << std::endl;
msg << "boxSize " << boxSize << std::endl;
msg << "spacingPerDimension " << spacingPerDimension << std::endl;
} else {
origin = _origin;
spacing = _spacing;
boxDimensions = _boxDimensions;
}
msg << getSettings() << std::endl;
if( _log ){
(void) fprintf( _log, "SimpleGrid %s\n", msg.str().c_str() );
}
errorFlag = setParticlesOnGrid( origin, boxDimensions, spacing, sfmt, _bondDistance, positions );
}
return errorFlag;
}
/**
* Place particles on a grid
*
* @param origin origin
* @param boxDimensions box dimensions
* @param spacing spacing
* @param array output vector of grid values
*
* @return -1 if particles will not fit on grid; 0 if they do
*/
int PositionGenerator::setParticlesOnGrid( const Vec3& origin, const Vec3& boxDimensions, const Vec3& spacing, OpenMM_SFMT::SFMT& sfmt,
double bondDistance, std::vector<Vec3>& array ) const {
const double pi = 3.14159265358979323846;
const double pi2 = 2.0*pi;
Vec3 start(origin);
if( array.size() != _numParticles ){
std::stringstream msg;
msg << "PositionGenerator::setParticlesOnGrid position vector size=" << array.size() << " != numParticles=" << _numParticles;
msg << getSettings();
throw OpenMMException( msg.str() );
}
// place molecule centers on grid
for( unsigned int ii = 0; ii < static_cast<unsigned int>(_numParticles); ii += _numParticlesPerMolecule ){
array[ii] = Vec3(start);
bool done = false;
for( unsigned int jj = 0; jj < 3 && !done; jj++ ){
start[jj] += spacing[jj];
if( (start[jj]+4.0*bondDistance) > boxDimensions[jj] ){
start[jj] = origin[jj];
} else {
done = true;
}
}
if( !done ){
std::stringstream msg;
msg << "PositionGenerator::setParticlesOnGrid error in grid settings at molecule=" << ii;
throw OpenMMException( msg.str() );
}
}
// add molecule atoms
for( unsigned int ii = 0; ii < static_cast<unsigned int>(_numMolecules); ii++ ){
int molecularIndex = ii*_numParticlesPerMolecule;
for( unsigned int jj = 1; jj < static_cast<unsigned int>(_numParticlesPerMolecule); jj++ ){
double theta = genrand_real2(sfmt)*pi2;
double phi = genrand_real2(sfmt)*pi;
array[molecularIndex+jj] = array[molecularIndex] + Vec3(_bondDistance*cos(theta)*cos(phi), _bondDistance*cos(theta)*sin(phi), _bondDistance*sin(theta) );
}
}
return 0;
}
/**
* Get enclosing box
*
* @param positions input vector of positions
* @param enclosingBox output Vec3[2] of minimum enclosing box ranges
*
*/
void PositionGenerator::getEnclosingBox( const std::vector<Vec3>& positions, Vec3 enclosingBox[2] ) const {
enclosingBox[0][0] = enclosingBox[1][0] = positions[0][0];
enclosingBox[0][1] = enclosingBox[1][1] = positions[0][1];
enclosingBox[0][2] = enclosingBox[1][2] = positions[0][2];
for( unsigned int ii = 1; ii < positions.size(); ii++ ){
if( enclosingBox[0][0] > positions[ii][0] ){
enclosingBox[0][0] = positions[ii][0];
}
if( enclosingBox[1][0] < positions[ii][0] ){
enclosingBox[1][0] = positions[ii][0];
}
if( enclosingBox[0][1] > positions[ii][1] ){
enclosingBox[0][1] = positions[ii][1];
}
if( enclosingBox[1][1] < positions[ii][1] ){
enclosingBox[1][1] = positions[ii][1];
}
if( enclosingBox[0][2] > positions[ii][2] ){
enclosingBox[0][2] = positions[ii][2];
}
if( enclosingBox[1][2] < positions[ii][2] ){
enclosingBox[1][2] = positions[ii][2];
}
}
return;
}
/**
* Show min/max distances between positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
* @param positionIndexVector list of entries to show min/max distances from
*
*/
void PositionGenerator::showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex,
const IntVector& positionIndexVector ){
if( !_log )return;
Vec3 box[2];
getEnclosingBox( positions, box );
(void) fprintf( _log, "Enclosing Box (in A): [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
box[0][0], box[1][0], box[0][1], box[1][1], box[0][2], box[1][2],
(box[1][0] - box[0][0]), (box[1][1] - box[0][1]), (box[1][2] - box[0][2]) );
for( unsigned int ii = 0; ii < positionIndexVector.size(); ii++ ){
if( positionIndexVector[ii] < static_cast<int>(positions.size()) ){
int positionIndex = positionIndexVector[ii];
IntDoublePairVector sortVector;
getSortedDistances( periodicBoundaryConditions, positionIndex, positions, sortVector );
(void) fprintf( _log, "Min/max distance from %6d:\n ", positionIndex );
for( unsigned int jj = 0; jj < sortVector.size() && jj < static_cast<unsigned int>(showIndex); jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( _log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( _log, "\n " );
for( unsigned int jj = (sortVector.size() - showIndex); jj < sortVector.size() && jj >= 0; jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( _log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( _log, "\n" );
}
}
return;
}
/**
* Show min/max distances between positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
*
*/
void PositionGenerator::showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex ){
if( !_log )return;
Vec3 box[2];
getEnclosingBox( positions, box );
(void) fprintf( _log, "Enclosing Box (in A): [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
box[0][0], box[1][0], box[0][1], box[1][1], box[0][2], box[1][2],
(box[1][0] - box[0][0]), (box[1][1] - box[0][1]), (box[1][2] - box[0][2]) );
IntDoublePairVector hitVector;
double minDistance = 1.0e+30;
double minDistanceCutoff = minDistance*1.1;
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
for( unsigned int jj = ii+1; jj < positions.size(); jj++ ){
double distance = periodicBoundaryConditions ? getPeriodicDistance( jj, ii, positions) :
getDistance( jj, ii, positions);
if( distance < minDistanceCutoff ){
if( distance < minDistance ){
minDistance = distance;
minDistanceCutoff = minDistance*1.1;
}
hitVector.push_back( IntDoublePair(ii*positions.size()+jj,distance ) );
}
}
}
std::sort( hitVector.begin(), hitVector.end(), TestIntDoublePair );
(void) fprintf( _log, "Min distances pbc=%d\n", periodicBoundaryConditions );
for( unsigned int jj = 0; jj < hitVector.size() && jj < static_cast<unsigned int>(showIndex); jj++ ){
IntDoublePair pair = hitVector[jj];
int index = pair.first;
int iIndex = static_cast<int>(index/positions.size());
int jIndex = index - iIndex*positions.size();
(void) fprintf( _log, " [%6d %6d %15.7e]\n", iIndex, jIndex, pair.second);
}
return;
}
/**
* Show particles within a specified distance of a given particle
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param particleIndex particle to check
* @param distanceToCheckFor distance to check for
* @param tolerance distance tolerance
*
*/
void PositionGenerator::showParticlesWithinDistance( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, unsigned int particleIndex,
double distanceToCheckFor, double tolerance){
if( !_log || particleIndex >= positions.size() )return;
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
double distance = periodicBoundaryConditions ? getPeriodicDistance( particleIndex, ii, positions) :
getDistance( particleIndex, ii, positions);
double delta = fabs( distanceToCheckFor - distance );
if( ii != particleIndex && delta < tolerance ){
(void) fprintf( _log, "Distance=%15.7e between particles %u %u.\n", distance, particleIndex, ii);
}
}
return;
}
/**
* Show distances between positions
*
* @param pairs particle indcies for which distance is to be reported
* @param positions input vector of positions
*
*/
void PositionGenerator::showDistances( const IntIntPairVector& pairs, const std::vector<Vec3>& positions ) const {
for( IntIntPairVectorCI ii = pairs.begin(); ii != pairs.end(); ii++ ){
if( ii->first < static_cast<int>(positions.size()) && ii->second < static_cast<int>(positions.size()) ){
double d = getDistance( ii->first, ii->second, positions );
(void) fprintf( _log, "Distance %6d %6d %15.7e d2=%15.7e\n", ii->first, ii->second, d, d*d );
}
}
return;
}
/**
* Get sorted distances from particular position
*
* @param periodicBoundaryConditions if set, apply PBC
* @param positionIndex input position index
* @param positions input vector of positions
* @param sortVector on output sorted IntDoublePairVector
*
*/
void PositionGenerator::getSortedDistances( int periodicBoundaryConditions, int positionIndex, const std::vector<Vec3>& positions,
IntDoublePairVector& sortVector ) const {
sortVector.resize( 0 );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
if( ii == positionIndex )continue;
double distance = periodicBoundaryConditions ? getPeriodicDistance( positionIndex, ii, positions) : getDistance( positionIndex, ii, positions);
sortVector.push_back( IntDoublePair( ii, distance ) );
}
std::sort( sortVector.begin(), sortVector.end(), TestIntDoublePair );
return;
}
/**---------------------------------------------------------------------------------------
*
* Set string field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setStringFromMap( MapStringString& argumentMap, std::string fieldToCheck, std::string& fieldToSet ){
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = (*check).second;
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set int field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setIntFromMap( MapStringString& argumentMap, std::string fieldToCheck, int& fieldToSet ){
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = atoi( (*check).second.c_str() );
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set int field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setIntFromMapStringToDouble( MapStringToDouble& argumentMap, std::string fieldToCheck, int& fieldToSet ){
// ---------------------------------------------------------------------------------------
MapStringToDoubleCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = static_cast<int>(check->second+0.0000001);
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
* Set float field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setFloatFromMap( MapStringString& argumentMap, std::string fieldToCheck, float& fieldToSet ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "setFloatFromMap";
// ---------------------------------------------------------------------------------------
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = static_cast<float>(atof( (*check).second.c_str() ));
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set double field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setDoubleFromMap( MapStringString& argumentMap, std::string fieldToCheck, double& fieldToSet ){
// ---------------------------------------------------------------------------------------
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = atof( (*check).second.c_str() );
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set double field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setDoubleFromMapStringToDouble( MapStringToDouble& argumentMap, std::string fieldToCheck, double& fieldToSet ){
// ---------------------------------------------------------------------------------------
MapStringToDoubleCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = check->second;
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Get relative difference between two forces
*
* @param f1 force1
* @param f2 force2
* @param forceNorm1 output norm of force1
* @param forceNorm2 output norm of force2
* @param relativeDiff output relative difference between force norms
* @param log if set, output forces
*
*
--------------------------------------------------------------------------------------- */
static void getForceRelativeDifference( const Vec3& f1, const Vec3& f2, double& forceNorm1, double& forceNorm2,
double& relativeDiff, FILE* log ) {
double diff = (f1[0] - f2[0])*(f1[0] - f2[0]) +
(f1[1] - f2[1])*(f1[1] - f2[1]) +
(f1[2] - f2[2])*(f1[2] - f2[2]);
forceNorm1 = sqrt( f1[0]*f1[0] + f1[1]*f1[1] + f1[2]*f1[2] );
forceNorm2 = sqrt( f2[0]*f2[0] + f2[1]*f2[1] + f2[2]*f2[2] );
if( forceNorm1 > 0.0 || forceNorm2 > 0.0 ){
relativeDiff = 2.0*sqrt( diff )/(forceNorm1+forceNorm2);
} else {
relativeDiff = 0.0;
}
return;
}
/**---------------------------------------------------------------------------------------
*
* Compare forces from two states
*
* @param state1 state1
* @param state2 state2
* @param relativeTolerance relative tolerance
* @param log if set, output forces
*
* @return number of entries with relative difference > tolerance
*
--------------------------------------------------------------------------------------- */
int compareForcesOfTwoStates( State& state1, State& state2, double relativeTolerance,
DoubleVector& stats, FILE* log ) {
int error = 0;
vector<Vec3> force1 = state1.getForces();
vector<Vec3> force2 = state2.getForces();
double averageRelativeDifference = 0.0;
double count = 0.0;
DoubleVector medians1( force1.size() );
DoubleVector medians2( force1.size() );
IntDoublePairVector relativeDifferences;
for( unsigned int ii = 0; ii < force1.size(); ii++ ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
getForceRelativeDifference( force1[ii], force2[ii], forceNorm1, forceNorm2, relativeDiff, log );
medians1[ii] = forceNorm1;
medians2[ii] = forceNorm2;
relativeDifferences.push_back( IntDoublePair(ii, relativeDiff ) );
averageRelativeDifference += relativeDiff;
count += 1.0;
if( relativeDiff > relativeTolerance ){
error++;
}
if( log ){
(void) fprintf( log, "F %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n", static_cast<unsigned int>(ii),
relativeDiff, force1[ii][0], force1[ii][1], force1[ii][2], force2[ii][0], force2[ii][1], force2[ii][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
// sort relative differences
std::sort( relativeDifferences.begin(), relativeDifferences.end(), TestIntDoublePair );
if( log ){
(void) fprintf( log, "\nEntries w/ largest relative differences.\n" );
for( unsigned int ii = relativeDifferences.size()-1; ii >= relativeDifferences.size()-10 && ii >= 0; ii-- ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
int index = relativeDifferences[ii].first;
getForceRelativeDifference( force1[index], force2[index], forceNorm1, forceNorm2, relativeDiff, log );
(void) fprintf( log, "Fs %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n",
static_cast<unsigned int>(index), relativeDiff,
force1[index][0], force1[index][1], force1[index][2],
force2[index][0], force2[index][1], force2[index][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
if( count > 0.0 ){
averageRelativeDifference /= count;
}
std::sort( medians1.begin(), medians1.end() );
std::sort( medians2.begin(), medians2.end() );
double median1 = medians1[medians1.size()/2];
double median2 = medians2[medians2.size()/2];
stats.resize( 4 );
stats[0] = averageRelativeDifference;
IntDoublePair pair = relativeDifferences[relativeDifferences.size()-1];
stats[1] = pair.second;
stats[2] = static_cast<double>(pair.first);
stats[3] = median1 < median2 ? median1 : median2;
return error;
}
/**
* Create nonbonded force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param reactionFieldDielectric reaction field dielectric
* @param parameterList list of parameters -- used via addParticle()
* @param bonds list of BondInfo_OpenMMTest containing info for exceptions
* @param log logging file (optional -- may be NULL)
*
*/
static NonbondedForce* getNonbondedForce( int nonbondedMethod, double cutoffDistance, double reactionFieldDielectric,
VectorOfDoubleVectors& parameterList, std::vector< BondInfo_OpenMMTest >& bonds, FILE* log ){
NonbondedForce* nonbondedForce = new NonbondedForce();
NonbondedForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = NonbondedForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = NonbondedForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = NonbondedForce::CutoffPeriodic;
break;
case Ewald_OpenMMTest:
method = NonbondedForce::Ewald;
break;
case PME_OpenMMTest:
method = NonbondedForce::PME;
break;
default:
method = NonbondedForce::NoCutoff;
}
nonbondedForce->setNonbondedMethod( method );
nonbondedForce->setCutoffDistance( cutoffDistance );
nonbondedForce->setReactionFieldDielectric( reactionFieldDielectric );
// load parameters
for( unsigned int ii = 0; ii < parameterList.size(); ii++ ){
DoubleVector parameters = parameterList[ii];
nonbondedForce->addParticle( parameters[ChargeIndex_OpenMMTest], parameters[SigmaIndex_OpenMMTest], parameters[EpsIndex_OpenMMTest] );
}
// add exceptions
for( unsigned int ii = 0; ii < bonds.size(); ii++ ){
BondInfo_OpenMMTest bond = bonds[ii];
nonbondedForce->addException( bond._particle1, bond._particle2, 0.0f, 1.0, 0.0f );
}
return nonbondedForce;
}
/**
* Create GBVI force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param useQuinticSpline if set use quintic spline for Born radii scaling; else use no scaling
* @param quinticLowerLimitFactor quintic lower limit factor
* @param quinticUpperBornRadiusLimit quintic upper Born radius limit
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param parameterList list of parameters -- used via addParticle()
* @param bonds list of BondInfo_OpenMMTest containing info for exceptions
* @param log logging file (optional -- may be NULL)
*
*/
static GBVIForce* getGBVIForce( int nonbondedMethod, double cutoffDistance, int useQuinticSpline,
double quinticLowerLimitFactor, double quinticUpperBornRadiusLimit,
double solventDielectric, double soluteDiecletric,
VectorOfDoubleVectors& parameterList, std::vector< BondInfo_OpenMMTest >& bonds, FILE* log ){
GBVIForce* gbviForce = new GBVIForce();
GBVIForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = GBVIForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBVIForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBVIForce::CutoffPeriodic;
break;
default:
method = GBVIForce::NoCutoff;
}
gbviForce->setNonbondedMethod( method );
gbviForce->setCutoffDistance( cutoffDistance );
gbviForce->setSolventDielectric( solventDielectric );
gbviForce->setSoluteDielectric( soluteDiecletric );
if( useQuinticSpline ){
gbviForce->setBornRadiusScalingMethod( GBVIForce::QuinticSpline );
gbviForce->setQuinticLowerLimitFactor( quinticLowerLimitFactor );
gbviForce->setQuinticUpperBornRadiusLimit( quinticUpperBornRadiusLimit );
} else {
gbviForce->setBornRadiusScalingMethod( GBVIForce::NoScaling );
}
// load parameters
for( unsigned int ii = 0; ii < parameterList.size(); ii++ ){
DoubleVector parameters = parameterList[ii];
gbviForce->addParticle( parameters[ChargeIndex_OpenMMTest], parameters[SigmaIndex_OpenMMTest], parameters[GammaIndex_OpenMMTest] );
}
// add exceptions
for( unsigned int ii = 0; ii < bonds.size(); ii++ ){
BondInfo_OpenMMTest bond = bonds[ii];
gbviForce->addBond( bond._particle1, bond._particle2, bond._distance);
}
return gbviForce;
}
/**
* Create OBC force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param parameterList list of parameters -- used via addParticle()
* @param log logging file (optional -- may be NULL)
*
*/
static GBSAOBCForce* getGBSAOBCForce( int nonbondedMethod, double cutoffDistance, double solventDielectric, double soluteDiecletric,
VectorOfDoubleVectors& parameterList, FILE* log ){
GBSAOBCForce* obcForce = new GBSAOBCForce();
GBSAOBCForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = GBSAOBCForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBSAOBCForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBSAOBCForce::CutoffPeriodic;
break;
default:
method = GBSAOBCForce::NoCutoff;
}
obcForce->setNonbondedMethod( method );
obcForce->setCutoffDistance( cutoffDistance );
obcForce->setSolventDielectric( solventDielectric );
obcForce->setSoluteDielectric( soluteDiecletric );
// load parameters
for( unsigned int ii = 0; ii < parameterList.size(); ii++ ){
DoubleVector parameters = parameterList[ii];
obcForce->addParticle( parameters[ChargeIndex_OpenMMTest], parameters[SigmaIndex_OpenMMTest], parameters[GammaIndex_OpenMMTest] );
}
return obcForce;
}
/**
* Create nonbonded softcore force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param reactionFieldDielectric reaction field dielectric
* @param log logging file (optional -- may be NULL)
*
*/
#ifdef USE_SOFTCORE
static NonbondedSoftcoreForce* getNonbondedSoftcoreForce( int nonbondedMethod, double cutoffDistance, double reactionFieldDielectric, FILE* log ){
NonbondedSoftcoreForce* nonbondedForce = new NonbondedSoftcoreForce();
NonbondedSoftcoreForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = NonbondedSoftcoreForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = NonbondedSoftcoreForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = NonbondedSoftcoreForce::CutoffPeriodic;
break;
default:
method = NonbondedSoftcoreForce::NoCutoff;
}
nonbondedForce->setNonbondedMethod( method );
nonbondedForce->setCutoffDistance( cutoffDistance );
nonbondedForce->setReactionFieldDielectric( reactionFieldDielectric );
return nonbondedForce;
}
/**
* Create GBVI softcore force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param useQuinticSpline if set use quintic spline for Born radii scaling; else use no scaling
* @param quinticLowerLimitFactor quintic lower limit factor
* @param quinticUpperBornRadiusLimit quintic upper Born radius limit
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param log logging file (optional -- may be NULL)
*
*/
static GBVISoftcoreForce* getGBVISoftcoreForce( int nonbondedMethod, double cutoffDistance, int useQuinticSpline,
double quinticLowerLimitFactor, double quinticUpperBornRadiusLimit,
double solventDielectric, double soluteDiecletric, FILE* log ){
GBVISoftcoreForce* gbviForce = new GBVISoftcoreForce();
GBVISoftcoreForce::NonbondedMethod method;
switch( gbviMethod ){
case NoCutoff_OpenMMTest:
method = GBVISoftcoreForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBVISoftcoreForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBVISoftcoreForce::CutoffPeriodic;
break;
default:
method = GBVISoftcoreForce::NoCutoff;
}
gbviForce->setNonbondedMethod( method );
gbviForce->setCutoffDistance( cutoffDistance );
gbviForce->setSolventDielectric( solventDielectric );
gbviForce->setSoluteDielectric( soluteDiecletric );
if( useQuinticSpline ){
gbviForce->setBornRadiusScalingMethod( GBVISoftcoreForce::QuinticSpline );
gbviForce->setQuinticLowerLimitFactor( quinticLowerLimitFactor );
gbviForce->setQuinticUpperBornRadiusLimit( quinticUpperBornRadiusLimit );
} else {
gbviForce->setBornRadiusScalingMethod( GBVISoftcoreForce::NoScaling );
}
return gbviForce;
}
/**
* Create OBC softcore force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param log logging file (optional -- may be NULL)
*
*/
static GBSAOBCSoftcoreForce* getGBSAOBCSoftcoreForce( int nonbondedMethod, double cutoffDistance,
double solventDielectric, double soluteDiecletric, FILE* log ){
GBSAOBCSoftcoreForce* obcForce = new GBSAOBCSoftcoreForce();
GBSAOBCSoftcoreForce::NonbondedMethod method;
switch( gbviMethod ){
case NoCutoff_OpenMMTest:
method = GBSAOBCSoftcoreForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBSAOBCSoftcoreForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBSAOBCSoftcoreForce::CutoffPeriodic;
break;
default:
method = GBSAOBCSoftcoreForce::NoCutoff;
}
obcForce->setNonbondedMethod( method );
obcForce->setCutoffDistance( cutoffDistance );
obcForce->setSolventDielectric( solventDielectric );
obcForce->setSoluteDielectric( soluteDiecletric );
return obcForce;
}
#endif
/**
* Get forces in system
*
* @param system system to serialize
* @param stringForceVector output stringForceVector[forceName] = force index
* @param log logging file (optional -- may be NULL)
*
*/
static void getStringForceMap( System& system, MapStringInt& stringForceVector, FILE* log ){
// print active forces and relevant parameters
for( int ii = 0; ii < system.getNumForces(); ii++ ) {
int hit = 0;
Force& force = system.getForce(ii);
if( !hit ){
try {
CMAPTorsionForce& castForce = dynamic_cast<CMAPTorsionForce&>(force);
stringForceVector["CMAPTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomAngleForce& castForce = dynamic_cast<CustomAngleForce&>(force);
stringForceVector["CustomAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomBondForce& castForce = dynamic_cast<CustomBondForce&>(force);
stringForceVector["CustomBond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomExternalForce& castForce = dynamic_cast<CustomExternalForce&>(force);
stringForceVector["CustomExternal"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomGBForce& castForce = dynamic_cast<CustomGBForce&>(force);
stringForceVector["CustomGB"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomHbondForce& castForce = dynamic_cast<CustomHbondForce&>(force);
stringForceVector["CustomHbond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomNonbondedForce& castForce = dynamic_cast<CustomNonbondedForce&>(force);
stringForceVector["CustomNonbonded"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomTorsionForce& castForce = dynamic_cast<CustomTorsionForce&>(force);
stringForceVector["CustomTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
GBSAOBCForce& castForce = dynamic_cast<GBSAOBCForce&>(force);
stringForceVector["GBSAOBC"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
GBVIForce& castForce = dynamic_cast<GBVIForce&>(force);
stringForceVector["GBVI"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
HarmonicAngleForce& castForce = dynamic_cast<HarmonicAngleForce&>(force);
stringForceVector["HarmonicAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
HarmonicBondForce& castForce = dynamic_cast<HarmonicBondForce&>(force);
stringForceVector["HarmonicBond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
NonbondedForce& castForce = dynamic_cast<NonbondedForce&>(force);
stringForceVector["Nonbonded"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
PeriodicTorsionForce& castForce = dynamic_cast<PeriodicTorsionForce&>(force);
stringForceVector["PeriodicTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
RBTorsionForce& castForce = dynamic_cast<RBTorsionForce&>(force);
stringForceVector["RBTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
MonteCarloBarostat& castForce = dynamic_cast<MonteCarloBarostat&>(force);
stringForceVector["MonteCarloBarostat"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AndersenThermostat& castForce = dynamic_cast<AndersenThermostat&>(force);
stringForceVector["AndersenThermostat"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
#ifdef USE_SOFTCORE
if( !hit ){
try {
GBSAOBCSoftcoreForce& castForce = dynamic_cast<GBSAOBCSoftcoreForce&>(force);
stringForceVector["GBSAOBCSoftcore"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
GBVISoftcoreForce& castForce = dynamic_cast<GBVISoftcoreForce&>(force);
stringForceVector["GBVISoftcore"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
NonbondedSoftcoreForce& castForce = dynamic_cast<NonbondedSoftcoreForce&>(force);
stringForceVector["NonbondedSoftcore"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
#endif
#ifdef INCLUDE_AMOEBA_FORCES
if( !hit ){
try {
AmoebaHarmonicBondForce& castForce = dynamic_cast<AmoebaHarmonicBondForce&>(force);
stringForceVector["AmoebaHarmonicBond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaHarmonicAngleForce& castForce = dynamic_cast<AmoebaHarmonicAngleForce&>(force);
stringForceVector["AmoebaHarmonicAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaHarmonicInPlaneAngleForce& castForce = dynamic_cast<AmoebaHarmonicInPlaneAngleForce&>(force);
stringForceVector["AmoebaHarmonicInPlaneAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaMultipoleForce& castForce = dynamic_cast<AmoebaMultipoleForce&>(force);
stringForceVector["AmoebaMultipole"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaOutOfPlaneBendForce& castForce = dynamic_cast<AmoebaOutOfPlaneBendForce&>(force);
stringForceVector["AmoebaOutOfPlaneBend"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaPiTorsionForce& castForce = dynamic_cast<AmoebaPiTorsionForce&>(force);
stringForceVector["AmoebaPiTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaStretchBendForce& castForce = dynamic_cast<AmoebaStretchBendForce&>(force);
stringForceVector["AmoebaStretchBend"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaTorsionForce& castForce = dynamic_cast<AmoebaTorsionForce&>(force);
stringForceVector["AmoebaTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaTorsionTorsionForce& castForce = dynamic_cast<AmoebaTorsionTorsionForce&>(force);
stringForceVector["AmoebaTorsionTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaUreyBradleyForce& castForce = dynamic_cast<AmoebaUreyBradleyForce&>(force);
stringForceVector["AmoebaUreyBradley"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaVdwForce& castForce = dynamic_cast<AmoebaVdwForce&>(force);
stringForceVector["AmoebaVdw"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaWcaDispersionForce& castForce = dynamic_cast<AmoebaWcaDispersionForce&>(force);
stringForceVector["AmoebaWcaDispersion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaGeneralizedKirkwoodForce& castForce = dynamic_cast<AmoebaGeneralizedKirkwoodForce&>(force);
stringForceVector["AmoebaGeneralizedKirkwood"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaTorsionTorsionForce& castForce = dynamic_cast<AmoebaTorsionTorsionForce&>(force);
stringForceVector["AmoebaTorsionTorsionForce"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
#endif
// COM
if( !hit ){
try {
CMMotionRemover& cMMotionRemover = dynamic_cast<CMMotionRemover&>(force);
hit++;
} catch( std::bad_cast ){
}
}
if( !hit && log ){
(void) fprintf( log, " entry=%2d force not recognized.\n", ii );
}
}
}
/**
* Get forces in system
*
* @param system system to serialize
* @param stringForceVector output stringForceVector[forceName] = force index
* @param log logging file (optional -- may be NULL)
*
*/
static Force* copyForce( const Force& force, FILE* log ){
// print active forces and relevant parameters
Force* forceCopy = NULL;
try {
const CMAPTorsionForce& castForce = dynamic_cast<const CMAPTorsionForce&>(force);
forceCopy = new CMAPTorsionForce( castForce );
} catch( std::bad_cast ){
}
if( forceCopy == NULL ){
try {
const CustomAngleForce& castForce = dynamic_cast<const CustomAngleForce&>(force);
forceCopy = new CustomAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomBondForce& castForce = dynamic_cast<const CustomBondForce&>(force);
forceCopy = new CustomBondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomExternalForce& castForce = dynamic_cast<const CustomExternalForce&>(force);
forceCopy = new CustomExternalForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomGBForce& castForce = dynamic_cast<const CustomGBForce&>(force);
forceCopy = new CustomGBForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomHbondForce& castForce = dynamic_cast<const CustomHbondForce&>(force);
forceCopy = new CustomHbondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomNonbondedForce& castForce = dynamic_cast<const CustomNonbondedForce&>(force);
forceCopy = new CustomNonbondedForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomTorsionForce& castForce = dynamic_cast<const CustomTorsionForce&>(force);
forceCopy = new CustomTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const GBSAOBCForce& castForce = dynamic_cast<const GBSAOBCForce&>(force);
forceCopy = new GBSAOBCForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const GBVIForce& castForce = dynamic_cast<const GBVIForce&>(force);
forceCopy = new GBVIForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const HarmonicAngleForce& castForce = dynamic_cast<const HarmonicAngleForce&>(force);
forceCopy = new HarmonicAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const HarmonicBondForce& castForce = dynamic_cast<const HarmonicBondForce&>(force);
forceCopy = new HarmonicBondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const NonbondedForce& castForce = dynamic_cast<const NonbondedForce&>(force);
forceCopy = new NonbondedForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const PeriodicTorsionForce& castForce = dynamic_cast<const PeriodicTorsionForce&>(force);
forceCopy = new PeriodicTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const RBTorsionForce& castForce = dynamic_cast<const RBTorsionForce&>(force);
forceCopy = new RBTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const MonteCarloBarostat& castForce = dynamic_cast<const MonteCarloBarostat&>(force);
forceCopy = new MonteCarloBarostat( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AndersenThermostat& castForce = dynamic_cast<const AndersenThermostat&>(force);
forceCopy = new AndersenThermostat( castForce );
} catch( std::bad_cast ){
}
}
#ifdef USE_SOFTCORE
if( forceCopy == NULL ){
try {
const GBSAOBCSoftcoreForce& castForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
forceCopy = new GBSAOBCSoftcoreForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const GBVISoftcoreForce& castForce = dynamic_cast<const GBVISoftcoreForce&>(force);
forceCopy = new GBVISoftcoreForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const NonbondedSoftcoreForce& castForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
forceCopy = new NonbondedSoftcoreForce( castForce );
} catch( std::bad_cast ){
}
}
#endif
#ifdef INCLUDE_AMOEBA_FORCES
if( forceCopy == NULL ){
try {
const AmoebaHarmonicBondForce& castForce = dynamic_cast<const AmoebaHarmonicBondForce&>(force);
forceCopy = new AmoebaHarmonicBondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaHarmonicAngleForce& castForce = dynamic_cast<const AmoebaHarmonicAngleForce&>(force);
forceCopy = new AmoebaHarmonicAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaHarmonicInPlaneAngleForce& castForce = dynamic_cast<const AmoebaHarmonicInPlaneAngleForce&>(force);
forceCopy = new AmoebaHarmonicInPlaneAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaMultipoleForce& castForce = dynamic_cast<const AmoebaMultipoleForce&>(force);
forceCopy = new AmoebaMultipoleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaOutOfPlaneBendForce& castForce = dynamic_cast<const AmoebaOutOfPlaneBendForce&>(force);
forceCopy = new AmoebaOutOfPlaneBendForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaPiTorsionForce& castForce = dynamic_cast<const AmoebaPiTorsionForce&>(force);
forceCopy = new AmoebaPiTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaStretchBendForce& castForce = dynamic_cast<const AmoebaStretchBendForce&>(force);
forceCopy = new AmoebaStretchBendForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaTorsionForce& castForce = dynamic_cast<const AmoebaTorsionForce&>(force);
forceCopy = new AmoebaTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaTorsionTorsionForce& castForce = dynamic_cast<const AmoebaTorsionTorsionForce&>(force);
forceCopy = new AmoebaTorsionTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaUreyBradleyForce& castForce = dynamic_cast<const AmoebaUreyBradleyForce&>(force);
forceCopy = new AmoebaUreyBradleyForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaVdwForce& castForce = dynamic_cast<const AmoebaVdwForce&>(force);
forceCopy = new AmoebaVdwForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaWcaDispersionForce& castForce = dynamic_cast<const AmoebaWcaDispersionForce&>(force);
forceCopy = new AmoebaWcaDispersionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaGeneralizedKirkwoodForce& castForce = dynamic_cast<const AmoebaGeneralizedKirkwoodForce&>(force);
forceCopy = new AmoebaGeneralizedKirkwoodForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaTorsionTorsionForce& castForce = dynamic_cast<const AmoebaTorsionTorsionForce&>(force);
forceCopy = new AmoebaTorsionTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
#endif
if( log && forceCopy == NULL ){
(void) fprintf( log, " force not recognized.\n" );
}
return forceCopy;
}
/**
* Return copy of system (but not forces)
*
* @param inputSystem system to copy
*
* @return system copy
*
*/
static void copySystem( const System& inputSystem, System& systemCopy, FILE* log ){
// add particle/mass
for( unsigned int ii = 0; ii < static_cast<unsigned int>(inputSystem.getNumParticles()); ii++ ){
systemCopy.addParticle( inputSystem.getParticleMass( static_cast<int>(ii) ) );
}
// box
Vec3 a;
Vec3 b;
Vec3 c;
inputSystem.getDefaultPeriodicBoxVectors( a, b, c );
systemCopy.setDefaultPeriodicBoxVectors( a, b, c );
// copy constraints
for( unsigned int ii = 0; ii < static_cast<unsigned int>(inputSystem.getNumConstraints()); ii++ ){
int particle1, particle2;
double distance;
inputSystem.getConstraintParameters( ii, particle1, particle2, distance);
systemCopy.addConstraint( particle1, particle2, distance);
}
// copy forces
for( unsigned int ii = 0; ii < static_cast<unsigned int>(inputSystem.getNumForces()); ii++ ){
systemCopy.addForce( copyForce( inputSystem.getForce(ii), log) );
}
return;
}
/**
* Randomize parameters
*
* @param parametersLowerBound vector of parameter lower bounds
* @param parametersUpperBound vector of parameter upper bounds
* @param sfmt SFMT random number generator
* @param parameters output vector of randomized parameter values
*
*/
static void randomizeParameters( const std::vector<double>& parametersLowerBound,
const std::vector<double>& parametersUpperBound,
OpenMM_SFMT::SFMT& sfmt, std::vector<double>& parameters ){
if( parametersLowerBound.size() != parametersUpperBound.size() ){
std::stringstream msg;
msg << " randomizeParameters parametersLowerBound size=" << parametersLowerBound.size() << " != parametersUpperBound size=" << parametersUpperBound.size();
throw OpenMMException( msg.str() );
}
if( parametersLowerBound.size() != parameters.size() ){
std::stringstream msg;
msg << " randomizeParameters parametersLowerBound size=" << parametersLowerBound.size() << " != parameter size=" << parameters.size();
throw OpenMMException( msg.str() );
}
for( unsigned int ii = 0; ii < parametersLowerBound.size(); ii++ ){
parameters[ii] = parametersLowerBound[ii] + (parametersUpperBound[ii] - parametersLowerBound[ii])*(genrand_real2(sfmt));
}
return;
}
/**
* Randomize Vec3 vector
*
* @param average mean value
* @param range +/- range
* @param sfmt SFMT random number generator
* @param array output vector of randomized values
*
*/
static void randomizeVec3( double average, double range,
OpenMM_SFMT::SFMT& sfmt, std::vector<Vec3>& array ){
range *= 2.0;
for( unsigned int ii = 0; ii < array.size(); ii++ ){
array[ii] = Vec3( average + range*(genrand_real2(sfmt) - 0.5),
average + range*(genrand_real2(sfmt) - 0.5),
average + range*(genrand_real2(sfmt) - 0.5) );
}
return;
}
/**
* Output contents of MapStringString
*
* @param inputArgumentMap map to output
* @param outputStream output stream
*
*/
static void streamArgumentMap( const MapStringString& inputArgumentMap, std::stringstream& outputStream ){
char buffer[2048];
for( MapStringStringCI ii = inputArgumentMap.begin(); ii != inputArgumentMap.end(); ii++ ){
std::string key = ii->first;
std::string value = ii->second;
(void) sprintf( buffer, " %30s %40s\n", key.c_str(), value.c_str() );
outputStream << buffer;
}
return;
}
/**
* Format argument/value
*
* @param buffer formatted output
* @param key argument name
* @param value argument value
* @param format format string
* @param call if call > 0, skip key name
* @param type type == 0, then use int value; else double
*
*/
static void formatArgument( char* buffer, const std::string& key, double value, const char* format, int call, int type ){
// if call > 0, skip key name
unsigned int index = 0;
while( index < key.size() ){
buffer[index] = call ? ' ' : key[index];
index++;
}
// add blank
buffer[index++] = ' ';
buffer[index] = static_cast<char>(NULL);
if( type == 0 ){
int valueInt = static_cast<int>(value+0.00001);
(void) sprintf( buffer + index, format, valueInt );
} else {
(void) sprintf( buffer + index, format, value );
}
return;
}
/**
* Output contents of MapStringString w/ all argument on one line
*
* @param inputArgumentMap map to output
* @param exclude map of keys to exclude from output
* @param outputStream output stream
*
*/
static void streamArgumentMapOneLine( const MapStringToDouble& inputArgumentMap, const MapStringToInt& exclude,
const StringVector& printFirst, int callId, std::stringstream& outputStream ){
char buffer[2048];
MapStringToInt excludeAll(exclude);
for( unsigned int ii = 0; ii < printFirst.size(); ii++ ){
MapStringToDoubleCI iter = inputArgumentMap.find( printFirst[ii] );
if( iter != inputArgumentMap.end() ){
std::string key = iter->first;
if( exclude.find( key ) == exclude.end() ){
double value = iter->second;
if( key == "numMolecules" ){
formatArgument( buffer, key, value, "%6d ", callId, 0 );
} else if( key == "nonbondedMethod" ){
formatArgument( buffer, key, value, "%1d ", callId, 0 );
} else if( key == "lambda1" || key == "lambda2" ){
formatArgument( buffer, key, value, "%4.2f ", callId, 1 );
} else if( key == "boxSize" ){
formatArgument( buffer, key, value, "%6.2f ", callId, 1 );
} else if( key == "cutoffDistance" ){
formatArgument( buffer, key, value, "%7.3f ", callId, 1 );
} else if( key == "relativeTolerance" ){
formatArgument( buffer, key, value, "%8.1e ", callId, 1 );
} else if( key == "positionPlacementMethod" || key == "applyAssert" || key == "serialize" ){
formatArgument( buffer, key, value, "%1d ", callId, 1 );
} else {
formatArgument( buffer, key, value, "%15.7e ", callId, 1 );
}
outputStream << buffer;
excludeAll[key] = 1;
}
}
}
for( MapStringToDoubleCI ii = inputArgumentMap.begin(); ii != inputArgumentMap.end(); ii++ ){
std::string key = ii->first;
if( excludeAll.find( key ) == excludeAll.end() ){
double value = ii->second;
int valueInt = static_cast<int>(value+0.00001);
double valueDouble = static_cast<double>(valueInt);
if( key == "numMolecules" ){
(void) sprintf( buffer, "%s=%6d ", key.c_str(), valueInt );
} else if( key == "nonbondedMethod" || key == "positionPlacementMethod" || key == "applyAssert" || key == "serialize" ){
(void) sprintf( buffer, "%s=%1d ", key.c_str(), valueInt );
} else if( key == "lambda1" || key == "lambda2" ){
(void) sprintf( buffer, "%s=%4.2f ", key.c_str(), value );
} else if( key == "boxSize" || key == "cutoffDistance" ){
(void) sprintf( buffer, "%s=%6.2f ", key.c_str(), value );
} else if( key == "relativeTolerance" ){
(void) sprintf( buffer, "%s=%8.1e ", key.c_str(), value );
} else if( valueDouble == value ){
(void) sprintf( buffer, "%s=%6d ", key.c_str(), valueInt );
} else {
(void) sprintf( buffer, "%s=%15.7e ", key.c_str(), value );
}
outputStream << buffer;
}
}
outputStream << std::endl;
return;
}
/**
* Get signature of a MapStringToDouble object
*
* @param inputArgumentMap map
* @return signature
*
*/
static double getMapStringToDoubleSignature( const MapStringToDouble& inputArgumentMap ){
double signature = 0.0;
double offset = 0.1;
for( MapStringToDoubleCI ii = inputArgumentMap.begin(); ii != inputArgumentMap.end(); ii++ ){
signature += (offset + ii->second);
offset += 0.1;
}
return signature;
}
/**
* Compare two MapStringToDouble to see if they have the same (key,value) pairs
*
* @param inputArgumentMap1 map 1
* @param inputArgumentMap2 map 2
*
* @return true if maps have same (key,value) pairs; otherwise false
*
*/
static bool compareMapStringToDoubles( const MapStringToDouble& inputArgumentMap1, const MapStringToDouble& inputArgumentMap2 ){
if( inputArgumentMap1.size() != inputArgumentMap1.size() ){
return false;
}
for( MapStringToDoubleCI ii = inputArgumentMap1.begin(); ii != inputArgumentMap1.end(); ii++ ){
MapStringToDoubleCI jj = inputArgumentMap2.find( (*ii).first );
if( jj == inputArgumentMap2.end() || jj->second != ii->second ){
return false;
}
}
return true;
}
/**
* Generate collection of inputArguments maps given
* list of DoubleVectors for each argument
*
* @param inputArguments map[argumentKey] = vector of double parameter values
* @param argumentMaps output vector of generated maps
*
*/
static void generateInputArgumentMapsFromStringVectors( const MapStringToDoubleVector& inputArguments,
VectorOfMapStringToDouble& argumentMaps ){
for( MapStringToDoubleVectorCI ii = inputArguments.begin(); ii != inputArguments.end(); ii++ ){
std::string argumentName = (*ii).first;
DoubleVector arguments = (*ii).second;
unsigned int initialArgumentMapSize = argumentMaps.size();
// generate signature map for each argument map
MapDoubleToInt signatures;
for( unsigned int kk = 0; kk < initialArgumentMapSize; kk++ ){
double signature = getMapStringToDoubleSignature( argumentMaps[kk] );
signatures[signature] = 1;
}
// for each current argumment map, add a new argument map w/ (key,value)
// check that no existing map has the same arguments before adding to the
// vector of argument maps
for( unsigned int kk = 0; kk < initialArgumentMapSize; kk++ ){
for( unsigned int jj = 0; jj < arguments.size(); jj++ ){
MapStringToDouble inputArgumentMap = MapStringToDouble(argumentMaps[kk]);
inputArgumentMap[argumentName] = arguments[jj];
double signature = getMapStringToDoubleSignature( inputArgumentMap );
if( signatures.find( signature ) == signatures.end() ){
argumentMaps.push_back( inputArgumentMap );
} else {
bool match = 0;
for( unsigned int mm = 0; mm < initialArgumentMapSize && !match; mm++ ){
match = compareMapStringToDoubles( inputArgumentMap, argumentMaps[mm] );
}
if( !match ){
argumentMaps.push_back( inputArgumentMap );
}
}
}
}
}
return;
}
/**
* Predicate for sorting map[string] = double
*
* @param d1 first MapStringToDouble to compare
* @param d2 second MapStringToDouble to compare
*
*/
bool TestMapSortPredicate( const MapStringToDouble& d1, const MapStringToDouble& d2 ){
StringVector sortOrder;
sortOrder.push_back( "numMolecules" );
sortOrder.push_back( "nonbondedMethod" );
sortOrder.push_back( "lambda2" );
sortOrder.push_back( "boxSize" );
for( unsigned int ii = 0; ii < sortOrder.size(); ii++ ){
if( d1.find( sortOrder[ii] ) != d1.end() &&
d2.find( sortOrder[ii] ) != d2.end() ){
MapStringToDoubleCI d1i = d1.find( sortOrder[ii] );
MapStringToDoubleCI d2i = d2.find( sortOrder[ii] );
if( d1i->second != d2i->second ){
return d1i->second < d2i->second;
}
}
}
return false;
}
#ifdef USE_SOFTCORE
static CustomNonbondedForce* buildCustomNonbondedSoftcoreForce( const NonbondedSoftcoreForce& nonbondedSoftcoreForce ){
CustomNonbondedForce* customNonbonded;
if( nonbondedSoftcoreForce.getNonbondedMethod() == NoCutoff ){
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customNonbonded->setNonbondedMethod( CustomNonbondedForce::NoCutoff );
} else {
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf);"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customNonbonded->setCutoffDistance( nonbondedSoftcoreForce.getCutoffDistance() );
if( nonbondedSoftcoreForce.getNonbondedMethod() == CutoffNonPeriodic ){
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffNonPeriodic );
} else {
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffPeriodic );
}
double cutoffDistance = nonbondedSoftcoreForce.getCutoffDistance();
double reactionFieldDielectric = nonbondedSoftcoreForce.getReactionFieldDielectric();
double eps2 = (reactionFieldDielectric - 1.0)/(2.0*reactionFieldDielectric+1.0);
double kValue = eps2/(cutoffDistance*cutoffDistance*cutoffDistance);
customNonbonded->addGlobalParameter("krf", kValue );
double cValue = (1.0/cutoffDistance)*(3.0*reactionFieldDielectric)/(2.0*reactionFieldDielectric + 1.0);
customNonbonded->addGlobalParameter("crf", cValue );
}
customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps");
customNonbonded->addPerParticleParameter("lambda");
vector<double> nonbondedParams(4);
for( unsigned int ii = 0; ii < nonbondedSoftcoreForce.getNumParticles(); ii++ ){
double charge;
double sigma;
double epsilon;
double softcoreLJLambda;
nonbondedSoftcoreForce.getParticleParameters(ii, charge, sigma, epsilon, softcoreLJLambda);
nonbondedParams[0] = charge;
nonbondedParams[1] = sigma;
nonbondedParams[2] = epsilon;
nonbondedParams[3] = softcoreLJLambda;
customNonbonded->addParticle( nonbondedParams );
}
return customNonbonded;
}
CustomBondForce* buildCustomBondForceForNonbondedExceptions( const NonbondedSoftcoreForce& nonbondedSoftcoreForce ){
CustomBondForce* customBond;
if( nonbondedSoftcoreForce.getNonbondedMethod() == NoCutoff ){
customBond = new CustomBondForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
} else {
customBond = new CustomBondForce("withinCutoff*(lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf));"
"withinCutoff=step(cutoff-r);"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
double cutoffDistance = nonbondedSoftcoreForce.getCutoffDistance();
double reactionFieldDielectric = nonbondedSoftcoreForce.getReactionFieldDielectric();
double eps2 = (reactionFieldDielectric - 1.0)/(2.0*reactionFieldDielectric+1.0);
double kValue = eps2/(cutoffDistance*cutoffDistance*cutoffDistance);
customBond->addGlobalParameter("krf", kValue );
double cValue = (1.0/cutoffDistance)*(3.0*reactionFieldDielectric)/(2.0*reactionFieldDielectric + 1.0);
customBond->addGlobalParameter("crf", cValue );
customBond->addGlobalParameter("cutoff", cutoffDistance );
}
customBond->addPerBondParameter("q");
customBond->addPerBondParameter("sigma");
customBond->addPerBondParameter("eps");
customBond->addPerBondParameter("lambda");
for( unsigned int ii = 0; ii < nonbondedSoftcoreForce.getNumExceptions(); ii++ ){
int particle1, particle2;
double chargeProd;
double sigma;
double epsilon;
double softcoreLJLambda;
nonbondedSoftcoreForce.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon, softcoreLJLambda );
vector<double> bondParams(4);
bondParams[0] = chargeProd;
bondParams[1] = sigma;
bondParams[2] = epsilon;
bondParams[3] = softcoreLJLambda;
customBond->addBond( particle1, particle2, bondParams );
}
return customBond;
}
#endif
/**
* Load plugins
*
* @param pluginDirectory plugin directory; if OPENMM_PLUGIN_DIR use ENV variable
* @param log logging file (optional -- may be NULL)
*
*/
static int loadPlugins( const std::string& pluginDirectory, std::vector<std::string>& loaded, FILE* log ){
const char* openmmPluginDirectory;
int envVariableIsSet = 0;
if( pluginDirectory.compare( "OPENMM_PLUGIN_DIR") == 0 ){
openmmPluginDirectory = getenv( "OPENMM_PLUGIN_DIR" );
} else {
openmmPluginDirectory = pluginDirectory.c_str();
}
try {
if( openmmPluginDirectory ){
envVariableIsSet = 1;
if( log ){
(void) fprintf( log, "openmmPluginDirectory=%s\n", openmmPluginDirectory );
(void) fflush( log );
}
loaded = Platform::loadPluginsFromDirectory( openmmPluginDirectory );
if( log ){
(void) fprintf( log, "\nLoaded following %u lib(s) from %s:\n", static_cast<unsigned int>(loaded.size()), openmmPluginDirectory ); (void) fflush( log );
for( unsigned int ii = 0; ii < loaded.size(); ii++ ){
(void) fprintf( log, " %s\n", loaded[ii].c_str() );
}
(void) fprintf( log, "\n" ); (void) fflush( log );
}
} else {
if( log && pluginDirectory == "OPENMM_PLUGIN_DIR" ){
(void) fprintf( log, "Env variable OPENMM_PLUGIN_DIR is not set.\n" );
(void) fflush( log );
}
}
} catch(const exception& e) {
(void) fprintf( log, "Exception: %s\n", e.what() );
(void) fflush( log );
}
return envVariableIsSet;
}
/**
* Set device id
*
* @param platform platform
* @param deviceId device id
* @param log logging file (optional -- may be NULL)
*
*/
static void setDeviceId( Platform& platform, int deviceId, FILE* log ){
std::stringstream deviceIdStr;
deviceIdStr << deviceId;
int wasSet = 0;
if( platform.getName().compare( "Cuda" ) == 0 ){
platform.setPropertyDefaultValue( "CudaDevice", deviceIdStr.str() );
wasSet = 1;
} else if( platform.getName().compare( "OpenCL" ) == 0 ){
platform.setPropertyDefaultValue( "OpenCLDeviceIndex", deviceIdStr.str());
wasSet = 1;
}
if( log && wasSet ){
(void) fprintf( log, "Set deviceId to %d\n", deviceId );
(void) fflush( log );
}
return;
}
/**
* Set device id
*
* @param platform platform
* @param deviceId device id
* @param log logging file (optional -- may be NULL)
*
*/
static void setDeviceIdUsingEnvVariable( Platform& platform, FILE* log ){
const char* deviceId = getenv( "GPU_DEVICE_ID" );
if( deviceId == NULL ){
return;
}
int wasSet = 0;
if( platform.getName().compare( "Cuda" ) == 0 ){
platform.setPropertyDefaultValue( "CudaDevice", deviceId );
wasSet = 1;
} else if( platform.getName().compare( "OpenCL" ) == 0 ){
platform.setPropertyDefaultValue( "OpenCLDeviceIndex", deviceId);
wasSet = 1;
}
if( log && wasSet ){
(void) fprintf( log, "Set deviceId to %s based on env variable GPU_DEVICE_ID setting.\n", deviceId );
(void) fflush( log );
}
return;
}
/**
* Get platform name
*
* @param platformId platformId( 0=Reference, 1=Cuda, 2=OpenCL)
* @param platformName output platform name
*
*/
static void getPlatformName( int platformId, std::string& platformName ){
switch( platformId ){
case Reference_OpenMMTest:
platformName = "Reference";
break;
case Cuda_OpenMMTest:
platformName = "Cuda";
break;
case OpenCL_OpenMMTest:
platformName = "OpenCL";
break;
default:
platformName = "NA";
break;
}
return;
}
/**
* Get lib name
*
* @param libPrefix lib prefix (lib or "")
* @param libSuffix lib suffix (.so, .dylib, .dll)
* @param baseName base name
*
* @return libname
*
*/
static std::string getLibName( const std::string& libPrefix, const std::string& libSuffix, const std::string& baseName ){
std::string fullName = libPrefix;
fullName.append( baseName );
fullName.append( libSuffix );
return fullName;
}
/**
* Get nonbonded method name
*
* @param nonbondedMethod nonbonded method flag
* @return nonbonded method name
*
*/
static std::string getNonbondedMethodName( int nonbondedMethod ){
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
return "NoCutoff";
case CutoffNonPeriodic_OpenMMTest:
return "CutoffNonPeriodic";
case CutoffPeriodic_OpenMMTest:
return "CutoffPeriodic";
case Ewald_OpenMMTest:
return "Ewald";
case PME_OpenMMTest:
return "PME";
default:
return "NA";
}
}
/**
* Check if required libs are available
*
* @param requiredLibs list of required libs
* @param loadedLibs list of available libs
* @param log optional logging reference
*
* @return 1 if all required libs are loaded; else 0
*
*/
static int checkRequiredLibsAreAvailable( const StringVector& requiredLibs, const StringVector& loadedLibs, FILE* log ){
unsigned int matchCount = 0;
for( unsigned int kk = 0; kk < requiredLibs.size(); kk++ ){
unsigned int match = 0;
for( unsigned int ii = 0; ii < loadedLibs.size() && match == 0; ii++ ){
if( loadedLibs[ii].compare( requiredLibs[kk] ) == 0 ){
match = 1;
}
}
if( log && !match ){
(void) fprintf( log, "Missing lib %s\n", requiredLibs[kk].c_str() );
}
matchCount += match;
}
int allPresent;
if( matchCount < requiredLibs.size() ){
allPresent = 0;
if( log ){
(void) fprintf( log, "Aborting tests due to missing libs.\n" );
}
} else {
allPresent = 1;
}
return allPresent;
}
/**
* Perform comparison of energies/forces for two systems
*
* @param system1 first system
* @param system2 second system
* @param platform1 first platform name (Reference, Cuda, OpenCL)
* @param platform2 second platform name (Reference, Cuda, OpenCL)
* @param positions positions
* @param inputArgumentMap arguments/flags (relativeTolerance, applyAssert, ...)
* @param idString id string
* @param log logging file (optional -- may be NULL)
*
*/
void runSystemComparisonTest( System& system1, System& system2,
const std::vector<Vec3>& positions, MapStringToDouble& inputArgumentMap,
const std::string& idString, const std::string& precision, FILE* log ){
int applyAssert = 0;
int platformId1 = 0;
int platformId2 = 0;
int deviceId1 = 0;
int deviceId2 = 0;
double relativeTolerance = 1.0e-04;
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
setIntFromMapStringToDouble( inputArgumentMap, "applyAssert", applyAssert ) ;
setIntFromMapStringToDouble( inputArgumentMap, "platformId1", platformId1 ) ;
setIntFromMapStringToDouble( inputArgumentMap, "platformId2", platformId2 ) ;
setIntFromMapStringToDouble( inputArgumentMap, "deviceId1", deviceId1 ) ;
setIntFromMapStringToDouble( inputArgumentMap, "deviceId2", deviceId2 ) ;
std::string platformName1;
std::string platformName2;
getPlatformName( platformId1, platformName1 );
getPlatformName( platformId2, platformName2 );
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
if( log ){
(void) fprintf( log, "System1: particles=%d forces=%d System2: particles=%d forces=%d\n",
system1.getNumParticles(), system1.getNumForces(),
system2.getNumParticles(), system2.getNumForces() );
(void) fprintf( log, "Positions=%u\n", static_cast<unsigned int>(positions.size()) );
(void) fprintf( log, "Platform1=%s Platform2=%s\n", platformName1.c_str(), platformName2.c_str() );
(void) fprintf( log, "deviceId1=%d deviceId2=%d\n", deviceId1, deviceId2 );
(void) fprintf( log, "relativeTolerance=%8.2e applyAssert=%d\n", relativeTolerance, applyAssert );
MapStringInt stringForceVector1;
MapStringInt stringForceVector2;
getStringForceMap( system1, stringForceVector1, log );
(void) fprintf( log, "Forces in system 1: [" );
for( MapStringIntCI ii = stringForceVector1.begin(); ii != stringForceVector1.end(); ii++ ){
(void) fprintf( log, " %s ", ii->first.c_str() );
}
getStringForceMap( system2, stringForceVector2, log );
(void) fprintf( log, "]\nForces in system 2: [" );
for( MapStringIntCI ii = stringForceVector2.begin(); ii != stringForceVector2.end(); ii++ ){
(void) fprintf( log, " %s ", ii->first.c_str() );
}
(void) fprintf( log, "]\n" );
}
if( system1.getNumParticles() != system2.getNumParticles() ){
std::stringstream msg;
msg << "Number of particles for systems to be compared are unequal: " << system1.getNumParticles() << " != " << system2.getNumParticles();
throw OpenMMException( msg.str() );
}
if( system1.getNumParticles() != static_cast<int>(positions.size()) ){
std::stringstream msg;
msg << "Number of particles for system does not equal size of position array: " << system1.getNumParticles() << " != " << positions.size();
throw OpenMMException( msg.str() );
}
#if TEST_PLATFORM == TEST_OPENCL_PLATFORM
ReferencePlatform platform1;
OpenCLPlatform platform2;
platform2.setPropertyDefaultValue("OpenCLPrecision", precision);
#elif TEST_PLATFORM == TEST_CUDA_PLATFORM
ReferencePlatform platform1;
CudaPlatform platform2;
platform2.setPropertyDefaultValue("CudaPrecision", precision);
#else
Platform& platform1 = Platform::getPlatformByName( platformName1 );
if( deviceId1 ){
setDeviceId( platform1, deviceId1, log );
}
setDeviceIdUsingEnvVariable( platform1, log );
Platform& platform2 = Platform::getPlatformByName( platformName2 );
if( deviceId2 ){
setDeviceId( platform2, deviceId2, log );
}
setDeviceIdUsingEnvVariable( platform2, log );
#endif
Context context1( system1, integrator1, platform1 );
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
Context context2( system2, integrator2, platform2 );
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
double energyDiff = 0.0;
if( fabs( state1.getPotentialEnergy() ) > 0.0 || fabs( state2.getPotentialEnergy()) > 0.0 ){
energyDiff = fabs( state1.getPotentialEnergy() - state2.getPotentialEnergy() )/( fabs( state1.getPotentialEnergy() ) + fabs( state2.getPotentialEnergy() ) );
}
if( log ){
DoubleVector stats;
compareForcesOfTwoStates( state1, state2, relativeTolerance, stats, log );
(void) fprintf( log, "%s %6d eDff=%15.7e fMx=%15.7e fAvg=%15.7e fMed=%15.7e eCd=%15.7e eRf=%15.7e mxFIdx=%d\n",
idString.c_str(), system1.getNumParticles(), energyDiff,
stats[1], stats[0], stats[3], state1.getPotentialEnergy(), state2.getPotentialEnergy(), static_cast<int>(stats[2]+0.0001));
(void) fflush( log );
}
if( applyAssert ){
ASSERT( energyDiff < relativeTolerance );
for( int ii = 0; ii < system1.getNumParticles(); ii++ ){
Vec3 f1 = state1.getForces()[ii];
Vec3 f2 = state2.getForces()[ii];
double f1N = sqrt( (f1[0]*f1[0]) + (f1[1]*f1[1]) + (f1[2]*f1[2]) );
double f2N = sqrt( (f2[0]*f2[0]) + (f2[1]*f2[1]) + (f2[2]*f2[2]) );
double diff = (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( f1N > 0.0 || f1N > 0.0 ){
diff = 2.0*sqrt( diff )/(f1N + f2N);
}
ASSERT( diff < relativeTolerance );
}
}
}
/**
* Serialize system
*
* @param system system to serialize
* @param serializeFileName file name for xml output
* @param log logging file (optional -- may be NULL)
*
*/
void serializeSystem( System& system, const std::string& serializeFileName, FILE* log ){
#ifdef OPENMM_SERIALIZE
//registerAmoebaSerializationProxies();
std::stringstream buffer;
XmlSerializer::serialize<System>(&system, "System", buffer);
FILE* filePtr = fopen( serializeFileName.c_str(), "w" );
if( filePtr == NULL ){
if( log ){
(void) fprintf( log, "Unable to open xml file %s\n", serializeFileName.c_str() );
return;
}
}
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
if( log ){
(void) fprintf( log, "Wrote system to xml file %s\n", serializeFileName.c_str() );
}
#endif
return;
}
/**
* Output vector of Vec3 to file
*
* @param positions system to serialize
* @param fileName file name for output
* @param log logging file (optional -- may be NULL)
*
*/
void serializeVectorOfVec3( const std::vector<Vec3>& positions, std::string fileName, FILE* log ){
#ifdef OPENMM_SERIALIZE
FILE* filePtr = fopen( fileName.c_str(), "w" );
if( filePtr == NULL ){
if( log ){
(void) fprintf( log, "Unable to open Vec3 file %s\n", fileName.c_str() );
return;
}
}
(void) fprintf( filePtr, "Positions %u\n", static_cast<unsigned int>(positions.size()) );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( filePtr, "%9u %17.10e %17.10e %17.10e\n", ii, positions[ii][0], positions[ii][1], positions[ii][2] );
}
(void) fclose( filePtr );
if( log ){
(void) fprintf( log, "Wrote to file %s\n", fileName.c_str() );
}
#endif
return;
}
/**
* Serialize system and positions
*
* @param system system to serialize
* @param positions positions to output
* @param baseFileName base file name for xml/txt output
* @param log logging file (optional -- may be NULL)
*
*/
void serializeSystemAndPositions( System& system, const std::vector<Vec3>& positions, const std::string& baseFileName, FILE* log ){
std::stringstream xmlfileName;
xmlfileName << baseFileName << ".xml";
serializeSystem( system, xmlfileName.str(), log );
std::stringstream posfileName;
posfileName << baseFileName << ".txt";
serializeVectorOfVec3( positions, posfileName.str(), log );
return;
}
void runTests( MapStringToDouble& inputArgumentMap, const std::string& precision, FILE* log ){
double lambda1 = 1.0;
double lambda2 = 1.0;
int nonbondedMethod = 0;
int numMolecules = 1;
int numParticlesPerMolecule = 2;
int useQuinticSpline = 1;
int applyAssert = 1;
int positionPlacementMethod = 0;
int serialize = 0;
double boxSize = 10.0;
double relativeTolerance = 1.0e-04;
double quinticLowerLimitFactor = 0.8;
double quinticUpperBornRadiusLimit = 2.0;
std::stringstream baseFileName;
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda1", lambda1 );
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda2", lambda2 );
setDoubleFromMapStringToDouble( inputArgumentMap, "boxSize", boxSize );
double cutoffDistance = boxSize*0.4;
setDoubleFromMapStringToDouble( inputArgumentMap, "cutoffDistance", cutoffDistance);
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
baseFileName << "Nb";
#if IMPLICIT_SOLVENT == TEST_GBVI
setDoubleFromMapStringToDouble( inputArgumentMap, "quinticLowerLimitFactor", quinticLowerLimitFactor );
setDoubleFromMapStringToDouble( inputArgumentMap, "quinticUpperBornRadiusLimit", quinticUpperBornRadiusLimit );
baseFileName << "Gbvi";
#endif
#if IMPLICIT_SOLVENT == TEST_OBC
baseFileName << "Obc";
#endif
setIntFromMapStringToDouble( inputArgumentMap, "positionPlacementMethod", positionPlacementMethod ) ;
setIntFromMapStringToDouble( inputArgumentMap, "nonbondedMethod", nonbondedMethod );
setIntFromMapStringToDouble( inputArgumentMap, "numMolecules", numMolecules );
setIntFromMapStringToDouble( inputArgumentMap, "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( inputArgumentMap, "serialize", serialize );
setIntFromMapStringToDouble( inputArgumentMap, "applyAssert", applyAssert );
double bondDistance = 0.1;
double minDistance = 0.1;
double cellSize = 2.0*bondDistance + minDistance;
double boxLength = cellSize*pow( static_cast<double>(numMolecules), 0.333333 );
if( positionPlacementMethod == 1 && boxLength > boxSize ){
boxSize = boxLength;
if( log ){
// (void) fprintf( log, "Updated box size: bL=%6.3f cell=%6.2e bond=%5.2f separation=%5.2f\n", boxLength, cellSize, bondDistance, minDistance );
}
}
if( nonbondedMethod >= 2 && cutoffDistance > boxSize*0.5 ){
cutoffDistance = boxSize*0.49;
}
int numParticles = numMolecules*numParticlesPerMolecule;
int includeGbvi = 1;
double reactionFieldDielectric = 80.0;
if( log ){
double particleDensity = static_cast<double>(numParticles)/(boxSize*boxSize*boxSize);
double particleCube = pow( particleDensity, (-1.0/3.0) );
(void) fprintf( log, "\n--------------------------------------------------------------------------------------\n" );
(void) fprintf( log, "Input arguments\n" );
(void) fflush( log );
//(void) fprintf( log, " includeGbvi %d\n", includeGbvi );
(void) fprintf( log, " nonbondedMethod %d\n", nonbondedMethod );
(void) fprintf( log, " numParticles %d\n", numParticles );
(void) fprintf( log, " numMolecules %d\n", numMolecules );
(void) fprintf( log, " numParticlesPerMolecule %d\n", numParticlesPerMolecule );
(void) fprintf( log, " positionPlacementMethod %d\n", positionPlacementMethod);
(void) fprintf( log, " boxSize %8.3f\n", boxSize );
(void) fprintf( log, " cutoffDistance %15.7e\n", cutoffDistance );
(void) fprintf( log, " reactionFieldDielectric %8.3f\n", reactionFieldDielectric );
#if IMPLICIT_SOLVENT == TEST_GBVI
(void) fprintf( log, " useQuinticSpline %d\n", useQuinticSpline );
(void) fprintf( log, " quinticLowerLimitFactor %8.3f\n", quinticLowerLimitFactor );
(void) fprintf( log, " quinticUpperBornRadiusLimit %8.3f\n", quinticUpperBornRadiusLimit );
#endif
#ifdef USE_SOFTCORE
(void) fprintf( log, " lambda1 %8.3f\n", lambda1 );
(void) fprintf( log, " lambda2 %8.3f\n", lambda2 );
#endif
(void) fprintf( log, " relativeTolerance %8.1e\n", relativeTolerance );
(void) fprintf( log, " particleDensity %8.2e\n", particleDensity );
(void) fprintf( log, " particleCube %8.2e\n", particleCube );
}
// Create two systems: one with GbviSoftcoreForce NonbondedSoftcoreForce forces, and one using a CustomNonbondedForce, CustomGBVI force to implement the same interaction.
System standardSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
double solventDielectric = 78.3; // 1.0 or 1.0e+10
double soluteDiecletric = 1.0;
std::vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
PositionGenerator positionGenerator( numMolecules, numParticlesPerMolecule, boxSize );
if( log ){
positionGenerator.setLog( log );
}
if( positionPlacementMethod == 1 ){
positionGenerator.setPositions( PositionGenerator::SimpleGrid, sfmt, positions );
} else {
positionGenerator.setBondDistance( 0.3 );
positionGenerator.setPositions( PositionGenerator::Random, sfmt, positions );
}
// show info on particle positions
if( log ){
int periodicBoundaryConditions = (nonbondedMethod > CutoffNonPeriodic_OpenMMTest) ? 1 : 0;
int showIndex = 5;
double distanceTolerance = 1.0e-04;
IntVector positionIndexVector;
positionIndexVector.push_back( 0 );
positionIndexVector.push_back( 5713 );
positionIndexVector.push_back( 6291 );
positionIndexVector.push_back( 3191 );
positionIndexVector.push_back( 3769 );
positionIndexVector.push_back( static_cast<int>(positions.size())-1 );
positionGenerator.showMinMaxDistances( positions, periodicBoundaryConditions, showIndex, positionIndexVector);
positionGenerator.showMinMaxDistances( positions, periodicBoundaryConditions, showIndex );
positionGenerator.showParticlesWithinDistance( positions, periodicBoundaryConditions, 5713, cutoffDistance, distanceTolerance );
positionGenerator.showParticlesWithinDistance( positions, periodicBoundaryConditions, 6291, cutoffDistance, distanceTolerance );
IntIntPairVector pairs;
pairs.push_back( IntIntPair( 5713, 6291 ) );
pairs.push_back( IntIntPair( 5713, 3191 ) );
pairs.push_back( IntIntPair( 5713, 3769 ) );
pairs.push_back( IntIntPair( 6291, 3191 ) );
positionGenerator.showDistances( pairs, positions );
}
const int numberOfParameters = 5;
std::vector<double> parameterLowerBound( numberOfParameters, 0.0 );
double fixedCharge = 0.1;
parameterLowerBound[ChargeIndex_OpenMMTest] = fixedCharge; // charge
parameterLowerBound[SigmaIndex_OpenMMTest] = 0.1; // sigma
parameterLowerBound[EpsIndex_OpenMMTest] = 0.5; // eps
parameterLowerBound[GammaIndex_OpenMMTest] = 0.1; // gamma
parameterLowerBound[LambdaIndex_OpenMMTest] = lambda1; // lambda
std::vector<double> parameterUpperBound( parameterLowerBound );
parameterUpperBound[ChargeIndex_OpenMMTest] = fixedCharge; // charge
parameterUpperBound[SigmaIndex_OpenMMTest] = 0.3; // sigma
parameterUpperBound[EpsIndex_OpenMMTest] = 40.0; // eps
parameterUpperBound[GammaIndex_OpenMMTest] = 40.0; // gamma
#if IMPLICIT_SOLVENT == TEST_OBC
parameterLowerBound[GammaIndex_OpenMMTest] = 0.1; // overlap factor
parameterUpperBound[GammaIndex_OpenMMTest] = 1.5;
#endif
std::vector<double> parameters( numberOfParameters );
VectorOfDoubleVectors parameterList;
std::vector< BondInfo_OpenMMTest > bonds;
double charge = fixedCharge;
for( int ii = 0; ii < numMolecules; ii++) {
charge *= -1.0;
double lambda = ii < (numMolecules/2) ? lambda1 : lambda2;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
parameters[ChargeIndex_OpenMMTest] = charge;
parameters[LambdaIndex_OpenMMTest] = lambda;
parameterList.push_back( parameters );
int baseParticleIndex = ii*numParticlesPerMolecule;
for( int jj = 1; jj < numParticlesPerMolecule; jj++) {
// alternate charges
charge *= -1.0;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
parameters[ChargeIndex_OpenMMTest] = charge;
parameters[LambdaIndex_OpenMMTest] = lambda;
parameterList.push_back( parameters );
double bondDistance = positionGenerator.getDistance( baseParticleIndex, baseParticleIndex+jj, positions );
bonds.push_back( BondInfo_OpenMMTest( baseParticleIndex, baseParticleIndex+jj, bondDistance ) );
}
// alternate charge if numParticlesPerMolecule is odd
if( (numParticlesPerMolecule % 2) ){
charge *= -1.0;
}
}
#ifdef USE_SOFTCORE
baseFileName << "Softcore";
NonbondedSoftcoreForce* nonbondedSoftcoreForce = getNonbondedSoftcoreForce( nonbondedMethod, cutoffDistance, reactionFieldDielectric,
parameterList, bonds, log);
#if IMPLICIT_SOLVENT == TEST_GBVI
GBVISoftcoreForce* gbviSoftcoreForce = getGBVISoftcoreForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
useQuinticSpline, quinticLowerLimitFactor, quinticUpperBornRadiusLimit,
solventDielectric, soluteDiecletric,
parameterList, bonds, log );
#endif
#if IMPLICIT_SOLVENT == TEST_OBC
GBSAOBCSoftcoreForce* gbviSoftcoreForce = getGBSAOBCSoftcoreForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
solventDielectric, soluteDiecletric,
parameterList, log );
#endif
#else
NonbondedForce* nonbondedSoftcoreForce = getNonbondedForce( nonbondedMethod, cutoffDistance, reactionFieldDielectric,
parameterList, bonds, log);
#if IMPLICIT_SOLVENT == TEST_GBVI
GBVIForce* gbviSoftcoreForce = getGBVIForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
useQuinticSpline, quinticLowerLimitFactor, quinticUpperBornRadiusLimit,
solventDielectric, soluteDiecletric, parameterList, bonds, log );
#endif
#if IMPLICIT_SOLVENT == TEST_OBC
GBSAOBCForce* gbviSoftcoreForce = getGBSAOBCForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
solventDielectric, soluteDiecletric, parameterList, log );
#endif
#endif
standardSystem.addForce(nonbondedSoftcoreForce);
#if IMPLICIT_SOLVENT > 0
if( includeGbvi ){
standardSystem.addForce(gbviSoftcoreForce);
}
#endif
// copy system and forces
System systemCopy;
copySystem( standardSystem, systemCopy, log );
// serialize
baseFileName << "_N" << positions.size();
baseFileName << "_Nb" << nonbondedMethod;
serializeSystemAndPositions( standardSystem, positions, baseFileName.str(), log);
// perform comparison
std::stringstream idString;
idString << "Nb " << nonbondedMethod << " l2 " << std::fixed << setprecision(2) << lambda2;
runSystemComparisonTest( standardSystem, systemCopy, positions, inputArgumentMap, idString.str(), precision, log );
}
int main(int argc, char* argv[]) {
std::string precision;
if (argc > 1)
precision = std::string(argv[1]);
else
precision = "single";
try {
#ifdef USE_SOFTCORE
registerFreeEnergyCudaKernelFactories( );
#endif
VectorOfMapStringToDouble vectorOfMapStringToDouble;
MapStringToDouble inputArgumentMap;
MapStringToDoubleVector generativeArgumentMaps;
//FILE* log = stderr;
FILE* log = NULL;
std::vector<std::string> loadedLibs;
int envVariableIsSet = loadPlugins( "OPENMM_PLUGIN_DIR", loadedLibs, log );
inputArgumentMap["platformId1"] = 0;
inputArgumentMap["platformId2"] = 1;
inputArgumentMap["deviceId1"] = 0;
inputArgumentMap["deviceId2"] = 0;
inputArgumentMap["lambda2"] = 1.0;
inputArgumentMap["nonbondedMethod"] = 0;
inputArgumentMap["numMolecules"] = 10;
inputArgumentMap["boxSize"] = 5.0;
inputArgumentMap["positionPlacementMethod"] = 1;
inputArgumentMap["cutoffDistance"] = 0.301*inputArgumentMap["boxSize"];
//inputArgumentMap["cutoffDistance"] = 1.0;
inputArgumentMap["relativeTolerance"] = 5.0e-04;
inputArgumentMap["applyAssert"] = 1;
inputArgumentMap["serialize"] = 1;
inputArgumentMap["numParticlesPerMolecule"] = 2;
#ifdef USE_SOFTCORE
DoubleVector lamda2;
lamda2.push_back( 1.0 );
lamda2.push_back( 0.5 );
lamda2.push_back( 0.0 );
if( lamda2.size() > 0 ){
generativeArgumentMaps["lambda2"] = lamda2;
inputArgumentMap["lambda2"] = lamda2[0];
}
#endif
DoubleVector numberOfMolecules;
//numberOfMolecules.push_back( 10 );
#if IMPLICIT_SOLVENT != TEST_NONBONDED
numberOfMolecules.push_back( 100 );
#endif
numberOfMolecules.push_back( 1000 );
numberOfMolecules.push_back( 2000 );
numberOfMolecules.push_back( 4000 );
//numberOfMolecules.push_back( 8000 );
if( numberOfMolecules.size() > 0 ){
generativeArgumentMaps["numMolecules"] = numberOfMolecules;
inputArgumentMap["numMolecules"] = numberOfMolecules[0];
}
DoubleVector nonbondedMethod;
nonbondedMethod.push_back( NoCutoff_OpenMMTest );
nonbondedMethod.push_back( CutoffNonPeriodic_OpenMMTest );
nonbondedMethod.push_back( CutoffPeriodic_OpenMMTest );
#if IMPLICIT_SOLVENT == TEST_NONBONDED
nonbondedMethod.push_back( Ewald_OpenMMTest );
nonbondedMethod.push_back( PME_OpenMMTest );
#endif
if( nonbondedMethod.size() > 0 ){
generativeArgumentMaps["nonbondedMethod"] = nonbondedMethod;
inputArgumentMap["nonbondedMethod"] = nonbondedMethod[0];
}
DoubleVector platformId2s;
#if TEST_PLATFORM == TEST_OPENCL_PLATFORM
platformId2s.push_back( OpenCL_OpenMMTest );
#elif TEST_PLATFORM == TEST_CUDA_PLATFORM
platformId2s.push_back( Cuda_OpenMMTest );
#else
platformId2s.push_back( Cuda_OpenMMTest );
#endif
// check that required libs are available for platform to be tested
// if unavailable, skip tests
std::string libPrefix = "lib";
std::string libSuffix = ".so";
#ifdef _MSC_VER
libPrefix = "";
libSuffix = ".dll";
#endif
#ifdef __APPLE__
libSuffix = ".dylib";
#endif
StringVector requiredLibs;
for( unsigned int kk = 0; kk < platformId2s.size(); kk++ ){
if( platformId2s[kk] == OpenCL_OpenMMTest ){
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMOpenCL" ) );
}
if( platformId2s[kk] == Cuda_OpenMMTest ){
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMCuda") );
#ifdef USE_SOFTCORE
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMFreeEnergy" ) );
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMFreeEnergyCuda" ) );
#endif
}
}
// if TEST_PLATFORM is not set, then check that required libs are available
#if TEST_PLATFORM != TEST_OPENCL_PLATFORM && TEST_PLATFORM != TEST_CUDA_PLATFORM
envVariableIsSet = checkRequiredLibsAreAvailable( requiredLibs, loadedLibs, log );
if( envVariableIsSet == 0 && log ){
(void) fprintf( log, "Aborting tests due to missing libs.\n" );
}
#else
// unit test path: force tests to run
envVariableIsSet = 1;
#endif
if( platformId2s.size() > 0 ){
generativeArgumentMaps["platformId2"] = platformId2s;
inputArgumentMap["platformId2"] = platformId2s[0];
}
vectorOfMapStringToDouble.push_back( inputArgumentMap );
generateInputArgumentMapsFromStringVectors( generativeArgumentMaps, vectorOfMapStringToDouble );
// modify relative tolerance for large systems
// case: Distance 433 669 1.5000001e+00 d2= 2.2500002e+00 w/ cutoff=1.500
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
int numMolecules = 0;
int numParticlesPerMolecule = 2;
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numMolecules", numMolecules );
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numParticlesPerMolecule", numParticlesPerMolecule );
if( numMolecules*numParticlesPerMolecule > 1000 ){
vectorOfMapStringToDouble[kk]["relativeTolerance"] = 6.0e-03;
}
}
if( log ){
MapStringToInt exclude;
exclude["lambda1"] = 1;
exclude["platformId1"] = 1;
exclude["platformId2"] = 1;
exclude["deviceId1"] = 1;
exclude["deviceId2"] = 1;
exclude["numParticlesPerMolecule"] = 1;
std::stringstream outputStream;
std::sort( vectorOfMapStringToDouble.begin(), vectorOfMapStringToDouble.end(), TestMapSortPredicate);
StringVector printOrder;
printOrder.push_back( "numMolecules" );
printOrder.push_back( "nonbondedMethod" );
printOrder.push_back( "lambda2" );
printOrder.push_back( "boxSize" );
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
streamArgumentMapOneLine( vectorOfMapStringToDouble[kk], exclude, printOrder, kk, outputStream );
}
(void) fprintf( log, "Initial argument maps: %u\n%s", static_cast<unsigned int>(vectorOfMapStringToDouble.size()), outputStream.str().c_str() );
}
// run tests
if( envVariableIsSet ){
int wasException = 0;
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size() && wasException < 3; kk++ ){
try {
runTests( vectorOfMapStringToDouble[kk], precision, log );
} catch(const exception& e) {
std::stringstream msg;
#if IMPLICIT_SOLVENT == TEST_NONBONDED
msg << "Nonbonded";
#elif IMPLICIT_SOLVENT == TEST_OBC
msg << "GBSAOBC";
#elif IMPLICIT_SOLVENT == TEST_GBVI
msg << "GBVI";
#endif
int numMolecules = 0;
int numParticlesPerMolecule = 0;
int nonbondedMethod = 0;
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numMolecules", numMolecules );
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "nonbondedMethod", nonbondedMethod);
msg << " test: system size=" << numMolecules*numParticlesPerMolecule << " nonbonded method=" << getNonbondedMethodName( nonbondedMethod );
msg << " exception: " << e.what() << endl;
// msg << "Note cases have been encountered for nonbonded methods with cutoffs where the error was due to particles being within 1.0e-05 of the cutoff." << endl;
cout << msg.str();
wasException += 1;
}
}
if( wasException ){
return 1;
}
}
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment