Skip to content

GitLab

Commit ce8d2142 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods

parent 9cfbfa4c
Hide whitespace changes
Inline Side-by-side
Showing with 725 additions and 381 deletions
platforms/brook/src/BrookBonded.cpp
View file @ ce8d2142
...@@ -826,7 +826,7 @@ int BrookBonded::_addPTorsions( int *nbondeds, int *particles, BrookOpenMMFloat* ...@@ -826,7 +826,7 @@ int BrookBonded::_addPTorsions( int *nbondeds, int *particles, BrookOpenMMFloat*
// note: parameters 0 & 2 switched // note: parameters 0 & 2 switched
PARAMS( ibonded, 1, 1 ) = (BrookOpenMMFloat) pTParameters[2]; PARAMS( ibonded, 1, 1 ) = (BrookOpenMMFloat) pTParameters[2];
PARAMS( ibonded, 1, 2 ) = (BrookOpenMMFloat) pTParameters[1]*DEG2RAD; PARAMS( ibonded, 1, 2 ) = (BrookOpenMMFloat) pTParameters[1];
PARAMS( ibonded, 1, 3 ) = (BrookOpenMMFloat) pTParameters[0]; PARAMS( ibonded, 1, 3 ) = (BrookOpenMMFloat) pTParameters[0];
if( debug ){ if( debug ){
......
platforms/brook/src/BrookFloatStreamInternal.cpp
View file @ ce8d2142
...@@ -529,37 +529,131 @@ const std::string BrookFloatStreamInternal::getContentsString( int level ) const ...@@ -529,37 +529,131 @@ const std::string BrookFloatStreamInternal::getContentsString( int level ) const
* *
* */ * */
int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log ){ int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log, int maxPrint ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookStreamInternal::_bodyPrintToFile"; //static const std::string methodName = "BrookStreamInternal::_bodyPrintToFile";
static const unsigned int MAX_LINE_CHARS = 256;
char value[MAX_LINE_CHARS];
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
assert( log ); assert( log );
void* dataArrayV = getData( 1 ); void* dataArrayV = getData( 1 );
//void* dataArrayV = getDataArray( );
//saveToArray( dataArrayV ); #ifdef WIN32
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );
#else
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );
#endif
int streamSize = getStreamSize(); int streamSize = getStreamSize();
int width = getWidth(); int width = getWidth();
int index = 0; int index = 0;
maxPrint = maxPrint < 0 ? streamSize : maxPrint;
maxPrint *= width;
const float* dataArray = (float*) dataArrayV; const float* dataArray = (float*) dataArrayV;
for( int ii = 0; ii < streamSize; ii++ ){ for( int ii = 0; ii < streamSize && ii < maxPrint; ii++ ){
std::stringstream message; std::stringstream message;
message.width( 15 ); (void) LOCAL_SPRINTF( value, "%6d ", ii );
message.precision( 8 ); message << value << " [ ";
message << ii << " [ ";
for( int jj = 0; jj < width; jj++ ){ for( unsigned int jj = 0; jj < width; jj++ ){
message << dataArray[index++] << " "; (void) LOCAL_SPRINTF( value, "%16.7e ", dataArray[index++] );
} message << value;
}
message << "]\n"; message << "]\n";
if( index == (_size+1)*width ){
(void) fprintf( log, "\n" );
}
(void) fprintf( log, "%s", message.str().c_str() ); (void) fprintf( log, "%s", message.str().c_str() );
} }
//delete[] dataArrayV; return DefaultReturnValue;
}
/*
* Get stats
*
* @return statistics vector
*
* */
int BrookFloatStreamInternal::getStatistics( std::vector<std::vector<double> >& statistics, int maxScan ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookStreamInternal::getStatistics";
static const int MinIndex = 4;
static const int MinIndexIndex = 5;
static const int MaxIndex = 6;
static const int MaxIndexIndex = 7;
static const int CountIndex = 8;
static const int vectorSize = CountIndex + 1;
static const double bigValue = 1.0e+10;
// ---------------------------------------------------------------------------------------
void* dataArrayV = getData( 1 );
statistics.resize( vectorSize );
int streamSize = getStreamSize();
int width = getWidth();
int index = 0;
const float* dataArray = (float*) dataArrayV;
for( int ii = 0; ii < vectorSize; ii++ ){
for( int jj = 0; jj < width; jj++ ){
if( ii == MinIndex ){
statistics[ii].push_back( bigValue );
} else if( ii == MaxIndex ){
statistics[ii].push_back( -bigValue );
} else {
statistics[ii].push_back( 0.0 );
}
}
}
for( int ii = 0; ii < streamSize && ii < maxScan; ii++ ){
for( int jj = 0; jj < width; jj++ ){
double value = (double) dataArray[index++];
statistics[0][jj] += value;
statistics[1][jj] += value*value;
statistics[2][jj] += fabs( value );
statistics[CountIndex][jj] += 1.0;
if( value < statistics[MinIndex][jj] ){
statistics[MinIndex][jj] = value;
statistics[MinIndexIndex][jj] = ii;
}
if( value > statistics[MaxIndex][jj] ){
statistics[MaxIndex][jj] = value;
statistics[MaxIndexIndex][jj] = ii;
}
}
}
for( int jj = 0; jj < width; jj++ ){
if( statistics[CountIndex][jj] > 0.0 ){
statistics[3][jj] = statistics[0][jj];
statistics[0][jj] /= statistics[CountIndex][jj];
statistics[2][jj] /= statistics[CountIndex][jj];
statistics[1][jj] = statistics[1][jj] - statistics[0][jj]*statistics[0][jj]*statistics[CountIndex][jj];
if( statistics[CountIndex][jj] > 1.0 ){
statistics[1][jj] = sqrt( statistics[1][jj]/( statistics[CountIndex][jj] - 1.0 ) );
}
}
}
return DefaultReturnValue; return DefaultReturnValue;
......
platforms/brook/src/BrookFloatStreamInternal.h
View file @ ce8d2142
...@@ -172,6 +172,15 @@ class BrookFloatStreamInternal : public BrookStreamInternal { ...@@ -172,6 +172,15 @@ class BrookFloatStreamInternal : public BrookStreamInternal {
int sumByDimension( int stopIndex, double* sum ); int sumByDimension( int stopIndex, double* sum );
/*
* Get stats
*
* @return statistics vector
*
* */
int getStatistics( std::vector<std::vector<double>>&, int maxScan );
private: private:
BrookOpenMMFloat _defaultDangleValue; BrookOpenMMFloat _defaultDangleValue;
...@@ -192,7 +201,7 @@ class BrookFloatStreamInternal : public BrookStreamInternal { ...@@ -192,7 +201,7 @@ class BrookFloatStreamInternal : public BrookStreamInternal {
* *
* */ * */
int _bodyPrintToFile( FILE* log ); int _bodyPrintToFile( FILE* log, int maxPrint );
}; };
} // namespace OpenMM } // namespace OpenMM
......
platforms/brook/src/BrookIntStreamInternal.cpp
View file @ ce8d2142
...@@ -256,7 +256,7 @@ void* BrookIntStreamInternal::getData( int readFromBoard ){ ...@@ -256,7 +256,7 @@ void* BrookIntStreamInternal::getData( int readFromBoard ){
* *
* */ * */
int BrookIntStreamInternal::_bodyPrintToFile( FILE* log ){ int BrookIntStreamInternal::_bodyPrintToFile( FILE* log, int maxPrint ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
platforms/brook/src/BrookIntStreamInternal.h
View file @ ce8d2142
...@@ -182,7 +182,7 @@ private: ...@@ -182,7 +182,7 @@ private:
* *
* */ * */
int _bodyPrintToFile( FILE* log ); int _bodyPrintToFile( FILE* log, int maxPrint );
}; };
} // namespace OpenMM } // namespace OpenMM
......
platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp
View file @ ce8d2142
...@@ -123,13 +123,17 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L ...@@ -123,13 +123,17 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const int printOn = 0; int printOn = 1;
static const std::string methodName = "BrookIntegrateLangevinStepKernel::initialize"; static const std::string methodName = "BrookIntegrateLangevinStepKernel::initialize";
FILE* log = NULL;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
FILE* log = getLog(); setLog( stderr );
if( printOn && log ){ printOn = (printOn && getLog()) ? printOn : 0;
if( printOn ){
log = getLog();
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fprintf( log, "%s\n", methodName.c_str() );
(void) fflush( log ); (void) fflush( log );
} }
...@@ -141,7 +145,7 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L ...@@ -141,7 +145,7 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
std::vector<double> masses; std::vector<double> masses;
masses.resize( numberOfParticles ); masses.resize( numberOfParticles );
if( printOn && log ){ if( printOn ){
(void) fprintf( log, "%s %d\n", methodName.c_str(), numberOfParticles ); (void) fprintf( log, "%s %d\n", methodName.c_str(), numberOfParticles );
(void) fflush( log ); (void) fflush( log );
} }
...@@ -154,7 +158,7 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L ...@@ -154,7 +158,7 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
int numberOfConstraints = system.getNumConstraints(); int numberOfConstraints = system.getNumConstraints();
if( printOn && log ){ if( printOn ){
(void) fprintf( log, "%s const=%d\n", methodName.c_str(), numberOfConstraints ); (void) fprintf( log, "%s const=%d\n", methodName.c_str(), numberOfConstraints );
(void) fflush( log ); (void) fflush( log );
} }
...@@ -188,15 +192,17 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L ...@@ -188,15 +192,17 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
BrookOpenMMFloat tolerance = static_cast<BrookOpenMMFloat>( integrator.getConstraintTolerance() ); BrookOpenMMFloat tolerance = static_cast<BrookOpenMMFloat>( integrator.getConstraintTolerance() );
_brookShakeAlgorithm->setShakeTolerance( tolerance ); _brookShakeAlgorithm->setShakeTolerance( tolerance );
_brookShakeAlgorithm->setMaxIterations( 30 ); _brookShakeAlgorithm->setMaxIterations( 40 );
_brookShakeAlgorithm->setLog( log ); if( log ){
_brookShakeAlgorithm->setLog( log );
}
// random number generator // random number generator
_brookRandomNumberGenerator = new BrookRandomNumberGenerator( ); _brookRandomNumberGenerator = new BrookRandomNumberGenerator( );
_brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() ); _brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
if( printOn && log ){ if( printOn ){
(void) fprintf( log, "%s done setup:\nBrookShakeAlgorithm:\n%s\nBrookRandomNumberGenerator:\n%s\n\n", methodName.c_str(), (void) fprintf( log, "%s done setup:\nBrookShakeAlgorithm:\n%s\nBrookRandomNumberGenerator:\n%s\n\n", methodName.c_str(),
_brookShakeAlgorithm->getContentsString().c_str(), _brookShakeAlgorithm->getContentsString().c_str(),
_brookRandomNumberGenerator->getContentsString().c_str() ); _brookRandomNumberGenerator->getContentsString().c_str() );
...@@ -234,10 +240,7 @@ void BrookIntegrateLangevinStepKernel::execute( OpenMMContextImpl& context, cons ...@@ -234,10 +240,7 @@ void BrookIntegrateLangevinStepKernel::execute( OpenMMContextImpl& context, cons
differences[1] = integrator.getFriction() - (double) _brookLangevinDynamics->getFriction(); differences[1] = integrator.getFriction() - (double) _brookLangevinDynamics->getFriction();
differences[2] = integrator.getStepSize() - (double) _brookLangevinDynamics->getStepSize(); differences[2] = integrator.getStepSize() - (double) _brookLangevinDynamics->getStepSize();
if( fabs( differences[0] ) > epsilon || fabs( differences[1] ) > epsilon || fabs( differences[2] ) > epsilon ){ if( fabs( differences[0] ) > epsilon || fabs( differences[1] ) > epsilon || fabs( differences[2] ) > epsilon ){
//printf( "%s calling updateParameters\n", methodName.c_str() );
_brookLangevinDynamics->updateParameters( integrator.getTemperature(), integrator.getFriction(), integrator.getStepSize() ); _brookLangevinDynamics->updateParameters( integrator.getTemperature(), integrator.getFriction(), integrator.getStepSize() );
} else {
//printf( "%s NOT calling updateParameters\n", methodName.c_str() );
} }
_brookLangevinDynamics->update( *(_openMMBrookInterface.getParticlePositions()), *(_openMMBrookInterface.getParticleVelocities()), _brookLangevinDynamics->update( *(_openMMBrookInterface.getParticlePositions()), *(_openMMBrookInterface.getParticleVelocities()),
......
platforms/brook/src/BrookLangevinDynamics.cpp
View file @ ce8d2142
...@@ -64,6 +64,7 @@ BrookLangevinDynamics::BrookLangevinDynamics( ){ ...@@ -64,6 +64,7 @@ BrookLangevinDynamics::BrookLangevinDynamics( ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_numberOfParticles = -1; _numberOfParticles = -1;
_internalStepCount = 0;
// mark stream dimension variables as unset // mark stream dimension variables as unset
...@@ -808,6 +809,46 @@ int BrookLangevinDynamics::setup( const std::vector<double>& masses, const Platf ...@@ -808,6 +809,46 @@ int BrookLangevinDynamics::setup( const std::vector<double>& masses, const Platf
return DefaultReturnValue; return DefaultReturnValue;
} }
/*
* Setup of LangevinDynamics parameters
*
* @param masses masses
* @param platform Brook platform
*
* @return nonzero value if error
*
* */
float BrookLangevinDynamics::getTemperature( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream,
BrookShakeAlgorithm& brookShakeAlgorithm ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookLangevinDynamics::getTemperature";
// ---------------------------------------------------------------------------------------
void* dataArrayV = velocities->getData( 1 );
float* velocitiesI = (float*) dataArrayV;
void* inverseMassStreamV = inverseMassStream->getData( 1 );
float* inverseMassStreamI = (float*) inverseMassStreamV;
float ke = 0.0f;
int index = 0;
for( int ii = 0; ii < getNumberOfParticles(); ii++ ){
ke += (velocitiesI[index]*velocitiesI[index] + velocitiesI[index+1]*velocitiesI[index+1] + velocitiesI[index+2]*velocitiesI[index+2] )/inverseMassStreamI[ii];
index += 3;
}
int degreesOfFreedom = 3*getNumberOfParticles() - brookShakeAlgorithm.getNumberOfConstraints();
ke /= ((float) BOLTZ)*((float) ( degreesOfFreedom ));
return ke;
}
/* /*
* Get contents of object * Get contents of object
* *
...@@ -887,6 +928,7 @@ std::string BrookLangevinDynamics::getContentsString( int level ) const { ...@@ -887,6 +928,7 @@ std::string BrookLangevinDynamics::getContentsString( int level ) const {
return message.str(); return message.str();
} }
/** /**
* Update * Update
* *
...@@ -908,19 +950,26 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -908,19 +950,26 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "\nBrookLangevinDynamics::update"; static const std::string methodName = "\nBrookLangevinDynamics::update";
static const int PrintOn = 0; int printOn = 1;
FILE* log;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_internalStepCount++;
setLog( stderr );
printOn = (printOn && getLog()) ? printOn : 0;
const BrookOpenMMFloat* derivedParameters = getDerivedParameters(); const BrookOpenMMFloat* derivedParameters = getDerivedParameters();
if( (0 || PrintOn) && getLog() ){ if( (0 || printOn) ){
log = getLog();
static int showAux = 1; static int showAux = 1;
if( PrintOn ){ if( printOn ){
(void) fprintf( getLog(), "%s shake=%d\n", methodName.c_str(), brookShakeAlgorithm.getNumberOfConstraints() ); (void) fprintf( log, "%s step=%d shake=%d\n", methodName.c_str(), _internalStepCount, brookShakeAlgorithm.getNumberOfConstraints() );
(void) fflush( getLog() ); (void) fflush( log );
} }
// show update // show update
...@@ -929,13 +978,12 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -929,13 +978,12 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
showAux = 0; showAux = 0;
std::string contents = brookRandomNumberGenerator.getContentsString( ); std::string contents = brookRandomNumberGenerator.getContentsString( );
(void) fprintf( getLog(), "%s RNG contents\n%s", methodName.c_str(), contents.c_str() ); (void) fprintf( log, "%s step=%d RNG contents\n%s", methodName.c_str(), _internalStepCount, contents.c_str() );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log );
contents = brookShakeAlgorithm.getContentsString( ); contents = brookShakeAlgorithm.getContentsString( );
(void) fprintf( getLog(), "%s Shake contents\n%s", methodName.c_str(), contents.c_str() ); (void) fprintf( log, "%s step=%d Shake contents\n%s", methodName.c_str(), _internalStepCount, contents.c_str() );
(void) fflush( log );
(void) fflush( getLog() );
} }
} }
...@@ -963,48 +1011,51 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -963,48 +1011,51 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
// diagnostics // diagnostics
if( PrintOn ){ if( 0 && printOn ){
(void) fprintf( getLog(), "\n%s Post kupdate_sd1_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d " (void) fprintf( log, "\n%s step=%d Post kupdate_sd1_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]", methodName.c_str(), "EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]", methodName.c_str(), _internalStepCount,
getLangevinDynamicsParticleStreamWidth(), getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(), brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(), brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] ); derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] );
(void) fprintf( getLog(), "\nSDPC1Stream\n" ); if( _internalStepCount == 1 ){
getSDPC1Stream()->printToFile( getLog() ); (void) fprintf( log, "\nSDPC1Stream %d\n", _internalStepCount );
getSDPC1Stream()->printToFile( log );
(void) fprintf( getLog(), "\nSD2XStream\n" );
getSD2XStream()->printToFile( getLog() ); (void) fprintf( log, "\nSD2XStream %d\n", _internalStepCount );
getSD2XStream()->printToFile( log );
(void) fprintf( getLog(), "\nInverseMassStream\n" );
getInverseMassStream()->printToFile( getLog() ); (void) fprintf( log, "\nInverseMassStream %d\n", _internalStepCount );
getInverseMassStream()->printToFile( log );
}
//StreamImpl& positionStreamImpl = positionStream.getImpl(); //StreamImpl& positionStreamImpl = positionStream.getImpl();
//const BrookStreamImpl brookPositions = dynamic_cast<BrookStreamImpl&> (positionStreamImpl); //const BrookStreamImpl brookPositions = dynamic_cast<BrookStreamImpl&> (positionStreamImpl);
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nForceStream\n" ); (void) fprintf( log, "\nForceStream %d\n", _internalStepCount );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() ); brookStreamInternalF->printToFile( log );
BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" ); (void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
brookStreamInternalV->printToFile( getLog() ); brookStreamInternalV->printToFile( log );
(void) fprintf( getLog(), "\nSD1VStream\n" ); (void) fprintf( log, "\nSD1VStream %d\n", _internalStepCount );
getSD1VStream()->printToFile( getLog() ); getSD1VStream()->printToFile( log );
(void) fprintf( getLog(), "\nVPrimeStream\n" ); (void) fprintf( log, "\nVPrimeStream %d\n", _internalStepCount );
getVPrimeStream()->printToFile( getLog() ); getVPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() ); (void) fprintf( log, "\nRvStreamIndex=%d step=%d\n", brookRandomNumberGenerator.getRvStreamIndex(), _internalStepCount );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log );
} }
// advance random number cursor // advance random number cursor
...@@ -1040,36 +1091,37 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -1040,36 +1091,37 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() ); getXPrimeStream()->getBrookStream() );
if( ( 0 || PrintOn) && getLog() ){ if( 0 && printOn ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d", (void) fprintf( log, "\n%s Post kshakeh_update2_fix1: particleStrW=%3d",
methodName.c_str(), getLangevinDynamicsParticleStreamWidth() ); methodName.c_str(), getLangevinDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeInverseMapStream\n" ); if( _internalStepCount == 1 ){
brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() ); (void) fprintf( log, "\nShakeInverseMapStream %d\n" );
brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( log );
}
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons0\n" ); (void) fprintf( log, "\nShakeXCons0 %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons1\n" ); (void) fprintf( log, "\nShakeXCons1 %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons2\n" ); (void) fprintf( log, "\nShakeXCons2 %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons3\n" ); (void) fprintf( log, "\nShakeXCons3 %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( log );
} }
} }
// second integration step // second integration step
...@@ -1093,39 +1145,39 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -1093,39 +1145,39 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
// diagnostics // diagnostics
if( 0 || PrintOn ){ if( 0 && printOn ){
(void) fprintf( getLog(), "\n%s Post kupdate_sd2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d " (void) fprintf( log, "\n%s step=%d Post kupdate_sd2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"Sd2pc[]=[%12.5e %12.5e]", methodName.c_str(), "Sd2pc[]=[%12.5e %12.5e]", methodName.c_str(), _internalStepCount,
getLangevinDynamicsParticleStreamWidth(), getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(), brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(), brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] ); derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
(void) fprintf( getLog(), "\nSDPC2Stream\n" ); (void) fprintf( log, "\nSDPC2Stream %d \n", _internalStepCount );
getSDPC2Stream()->printToFile( getLog() ); getSDPC2Stream()->printToFile( log );
(void) fprintf( getLog(), "\nSD2XStream\n" ); (void) fprintf( log, "\nSD2XStream %d \n", _internalStepCount );
getSD1VStream()->printToFile( getLog() ); getSD1VStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d \n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nVPrimeStream\n" ); (void) fprintf( log, "\nVPrimeStream %d \n", _internalStepCount );
getVPrimeStream()->printToFile( getLog() ); getVPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d \n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\ngetSD2XStream\n" ); (void) fprintf( log, "\ngetSD2XStream %d \n", _internalStepCount );
getSD2XStream()->printToFile( getLog() ); getSD2XStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalVel = velocityStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalVel = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" ); (void) fprintf( log, "\nVelocityStream %d \n", _internalStepCount );
brookStreamInternalVel->printToFile( getLog() ); brookStreamInternalVel->printToFile( log );
(void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() ); (void) fprintf( log, "\nRvStreamIndex=%d step=%d\n", brookRandomNumberGenerator.getRvStreamIndex(), _internalStepCount );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log );
} }
// advance random number cursor // advance random number cursor
...@@ -1164,41 +1216,58 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -1164,41 +1216,58 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
// diagnostics // diagnostics
if( PrintOn ){ if( printOn ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d " (void) fprintf( log, "\n%s step=%d Post kshakeh_update2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"Sd2pc[]=[%12.5e %12.5e]", methodName.c_str(), "Sd2pc[]=[%12.5e %12.5e]\n", methodName.c_str(), _internalStepCount,
getLangevinDynamicsParticleStreamWidth(), getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(), brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(), brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] ); derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, log, 0.0001f );
brookStreamInternalPos->printToFile( getLog() ); (void) fprintf( log, "\nPositionStream\n" );
brookStreamInternalPos->printToFile( log );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
std::vector<std::vector<double> > forceStatistics;
brookStreamInternalF->getStatistics( forceStatistics, getNumberOfParticles() );
std::stringstream tag;
tag << _internalStepCount << " Fxx ";
std::string stats = brookStreamInternalF->printStatistics( tag.str(), forceStatistics );
(void) fprintf( log, "\nStep %d Force stats:\n%s", _internalStepCount, stats.c_str() );
brookStreamInternalF->printToFile( log );
brookStreamInternalPos = velocityStream.getBrookStreamImpl(); brookStreamInternalPos = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" ); std::vector<std::vector<double> > velocityStatistics;
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->getStatistics( velocityStatistics, getNumberOfParticles() );
std::stringstream tagV;
(void) fprintf( getLog(), "\nXPrimeStream\n" ); tagV << _internalStepCount << " Vxx ";
getXPrimeStream()->printToFile( getLog() ); stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
(void) fprintf( log, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
float temperature = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm );
(void) fprintf( log, "\nVelocityStream %d T=%.3f\n", _internalStepCount, temperature );
brookStreamInternalPos->printToFile( log );
(void) fprintf( log, "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( log );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log );
} }
} else { } else {
// no constraints // no constraints
if( PrintOn ){ if( printOn ){
(void) fprintf( getLog(), "\n%s Pre ksetStr3 (no constraints)", methodName.c_str() ); (void) fprintf( log, "\n%s Pre ksetStr3 (no constraints)", methodName.c_str() );
// (void) fprintf( getLog(), "\nXPrimeStream\n" ); // (void) fprintf( log, "\nXPrimeStream\n" );
// getXPrimeStream()->printToFile( getLog() ); // getXPrimeStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
} }
kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream(), positionStream.getBrookStream() ); kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream(), positionStream.getBrookStream() );
...@@ -1206,17 +1275,17 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI ...@@ -1206,17 +1275,17 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
// diagnostics // diagnostics
if( PrintOn ){ if( printOn ){
(void) fprintf( getLog(), "\n%s Post ksetStr3 (no constraints)", methodName.c_str() ); (void) fprintf( log, "\n%s Post ksetStr3 (no constraints)", methodName.c_str() );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log );
} }
} }
......
platforms/brook/src/BrookLangevinDynamics.h
View file @ ce8d2142
...@@ -258,6 +258,17 @@ class BrookLangevinDynamics : public BrookCommon { ...@@ -258,6 +258,17 @@ class BrookLangevinDynamics : public BrookCommon {
BrookFloatStreamInternal* getInverseMassStream( void ) const; BrookFloatStreamInternal* getInverseMassStream( void ) const;
/**
* Get T
*
* @return T
*
*/
float getTemperature( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream,
BrookShakeAlgorithm& brookShakeAlgorithm ) const;
private: private:
enum DerivedParameters { GDT, EPH, EMH, EP, EM, B, C, D, V, X, Yv, Yx, enum DerivedParameters { GDT, EPH, EMH, EP, EM, B, C, D, V, X, Yv, Yx,
...@@ -286,6 +297,10 @@ class BrookLangevinDynamics : public BrookCommon { ...@@ -286,6 +297,10 @@ class BrookLangevinDynamics : public BrookCommon {
BrookOpenMMFloat _temperature; BrookOpenMMFloat _temperature;
BrookOpenMMFloat _stepSize; BrookOpenMMFloat _stepSize;
// internal step count
int _internalStepCount;
// Particle stream dimensions // Particle stream dimensions
int _sdParticleStreamWidth; int _sdParticleStreamWidth;
......
platforms/brook/src/BrookNonBonded.cpp
View file @ ce8d2142
...@@ -1279,17 +1279,17 @@ void BrookNonBonded::computeForces( BrookStreamImpl& positionStream, BrookStream ...@@ -1279,17 +1279,17 @@ void BrookNonBonded::computeForces( BrookStreamImpl& positionStream, BrookStream
static const std::string methodName = "BrookNonBonded::computeForces"; static const std::string methodName = "BrookNonBonded::computeForces";
static int PrintOn = 1; static int printOn = 0;
static const int MaxErrorMessages = 2; static const int MaxErrorMessages = 2;
static int ErrorMessages = 0; static int ErrorMessages = 0;
static const float4 dummyParameters( 0.0, 0.0, 0.0, 0.0 ); FILE* log;
// static const int debug = 1; // static const int debug = 1;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
PrintOn = (PrintOn && getLog()) ? 1 : 0; printOn = (printOn && getLog()) ? printOn : 0;
// nonbonded forces // nonbonded forces
...@@ -1305,7 +1305,6 @@ void BrookNonBonded::computeForces( BrookStreamImpl& positionStream, BrookStream ...@@ -1305,7 +1305,6 @@ void BrookNonBonded::computeForces( BrookStreamImpl& positionStream, BrookStream
(float) getJStreamWidth( ), (float) getJStreamWidth( ),
(float) getPartialForceStreamWidth( ), (float) getPartialForceStreamWidth( ),
epsfac, epsfac,
dummyParameters,
positionStream.getBrookStream(), positionStream.getBrookStream(),
getChargeStream()->getBrookStream(), getChargeStream()->getBrookStream(),
getOuterVdwStream()->getBrookStream(), getOuterVdwStream()->getBrookStream(),
...@@ -1328,7 +1327,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -1328,7 +1327,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
// diagnostics // diagnostics
//if( 1 && PrintOn ){ //if( 1 && printOn ){
static int step = 0; static int step = 0;
if( step++ < 1 ){ if( step++ < 1 ){
//FILE* log = getLog(); //FILE* log = getLog();
...@@ -1387,15 +1386,12 @@ static int step = 0; ...@@ -1387,15 +1386,12 @@ static int step = 0;
// diagnostics // diagnostics
//if( PrintOn ){ if( printOn ){
if( 1 ){
//FILE* log = getLog(); (void) fprintf( log, "\n%s NB forces\n", methodName.c_str() );
FILE* log = stderr;
(void) fprintf( log, "\n%s NB forces", methodName.c_str() );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( log ); brookStreamInternalF->printToFile( log );
(void) fprintf( log, "\n%s Done", methodName.c_str() ); (void) fflush( log ); (void) fprintf( log, "\n%s Done\n", methodName.c_str() ); (void) fflush( log );
} }
......
platforms/brook/src/BrookShakeAlgorithm.cpp
View file @ ce8d2142
...@@ -29,48 +29,13 @@ ...@@ -29,48 +29,13 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
#include <sstream>
#include "BrookShakeAlgorithm.h" #include "BrookShakeAlgorithm.h"
#include "BrookPlatform.h" #include "BrookPlatform.h"
#include "OpenMMException.h"
#include "BrookStreamImpl.h" #include "BrookStreamImpl.h"
using namespace OpenMM; using namespace OpenMM;
using namespace std; using namespace std;
struct ShakeCluster {
int _centralID;
int _peripheralID[3];
int _size;
float _distance;
float _centralInvMass, _peripheralInvMass;
ShakeCluster(){};
ShakeCluster( int centralID, float invMass) : _centralID(centralID), _centralInvMass(invMass), _size(0) {};
void addAtom( int id, float dist, float invMass){
if( _size == 3 ){
std::stringstream message;
message << "ShakeCluster::addAtom: " << "atom " << id << " has more than 3 constraints!." << std::endl;
throw OpenMMException( message.str() );
}
if( _size > 0 && dist != _distance ){
std::stringstream message;
message << "ShakeCluster::addAtom: " << "atom " << id << " has different constraint distances: " << dist << " and " << _distance << std::endl;
throw OpenMMException( message.str() );
}
if( _size > 0 && invMass != _peripheralInvMass ){
std::stringstream message;
message << "ShakeCluster::addAtom: " << " constrainted atoms associated w/ atom " << id << " have different masses: " << invMass << " and " << _peripheralInvMass << std::endl;
throw OpenMMException( message.str() );
}
_peripheralID[_size++] = id;
_distance = dist;
_peripheralInvMass = invMass;
}
};
/** /**
* *
* Constructor * Constructor
...@@ -496,7 +461,6 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co ...@@ -496,7 +461,6 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
// Find clusters consisting of a central atom with up to three peripheral atoms. // Find clusters consisting of a central atom with up to three peripheral atoms.
map<int, ShakeCluster> clusters;
particleIterator = constraintIndices.begin(); particleIterator = constraintIndices.begin();
std::vector<double>::const_iterator distanceIterator = constraintLengths.begin(); std::vector<double>::const_iterator distanceIterator = constraintLengths.begin();
while( particleIterator != constraintIndices.end() ){ while( particleIterator != constraintIndices.end() ){
...@@ -537,10 +501,10 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co ...@@ -537,10 +501,10 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
// Add it to the cluster // Add it to the cluster
if( clusters.find(centralID) == clusters.end() ){ if( _clusters.find(centralID) == _clusters.end() ){
clusters[centralID] = ShakeCluster( centralID, centralInvMass ); _clusters[centralID] = ShakeCluster( centralID, 1.0f/centralInvMass );
} }
clusters[centralID].addAtom( peripheralID, distance, peripheralInvMass ); _clusters[centralID].addAtom( peripheralID, distance, 1.0f/peripheralInvMass );
particleIterator++; particleIterator++;
distanceIterator++; distanceIterator++;
...@@ -548,7 +512,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co ...@@ -548,7 +512,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
// allocate space for streams // allocate space for streams
_initializeStreamSizes( (int) masses.size(), (int) clusters.size(), platform ); _initializeStreamSizes( (int) masses.size(), (int) _clusters.size(), platform );
_initializeStreams( platform ); _initializeStreams( platform );
int shakeParticleStreamSize = getShakeParticleStreamSize(); int shakeParticleStreamSize = getShakeParticleStreamSize();
...@@ -567,8 +531,8 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co ...@@ -567,8 +531,8 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
// load indices & parameters // load indices & parameters
int constraintIndex = 0; int constraintIndex = 0;
_numberOfConstraints = static_cast<int>(clusters.size()); _numberOfConstraints = static_cast<int>(_clusters.size());
for( map<int, ShakeCluster>::const_iterator iter = clusters.begin(); iter != clusters.end(); ++iter ){ for( map<int, ShakeCluster>::const_iterator iter = _clusters.begin(); iter != _clusters.end(); ++iter ){
const ShakeCluster& cluster = iter->second; const ShakeCluster& cluster = iter->second;
...@@ -577,10 +541,10 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co ...@@ -577,10 +541,10 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
particleIndices[constraintIndex+ii+1] = static_cast<BrookOpenMMFloat>( cluster._peripheralID[ii] ); particleIndices[constraintIndex+ii+1] = static_cast<BrookOpenMMFloat>( cluster._peripheralID[ii] );
} }
shakeParameters[constraintIndex] = one/( static_cast<BrookOpenMMFloat>( cluster._centralInvMass ) ); shakeParameters[constraintIndex] = static_cast<BrookOpenMMFloat>( cluster._centralInvMass );
shakeParameters[constraintIndex+1] = half/( static_cast<BrookOpenMMFloat>(cluster._centralInvMass + cluster._peripheralInvMass) ); shakeParameters[constraintIndex+1] = half/( static_cast<BrookOpenMMFloat>( cluster._centralInvMass + cluster._peripheralInvMass) );
shakeParameters[constraintIndex+2] = static_cast<BrookOpenMMFloat>( cluster._distance*cluster._distance ); shakeParameters[constraintIndex+2] = static_cast<BrookOpenMMFloat>( cluster._distance*cluster._distance );
shakeParameters[constraintIndex+3] = one/( static_cast<BrookOpenMMFloat>( cluster._peripheralInvMass ) ); shakeParameters[constraintIndex+3] = static_cast<BrookOpenMMFloat>( cluster._peripheralInvMass );
constraintIndex += 4; constraintIndex += 4;
} }
...@@ -650,6 +614,106 @@ int BrookShakeAlgorithm::setup( const std::vector<double>& masses, const std::ve ...@@ -650,6 +614,106 @@ int BrookShakeAlgorithm::setup( const std::vector<double>& masses, const std::ve
return DefaultReturnValue; return DefaultReturnValue;
} }
/*
* Check constraints
*
* @param positions atom positions
* @param log file to print to (can be NULL)
* @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 {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookShakeAlgorithm::checkConstraints";
static const int maxErrorToPrint = -10;
// ---------------------------------------------------------------------------------------
void* posArrayV = positions->getData( 1 );
const float* posArray = (float*) posArrayV;
/*
fprintf( log, "ShakeCoords\n" );
for( int ii = 0; ii < 30; ii++ ){
if( !(ii % 3 ) )fprintf( log, "%d [", ii );
fprintf( log, "%16.7e ", posArray[ii] );
if( !(ii % 2 ) )fprintf( log, "]\n", posArray[ii] );
} */
// loop over constraints checking if distances are satisfied
// build up message string if errors are detected
// also track max difference in case all constraints satisfied
// to within specified tolerance
std::stringstream message;
int errors = 0;
float maxDiff = -1.0f;
int maxDiffCentralIndex = -1;
int maxDiffPeripheralIndex = -1;
tolerance = tolerance > 0.0f ? tolerance : _shakeTolerance;
for( map<int, ShakeCluster>::const_iterator iter = _clusters.begin(); iter != _clusters.end(); ++iter ){
const ShakeCluster& cluster = iter->second;
// multiply by 3 to get index into 1d array
int centralAtomIndex = 3*cluster._centralID;
for( int ii = 0; ii < cluster._size; ii++ ){
int peripheralAtomIndex = 3*cluster._peripheralID[ii];
if( peripheralAtomIndex > -1 && centralAtomIndex > -1 ){
float distance = sqrtf(
(posArray[centralAtomIndex] - posArray[peripheralAtomIndex] )*(posArray[centralAtomIndex] - posArray[peripheralAtomIndex] ) +
(posArray[centralAtomIndex+1] - posArray[peripheralAtomIndex+1])*(posArray[centralAtomIndex+1] - posArray[peripheralAtomIndex+1]) +
(posArray[centralAtomIndex+2] - posArray[peripheralAtomIndex+2])*(posArray[centralAtomIndex+2] - posArray[peripheralAtomIndex+2]) );
float diff = fabsf( distance - cluster._distance );
if( diff > tolerance && (errors++ < maxErrorToPrint || maxErrorToPrint < 0) ){
char value[1024];
#ifdef WIN32
sprintf_s( value, 1024, "Error: Atom [%6d %6d] d[%16.7e %16.7e] diff=%16.7e [%16.7e %16.7e %16.7e] [%16.7e %16.7e %16.7e]\n",
centralAtomIndex/3, peripheralAtomIndex/3, distance, cluster._distance, diff,
posArray[centralAtomIndex], posArray[centralAtomIndex+1], posArray[centralAtomIndex+2],
posArray[peripheralAtomIndex], posArray[peripheralAtomIndex+1], posArray[peripheralAtomIndex+2] );
#else
sprintf( value, "Error: Atom [%6d %6d] diff=%16.7e d[%16.7e %16.7e] [%16.7e %16.7e %16.7e] [%16.7e %16.7e %16.7e]\n",
centralAtomIndex/3, peripheralAtomIndex/3, diff, distance, cluster._distance,
posArray[centralAtomIndex][0], posArray[centralAtomIndex][1], posArray[centralAtomIndex][2],
posArray[peripheralAtomIndex][0], posArray[peripheralAtomIndex][1], posArray[peripheralAtomIndex][2] );
#endif
message << value;
}
if( diff > maxDiff ){
maxDiffCentralIndex = centralAtomIndex/3;
maxDiffPeripheralIndex = peripheralAtomIndex/3;
maxDiff = diff;
}
}
}
}
// report findings
if( errors && log ){
(void) fprintf( log, "Shake errors=%d tolerance=%.3e max diff=%.3e atoms[%d %d]", errors, tolerance, maxDiff, maxDiffCentralIndex, maxDiffPeripheralIndex );
if( errors >= maxErrorToPrint ){
(void) fprintf( log, " only printing first %d errors", maxErrorToPrint );
}
(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 );
}
return errors;
}
/* /*
* Get contents of object * Get contents of object
* *
...@@ -710,7 +774,23 @@ std::string BrookShakeAlgorithm::getContentsString( int level ) const { ...@@ -710,7 +774,23 @@ std::string BrookShakeAlgorithm::getContentsString( int level ) const {
(void) LOCAL_SPRINTF( value, "%d", getShakeConstraintStreamSize() ); (void) LOCAL_SPRINTF( value, "%d", getShakeConstraintStreamSize() );
message << _getLine( tab, "Constraint stream size:", value ); message << _getLine( tab, "Constraint stream size:", value );
message << _getLine( tab, "Constraints:", " " );
int index = 0;
for( map<int, ShakeCluster>::const_iterator iter = _clusters.begin(); iter != _clusters.end(); ++iter ){
char buffer[1024];
const ShakeCluster& cluster = iter->second;
#ifdef WIN32
sprintf_s( buffer, 1024, "%5d %6d [%6d %6d %6d] d=%.3f m[%12.5f %12.5f]", index++, cluster._centralID,
cluster._peripheralID[0], cluster._peripheralID[1], cluster._peripheralID[2],
cluster._distance, cluster._centralInvMass, cluster._peripheralInvMass );
#endif
message << _getLine( tab, " ", buffer );
}
message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) ); message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) );
message << _getLine( tab, "ParticleIndices:", (getShakeParticleIndicesStream() ? Set : NotSet) ); message << _getLine( tab, "ParticleIndices:", (getShakeParticleIndicesStream() ? Set : NotSet) );
......
platforms/brook/src/BrookShakeAlgorithm.h
View file @ ce8d2142
...@@ -34,13 +34,50 @@ ...@@ -34,13 +34,50 @@
#include <vector> #include <vector>
#include <set> #include <set>
#include <sstream>
#include "BrookFloatStreamInternal.h" #include "BrookFloatStreamInternal.h"
#include "BrookPlatform.h" #include "BrookPlatform.h"
#include "BrookCommon.h" #include "BrookCommon.h"
#include "OpenMMException.h"
namespace OpenMM { namespace OpenMM {
struct ShakeCluster {
int _centralID;
int _peripheralID[3];
int _size;
float _distance;
float _centralInvMass, _peripheralInvMass;
ShakeCluster(){};
ShakeCluster( int centralID, float invMass) : _centralID(centralID), _centralInvMass(invMass), _size(0) {
_peripheralID[0] = _peripheralID[1] = _peripheralID[2] = -1;
};
void addAtom( int id, float dist, float invMass){
if( _size == 3 ){
std::stringstream message;
message << "ShakeCluster::addAtom: " << "atom " << id << " has more than 3 constraints!." << std::endl;
throw OpenMMException( message.str() );
}
if( _size > 0 && dist != _distance ){
std::stringstream message;
message << "ShakeCluster::addAtom: " << "atom " << id << " has different constraint distances: " << dist << " and " << _distance << std::endl;
throw OpenMMException( message.str() );
}
if( _size > 0 && invMass != _peripheralInvMass ){
std::stringstream message;
message << "ShakeCluster::addAtom: " << " constrainted atoms associated w/ atom " << id << " have different masses: " << invMass << " and " << _peripheralInvMass << std::endl;
throw OpenMMException( message.str() );
}
_peripheralID[_size++] = id;
_distance = dist;
_peripheralInvMass = invMass;
}
};
/** /**
* *
* Encapsulates stochastic dynamics algorithm * Encapsulates stochastic dynamics algorithm
...@@ -260,6 +297,19 @@ class BrookShakeAlgorithm : public BrookCommon { ...@@ -260,6 +297,19 @@ class BrookShakeAlgorithm : public BrookCommon {
BrookFloatStreamInternal* getShakeInverseMapStream( void ) const; BrookFloatStreamInternal* getShakeInverseMapStream( void ) const;
/*
* Check constraints
*
* @param positions atom positions
* @param log file to print to (can be NULL)
* @param tolerance tolerance to compare (if < 0, then use algorithm tolerance
*
* @return number of errors
*
*/
int checkConstraints( BrookStreamInternal* positions, FILE* log, float tolerance ) const;
private: private:
// streams indices // streams indices
...@@ -307,6 +357,10 @@ class BrookShakeAlgorithm : public BrookCommon { ...@@ -307,6 +357,10 @@ class BrookShakeAlgorithm : public BrookCommon {
BrookFloatStreamInternal* _shakeStreams[LastStreamIndex]; BrookFloatStreamInternal* _shakeStreams[LastStreamIndex];
// Shake cluster (used for debugging)
std::map<int, ShakeCluster> _clusters;
/* /*
* Setup of stream dimensions * Setup of stream dimensions
* *
......
platforms/brook/src/BrookStreamInternal.cpp
View file @ ce8d2142
...@@ -293,7 +293,7 @@ std::string BrookStreamInternal::_getLine( const std::string& tab, ...@@ -293,7 +293,7 @@ std::string BrookStreamInternal::_getLine( const std::string& tab,
* *
* */ * */
int BrookStreamInternal::printToFile( FILE* log ){ int BrookStreamInternal::printToFile( FILE* log, int maxPrint ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -307,7 +307,7 @@ int BrookStreamInternal::printToFile( FILE* log ){ ...@@ -307,7 +307,7 @@ int BrookStreamInternal::printToFile( FILE* log ){
std::string contents = getContentsString(); std::string contents = getContentsString();
(void) fprintf( log, "%s\n", contents.c_str() ); (void) fprintf( log, "%s\n", contents.c_str() );
_bodyPrintToFile( log ); _bodyPrintToFile( log, maxPrint );
return DefaultReturnValue; return DefaultReturnValue;
...@@ -326,7 +326,7 @@ const std::string BrookStreamInternal::getContentsString( int level ) const { ...@@ -326,7 +326,7 @@ const std::string BrookStreamInternal::getContentsString( int level ) const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const std::string methodName = "BrookIntStreamInternal::getContentsString"; // static const std::string methodName = "BrookStreamInternal::getContentsString";
static const unsigned int MAX_LINE_CHARS = 256; static const unsigned int MAX_LINE_CHARS = 256;
char value[MAX_LINE_CHARS]; char value[MAX_LINE_CHARS];
...@@ -363,3 +363,88 @@ const std::string BrookStreamInternal::getContentsString( int level ) const { ...@@ -363,3 +363,88 @@ const std::string BrookStreamInternal::getContentsString( int level ) const {
return message.str(); return message.str();
} }
/*
* Get stats
*
* @return statistics vector
*
* */
int BrookStreamInternal::getStatistics( std::vector<std::vector<double> >& statistics, int maxScan ){
return 0;
}
/*
* Get stat string
*
* @param tag id tag
* @param statistics stat vector
* @return stat string
*
* */
std::string BrookStreamInternal::printStatistics( std::string tag, std::vector<std::vector<double> >& statistics ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookStreamInternal::getContentsString";
static const unsigned int MAX_LINE_CHARS = 256;
char value[MAX_LINE_CHARS];
std::string tab = " ";
static int initialized = 0;
static std::vector<std::string> header;
// ---------------------------------------------------------------------------------------
// row headers
if( !initialized ){
initialized = 1;
header.push_back( "Average" );
header.push_back( "StdDev" );
header.push_back( "|Average|" );
header.push_back( "RawSum" );
header.push_back( "Min" );
header.push_back( "MinIndex" );
header.push_back( "Max" );
header.push_back( "MaxIndex" );
header.push_back( "Count" );
}
#ifdef WIN32
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );
#else
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );
#endif
std::stringstream message;
// loop over rows (Average, StdDev, ... )
// building message string: tag + header + tab + values
for( unsigned int ii = 0; ii < statistics.size(); ii++ ){
std::vector<double> row = statistics[ii];
std::stringstream head;
std::stringstream valueString;
head << tag << " ";
if( ii < header.size() ){
head << header[ii];
}
head << " ";
valueString << "[ ";
for( unsigned int jj = 0; jj < row.size(); jj++ ){
(void) LOCAL_SPRINTF( value, "%16.7e ", row[jj] );
valueString << value;
}
valueString << "]";
message << _getLine( tab, head.str(), valueString.str() );
}
return message.str();
}
platforms/brook/src/BrookStreamInternal.h
View file @ ce8d2142
...@@ -231,7 +231,7 @@ class BrookStreamInternal { ...@@ -231,7 +231,7 @@ class BrookStreamInternal {
* *
* */ * */
int printToFile( FILE* log ); int printToFile( FILE* log, int maxPrint = -1 );
/** /**
* Sum over stream dimensions * Sum over stream dimensions
...@@ -245,6 +245,27 @@ class BrookStreamInternal { ...@@ -245,6 +245,27 @@ class BrookStreamInternal {
*/ */
int sumByDimension( int stopIndex, double* sum ); int sumByDimension( int stopIndex, double* sum );
/*
* Get stats
*
* @return statistics vector
*
* */
virtual int getStatistics( std::vector<std::vector<double>>&, int maxScan );
/*
* Get stat string
*
* @param tag id tag
* @param statistics stat vector
*
* @return stat string
*
* */
std::string printStatistics( std::string tag, std::vector<std::vector<double> >& statistics ) const;
protected: protected:
...@@ -284,7 +305,7 @@ class BrookStreamInternal { ...@@ -284,7 +305,7 @@ class BrookStreamInternal {
* *
* */ * */
virtual int _bodyPrintToFile( FILE* log ) = 0; virtual int _bodyPrintToFile( FILE* log, int maxPrint ) = 0;
}; };
} // namespace OpenMM } // namespace OpenMM
......
platforms/brook/src/BrookVerletDynamics.cpp
View file @ ce8d2142
...@@ -60,6 +60,7 @@ BrookVerletDynamics::BrookVerletDynamics( ){ ...@@ -60,6 +60,7 @@ BrookVerletDynamics::BrookVerletDynamics( ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_numberOfParticles = -1; _numberOfParticles = -1;
_internalStepCount = 0;
// mark stream dimension variables as unset // mark stream dimension variables as unset
...@@ -75,7 +76,7 @@ BrookVerletDynamics::BrookVerletDynamics( ){ ...@@ -75,7 +76,7 @@ BrookVerletDynamics::BrookVerletDynamics( ){
// setup inverse sqrt masses // setup inverse sqrt masses
_inverseMasses = NULL; _inverseMasses = NULL;
} }
...@@ -466,25 +467,31 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -466,25 +467,31 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static std::string methodName = "\nBrookVerletDynamics::update"; static std::string methodName = "\nBrookVerletDynamics::update";
static const int PrintOn = 0; static int printOn = 1;
FILE* log;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_internalStepCount++;
setLog( stderr );
printOn = (printOn && getLog()) ? printOn : 0;
BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC); BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
if( (1 || PrintOn) && getLog() ){ if( (1 || printOn) ){
static int showAux = 1; static int showAux = 1;
log = getLog();
if( showAux ){ if( showAux ){
showAux = 0; showAux = 0;
/* /*
std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( ); std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
(void) fprintf( getLog(), "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() ); (void) fprintf( log, "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
*/ */
(void) fprintf( getLog(), "%s Shake contents\n%s", methodName.c_str(), brookShakeAlgorithm.getContentsString().c_str() ); (void) fprintf( log, "%s step=%d Shake contents\n%s", methodName.c_str(), _internalStepCount, brookShakeAlgorithm.getContentsString().c_str() );
(void) fflush( getLog() ); (void) fflush( log );
} }
} }
...@@ -505,27 +512,41 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -505,27 +512,41 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
// diagnostics // diagnostics
if( PrintOn && getLog() ){ if( printOn ){
(void) fprintf( getLog(), "\n%s Post kupdate_md1: particleStrW=%3d step=%.5f", (void) fprintf( log, "\n%s step=%d Post kupdate_md1: particleStrW=%3d step=%.5f",
methodName.c_str(), getVerletDynamicsParticleStreamWidth(), getStepSize() ); methodName.c_str(), _internalStepCount, getVerletDynamicsParticleStreamWidth(), getStepSize() );
(void) fprintf( getLog(), "\nInverseMassStream\n" ); (void) fprintf( log, "\nInverseMassStream %d\n", _internalStepCount );
getInverseMassStream()->printToFile( getLog() ); getInverseMassStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nForceStream\n" ); (void) fprintf( log, "\nForceStream %d\n", _internalStepCount );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() ); std::vector<std::vector<double> > forceStatistics;
brookStreamInternalF->getStatistics( forceStatistics, getNumberOfParticles() );
std::stringstream tag;
tag << _internalStepCount << " Fxx ";
std::string stats = brookStreamInternalF->printStatistics( tag.str(), forceStatistics );
(void) fprintf( log, "\nStep %d Force stats:\n%s", _internalStepCount, stats.c_str() );
brookStreamInternalF->printToFile( log );
BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" ); std::vector<std::vector<double> > velocityStatistics;
brookStreamInternalV->printToFile( getLog() ); brookStreamInternalV->getStatistics( velocityStatistics, getNumberOfParticles() );
std::stringstream tagV;
tagV << _internalStepCount << " Vxx ";
stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
(void) fprintf( log, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
(void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
brookStreamInternalV->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
} }
// Shake // Shake
...@@ -543,34 +564,34 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -543,34 +564,34 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() ); brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
if( (0|| PrintOn) && getLog() ){ if( printOn ){
(void) fprintf( getLog(), "\n%s Post kshakeh_fix1: particleStrW=%3d", methodName.c_str(), getVerletDynamicsParticleStreamWidth() ); (void) fprintf( log, "\n%s Post kshakeh_fix1: particleStrW=%3d", methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeParticleIndicesStream\n" ); (void) fprintf( log, "\nShakeParticleIndicesStream %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeParticleIndicesStream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeParticleIndicesStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeParticleParameterStream\n" ); (void) fprintf( log, "\nShakeParticleParameterStream %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeParticleParameterStream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeParticleParameterStream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons0\n" ); (void) fprintf( log, "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons1\n" ); (void) fprintf( log, "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons2\n" ); (void) fprintf( log, "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons3\n" ); (void) fprintf( log, "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( log );
} }
...@@ -587,32 +608,32 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -587,32 +608,32 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
getXPrimeStream()->getBrookStream() ); getXPrimeStream()->getBrookStream() );
//positionStream.getBrookStream() ); //positionStream.getBrookStream() );
if( ( 0 || PrintOn) && getLog() ){ if( 0 && printOn ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d", (void) fprintf( log, "\n%s Post kshakeh_update2_fix1: particleStrW=%3d",
methodName.c_str(), getVerletDynamicsParticleStreamWidth() ); methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeInverseMapStream\n" ); (void) fprintf( log, "\nShakeInverseMapStream %d\n", _internalStepCount );
brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( log );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons0\n" ); (void) fprintf( log, "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons1\n" ); (void) fprintf( log, "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons2\n" ); (void) fprintf( log, "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( log );
(void) fprintf( getLog(), "\nShakeXCons3\n" ); (void) fprintf( log, "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() ); brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( log );
} }
...@@ -625,24 +646,25 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -625,24 +646,25 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
velocityStream.getBrookStream(), velocityStream.getBrookStream(),
positionStream.getBrookStream() ); positionStream.getBrookStream() );
if( ( 0 || PrintOn) && getLog() ){ if( printOn ){
(void) fprintf( getLog(), "\n%s Post kupdate_md2: inverseStepSize=%3e", (void) fprintf( log, "\n%s step=%d Post kupdate_md2: inverseStepSize=%3e",
methodName.c_str(), inverseStepSize ); methodName.c_str(), _internalStepCount, inverseStepSize );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" ); (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
(void) fprintf( getLog(), "\nXPrimeStream\n" ); (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
getXPrimeStream()->printToFile( getLog() ); getXPrimeStream()->printToFile( log );
brookStreamInternalPos = velocityStream.getBrookStreamImpl(); brookStreamInternalPos = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" ); (void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
brookStreamInternalPos->printToFile( getLog() ); brookStreamInternalPos->printToFile( log );
} }
} else { } else {
kupdateMdNoShake( getStepSize(), kupdateMdNoShake( getStepSize(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
velocityStream.getBrookStream(), velocityStream.getBrookStream(),
......
platforms/brook/src/BrookVerletDynamics.h
View file @ ce8d2142
...@@ -204,6 +204,10 @@ class BrookVerletDynamics : public BrookCommon { ...@@ -204,6 +204,10 @@ class BrookVerletDynamics : public BrookCommon {
LastStreamIndex LastStreamIndex
}; };
// track step (diagnostics)
int _internalStepCount;
BrookOpenMMFloat _stepSize; BrookOpenMMFloat _stepSize;
// Particle stream dimensions // Particle stream dimensions
......
platforms/brook/src/OpenMMBrookInterface.cpp
View file @ ce8d2142
...@@ -538,7 +538,11 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){ ...@@ -538,7 +538,11 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
//setLog( stderr ); (void) fprintf( stderr, "%s done nonbonded & bonded=%d completed=%d\n", methodName.c_str(),
_brookBonded.isActive(), _brookBonded.isSetupCompleted() );
(void) fflush( stderr );
setLog( stderr );
printOn = (printOn && getLog()) ? printOn : 0; printOn = (printOn && getLog()) ? printOn : 0;
BrookStreamImpl* positions = getParticlePositions(); BrookStreamImpl* positions = getParticlePositions();
...@@ -547,10 +551,14 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){ ...@@ -547,10 +551,14 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
// info // info
//if( printOn > 1 ){ //if( printOn > 1 ){
if( 1 ){ if( 0 ){
printForcesToFile( context ); printForcesToFile( context );
} }
(void) fprintf( stderr, "%s x1 done nonbonded & bonded=%d completed=%d\n", methodName.c_str(),
_brookBonded.isActive(), _brookBonded.isSetupCompleted() );
(void) fflush( stderr );
// nonbonded forces // nonbonded forces
if( _brookNonBonded.isActive() ){ if( _brookNonBonded.isActive() ){
...@@ -579,11 +587,9 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){ ...@@ -579,11 +587,9 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
getParticleStreamWidth(), getParticleStreamSize() ); getParticleStreamWidth(), getParticleStreamSize() );
} }
/* (void) fprintf( stderr, "%s done nonbonded & bonded=%d completed=%d\n", methodName.c_str(),
(void) fprintf( getLog(), "%s done nonbonded: bonded=%d completed=%d\n", methodName.c_str(),
_brookBonded.isActive(), _brookBonded.isSetupCompleted() ); _brookBonded.isActive(), _brookBonded.isSetupCompleted() );
(void) fflush( getLog() ); (void) fflush( stderr );
*/
_brookBonded.computeForces( *positions, *forces ); _brookBonded.computeForces( *positions, *forces );
...@@ -622,6 +628,9 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){ ...@@ -622,6 +628,9 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
_brookGbsa.computeForces( *positions, *forces ); _brookGbsa.computeForces( *positions, *forces );
} }
(void) fprintf( stderr, "%s done computeForces\n", methodName.c_str() );
(void) fflush( stderr );
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
} }
...@@ -644,7 +653,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){ ...@@ -644,7 +653,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
// first step only? // first step only?
if( step > 0 ){ if( step > 20 ){
return; return;
} }
...@@ -734,6 +743,8 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){ ...@@ -734,6 +743,8 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
printStreamsToFile( fileName, streams ); printStreamsToFile( fileName, streams );
} }
(void) fprintf( stderr, "%s done computeForces\n", methodName.c_str() );
(void) fflush( stderr );
forces->fillWithValue( &zero ); forces->fillWithValue( &zero );
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -822,6 +833,7 @@ int OpenMMBrookInterface::printStreamsToFile( std::string fileName, std::vector< ...@@ -822,6 +833,7 @@ int OpenMMBrookInterface::printStreamsToFile( std::string fileName, std::vector<
minIndex = stream->getStreamSize(); minIndex = stream->getStreamSize();
} }
} }
minIndex = minIndex > getNumberOfParticles() ? getNumberOfParticles() : minIndex;
// sum columns // sum columns
......
platforms/brook/src/gpu/kforce.h
View file @ ce8d2142
...@@ -60,23 +60,6 @@ void knbforce_CDLJ (const float natoms, ...@@ -60,23 +60,6 @@ void knbforce_CDLJ (const float natoms,
const float jstrwidth, const float jstrwidth,
const float fstrwidth, const float fstrwidth,
const float epsfac, const float epsfac,
const float4 params,
::brook::stream posq,
::brook::stream isigeps,
::brook::stream sigma,
::brook::stream epsilon,
::brook::stream excl,
::brook::stream outforce1,
::brook::stream outforce2);
void knbforce_CDLJ2(const float natoms,
const float dupfac,
const float strheight,
const float strwidth,
const float jstrwidth,
const float fstrwidth,
const float epsfac,
const float4 params,
::brook::stream posq, ::brook::stream posq,
::brook::stream isigeps, ::brook::stream isigeps,
::brook::stream sigma, ::brook::stream sigma,
...@@ -94,7 +77,6 @@ void knbforce_CDLJ4( ...@@ -94,7 +77,6 @@ void knbforce_CDLJ4(
const float jstrwidth, const float jstrwidth,
const float fstrwidth, const float fstrwidth,
const float epsfac, const float epsfac,
const float4 params,
::brook::stream posq, ::brook::stream posq,
::brook::stream charge, ::brook::stream charge,
::brook::stream isigeps, ::brook::stream isigeps,
...@@ -106,108 +88,6 @@ void knbforce_CDLJ4( ...@@ -106,108 +88,6 @@ void knbforce_CDLJ4(
::brook::stream outforce3, ::brook::stream outforce3,
::brook::stream outforce4); ::brook::stream outforce4);
void knbforce_CDLJ4Debug(
const float natoms,
const float nAtomsCeiling,
const float dupfac,
const float strheight,
const float strwidth,
const float jstrwidth,
const float fstrwidth,
const float epsfac,
const float4 params,
::brook::stream posq,
::brook::stream isigeps,
::brook::stream sigma,
::brook::stream epsilon,
::brook::stream excl,
::brook::stream outforce1,
::brook::stream outforce2,
::brook::stream outforce3,
::brook::stream outforce4);
void knbforce_CDLJ_1(
const float natoms,
const float nAtomsCeiling,
const float strwidth,
const float epsfac,
::brook::stream posq,
::brook::stream isigeps,
::brook::stream excl,
::brook::stream outforce );
void knbforce_CDLJ4NoEx(
const float natoms,
const float nAtomsCeiling,
const float dupfac,
const float strheight,
const float strwidth,
const float jstrwidth,
const float fstrwidth,
const float epsfac,
const float4 params,
::brook::stream posq,
::brook::stream charge,
::brook::stream isigeps,
::brook::stream sigma,
::brook::stream epsilon,
::brook::stream excl,
::brook::stream outforce1,
::brook::stream outforce2,
::brook::stream outforce3,
::brook::stream outforce4);
void knbforce_LDLJ (const float natoms,
const float dupfac,
const float strheight,
const float strwidth,
const float jstrwidth,
const float fstrwidth,
const float epsfac,
const float4 params,
const __BRTIter& a_iatom2,
::brook::stream posq,
::brook::stream isigeps,
::brook::stream sigma,
::brook::stream epsilon,
::brook::stream excl,
::brook::stream outforce1,
::brook::stream outforce2);
void knbforce_SFDLJ (const float natoms,
const float dupfac,
const float strheight,
const float strwidth,
const float jstrwidth,
const float fstrwidth,
const float epsfac,
const float4 params,
const __BRTIter& a_iatom2,
::brook::stream posq,
::brook::stream isigeps,
::brook::stream sigma,
::brook::stream epsilon,
::brook::stream excl,
::brook::stream outforce1,
::brook::stream outforce2);
void knbforce_SHEFFIELD (const float natoms,
const float dupfac,
const float strheight,
const float strwidth,
const float jstrwidth,
const float fstrwidth,
const float epsfac,
const float4 params,
const __BRTIter& a_iatom2,
::brook::stream posq,
::brook::stream isigeps,
::brook::stream sigma,
::brook::stream epsilon,
::brook::stream excl,
::brook::stream outforce1,
::brook::stream outforce2);
void kmerge_partial_forces (const float repfac, void kmerge_partial_forces (const float repfac,
const float atomStrWidth, const float atomStrWidth,
const float pforceStrWidth, const float pforceStrWidth,
......
platforms/brook/src/gpu/kforce_CDLJ.br
View file @ ce8d2142
...@@ -29,7 +29,7 @@ ...@@ -29,7 +29,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
kernel float4 scalar_force_CDLJ( float4 qq, float epsfac, float4 sig, float4 eps, float4 r2, float4 params ) { kernel float4 scalar_force_CDLJ( float4 qq, float epsfac, float4 sig, float4 eps, float4 r2 ) {
float4 invr, invrsig2, invrsig6; float4 invr, invrsig2, invrsig6;
float4 f; float4 f;
...@@ -196,7 +196,7 @@ kernel void knbforce_CDLJ ( ...@@ -196,7 +196,7 @@ kernel void knbforce_CDLJ (
r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f; r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params ); fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2 );
outforce1 += fs.x * d1; outforce1 += fs.x * d1;
outforce1 += fs.y * d2; outforce1 += fs.y * d2;
...@@ -233,7 +233,7 @@ if( exclmask.w > 0.5f ){ ...@@ -233,7 +233,7 @@ if( exclmask.w > 0.5f ){
d4 = ipos2 - jpos4; d4 = ipos2 - jpos4;
r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f; r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params ); fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2 );
outforce2 += fs.x * d1; outforce2 += fs.x * d1;
outforce2 += fs.y * d2; outforce2 += fs.y * d2;
...@@ -280,7 +280,6 @@ kernel void knbforce_CDLJ4( ...@@ -280,7 +280,6 @@ kernel void knbforce_CDLJ4(
float jstrwidth, float jstrwidth,
float fstrwidth, float fstrwidth,
float epsfac, float epsfac,
float4 params,
float3 posq[][], float3 posq[][],
float charge[][], float charge[][],
float4 isigeps[][], float4 isigeps[][],
...@@ -446,7 +445,7 @@ kernel void knbforce_CDLJ4( ...@@ -446,7 +445,7 @@ kernel void knbforce_CDLJ4(
d4 = ipos1 - jpos4; d4 = ipos1 - jpos4;
r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier); r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params ); fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2 );
outforce1 += fs.x * d1; outforce1 += fs.x * d1;
outforce1 += fs.y * d2; outforce1 += fs.y * d2;
...@@ -470,7 +469,7 @@ kernel void knbforce_CDLJ4( ...@@ -470,7 +469,7 @@ kernel void knbforce_CDLJ4(
d4 = ipos2 - jpos4; d4 = ipos2 - jpos4;
r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier); r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params ); fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2 );
outforce2 += fs.x * d1; outforce2 += fs.x * d1;
outforce2 += fs.y * d2; outforce2 += fs.y * d2;
...@@ -494,7 +493,7 @@ kernel void knbforce_CDLJ4( ...@@ -494,7 +493,7 @@ kernel void knbforce_CDLJ4(
d4 = ipos3 - jpos4; d4 = ipos3 - jpos4;
r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier); r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params ); fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2 );
outforce3 += fs.x * d1; outforce3 += fs.x * d1;
outforce3 += fs.y * d2; outforce3 += fs.y * d2;
...@@ -518,7 +517,7 @@ kernel void knbforce_CDLJ4( ...@@ -518,7 +517,7 @@ kernel void knbforce_CDLJ4(
d4 = ipos4 - jpos4; d4 = ipos4 - jpos4;
r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier); r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params ); fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2 );
outforce4 += fs.x * d1; outforce4 += fs.x * d1;
outforce4 += fs.y * d2; outforce4 += fs.y * d2;
......
platforms/brook/src/gpu/kshakeh.br
View file @ ce8d2142
...@@ -134,7 +134,8 @@ kshakeh_fix1( ...@@ -134,7 +134,8 @@ kshakeh_fix1(
i = 0.0f; i = 0.0f;
converged = 1.0f; converged = 1.0f;
tolerance = 2.0f*inputTolerance; //tolerance = 2.0f*inputTolerance;
tolerance = inputTolerance;
while( i < maxIterations && converged > 0.0f ){ while( i < maxIterations && converged > 0.0f ){
//First hydrogen //First hydrogen
......
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