Commit 180c245d authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods

parent e3a1a218
#----------------------------------------------------
# OpenMM Brook Platform
#
# Creates OpenMM library, base name=OpenMM_Brook.
# Default libraries are shared & optimized. Variants
# are created for static (_static) and debug (_d).
#
# Windows:
# OpenMM_Brook[_d].dll
# OpenMM_Brook[_d].lib
# OpenMM_BROOK_static[_d].lib
# Unix:
# libOpenMM_Brook[_d].so
# libOpenMM_BROOK_static[_d].a
#----------------------------------------------------
SUBDIRS (tests)
# logging
SET(LOG TRUE)
IF(LOG)
SET(LOG_FILE "CMakeLog.txt" )
FILE( WRITE ${LOG_FILE} "In Brook Cmake\n")
ENDIF(LOG)
# The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
SET(OPENMM_SOURCE_SUBDIRS .)
# Collect up information about the version of the OpenMM library we're building
# and make it available to the code so it can be built into the binaries.
SET(OpenMM_BROOK_LIBRARY_NAME OpenMM_Brook)
SET(SHARED_BROOK_TARGET ${OpenMM_BROOK_LIBRARY_NAME})
SET(STATIC_BROOK_TARGET ${OpenMM_BROOK_LIBRARY_NAME}_static)
# Ensure that debug libraries have "_d" appended to their names.
# CMake gets this right on Windows automatically with this definition.
IF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
SET(CMAKE_DEBUG_POSTFIX "_d" CACHE INTERNAL "" FORCE)
ENDIF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
# But on Unix or Cygwin we have to add the suffix manually
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_BROOK_TARGET ${SHARED_BROOK_TARGET}_d)
SET(STATIC_BROOK_TARGET ${STATIC_BROOK_TARGET}_d)
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
# These are all the places to search for header files which are
# to be part of the API.
SET(API_INCLUDE_DIRS) # start empty
FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
# append
SET(API_INCLUDE_DIRS ${API_INCLUDE_DIRS}
${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include
${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include/internal)
ENDFOREACH(subdir)
# We'll need both *relative* path names, starting with their API_INCLUDE_DIRS,
# and absolute pathnames.
SET(API_REL_INCLUDE_FILES) # start these out empty
SET(API_ABS_INCLUDE_FILES)
FOREACH(dir ${API_INCLUDE_DIRS})
FILE(GLOB fullpaths ${dir}/*.h) # returns full pathnames
SET(API_ABS_INCLUDE_FILES ${API_ABS_INCLUDE_FILES} ${fullpaths})
FOREACH(pathname ${fullpaths})
GET_FILENAME_COMPONENT(filename ${pathname} NAME)
SET(API_REL_INCLUDE_FILES ${API_REL_INCLUDE_FILES} ${dir}/${filename})
ENDFOREACH(pathname)
ENDFOREACH(dir)
# ----------------------------------------------------------------------------
IF(LOG)
LOG_DIR( ${LOG_FILE} "API_ABS_INCLUDE_FILES" ${API_ABS_INCLUDE_FILES} )
LOG_DIR( ${LOG_FILE} "API_REL_INCLUDE_FILES" ${API_REL_INCLUDE_FILES} )
LOG_DIR( ${LOG_FILE} "CMAKE_CURRENT_SOURCE_DIR" ${CMAKE_CURRENT_SOURCE_DIR} )
ENDIF(LOG)
# ----------------------------------------------------------------------------
# collect cpp source files
SET(SOURCE_FILES) # empty
SET(SOURCE_INCLUDE_FILES)
# SET( CMAKE_CURRENT_SOURCE_DIR /home/friedrim/src/openmm/trunk/OpenMM/platforms/brook )
FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
FILE(GLOB src_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*/*.cpp)
FILE(GLOB incl_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.h)
SET(SOURCE_FILES ${SOURCE_FILES} ${src_files}) #append
SET(SOURCE_INCLUDE_FILES ${SOURCE_INCLUDE_FILES} ${incl_files})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
ENDFOREACH(subdir)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
# ----------------------------------------------------------------------------
IF(LOG)
LOG_DIR( ${LOG_FILE} "BROOK_SOURCE_FILES" ${SOURCE_FILES} )
LOG_DIR( ${LOG_FILE} "SOURCE_INCLUDE_FILES" ${SOURCE_INCLUDE_FILES} )
ENDIF(LOG)
# ----------------------------------------------------------------------------
# Brook setup
# MESSAGE("CMAKE_CURRENT_SOURCE_DIR=${CMAKE_CURRENT_SOURCE_DIR}" )
INCLUDE(${CMAKE_CURRENT_SOURCE_DIR}/brook-cmake/FindBrook.cmake)
INCLUDE_DIRECTORIES(${BROOK_INCLUDE_DIR})
LINK_DIRECTORIES(${${BROOK_brook_LIBRARY}})
# get *br files
FILE(GLOB BROOK_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/src/gpu/*.br)
FILE(GLOB BROOK_INCLUDE_FILES ${CMAKE_CURRENT_SOURCE_DIR}/src/gpu/*.h)
FILE( APPEND ${LOG_FILE} "BROOK_SRC_FILES=${BROOK_SRC_FILES}\n" )
# ----------------------------------------------------------------------------
IF(LOG)
LOG_DIR( ${LOG_FILE} "Brook src" ${BROOK_SRC_FILES} )
LOG_DIR( ${LOG_FILE} "Brook include" ${BROOK_INCLUDE_FILES} )
ENDIF(LOG)
# ----------------------------------------------------------------------------
# create Brook custom rules
SET(BROOK_CPP_FILES)
FOREACH(brookFile ${BROOK_SRC_FILES})
BROOK_FILE( ${brookFile} )
ENDFOREACH(brookFile)
# ----------------------------------------------------------------------------
IF(LOG)
LOG_DIR( ${LOG_FILE} "Brook cpp" ${BROOK_CPP_FILES} )
ENDIF(LOG)
# ----------------------------------------------------------------------------
# BROOK_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
ADD_LIBRARY(${SHARED_BROOK_TARGET} SHARED ${BROOK_CPP_FILES} ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} )
ADD_LIBRARY(${STATIC_BROOK_TARGET} STATIC ${BROOK_CPP_FILES} ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} )
# ----------------------------------------------------------------------------
IF(LOG)
FILE( APPEND ${LOG_FILE} "\nSTATIC_BROOK_TARGET=${STATIC_BROOK_TARGET} OPENMM_LIBRARY_NAME=${OPENMM_LIBRARY_NAME}\n" )
FILE( APPEND ${LOG_FILE} "SHARED_BROOK_TARGET=${SHARED_BROOK_TARGET} SHARED_TARGET=${SHARED_TARGET} STATIC_TARGET=${STATIC_TARGET}\n" )
FILE( APPEND ${LOG_FILE} "PROJECT_BINARY_DIR=${PROJECT_BINARY_DIR} BROOK_LIB_PATH=${BROOK_LIB_PATH}\n" )
ENDIF(LOG)
# ----------------------------------------------------------------------------
TARGET_LINK_LIBRARIES(${SHARED_BROOK_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME})
TARGET_LINK_LIBRARIES(${SHARED_BROOK_TARGET} debug brook_d optimized brook )
TARGET_LINK_LIBRARIES(${SHARED_BROOK_TARGET} ${SHARED_TARGET} )
LINK_DIRECTORIES(${SHARED_BROOK_TARGET} ${BROOK_LIB_PATH})
LINK_DIRECTORIES(${SHARED_BROOK_TARGET} ${PROJECT_BINARY_DIR})
TARGET_LINK_LIBRARIES(${STATIC_BROOK_TARGET} debug ${OPENMM_LIBRARY_NAME}_static_d optimized ${OPENMM_LIBRARY_NAME}_static)
TARGET_LINK_LIBRARIES(${STATIC_BROOK_TARGET} debug brook_d optimized brook)
TARGET_LINK_LIBRARIES(${STATIC_BROOK_TARGET} ${STATIC_TARGET})
LINK_DIRECTORIES(${STATIC_BROOK_TARGET} ${BROOK_LIB_PATH})
#---------------------------------------------------
......@@ -40,7 +40,8 @@ FIND_PROGRAM(BROOK_CC brcc
FIND_LIBRARY(BROOK_brook_LIBRARY
NAMES
${sub_lib}
brook
brook_d
PATHS
$ENV{BROOKDIR}/lib
$ENV{BROOKDIR}/bin
......@@ -67,6 +68,8 @@ IF(LOG)
FILE( APPEND ${LOG_FILE} "BROOK_INCLUDE_DIR=${BROOK_INCLUDE_DIR}\n" )
FILE( APPEND ${LOG_FILE} "BROOK_CC=${BROOK_CC}\n" )
FILE( APPEND ${LOG_FILE} "BROOK_brook_LIBRARY=${BROOK_brook_LIBRARY}\n" )
FILE( APPEND ${LOG_FILE} "BROOKROOT=<$ENV{BROOKROOT}>\n" )
FILE( APPEND ${LOG_FILE} "sub_lib=<${sub_lib}>\n" )
ENDIF(LOG)
# ----------------------------------------------------------------------------
......@@ -138,7 +141,8 @@ IF (BROOK_INCLUDE_DIR AND BROOK_brook_LIBRARY AND BROOK_CC)
SET(BROOK_CPP_FILES ${BROOK_CPP_FILES} ${OUTFILE})
ENDMACRO(BROOK_FILE)
ELSE (BROOK_INCLUDE_DIR AND BROOK_brook_LIBRARY AND BROOK_CC)
SET(BROOK_FOUND FALSE)
ENDIF (BROOK_INCLUDE_DIR AND BROOK_brook_LIBRARY AND BROOK_CC)
# Some verbosity
......@@ -151,7 +155,7 @@ ENDIF (NOT BROOK_FOUND)
# ----------------------------------------------------------------------------
IF(LOG)
FILE( APPEND ${LOG_FILE} "BROOK_FOUND=${BROOK_FOUND}\n" )
FILE( APPEND ${LOG_FILE} "BROOK_FOUND=<${BROOK_FOUND}>\n" )
FILE( APPEND ${LOG_FILE} "\nLeaving FindBrook.cmake\n" )
ENDIF(LOG)
# ----------------------------------------------------------------------------
......@@ -41,11 +41,116 @@
namespace OpenMM {
class KernelImpl;
class OpenMMBrookInterface;
class BrookPlatformData {
public:
/**
* Constructor
*
*/
BrookPlatformData( void );
/**
* Destructor
*
*/
~BrookPlatformData( void );
/**
* Get _removeCOM flag
*
* @return _removeCOM
*
*/
int removeCOM( void ) const;
/**
* Get _useOBC flag
*
* @return _useOBC
*
*/
int useOBC( void ) const;
/**
* Get _hasBonds flag
*
* @return _hasBonds
*
*/
int hasBonds( void ) const;
/**
* Get _hasAngles
*
* @return _hasAngles
*
*/
int hasAngles( void ) const;
/**
* Get _hasPeriodicTorsions
*
* @return _hasPeriodicTorsions
*
*/
int hasPeriodicTorsions( void ) const;
/**
* Get _hasRB
*
* @return _hasRB
*
*/
int hasRB( void ) const;
/**
* Get _hasNonbonded
*
* @return _hasNonbonded
*
*/
int hasNonbonded( void ) const;
/**
* Get _cmMotionFrequency
*
* @return _cmMotionFrequency
*
*/
int cmMotionFrequency( void ) const;
private:
int _removeCOM;
int _useOBC;
int _hasBonds;
int _hasAngles;
int _hasPeriodicTorsions;
int _hasRB;
int _hasNonbonded;
int _cmMotionFrequency;
};
/**
* This Platform subclass uses the Brook implementations of all the OpenMM kernels.
*/
class BrookPlatform : public Platform {
class OPENMM_EXPORT BrookPlatform : public Platform {
public:
......@@ -53,6 +158,7 @@ class BrookPlatform : public Platform {
static const int DefaultReturnValue = 0;
static const int DefaultErrorValue = -1;
/**
* BrookPlatform constructor
*
......@@ -63,13 +169,13 @@ class BrookPlatform : public Platform {
/**
* BrookPlatform constructor
*
* @param defaultAtomStreamWidth stream width
* @param defaultParticleStreamWidth stream width
* @param runtime Brook runtime (cal/cpu)
* @param log log file reference
*
*/
BrookPlatform( int atomStreamWdith, const std::string& runtime, FILE* log = NULL );
BrookPlatform( int particleStreamWdith, const std::string& runtime, FILE* log = NULL );
/**
* BrookPlatform destructor
......@@ -102,6 +208,12 @@ class BrookPlatform : public Platform {
bool supportsDoublePrecision( void ) const;
/**
* Return default Brook stream factory
*
* @return Brook stream factory
*/
const StreamFactory& getDefaultStreamFactory( void ) const;
/**
......@@ -130,7 +242,7 @@ class BrookPlatform : public Platform {
* @return default stream width
*/
int getAtomStreamWidth( void ) const;
int getParticleStreamWidth( void ) const;
/**
* Set log file reference
......@@ -163,6 +275,19 @@ class BrookPlatform : public Platform {
FILE* getLog( void ) const;
/**
* This is called whenever a new OpenMMContext is created. It gives the Platform a chance to initialize
* the context and store platform-specific data in it.
*/
void contextCreated( OpenMMContextImpl& context ) const;
/**
* This is called whenever an OpenMMContext is deleted. It gives the Platform a chance to clean up
* any platform-specific data that was stored in it.
*/
void contextDestroyed( OpenMMContextImpl& context ) const;
private:
// log file reference
......@@ -173,11 +298,11 @@ class BrookPlatform : public Platform {
// default stream width
static const int DefaultAtomStreamWidth = 32;
static const int DefaultParticleStreamWidth = 32;
// atom streamwidth
// particle streamwidth
int _atomStreamWidth;
int _particleStreamWidth;
// Brook runtime
......@@ -202,6 +327,7 @@ class BrookPlatform : public Platform {
};
} // namespace OpenMM
#endif /*OPENMM_BROOKPLATFORM_H_*/
......@@ -50,9 +50,9 @@ class BrookStreamFactory : public StreamFactory {
// 'external' streams
static const std::string AtomPositions;
static const std::string AtomVelocities;
static const std::string AtomForces;
static const std::string ParticlePositions;
static const std::string ParticleVelocities;
static const std::string ParticleForces;
/**
* Create StreamImpl
......@@ -69,27 +69,27 @@ class BrookStreamFactory : public StreamFactory {
StreamImpl* createStreamImpl( std::string name, int size, Stream::DataType type, const Platform& platform, OpenMMContextImpl& context ) const;
/**
* Get atom stream width
* Get particle stream width
*
* @return atom stream width
* @return particle stream width
*
*
*/
int getDefaultAtomStreamWidth( void ) const;
int getDefaultParticleStreamWidth( void ) const;
/**
* Set atom stream width
* Set particle stream width
*
* @param atomStreamWidth atom stream width
* @param particleStreamWidth particle stream width
*
* @return DefaultReturnValue
*
* @throw OpenMMException if atomStreamWidth < 1
* @throw OpenMMException if particleStreamWidth < 1
*
*/
int setDefaultAtomStreamWidth( int atomStreamWidth );
int setDefaultParticleStreamWidth( int particleStreamWidth );
/**
* Get randomNumber stream width
......@@ -159,17 +159,17 @@ class BrookStreamFactory : public StreamFactory {
private:
static const int DefaultStreamAtomWidth = 32;
static const int DefaultStreamParticleWidth = 32;
static const int DefaultStreamRandomNumberWidth = 32;
static const int DefaultStreamRandomNumberSize = 1024;
static const int DefaultStreamRandomNumberWidth = 1024;
static const int DefaultStreamRandomNumberSize = 1024*1024;
static const double DefaultDangleValue;
static const int DefaultReturnValue = 0;
static const int ErrorReturnValue = -1;
int _defaultAtomStreamWidth;
int _defaultParticleStreamWidth;
int _defaultStreamRandomNumberWidth;
int _defaultStreamRandomNumberSize;
......
......@@ -39,6 +39,7 @@ using namespace std;
/**
* BrookBondParameters constructor
*
* @param bondName bond name
* @param numberOfParticlesInBond no. of particles in each bond
* @param numberOfParametersInBond no. of parameters in each bond
* @param numberOfBonds no. of bonds
......@@ -189,9 +190,6 @@ int BrookBondParameters::setBond( int bondIndex, int* particleIndices, double* b
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookBondParameters::setBond";
// ---------------------------------------------------------------------------------------
FILE* log = getLog();
// ---------------------------------------------------------------------------------------
......@@ -229,7 +227,7 @@ int BrookBondParameters::setBond( int bondIndex, int* particleIndices, double* b
}
/*
* Get contents of object
* Format line
*
* @param tab tab
* @param description description
......@@ -246,7 +244,6 @@ std::string BrookBondParameters::_getLine( const std::string& tab,
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookStreamInternal::_getLine";
static const unsigned int MAX_LINE_CHARS = 256;
char line[MAX_LINE_CHARS];
......@@ -268,7 +265,6 @@ std::string BrookBondParameters::_getLine( const std::string& tab,
/*
* Get contents of object
*
*
* @param level level of dump
*
* @return string containing contents
......@@ -311,8 +307,6 @@ std::string BrookBondParameters::getContentsString( int level ) const {
(void) LOCAL_SPRINTF( value, "%d", getNumberOfParametersInBond() );
message << _getLine( tab, "Parameters/bond:", value );
//(void) fprintf( getLog(), "%s %s QQQ1\n", methodName.c_str(), message.str().c_str() ); (void) fflush( getLog() );
message << "Bonds:" << std::endl;
for( int ii = 0; ii < getNumberOfBonds(); ii++ ){
const static size_t descriptionSz = 1024;
......@@ -364,4 +358,3 @@ std::string BrookBondParameters::getContentsString( int level ) const {
return message.str();
}
......@@ -48,6 +48,17 @@ class BrookBondParameters {
static const int DefaultReturnValue = 0;
static const int ErrorReturnValue = -1;
/**
* BrookBondParameters constructor
*
* @param bondName bond name
* @param numberOfParticlesInBond no. of particles in each bond
* @param numberOfParametersInBond no. of parameters in each bond
* @param numberOfBonds no. of bonds
* @param log optional log reference
*
*/
BrookBondParameters( std::string bondName, int numberOfParticlesInBond, int numberOfParametersInBond, int numberOfBonds, FILE* log );
~BrookBondParameters();
......
......@@ -568,7 +568,7 @@ int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log, int maxPrint ){
(void) LOCAL_SPRINTF( value, "%6d ", ii );
message << value << " [ ";
for( unsigned int jj = 0; jj < width; jj++ ){
for( int jj = 0; jj < width; jj++ ){
(void) LOCAL_SPRINTF( value, "%16.7e ", dataArray[index++] );
message << value;
}
......
......@@ -136,8 +136,8 @@ BrookOpenMMFloat BrookLangevinDynamics::getTau( void ) const {
*/
BrookOpenMMFloat BrookLangevinDynamics::getFriction( void ) const {
static const BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
static const BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0;
static const BrookOpenMMFloat zero = static_cast<BrookOpenMMFloat>( 0.0 );
static const BrookOpenMMFloat one = static_cast<BrookOpenMMFloat>( 1.0 );
return ( (_tau == zero) ? zero : (one/_tau) );
}
......@@ -187,7 +187,7 @@ int BrookLangevinDynamics::_setTau( BrookOpenMMFloat tau ){
*/
int BrookLangevinDynamics::_setFriction( BrookOpenMMFloat friction ){
_tau = (BrookOpenMMFloat) ( (friction != 0.0) ? 1.0/friction : 0.0);
_tau = static_cast<BrookOpenMMFloat>( (friction != 0.0) ? 1.0/friction : 0.0);
return DefaultReturnValue;
}
......@@ -224,7 +224,7 @@ int BrookLangevinDynamics::_setStepSize( BrookOpenMMFloat stepSize ){
*
* @return DefaultReturnValue
*
* @throw if tau too small
* @throw OpenMMException if tau too small
*
*/
......@@ -263,7 +263,7 @@ int BrookLangevinDynamics::_updateDerivedParameters( void ){
_derivedParameters[EM] = EXP( -_derivedParameters[GDT] );
_derivedParameters[EP] = EXP( _derivedParameters[GDT] );
if( _derivedParameters[GDT] >= (BrookOpenMMFloat) 0.1 ){
if( _derivedParameters[GDT] >= static_cast<BrookOpenMMFloat>( 0.1 ) ){
BrookOpenMMFloat term1 = _derivedParameters[EPH] - one;
term1 *= term1;
......@@ -278,20 +278,20 @@ int BrookLangevinDynamics::_updateDerivedParameters( void ){
BrookOpenMMFloat term2 = term1*term1;
BrookOpenMMFloat term4 = term2*term2;
BrookOpenMMFloat third = (BrookOpenMMFloat) ( 1.0/3.0 );
BrookOpenMMFloat o7_9 = (BrookOpenMMFloat) ( 7.0/9.0 );
BrookOpenMMFloat o1_12 = (BrookOpenMMFloat) ( 1.0/12.0 );
BrookOpenMMFloat o17_90 = (BrookOpenMMFloat) ( 17.0/90.0 );
BrookOpenMMFloat o7_30 = (BrookOpenMMFloat) ( 7.0/30.0 );
BrookOpenMMFloat o31_1260 = (BrookOpenMMFloat) ( 31.0/1260.0 );
BrookOpenMMFloat o_360 = (BrookOpenMMFloat) ( 1.0/360.0 );
BrookOpenMMFloat third = static_cast<BrookOpenMMFloat>( ( 1.0/3.0 ) );
BrookOpenMMFloat o7_9 = static_cast<BrookOpenMMFloat>( ( 7.0/9.0 ) );
BrookOpenMMFloat o1_12 = static_cast<BrookOpenMMFloat>( ( 1.0/12.0 ) );
BrookOpenMMFloat o17_90 = static_cast<BrookOpenMMFloat>( ( 17.0/90.0 ) );
BrookOpenMMFloat o7_30 = static_cast<BrookOpenMMFloat>( ( 7.0/30.0 ) );
BrookOpenMMFloat o31_1260 = static_cast<BrookOpenMMFloat>( ( 31.0/1260.0 ) );
BrookOpenMMFloat o_360 = static_cast<BrookOpenMMFloat>( ( 1.0/360.0 ) );
_derivedParameters[B] = term4*( third + term1*( third + term1*( o17_90 + term1*o7_9 )));
_derivedParameters[C] = term2*term1*( two*third + term1*( -half + term1*( o7_30 + term1*(-o1_12 + term1*o31_1260 ))));
_derivedParameters[D] = term2*( -one + term2*(-o1_12 - term2*o_360));
}
BrookOpenMMFloat kT = ((BrookOpenMMFloat) BOLTZ)*temperature;
BrookOpenMMFloat kT = static_cast<BrookOpenMMFloat>( BOLTZ )*temperature;
_derivedParameters[V] = SQRT( kT*( one - _derivedParameters[EM]) );
_derivedParameters[X] = tau*SQRT( kT*_derivedParameters[C] );
......@@ -362,8 +362,7 @@ int BrookLangevinDynamics::updateParameters( double temperature, double friction
return DefaultReturnValue;
};
}
/**
*
......@@ -687,6 +686,8 @@ int BrookLangevinDynamics::_updateSdStreams( void ){
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize();
// create and initialize sdpc streams
BrookOpenMMFloat* sdpc[2];
for( int ii = 0; ii < 2; ii++ ){
sdpc[ii] = new BrookOpenMMFloat[2*sdParticleStreamSize];
......@@ -775,8 +776,8 @@ int BrookLangevinDynamics::_setInverseSqrtMasses( const std::vector<double>& mas
//static const std::string methodName = "BrookLangevinDynamics::_setInverseSqrtMasses";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0;
BrookOpenMMFloat zero = static_cast<BrookOpenMMFloat>( 0.0 );
BrookOpenMMFloat one = static_cast<BrookOpenMMFloat>( 1.0 );
// ---------------------------------------------------------------------------------------
......@@ -861,13 +862,6 @@ float BrookLangevinDynamics::getTemperature( BrookStreamInternal* velocities, Br
int numberOfParticles = getNumberOfParticles();
for( int ii = 0; ii < numberOfParticles; ii++, index += 3 ){
ke += (velocitiesI[index]*velocitiesI[index] + velocitiesI[index+1]*velocitiesI[index+1] + velocitiesI[index+2]*velocitiesI[index+2] )/inverseMassStreamI[ii];
/*
if( ii < 3 || ii >= (numberOfParticles - 3) ){
(void) fprintf( stderr, "%s %6d m=%14.5e v[%14.5e %14.5e %14.5e]\n", methodName.c_str(), ii, 1.0/inverseMassStreamI[ii],
velocitiesI[index], velocitiesI[index+1], velocitiesI[index+2] );
} */
}
int degreesOfFreedom = 3*getNumberOfParticles() - numberOfConstraints;
......@@ -879,7 +873,7 @@ if( ii < 3 || ii >= (numberOfParticles - 3) ){
}
/**
* Remove velocity com
* Remove velocity com (diagnostics)
*
* @param velocities velocities
* @param inverseMassStream inverse masses
......@@ -895,8 +889,6 @@ int BrookLangevinDynamics::removeCom( BrookStreamInternal* velocities, BrookFloa
// ---------------------------------------------------------------------------------------
// (void) fprintf( stderr, "%s\n", methodName.c_str() ); fflush( stderr );
void* dataArrayV = velocities->getData( 1 );
float* velocitiesI = (float*) dataArrayV;
......@@ -907,7 +899,6 @@ int BrookLangevinDynamics::removeCom( BrookStreamInternal* velocities, BrookFloa
float com[3] = { 0.0f, 0.0f, 0.0f };
int index = 0;
//(void) fprintf( stderr, "%s strm %d %d\n", methodName.c_str(), velocities->getStreamSize(),velocities->getWidth() ); fflush( stderr );
for( int ii = 0; ii < getNumberOfParticles(); ii++ ){
float mass = 1.0f/inverseMassStreamI[ii];
totalMass += mass;
......@@ -931,13 +922,6 @@ int BrookLangevinDynamics::removeCom( BrookStreamInternal* velocities, BrookFloa
index += 3;
}
/*
(void) fprintf( stderr, "%s com[%14.5e %14.5e %14.5e]\n", methodName.c_str(), com[0], com[1], com[2] );
for( int ii = 0; ii < velocities->getStreamSize()*3; ii += 3 ){
(void) fprintf( stderr, "%s %d newVelocities[%14.5e %14.5e %14.5e]\n", methodName.c_str(), ii/3, newVelocities[ii], newVelocities[ii+1], newVelocities[ii+2] );
}
*/
velocities->loadFromArray( newVelocities );
dataArrayV = velocities->getData( 1 );
......@@ -973,12 +957,12 @@ int BrookLangevinDynamics::resetVelocities( BrookStreamInternal* velocities ) co
// ---------------------------------------------------------------------------------------
(void) fprintf( stderr, "%s\n", methodName.c_str() ); fflush( stderr );
// resset v to determinisitic values
// reset velocities to determinisitic values
// note use of double instead of float for the load array
double* newVelocities = new double[velocities->getStreamSize()*velocities->getWidth()];
memset( newVelocities, 0, sizeof( double )*velocities->getStreamSize()*velocities->getWidth() );
for( int ii = 1; ii <= 3*getNumberOfParticles(); ii++ ){
int jj = ii % 10;
double sign = jj % 2 ? 0.1 : -0.1;
......@@ -1137,6 +1121,8 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
}
}
// diagnostics
if( 0 && _internalStepCount == 1 ){
//resetVelocities( velocityStream.getBrookStreamInternal() );
std::string velocityFileName = "kupdate_sd1_strV.in";
......@@ -1149,6 +1135,8 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
getSD2XStream()->loadStreamGivenFileName( sd2xName );
}
// more diagnostics
if( 0 && (_internalStepCount % 10) == 0 ){
FILE* log1 = stderr;
(void) fprintf( log1, "\nVelocityStream %d XX\n", _internalStepCount ); fflush( log1 );
......@@ -1429,8 +1417,9 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamInternal();
brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, log, 0.0001f );
(void) fprintf( log, "\nPositionStream output sd shake2 step=%d\n", _internalStepCount );
std::string violationString;
brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, violationString, 0.0001f );
(void) fprintf( log, "\nPositionStream output sd shake2 step=%d %s\n", _internalStepCount, violationString.c_str() );
brookStreamInternalPos->printToFile( log );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamInternal();
......@@ -1466,16 +1455,12 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
if( printOn ){
(void) fprintf( log, "\n%s Pre ksetStr3 (no constraints)", methodName.c_str() );
// (void) fprintf( log, "\nXPrimeStream\n" );
// getXPrimeStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamInternal();
(void) fprintf( log, "\nPositionStream\n" );
brookStreamInternalPos->printToFile( log );
}
kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream(), positionStream.getBrookStream() );
//ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
// diagnostics
......@@ -1493,9 +1478,69 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
}
}
if( (_internalStepCount % 400) == 0 ){
// diagnostics
if( (_internalStepCount % 10) == 0 ){
FILE* log1 = stderr;
BrookStreamInternal* brookStreamInternalPos = velocityStream.getBrookStreamInternal();
float epsilon = 1.0e-01f;
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamInternal();
BrookStreamInternal* brookStreamInternalVel = velocityStream.getBrookStreamInternal();
BrookStreamInternal* brookStreamInternalFrc = forceStream.getBrookStreamInternal();
// check for nan and infinities
int coordinateNans = brookStreamInternalPos->checkForNans( );
int velocityNans = brookStreamInternalVel->checkForNans( );
int forceNans = brookStreamInternalFrc->checkForNans( );
int abort = abs( coordinateNans ) + abs( velocityNans ) + abs( forceNans );
// Shake violations
std::string violationString;
int constraintViolations = brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, violationString, 0.0001f );
abort += abs( constraintViolations );
// check T consistent w/ specified value
float temperature = getTemperature( brookStreamInternalVel, getInverseMassStream(), brookShakeAlgorithm.getNumberOfConstraints() );
if( fabsf( temperature - getTemperature() ) > 2.0f*temperature ){
abort++;
}
// force sums ~ 0?
std::vector<float> sums;
brookStreamInternalFrc->sumColumns( sums );
// check if should abort
(void) fprintf( log1, "%d T=%.3f Nans: x=%d v=%d f=%d ", _internalStepCount, temperature, coordinateNans, velocityNans, forceNans );
(void) fprintf( log1, " Fsum[" );
for( int ii = 0; ii < 3; ii++ ){
if( fabsf( sums[ii] ) > epsilon ){
abort++;
}
(void) fprintf( log1, "%12.4e ", sums[ii] );
}
(void) fprintf( log1, "] %s abort=%d\n", violationString.c_str(), abort );
(void) fflush( log1 );
if( abort ){
int nans[2];
nans[0] = brookRandomNumberGenerator.getRandomNumberStream( 0 )->checkForNans();
nans[1] = brookRandomNumberGenerator.getRandomNumberStream( 1 )->checkForNans();
(void) fprintf( log1, "Aborting: Nans rng: active index=%d %d %d\n", brookRandomNumberGenerator.getRvStreamIndex(), nans[0], nans[1] );
brookStreamInternalPos->printToFile( log );
brookStreamInternalVel->printToFile( log );
brookStreamInternalFrc->printToFile( log );
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log );
exit(1);
}
/*
std::vector<std::vector<double> > velocityStatistics;
brookStreamInternalPos->getStatistics( velocityStatistics, getNumberOfParticles() );
......@@ -1504,9 +1549,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
std::string stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
(void) fprintf( log1, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
*/
float temperature = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm.getNumberOfConstraints() );
//removeCom( brookStreamInternalPos, getInverseMassStream() );
(void) fprintf( log1, "\nVelocityStream %d T=%.3f\n", _internalStepCount, temperature );
}
return DefaultReturnValue;
......
......@@ -76,7 +76,7 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
_auxiliaryStreams[ii] = NULL;
}
// set randomNumber seed
// set randomNumber seed & generator
_randomNumberSeed = 1393;
//_randomNumberGenerator = Mersenne;
......@@ -391,14 +391,6 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
methodName.c_str(), stateInitialized, reseed, state[0], state[1], state[2], state[3] );
(void) fflush( log );
}
/*
state[0] = 9578;
state[1] = 29245;
state[2] = 16266;
state[3] = 27587;
*/
}
stateInitialized++;
......@@ -435,7 +427,6 @@ state[3] = 27587;
return DefaultReturnValue;
}
/**
* Load random number streams using Mersenne algorithm
*
......@@ -708,7 +699,7 @@ int BrookRandomNumberGenerator::_shuffleGVStreams( void ){
int numberOfRvStreams = getNumberOfRandomNumberStreams();
for( int ii = 0; ii < numberOfRvStreams - 1; ii++ ){
for( int ii = 0; ii < (numberOfRvStreams - 1); ii++ ){
kpermute_vectors( (float) getRandomNumberStreamWidth(),
_getShuffleStream()->getBrookStream(),
getRandomNumberStream( ii + 1 )->getBrookStream(),
......
......@@ -618,14 +618,14 @@ int BrookShakeAlgorithm::setup( const std::vector<double>& masses, const std::ve
* Check constraints
*
* @param positions atom positions
* @param log file to print to (can be NULL)
* @param outputString output message
* @param tolerance tolerance to compare (if < 0, then use algorithm tolerance
*
* @return number of errors
*
*/
int BrookShakeAlgorithm::checkConstraints( BrookStreamInternal* positions, FILE* log, float tolerance ) const {
int BrookShakeAlgorithm::checkConstraints( BrookStreamInternal* positions, std::string& outputString, float tolerance ) const {
// ---------------------------------------------------------------------------------------
......@@ -701,16 +701,42 @@ if( !(ii % 2 ) )fprintf( log, "]\n", posArray[ii] );
// report findings
if( errors && log ){
(void) fprintf( log, "Shake errors=%d tolerance=%.3e max diff=%.3e atoms[%d %d]", errors, tolerance, maxDiff, maxDiffCentralIndex, maxDiffPeripheralIndex );
std::stringstream outputMessage;
char text[1024];
if( errors ){
#ifdef WIN32
(void) sprintf_s( text, 1024, "Shake errors=%d tol=%.3e mxDff=%.3e atoms[%d %d]", errors, tolerance, maxDiff, maxDiffCentralIndex, maxDiffPeripheralIndex );
#else
(void) sprintf( text, "Shake errors=%d tol=%.3e mxDff=%.3e atoms[%d %d]", errors, tolerance, maxDiff, maxDiffCentralIndex, maxDiffPeripheralIndex );
#endif
outputMessage << text;
if( errors >= maxErrorToPrint ){
(void) fprintf( log, " only printing first %d errors", maxErrorToPrint );
#ifdef WIN32
(void) sprintf_s( text, 1024, " only printing first %d errors", maxErrorToPrint );
#else
(void) sprintf( text, " only printing first %d errors", maxErrorToPrint );
#endif
outputMessage << text;
}
(void) fprintf( log, "\n%s", message.str().c_str() );
} else if( log ){
(void) fprintf( log, "Shake no errors: tolerance=%.3e max diff=%.3e", tolerance, maxDiff );
outputMessage << message.str();
} else {
#ifdef WIN32
(void) sprintf_s( text, 1024, "Shake no errors: tol=%.3e mxDff=%.3e", tolerance, maxDiff );
#else
(void) sprintf( text, "Shake no errors: tol=%.3e mxDff=%.3e", tolerance, maxDiff );
#endif
outputMessage << text;
}
outputString = outputMessage.str();
return errors;
}
......
......@@ -301,14 +301,14 @@ class BrookShakeAlgorithm : public BrookCommon {
* Check constraints
*
* @param positions atom positions
* @param log file to print to (can be NULL)
* @param outputString output message
* @param tolerance tolerance to compare (if < 0, then use algorithm tolerance
*
* @return number of errors
*
*/
int checkConstraints( BrookStreamInternal* positions, FILE* log, float tolerance ) const;
int checkConstraints( BrookStreamInternal* positions, std::string& outputString, float tolerance ) const;
private:
......
......@@ -33,6 +33,15 @@
#include "OpenMMException.h"
#include "BrookStreamInternal.h"
#ifdef _WIN32
#include <float.h>
#define isnan _isnan
#define isinf !_finite
#endif
using namespace OpenMM;
using namespace std;
......@@ -288,6 +297,7 @@ std::string BrookStreamInternal::_getLine( const std::string& tab,
* Print contents of object to file
*
* @param log file to print to
* @param maxPrint max values to print; if < 0, then all values printed; default value is -1
*
* @return DefaultReturnValue
*
......@@ -304,6 +314,7 @@ int BrookStreamInternal::printToFile( FILE* log, int maxPrint ){
if( log == NULL ){
log = stderr;
}
std::string contents = getContentsString();
(void) fprintf( log, "%s\n", contents.c_str() );
......@@ -364,9 +375,12 @@ const std::string BrookStreamInternal::getContentsString( int level ) const {
}
/*
* Get stats
* Get stats (virtual method -- only BrookStreamInternalFloat implemented for now)
*
* @param statistics output vector of stats
* @param maxScan number of points to use in computing stats
*
* @return statistics vector
* @return 0
*
* */
......@@ -379,6 +393,7 @@ int BrookStreamInternal::getStatistics( std::vector<std::vector<double> >& stati
*
* @param tag id tag
* @param statistics stat vector
*
* @return stat string
*
* */
......@@ -457,7 +472,7 @@ std::string BrookStreamInternal::printStatistics( std::string tag, std::vector<s
*
* @return DefaultReturnValue
*
* */
**/
int BrookStreamInternal::printStreamsToFile( std::string fileName, std::vector<BrookStreamInternal*>& streams ){
......@@ -558,6 +573,15 @@ int BrookStreamInternal::printStreamsToFile( std::string fileName, std::vector<B
}
/*
* Load data into stream from file
*
* @param fileName file name
*
* @return DefaultReturnValue
*
**/
typedef struct {
unsigned int type;
unsigned int dimensions;
......@@ -667,3 +691,69 @@ int BrookStreamInternal::loadStreamGivenFileName( std::string& filename ){
return DefaultReturnValue;
}
/*
* Check for NANs
*
* @return number of Nans found
*
**/
int BrookStreamInternal::checkForNans( void ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookStreamInternal::checkForNans";
// ---------------------------------------------------------------------------------------
int numberOfNans = 0;
void* dataArrayV = getData( 1 );
float* array = static_cast<float*>( dataArrayV );
int width = getWidth();
int streamSize = getSize();
for( int ii = 0; ii < width*streamSize; ii++ ){
if( isnan( array[ii] ) || isinf( array[ii] ) ){
numberOfNans++;
}
}
return numberOfNans;
}
/*
* Sum columns
*
* @param sums output vector of column sums
*
* @return DefaultReturnValue
*
**/
int BrookStreamInternal::sumColumns( std::vector<float>& sums ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookStreamInternal::sumColumns";
// ---------------------------------------------------------------------------------------
void* dataArrayV = getData( 1 );
float* array = static_cast<float*>( dataArrayV );
int width = getWidth();
sums.resize( width );
for( int ii = 0; ii < width; ii++ ){
sums[ii] = 0.0f;
}
int streamSize = getSize();
int index = 0;
for( int ii = 0; ii < streamSize; ii++ ){
for( int jj = 0; jj < width; jj++ ){
sums[jj] += array[index++];
}
}
return DefaultReturnValue;
}
......@@ -226,6 +226,7 @@ class BrookStreamInternal {
* Print to file
*
* @param log log file
* @param maxPrint max values to print; if < 0, then all values printed; default value is -1
*
* @return DefaultReturnValue
*
......@@ -249,11 +250,14 @@ class BrookStreamInternal {
/*
* Get stats
*
* @param statistics output vector of stats
* @param maxScan number of points to use in computing stats
*
* @return statistics vector
*
* */
virtual int getStatistics( std::vector<std::vector<double>>&, int maxScan );
virtual int getStatistics( std::vector<std::vector<double>>& statVector, int maxScan );
/*
* Get stat string
......@@ -290,6 +294,27 @@ class BrookStreamInternal {
static int printStreamsToFile( std::string fileName, std::vector<BrookStreamInternal*>& streams );
/*
* Check for NANs
*
* @return number of Nans found
*
**/
int checkForNans( void );
/*
* Sum columns
*
* @param sums output vector of column sums
*
* @return DefaultReturnValue
*
**/
int sumColumns( std::vector<float>& sums );
protected:
std::string _name;
......
......@@ -51,7 +51,7 @@ BrookVelocityCenterOfMassRemoval::BrookVelocityCenterOfMassRemoval( ){
//static const std::string methodName = "BrookVelocityCenterOfMassRemoval::BrookVelocityCenterOfMassRemoval";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
BrookOpenMMFloat zero = static_cast<BrookOpenMMFloat>( 0.0 );
// ---------------------------------------------------------------------------------------
......@@ -150,7 +150,7 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
// static const std::string methodName = "\nBrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
BrookOpenMMFloat zero = static_cast<BrookOpenMMFloat>( 0.0 );
// ---------------------------------------------------------------------------------------
......@@ -158,13 +158,13 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
// subtract it (/totalMass) from velocities
BrookStreamInternal* velocityStream = dynamic_cast<BrookStreamInternal*> (vStream.getBrookStreamInternal());
*ke = static_cast<BrookOpenMMFloat>(0.0);
*ke = zero;
void* velV = velocityStream->getData( 1 );
const float* vArray = (float*) velV;
const float* vArray = static_cast<float*>( velV );
void* massV = getMassStream()->getData( 1 );
const float* mArray = (float*) massV;
const float* mArray = static_cast<float*>( massV );
int numberOfParticles = getNumberOfParticles();
int index = 0;
......@@ -176,13 +176,11 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
velocityCom[1] += mArray[ii]*vArray[index+1];
velocityCom[2] += mArray[ii]*vArray[index+2];
*ke += mArray[ii]*( vArray[index]*vArray[index] + vArray[index+1]*vArray[index+1] + vArray[index+2]*vArray[index+2] );
if( ii < 3 ){
if( 0 && ii < 3 ){
fprintf( stderr, "getVelocityCenterOfMass %5d %5d m=%15.6e v[%15.6e %15.6e %15.6e ] %d\n", ii, index, mArray[ii], vArray[index], vArray[index+1], vArray[index+2], numberOfParticles );
}
}
*ke *= static_cast<BrookOpenMMFloat>(0.5);
return DefaultReturnValue;
}
......@@ -258,7 +256,7 @@ int BrookVelocityCenterOfMassRemoval::_initializeStreams( const Platform& platfo
//static const std::string methodName = "BrookVelocityCenterOfMassRemoval::_initializeStreams";
BrookOpenMMFloat dangleValue = (BrookOpenMMFloat) 0.0;
BrookOpenMMFloat dangleValue = static_cast<BrookOpenMMFloat>( 0.0 );
// ---------------------------------------------------------------------------------------
......@@ -295,8 +293,8 @@ int BrookVelocityCenterOfMassRemoval::_setMasses( const std::vector<double>& mas
static const std::string methodName = "BrookVelocityCenterOfMassRemoval::_setMasses";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0;
BrookOpenMMFloat zero = static_cast<BrookOpenMMFloat>( 0.0 );
BrookOpenMMFloat one = static_cast<BrookOpenMMFloat>( 1.0 );
// ---------------------------------------------------------------------------------------
......@@ -314,7 +312,7 @@ int BrookVelocityCenterOfMassRemoval::_setMasses( const std::vector<double>& mas
memset( localMasses, 0, sizeof( BrookOpenMMFloat )*getComParticleStreamSize() );
int index = 0;
_totalInverseMass = (BrookOpenMMFloat) 0.0;
_totalInverseMass = zero;
for( std::vector<double>::const_iterator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
if( *ii != 0.0 ){
BrookOpenMMFloat value = static_cast<BrookOpenMMFloat>(*ii);
......@@ -369,9 +367,8 @@ int BrookVelocityCenterOfMassRemoval::setup( const std::vector<double>& masses,
_setMasses( masses );
//if( 1 && getLog() ){
if( 1 ){
FILE* log = stderr;
if( 1 && getLog() ){
FILE* log = getLog();
std::string contents = getContentsString( 0 );
(void) fprintf( log, "%s contents:\n%s\n", methodName.c_str(), contents.c_str() );
(void) fflush( log );
......@@ -457,31 +454,42 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
static const int debug = 1;
static int printOn = 0;
static const double boltz = 8.3145112119486e-03;
FILE* log;
BrookOpenMMFloat ke;
// ---------------------------------------------------------------------------------------
setLog( stderr );
if( debug && getLog() ){
//setLog( stderr );
if( printOn && getLog() ){
log = getLog();
} else {
printOn = 0;
}
// diagnostics
if( printOn ){
BrookOpenMMFloat com[3];
getVelocityCenterOfMass( velocityStream, com, &ke );
BrookOpenMMFloat denomiator = (( (float) 3*getNumberOfParticles()))*boltz;
(void) fprintf( getLog(), "%s Pre removal com: [%12.5e %12.5e %12.5e] ke=%.3f ~T=%.3f\n",
BrookOpenMMFloat denomiator = static_cast<float>( 3*getNumberOfParticles())*static_cast<float>( boltz );
(void) fprintf( log, "%s Pre removal com: [%12.5e %12.5e %12.5e] ke=%.3f ~T=%.3f\n",
methodName.c_str(), com[0], com[1], com[2], ke, ke/denomiator );
BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamInternal());
void* velV = outputVelocityStream->getData( 1 );
const float* vArray = (float*) velV;
const float* vArray = static_cast<float*>( velV );
int index = 0;
if( debug > 1 ){
if( printOn > 1 ){
for( int ii = 0; ii < getNumberOfParticles(); ii++, index += 3 ){
(void) fprintf( getLog(), "V %d [%12.5e %12.5e %12.5e]\n", ii, vArray[index], vArray[index+1], vArray[index+2] );
(void) fprintf( log, "V %d [%12.5e %12.5e %12.5e]\n", ii, vArray[index], vArray[index+1], vArray[index+2] );
}
}
(void) fflush( getLog() );
(void) fflush( log );
}
// calculate linear momentum via reduction
......@@ -492,37 +500,41 @@ setLog( stderr );
getWorkStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
kRemoveLinearMomentum( getLinearMomentumStream()->getBrookStream(), velocityStream.getBrookStream(), velocityStream.getBrookStream() );
if( (0 || debug) && getLog() ){
// diagnostics
if( printOn ){
BrookOpenMMFloat com[3];
getVelocityCenterOfMass( velocityStream, com, &ke );
BrookOpenMMFloat denomiator = (( (float) 3*getNumberOfParticles()))*boltz;
(void) fprintf( getLog(), "%s strW=%d iatm=%d invMass=%.4e\n Post removal com: [%12.5e %12.5e %12.5e] ke=%.3f ~T=%.3f\n",
BrookOpenMMFloat denomiator = static_cast<float>( 3*getNumberOfParticles() )*static_cast<float>( boltz );
(void) fprintf( log, "%s strW=%d iatm=%d invMass=%.4e\n Post removal com: [%12.5e %12.5e %12.5e] ke=%.3f ~T=%.3f\n",
methodName.c_str(),
getComParticleStreamWidth(), getNumberOfParticles(), getTotalInverseMass(), com[0], com[1], com[2],
ke, ke/denomiator );
void* linMoV = getLinearMomentumStream()->getData( 1 );
float* linMo = (float*) linMoV;
(void) fprintf( getLog(), " LM [%12.5e %12.5e %12.5e]\n", linMo[0], linMo[1], linMo[2] );
float* linMo = static_cast<float*>( linMoV );
(void) fprintf( log, " LM [%12.5e %12.5e %12.5e]\n", linMo[0], linMo[1], linMo[2] );
BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamInternal());
void* velV = outputVelocityStream->getData( 1 );
const float* vArray = (float*) velV;
const float* vArray = static_cast<float*>( velV );
void* w1 = getWorkStream()->getData( 1 );
const float* w2 = (float*) w1;
const float* w2 = static_cast<float*>( w1 );
if( debug > 1 ){
if( printOn > 1 ){
int index = 0;
for( int ii = 0; ii < getNumberOfParticles(); ii++, index += 3 ){
(void) fprintf( getLog(), "V %d [%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]\n", ii,
(void) fprintf( log, "V %d [%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]\n", ii,
vArray[index], vArray[index+1], vArray[index+2], w2[index], w2[index+1], w2[index+2] );
}
}
(void) fflush( getLog() );
(void) fflush( log );
//exit(0);
}
......
......@@ -652,7 +652,9 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
methodName.c_str(), _internalStepCount, inverseStepSize );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamInternal();
brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, log, 0.0001f );
std::string violationString;
brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, violationString, 0.0001f );
(void) fprintf( log, "Shake: %s\n", violationString.c_str() );
(void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( log );
......
......@@ -29,26 +29,37 @@ SET(STATIC_BROOK_TARGET ${OpenMM_BROOK_LIBRARY_NAME}_static)
# Automatically create tests using files named "Test*.cpp"
FILE(GLOB TEST_PROGS "*Test*.cpp")
FOREACH(TEST_PROG ${TEST_PROGS})
GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
# Link with shared library
#ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
#TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
# ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
#ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
# ----------------------------------------------------------------------------
IF(LOG)
FILE( APPEND ${LOG_FILE} "TARGET_LINK_LIBRARIES: ${TEST_PROG} TARGET=${SHARED_TARGET} BROOK_TARGET=${ROOK_TARGET} BROOK_LIB=${BROOK_LIB}\n")
ENDIF(LOG)
# ----------------------------------------------------------------------------
SET( CMAKE_EXE_LINKER_FLAGS "/NODEFAULTLIB:\"LIBCMT.lib\"")
SET( CMAKE_EXE_LINKER_FLAGS_DEBUG "/NODEFAULTLIB:\"LIBCMTD.lib\"")
# SET( CMAKE_EXE_LINKER_FLAGS "/NODEFAULTLIB:\"LIBCMT.lib\"")
ADD_DEFINITIONS(-D_WIN32 )
ADD_DEFINITIONS( -D_WIN32 )
# Link with static library
SET(TEST_STATIC ${TEST_ROOT}Static)
ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
# SET_TARGET_PROPERTIES(${TEST_STATIC}
# PROPERTIES
# COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
# )
TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET} ${STATIC_BROOK_TARGET} ${BROOK_LIB})
# ADD_TEST(${TEST_STATIC} ${EXECUTABLE_OUTPUT_PATH}/${TEST_STATIC})
......
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