"platforms/vscode:/vscode.git/clone" did not exist on "637de3d3a6ccd61a0c207403b23e1016466c3f59"
Commit ff8bd4e1 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Brook SD code

parent a475f0a6
...@@ -89,6 +89,10 @@ const std::string BrookCommon::ShakeXCons2Stream ...@@ -89,6 +89,10 @@ const std::string BrookCommon::ShakeXCons2Stream
const std::string BrookCommon::ShakeXCons3Stream = "ShakeXCons3Stream"; const std::string BrookCommon::ShakeXCons3Stream = "ShakeXCons3Stream";
const std::string BrookCommon::ShakeInverseMapStream = "ShakeInverseMapStream"; const std::string BrookCommon::ShakeInverseMapStream = "ShakeInverseMapStream";
// Random number streams
const std::string BrookCommon::ShuffleStream = "ShuffleStream";
/** /**
* Constructor * Constructor
* *
......
...@@ -104,6 +104,10 @@ class BrookCommon { ...@@ -104,6 +104,10 @@ class BrookCommon {
static const std::string ShakeXCons3Stream; static const std::string ShakeXCons3Stream;
static const std::string ShakeInverseMapStream; static const std::string ShakeInverseMapStream;
// Random number generator streams
static const std::string ShuffleStream;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
/** /**
...@@ -224,7 +228,7 @@ class BrookCommon { ...@@ -224,7 +228,7 @@ class BrookCommon {
FILE* getLog( void ) const; FILE* getLog( void ) const;
/** /**
* Be verbose flag * Get verbose flag
* *
* @return verbosity flag * @return verbosity flag
* *
...@@ -239,7 +243,7 @@ class BrookCommon { ...@@ -239,7 +243,7 @@ class BrookCommon {
* *
*/ */
int setVerbosity( int verbositym); int setVerbosity( int verbosity );
/* /*
* Given number of stream elements and width, returns the appropriate * Given number of stream elements and width, returns the appropriate
......
...@@ -46,6 +46,7 @@ BrookIntegrateLangevinStepKernel::BrookIntegrateLangevinStepKernel( std::string ...@@ -46,6 +46,7 @@ BrookIntegrateLangevinStepKernel::BrookIntegrateLangevinStepKernel( std::string
_brookStochasticDynamics = NULL; _brookStochasticDynamics = NULL;
_brookShakeAlgorithm = NULL; _brookShakeAlgorithm = NULL;
_brookRandomNumberGenerator = NULL;
} }
...@@ -59,6 +60,7 @@ BrookIntegrateLangevinStepKernel::~BrookIntegrateLangevinStepKernel( ){ ...@@ -59,6 +60,7 @@ BrookIntegrateLangevinStepKernel::~BrookIntegrateLangevinStepKernel( ){
delete _brookStochasticDynamics; delete _brookStochasticDynamics;
delete _brookShakeAlgorithm; delete _brookShakeAlgorithm;
delete _brookRandomNumberGenerator;
} }
...@@ -77,6 +79,9 @@ void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses, ...@@ -77,6 +79,9 @@ void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses,
_brookShakeAlgorithm = new BrookShakeAlgorithm( ); _brookShakeAlgorithm = new BrookShakeAlgorithm( );
_brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() ); _brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() );
_brookRandomNumberGenerator = new BrookRandomNumberGenerator( );
_brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
} }
/** /**
...@@ -116,6 +121,7 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc ...@@ -116,6 +121,7 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
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 ){
_brookStochasticDynamics->updateParameters( temperature, friction, stepSize ); _brookStochasticDynamics->updateParameters( temperature, friction, stepSize );
} }
_brookStochasticDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm ); _brookStochasticDynamics->update( positions, velocities, forces, (_brookShakeAlgorithm ? *_brookShakeAlgorithm : NULL),
*_brookRandomNumberGenerator );
} }
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include "kernels.h" #include "kernels.h"
#include "BrookStochasticDynamics.h" #include "BrookStochasticDynamics.h"
#include "BrookShakeAlgorithm.h" #include "BrookShakeAlgorithm.h"
#include "BrookRandomNumberGenerator.h"
namespace OpenMM { namespace OpenMM {
...@@ -91,6 +92,7 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel { ...@@ -91,6 +92,7 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
BrookStochasticDynamics* _brookStochasticDynamics; BrookStochasticDynamics* _brookStochasticDynamics;
BrookShakeAlgorithm* _brookShakeAlgorithm; BrookShakeAlgorithm* _brookShakeAlgorithm;
BrookRandomNumberGenerator* _brookRandomNumberGenerator;
}; };
......
...@@ -145,8 +145,8 @@ void BrookPlatform::_initializeKernelFactory( void ){ ...@@ -145,8 +145,8 @@ void BrookPlatform::_initializeKernelFactory( void ){
registerKernelFactory( CalcStandardMMForceFieldKernel::Name(), factory); registerKernelFactory( CalcStandardMMForceFieldKernel::Name(), factory);
// registerKernelFactory( CalcGBSAOBCForceFieldKernel::Name(), factory); // registerKernelFactory( CalcGBSAOBCForceFieldKernel::Name(), factory);
registerKernelFactory( IntegrateVerletStepKernel::Name(), factory); //registerKernelFactory( IntegrateVerletStepKernel::Name(), factory);
//registerKernelFactory( IntegrateLangevinStepKernel::Name(), factory); registerKernelFactory( IntegrateLangevinStepKernel::Name(), factory);
//registerKernelFactory( IntegrateBrownianStepKernel::Name(), factory); //registerKernelFactory( IntegrateBrownianStepKernel::Name(), factory);
//registerKernelFactory( ApplyAndersenThermostatKernel::Name(), factory); //registerKernelFactory( ApplyAndersenThermostatKernel::Name(), factory);
registerKernelFactory( CalcKineticEnergyKernel::Name(), factory); registerKernelFactory( CalcKineticEnergyKernel::Name(), factory);
...@@ -180,7 +180,7 @@ void BrookPlatform::_setBrookRuntime( const std::string& runtime ){ ...@@ -180,7 +180,7 @@ void BrookPlatform::_setBrookRuntime( const std::string& runtime ){
} }
if( getLog() ){ if( getLog() ){
(void) fprintf( getLog(), "%s Brook initializing to runtime=<%s>", methodName.c_str(), _runtime.c_str() ); (void) fprintf( getLog(), "%s Brook initializing to runtime=<%s>\n", methodName.c_str(), _runtime.c_str() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
......
working on shuffleGVStreams
/* -------------------------------------------------------------------------- * /* -------------------------------------------------------------------------- *
* OpenMM * * OpenMM *
* -------------------------------------------------------------------------- * * -------------------------------------------------------------------------- *
...@@ -32,9 +31,8 @@ working on shuffleGVStreams ...@@ -32,9 +31,8 @@ working on shuffleGVStreams
#include <sstream> #include <sstream>
#include "BrookRandomNumberGenerator.h" #include "BrookRandomNumberGenerator.h"
#include "BrookPlatform.h"
#include "OpenMMException.h" #include "OpenMMException.h"
#include "BrookStreamImpl.h" #include "kupdatesd.h"
// use random number generator // use random number generator
...@@ -60,6 +58,7 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){ ...@@ -60,6 +58,7 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
// fixed for now // fixed for now
_numberOfRandomNumberStreams = 2; _numberOfRandomNumberStreams = 2;
_randomNumberGeneratorStreams = NULL;
// mark stream dimension variables as unset // mark stream dimension variables as unset
...@@ -67,8 +66,16 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){ ...@@ -67,8 +66,16 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
_randomNumberStreamHeight = -1; _randomNumberStreamHeight = -1;
_randomNumberStreamSize = -1; _randomNumberStreamSize = -1;
_rvStreamIndex = 0;
_rvStreamOffset = 0;
_numberOfShuffles = 0;
_maxShuffles = 100;
_loadBuffer = NULL;
_shuffleIndices = NULL;
for( int ii = 0; ii < LastStreamIndex; ii++ ){ for( int ii = 0; ii < LastStreamIndex; ii++ ){
_randomNumberStreams[ii] = NULL; _auxiliaryStreams[ii] = NULL;
} }
// set randomNumber seed // set randomNumber seed
...@@ -93,9 +100,13 @@ BrookRandomNumberGenerator::~BrookRandomNumberGenerator( ){ ...@@ -93,9 +100,13 @@ BrookRandomNumberGenerator::~BrookRandomNumberGenerator( ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
for( int ii = 0; ii < LastStreamIndex; ii++ ){ for( int ii = 0; ii < LastStreamIndex; ii++ ){
delete _sdStreams[ii]; delete _auxiliaryStreams[ii];
} }
delete[] _randomNumberGeneratorStreams;
delete[] _loadBuffer;
delete[] _shuffleIndices;
} }
/** /**
...@@ -110,15 +121,27 @@ int BrookRandomNumberGenerator::getNumberOfRandomNumberStreams( void ) const { ...@@ -110,15 +121,27 @@ int BrookRandomNumberGenerator::getNumberOfRandomNumberStreams( void ) const {
} }
/** /**
* Get number of random number streams * Get random number stream width
* *
* @return number of random number streams * @return nndom number stream width
* *
*/ */
int BrookRandomNumberGenerator::getNumberOfRandomNumberStreams( void ) const { int BrookRandomNumberGenerator::getRandomNumberStreamWidth( void ) const {
return _numberOfRandomNumberStreams; return _randomNumberStreamWidth;
} }
/**
* Get random number stream height
*
* @return nndom number stream height
*
*/
int BrookRandomNumberGenerator::getRandomNumberStreamHeight( void ) const {
return _randomNumberStreamHeight;
}
/** /**
* Get random number seed * Get random number seed
* *
...@@ -150,7 +173,7 @@ unsigned long int BrookRandomNumberGenerator::incrementRandomNumberSeed( unsigne ...@@ -150,7 +173,7 @@ unsigned long int BrookRandomNumberGenerator::incrementRandomNumberSeed( unsigne
* @return random number seed * @return random number seed
*/ */
unsigned long int setRandomNumberSeed( unsigned long int seed = 1 ); unsigned long int BrookRandomNumberGenerator::setRandomNumberSeed( unsigned long int seed ){
_randomNumberSeed = seed; _randomNumberSeed = seed;
return _randomNumberSeed; return _randomNumberSeed;
} }
...@@ -164,11 +187,11 @@ unsigned long int setRandomNumberSeed( unsigned long int seed = 1 ); ...@@ -164,11 +187,11 @@ unsigned long int setRandomNumberSeed( unsigned long int seed = 1 );
* *
*/ */
BrookOpenMMFloat BrookRandomNumberGenerator::generateGromacsRandomNumber( int* ig ){ BrookOpenMMFloat BrookRandomNumberGenerator::_generateGromacsRandomNumber( unsigned long int* ig ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::generateGromacsRandomNumber"; // static const char* methodName = "\nBrookRandomNumberGenerator::_generateGromacsRandomNumber";
int irand; int irand;
...@@ -182,7 +205,8 @@ BrookOpenMMFloat BrookRandomNumberGenerator::generateGromacsRandomNumber( int* i ...@@ -182,7 +205,8 @@ BrookOpenMMFloat BrookRandomNumberGenerator::generateGromacsRandomNumber( int* i
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
irand = abs(*ig) % m; unsigned long int igg = (*ig > 0) ? *ig : -1*(*ig);
irand = igg % m;
/* multiply irand by mult, but take into account that overflow /* multiply irand by mult, but take into account that overflow
* must be discarded, and do not generate an error. * must be discarded, and do not generate an error.
...@@ -209,7 +233,7 @@ BrookOpenMMFloat BrookRandomNumberGenerator::generateGromacsRandomNumber( int* i ...@@ -209,7 +233,7 @@ BrookOpenMMFloat BrookRandomNumberGenerator::generateGromacsRandomNumber( int* i
inline int MaxInt( unsigned int x, unsigned int y ){ return x > y ? x : y; } inline int MaxInt( unsigned int x, unsigned int y ){ return x > y ? x : y; }
/** /**
* Generate a random number using algorithm in Nvidia code * Generate a random number using algorithm in Kiss code
* http://www.helsbreth.org/random/rng_kiss.html * http://www.helsbreth.org/random/rng_kiss.html
* *
* @param randomV1 output random value * @param randomV1 output random value
...@@ -219,12 +243,12 @@ inline int MaxInt( unsigned int x, unsigned int y ){ return x > y ? x : y; } ...@@ -219,12 +243,12 @@ inline int MaxInt( unsigned int x, unsigned int y ){ return x > y ? x : y; }
* *
*/ */
void BrookRandomNumberGenerator::generateRandomsAlaNvidia( float* randomV1, float* randomV2, float* randomV3, void BrookRandomNumberGenerator::_generateRandomsKiss( float* randomV1, float* randomV2, float* randomV3,
unsigned int state[4] ){ unsigned int state[4] ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::generateRandomsAlaNvidia"; // static const char* methodName = "\nBrookRandomNumberGenerator::_generateRandomsKiss";
unsigned int carry = 0; unsigned int carry = 0;
...@@ -327,22 +351,21 @@ if( x1 < 0.0f ){ ...@@ -327,22 +351,21 @@ if( x1 < 0.0f ){
} }
/** /**
* Load random number streams using Nvidia algorithm * Load random number streams using Kiss algorithm
* *
* *
* @return DefaultReturnValue; * @return DefaultReturnValue;
*/ */
int BrookRandomNumberGenerator::loadRandomNumberStreamsNvidia( void ){ int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static float *buf = NULL;
static unsigned int state[4]; static unsigned int state[4];
static int stateInitialized = 0; static int stateInitialized = 0;
static const int reseed = 5; static const int reseed = 5;
// static const char* methodName = "\nBrookRandomNumberGenerator::loadGVStreamsNvidia"; // static const char* methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsKiss";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -356,7 +379,7 @@ int BrookRandomNumberGenerator::loadRandomNumberStreamsNvidia( void ){ ...@@ -356,7 +379,7 @@ int BrookRandomNumberGenerator::loadRandomNumberStreamsNvidia( void ){
state[3] = rand(); state[3] = rand();
if( getVerbosity() && getLog() ){ if( getVerbosity() && getLog() ){
(void) fprintf( getLog(), "LoadGVStreamsNvidia: reset state seeds stateInitialized=%d reseed=%d\n", (void) fprintf( getLog(), "LoadGVStreamsKiss: reset state seeds stateInitialized=%d reseed=%d\n",
stateInitialized, reseed ); stateInitialized, reseed );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -373,28 +396,26 @@ state[3] = 27587; ...@@ -373,28 +396,26 @@ state[3] = 27587;
// allocate memory once for download of random nos. // allocate memory once for download of random nos.
if( buf == NULL ){ float* loadBuffer = _getLoadBuffer();
buf = (float*) malloc( sizeof(float) * 3 * getRandomNumberStreamSize() );
}
if( getVerbosity() && getLog() ){ if( getVerbosity() && getLog() ){
static float count = 0.0f; static float count = 0.0f;
float block = (float) (3*sdp->gvSize); float block = (float) (3*getRandomNumberStreamSize() );
count += 1.0f; count += 1.0f;
(void) fprintf( getLog(), "LoadGVStreamsNvidia: count=%.1f ttl=%.3e no./count=%.1f %d %d\n", (void) fprintf( getLog(), "LoadGVStreamsKiss: count=%.1f ttl=%.3e no./count=%.1f %d %d\n",
count, block*count, block, sdp->gvSize, NGVSTREAMS ); count, block*count, block, getRandomNumberStreamSize(), getNumberOfRandomNumberStreams() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){ for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){ for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){
float v1,v2,v3; float v1,v2,v3;
generateRandomsAlaNvidia( &v1, &v2, &v3, state, NULL ); _generateRandomsKiss( &v1, &v2, &v3, state );
buf[ii] = v1; loadBuffer[ii] = v1;
buf[ii+1] = v2; loadBuffer[ii+1] = v2;
buf[ii+2] = v3; loadBuffer[ii+2] = v3;
} }
getRandomNumberStream( jj )->loadFromArray( buf ); getRandomNumberStream( jj )->loadFromArray( loadBuffer );
} }
return DefaultReturnValue; return DefaultReturnValue;
...@@ -407,30 +428,26 @@ state[3] = 27587; ...@@ -407,30 +428,26 @@ state[3] = 27587;
* @return DefaultReturnValue; * @return DefaultReturnValue;
*/ */
int BrookRandomNumberGenerator::loadGVStreamsOriginal( void ){ int BrookRandomNumberGenerator::_loadGVStreamsOriginal( void ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
int i, j;
static float *buf = NULL;
unsigned long int jran; unsigned long int jran;
// static const char* methodName = "\nBrookRandomNumberGenerator::LoadGVStreamsOriginal"; // static const char* methodName = "\nBrookRandomNumberGenerator::_loadGVStreamsOriginal";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
if( buf == NULL ){ float* loadBuffer = _getLoadBuffer();
buf = (float*) malloc( sizeof(float) * 3 * sdp->gvSize );
}
jran = getRandomNumberSeed(); jran = getRandomNumberSeed();
for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){ for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){ for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii++ ){
//buf[i] = sdp->fgauss( &jran ); //loadBuffer[i] = sdp->fgauss( &jran );
buf[i] = generateGromacsRandomNumber( &jran ); loadBuffer[ii] = _generateGromacsRandomNumber( &jran );
} }
getRandomNumberStream( jj )->loadFromArray( buf ); getRandomNumberStream( jj )->loadFromArray( loadBuffer );
} }
incrementRandomNumberSeed( 1 ); incrementRandomNumberSeed( 1 );
...@@ -438,6 +455,70 @@ int BrookRandomNumberGenerator::loadGVStreamsOriginal( void ){ ...@@ -438,6 +455,70 @@ int BrookRandomNumberGenerator::loadGVStreamsOriginal( void ){
return DefaultReturnValue; return DefaultReturnValue;
} }
/**
* Get load buffer
*
* @return ptr to load buffer
*
* @throw OpenMMException if rv stream size is < 1
*
**/
float* BrookRandomNumberGenerator::_getLoadBuffer( void ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookRandomNumberGenerator::_getLoadBuffer";
// ---------------------------------------------------------------------------------------
if( getRandomNumberStreamSize() < 1 ){
std::stringstream message;
message << methodName << " rv stream size=" << getRandomNumberStreamSize() << " is less than 1.";
throw OpenMMException( message.str() );
return NULL;
}
if( _loadBuffer == NULL ){
_loadBuffer = (float*) malloc( sizeof(float)*3*getRandomNumberStreamSize() );
}
return _loadBuffer;
}
/**
* Get ptr to shuffle indices
*
* @return ptr to shuffle indices
*
* @throw OpenMMException if size is < 1
*
**/
int* BrookRandomNumberGenerator::_getShuffleIndices( int size ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookRandomNumberGenerator::_getShuffleIndices";
// ---------------------------------------------------------------------------------------
if( size < 1 ){
std::stringstream message;
message << methodName << " size=" << size << " is less than 1.";
throw OpenMMException( message.str() );
return NULL;
}
if( _shuffleIndices == NULL ){
_shuffleIndices = (int*) malloc( sizeof(int)*size );
}
return _shuffleIndices;
}
/** /**
* Loads a permutation of indices from 0 to gvSize-1 in * Loads a permutation of indices from 0 to gvSize-1 in
* sdp->strShuffle. To make sure that the order of the * sdp->strShuffle. To make sure that the order of the
...@@ -456,7 +537,7 @@ int BrookRandomNumberGenerator::loadGVStreamsOriginal( void ){ ...@@ -456,7 +537,7 @@ int BrookRandomNumberGenerator::loadGVStreamsOriginal( void ){
* @return DefaultReturnValue; * @return DefaultReturnValue;
**/ **/
int BrookRandomNumberGenerator::loadGVShuffle( void ){ int BrookRandomNumberGenerator::_loadGVShuffle( void ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -464,47 +545,46 @@ int BrookRandomNumberGenerator::loadGVShuffle( void ){ ...@@ -464,47 +545,46 @@ int BrookRandomNumberGenerator::loadGVShuffle( void ){
const int np = sizeof(p) / sizeof(p[0]); const int np = sizeof(p) / sizeof(p[0]);
const int pmax = p[np-1]; const int pmax = p[np-1];
static float *buf = 0;
int iter, i, j;
static int *indices;
float tmp;
// static const char* methodName = "\nBrookRandomNumberGenerator::loadGVShuffle"; // static const char* methodName = "\nBrookRandomNumberGenerator::loadGVShuffle";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
if( buf == NULL ){ float* loadBuffer = _getLoadBuffer();
buf = (float*) malloc( sizeof(float) * sdp->gvSize ); int* indices = _getShuffleIndices( pmax );
indices = (int*) malloc( sizeof(int) * pmax );
}
int rvSize = getRandomNumberStreamSize(); int rvSize = getRandomNumberStreamSize();
for ( i = 0; i < rvSize; i++ ) { for ( int ii = 0; ii < rvSize; ii++ ) {
buf[i] = (float) i; loadBuffer[ii] = (float) ii;
} }
//How to come up with this number here? //How to come up with this number here?
unsigned long int seed = getRandomNumberSeed(); unsigned long int seed = getRandomNumberSeed();
for ( iter = 0; iter < 1000000; iter++ ) { for( int iter = 0; iter < 1000000; iter++ ) {
//for each p //for each p
for ( i = 0; i < np; i++ ){
for( int ii = 0; ii < np; ii++ ){
//Come up with p random indices //Come up with p random indices
//Not checking that they are distinct //Not checking that they are distinct
//because that should be fairly rare //because that should be fairly rare
for ( j = 0; j < p[i]; j++ ) {
for ( int jj = 0; jj < p[ii]; jj++ ) {
//indices[j] = (int) ( gmx_rando( &sdp->seed ) * sdp->gvSize ); //indices[j] = (int) ( gmx_rando( &sdp->seed ) * sdp->gvSize );
indices[j] = (int) ( generateGromacsRandomNumber( &seed )*rvSize ); indices[jj] = (int) ( _generateGromacsRandomNumber( &seed )*rvSize );
} }
//do a p-cycle
tmp = buf[ indices[0] ]; // do a p-cycle
for ( j = 0; j < p[i]-1; j++ ) {
buf[ indices[j] ] = buf[ indices[j+1] ]; float tmp = loadBuffer[ indices[0] ];
for ( int jj = 0; jj < p[ii]-1; jj++ ) {
loadBuffer[ indices[jj] ] = loadBuffer[ indices[jj+1] ];
} }
buf[ indices[ p[i]-1 ] ] = tmp; loadBuffer[ indices[ p[ii]-1 ] ] = tmp;
} }
} }
getShuffleStream()->loadFromArray( buf ); _getShuffleStream()->loadFromArray( loadBuffer );
return DefaultReturnValue; return DefaultReturnValue;
} }
...@@ -515,455 +595,197 @@ int BrookRandomNumberGenerator::loadGVShuffle( void ){ ...@@ -515,455 +595,197 @@ int BrookRandomNumberGenerator::loadGVShuffle( void ){
* @return DefaultReturnValue; * @return DefaultReturnValue;
*/ */
int BrookRandomNumberGenerator::shuffleGVStreams( void ){ int BrookRandomNumberGenerator::_shuffleGVStreams( void ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::shuffleGVStreams"; // static const char* methodName = "\nBrookRandomNumberGenerator::_shuffleGVStreams";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
int numberOfRvStreams = getNumberOfRandomNumberStreams(); int numberOfRvStreams = getNumberOfRandomNumberStreams();
for( int ii = 0; ii < numberOfRvStreams - 1; ii++ ){ for( int ii = 0; ii < numberOfRvStreams - 1; ii++ ){
kpermute_vectors( (float) , sdp->strShuffle, sdp->strVGauss[i+1], kpermute_vectors( (float) getRandomNumberStreamWidth(),
sdp->strVGauss[i] ); _getShuffleStream()->getBrookStream(),
getRandomNumberStream( ii + 1 )->getBrookStream(),
getRandomNumberStream( ii )->getBrookStream() );
} }
for( int jj = 0; jj < numberOfRvStreams - 1; jj++ ){ kpermute_vectors( (float) getRandomNumberStreamWidth(),
kpermute_vectors( (float) sdp->gvWidth, sdp->strShuffle, sdp->strVGauss[0], _getShuffleStream()->getBrookStream(),
sdp->strVGauss[ NGVSTREAMS-1 ] ); getRandomNumberStream( 0 )->getBrookStream(),
} getRandomNumberStream( numberOfRvStreams - 1 )->getBrookStream() );
return DefaultReturnValue; return DefaultReturnValue;
} }
static int /**
PrintGVStreams( gpuContext gpu ) * Advances the current position by 2*gpu->natoms
{ * If we run out of rand numbers, we shuffle and
int ii; * reuse a few times before reloading from the cpu
*
* @param numberOfRandomValuesConsumedPerIteration number of random values consumed/iteration
*
* @return DefaultReturnValue;
*/
gpuSDParams sdp = gpu->sdparams; int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsumedPerIteration ){
char testName[128];
char fileName[128];
/* dump inputs */ // ---------------------------------------------------------------------------------------
(void) sprintf( testName, "ShuffleGVStreams" ); // static const char* methodName = "\nBrookRandomNumberGenerator::advanceGVCursor";
(void) sprintf( fileName, "%s.32bits.in", testName );
Save32Bits( fileName, "iiiii", sdp->gvWidth, sdp->seed, sdp->gvOffset, sdp->gvCurStream, sdp->gvNumShuffles );
(void) sprintf( fileName, "%s.strShuffle.in", testName ); // ---------------------------------------------------------------------------------------
SaveStream( fileName, sdp->strShuffle );
for ( ii = 0; ii < NGVSTREAMS; ii++ ){ int rvStreamSize = getRandomNumberStreamSize();
(void) sprintf( testName, "ShuffleGVStreams%d", ii );
(void) sprintf( fileName, "%s.strVGaussIn.in", testName );
SaveStream( fileName, sdp->strVGauss[ii] );
} // use 2 random values per sd
return 1; _rvStreamOffset += numberOfRandomValuesConsumedPerIteration;
}
/* Advances the current position by 2*gpu->natoms //Check if we've used up this texture
* If we run out of rand numbers, we shuffle and
* reuse a few times before reloading from the cpu
* */
static int if ( _rvStreamOffset > rvStreamSize - numberOfRandomValuesConsumedPerIteration ){
AdvanceGVCursor( gpuSDParams sdp, gpuContext gpu, int step )
{
static clock_t rvTime = 0;
static clock_t shuffleTime = 0;
static const int trackTime = 0;
clock_t time;
if( trackTime ){ // next one if available
time = clock();
}
sdp->gvOffset += 2*gpu->natoms; _rvStreamOffset = 0;
//Check if we've used up this texture
if ( sdp->gvOffset > sdp->gvSize - 2*gpu->natoms ) { if ( _rvStreamIndex < _numberOfRandomNumberStreams - 1 ){
//Next one if available _rvStreamIndex++;
sdp->gvOffset = 0;
if ( sdp->gvCurStream < NGVSTREAMS - 1 ) {
sdp->gvCurStream++;
//printf( "Using stream %d\n", sdp->gvCurStream );
//fflush( stdout );
} else { } else {
//No more textures, need to shuffle //No more textures, need to shuffle
if ( sdp->gvNumShuffles < NGVSHUFFLE ) {
//printf( "Shuffling streams\n" );
//fflush( stdout );
clock_t localTime; if( _numberOfShuffles < _maxShuffles ){
if( trackTime ){
localTime = clock();
}
ShuffleGVStreams( sdp, gpu ); _shuffleGVStreams( );
_numberOfShuffles++;
if( trackTime ){
shuffleTime += clock() - localTime;
}
sdp->gvNumShuffles++;
} else { //Need to refresh random numbers from cpu } else { //Need to refresh random numbers from cpu
static int reducePrint = 0;
fflush( gpu->log );
if( UseOriginalRng ){ if( UseOriginalRng ){
LoadGVStreamsOriginal( sdp ); _loadGVStreamsOriginal( );
} else { } else {
LoadGVStreamsNvidia( sdp, gpu->log ); _loadRandomNumberStreamsKiss( );
} }
sdp->gvNumShuffles = 0; _numberOfShuffles = 0;
if( trackTime && !(reducePrint++ % 20) ){
rvTime += clock() - time;
(void) fprintf( gpu->log, "Reloading random number streams from cpu: seed=%d at step=%d", sdp->seed, step );
(void) fprintf( gpu->log, " cmltv times: RValue=%.5f Shuffle=%.5f\n", ((double) rvTime/(double)CLOCKS_PER_SEC), ((double) shuffleTime/(double)CLOCKS_PER_SEC) );
(void) fflush( gpu->log );
}
} }
sdp->gvCurStream = 0; _rvStreamIndex = 0;
} }
} }
/*
if( trackTime ){
rvTime += clock() - time;
if( (!sdp->gvNumShuffles && rvTime > 1) || !(sdp->gvNumShuffles %20) && (rvTime > 2*CLOCKS_PER_SEC) ){
(void) fprintf( gpu->log, "\n Total RValue=%.5f Shuffle=%.5f", ((double) rvTime/(double)CLOCKS_PER_SEC), ((double) shuffleTime/(double)CLOCKS_PER_SEC) );
}
}
*/
return 1;
}
extern "C" int
gpuResetSdSeed( gpuContext gpu, int step, int seed )
{
gpuSDParams sdp = gpu->sdparams;
sdp->seed = seed;
(void) fprintf( gpu->log, "gpuResetSdSeed: reset SD seed to %d at step=%d.\n", seed, step );
(void) fflush( gpu->log );
return 0;
}
/**
* Update parameters -- only way parameters can be set
*
* @param temperature temperature
* @param friction friction
* @param step size step size
*
* @return solute dielectric
*
*/
int BrookRandomNumberGenerator::updateParameters( double temperature, double friction, double stepSize ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::updateParameters";
// ---------------------------------------------------------------------------------------
_setStepSize( (BrookOpenMMFloat) stepSize );
_setFriction( (BrookOpenMMFloat) friction );
_setTemperature( (BrookOpenMMFloat) temperature );
_updateDerivedParameters( );
_updateSdStreams( );
return DefaultReturnValue; return DefaultReturnValue;
}
};
/** /**
* Update * Get random number stream size
*
* @param positions atom positions
* @param velocities atom velocities
* @param forces atom forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
* *
* @return DefaultReturnValue * @return random number stream size
* *
*/ */
int BrookRandomNumberGenerator::update( Stream& positions, Stream& velocities, int BrookRandomNumberGenerator::getRandomNumberStreamSize( void ) const {
const Stream& forces, return _randomNumberStreamSize;
BrookShakeAlgorithm& brookShakeAlgorithm ){ }
// ---------------------------------------------------------------------------------------
// unused Shake parameter
float omega = 1.0f;
// static const char* methodName = "\nBrookRandomNumberGenerator::update";
// ---------------------------------------------------------------------------------------
BrookOpenMMFloat* derivedParameters = getDerivedParameters();
BrookStreamImpl& positionStream = dynamic_cast<BrookStreamImpl&> (positions.getImpl());
BrookStreamImpl& velocityStream = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
BrookStreamImpl& forceStreamC = dynamic_cast<BrookStreamImpl&> (forces.getImpl());
const BrookStreamImpl& forceStream = dynamic_cast<const BrookStreamImpl&> (forceStream.getImpl());
// first integration step
kupdate_sd1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(),
(float) sdp->gvWidth,
(float) sdp->gvOffset,
derivedParameters[EM],
derivedParameters[Sd1pc1],
derivedParameters[Sd1pc2],
derivedParameters[Sd1pc3],
getSDPC1Stream()->getBrookStream(),
sdp->strVGauss[ sdp->gvCurStream ],
getSD2XStream()->getBrookStream(),
positionStream.getBrookStream(),
forceStream.getBrookStream(),
velocityStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
getSD1VStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getXPrimeStream()->getBrookStream()
);
AdvanceGVCursor( sdp, gpu, step );
// first Shake step
kshakeh_fix1(
10.0f,
(float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm->getInverseHydrogenMass(),
omega,
brookShakeAlgorithm->getShakeAtomIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeAtomParameterStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream() );
// first Shake gather
kshakeh_update1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(),
derivedParameters[Sd2pc1],
brookShakeAlgorithm->getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() );
// second integration step
kupdate_sd2_fix1(
(float) getStochasticDynamicsAtomStreamWidth(),
(float) sdp->gvWidth,
(float) sdp->gvOffset,
derivedParameters[Sd2pc1],
derivedParameters[Sd2pc2],
getSDPC2Stream()->getBrookStream(),
sdp->strVGauss[ sdp->gvCurStream ],
getSD1VStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getS2XStream()->getBrookStream(),
velocityStream.getBrookStream(),
getXPrimeStream()->getBrookStream()
);
AdvanceGVCursor( sdp, gpu, step );
// second Shake step
kshakeh_fix1(
10.0f,
(float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm->getInverseHydrogenMass(),
omega,
brookShakeAlgorithm->getShakeAtomIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeAtomParameterStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream() );
// second Shake gather
kshakeh_update1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(),
derivedParameters[Sd2pc1],
brookShakeAlgorithm->getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() );
kshakeh_update2_fix1(
(float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm->getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream(),
positionStream.getBrookStream() );
return DefaultReturnValue;
};
/** /**
* Get Shuffle stream
* *
* Get array of derived parameters indexed by 'DerivedParameters' enums * @return shuffle stream
*
* @return array
* *
*/ */
const BrookOpenMMFloat* BrookRandomNumberGenerator::getDerivedParameters( void ) const { BrookFloatStreamInternal* BrookRandomNumberGenerator::_getShuffleStream( void ) const {
return _auxiliaryStreams[ShuffleStream];
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::getDerivedParameters";
// ---------------------------------------------------------------------------------------
return _derivedParameters;
} }
/** /**
* Get Atom stream size * Get random number stream
* *
* @return Atom stream size * @param index random number stream index
* *
*/ * @return random number stream
int BrookRandomNumberGenerator::getStochasticDynamicsAtomStreamSize( void ) const {
return _sdAtomStreamSize;
}
/**
* Get atom stream width
* *
* @return atom stream width * @throw OpenMMException if index is invalid or _randomNumberGeneratorStreams not set
* *
*/ */
int BrookRandomNumberGenerator::getStochasticDynamicsAtomStreamWidth( void ) const { BrookFloatStreamInternal* BrookRandomNumberGenerator::getRandomNumberStream( int index ) const {
return _sdAtomStreamWidth;
}
/** // ---------------------------------------------------------------------------------------
* Get atom stream height
*
* @return atom stream height
*/
int BrookRandomNumberGenerator::getStochasticDynamicsAtomStreamHeight( void ) const {
return _sdAtomStreamHeight;
}
/** static const char* methodName = "\nBrookRandomNumberGenerator::getRandomNumberStream";
* Get SDPC1 stream
*
* @return SDPC1 stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSDPC1Stream( void ) const { // ---------------------------------------------------------------------------------------
return _sdStreams[SDPC1Stream];
}
/** if( index < 0 || index >= _numberOfRandomNumberStreams ){
* Get SDPC2 stream std::stringstream message;
* message << methodName << " index=" << index << " is out of range: [0," << _numberOfRandomNumberStreams;
* @return SDPC2 stream throw OpenMMException( message.str() );
* return NULL;
*/ }
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSDPC2Stream( void ) const { if( index < 0 || index >= _numberOfRandomNumberStreams ){
return _sdStreams[SDPC2Stream]; std::stringstream message;
} message << methodName << " index=" << index << " is out of range: [0," << _numberOfRandomNumberStreams;
throw OpenMMException( message.str() );
return NULL;
}
/** if( _randomNumberGeneratorStreams == NULL ){
* Get SD2X stream std::stringstream message;
* message << methodName << " rv streams not initialized; input index=" << index;
* @return SD2X stream throw OpenMMException( message.str() );
* return NULL;
*/ }
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSD2XStream( void ) const { return _randomNumberGeneratorStreams[index];
return _sdStreams[SD2XStream];
} }
/** /**
* Get SD1V stream * Get random value stream index
* *
* @return SD1V stream * @return random value stream index
* *
*/ */
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSD1VStream( void ) const { int BrookRandomNumberGenerator::getRvStreamIndex( void ) const {
return _sdStreams[SD1VStream]; return _rvStreamIndex;
} }
/** /**
* Get VPrime stream * Get random value stream offset
* *
* @return Vprime stream * @return random value stream offset
* *
*/ */
BrookFloatStreamInternal* BrookRandomNumberGenerator::getVPrimeStream( void ) const { int BrookRandomNumberGenerator::getRvStreamOffset( void ) const {
return _sdStreams[VPrimeStream]; return _rvStreamOffset;
} }
/** /**
* Get XPrime stream * Get max shuffles
* *
* @return Xprime stream * @return max shuffles
* *
*/ */
BrookFloatStreamInternal* BrookRandomNumberGenerator::getXPrimeStream( void ) const { int BrookRandomNumberGenerator::getMaxShuffles( void ) const {
return _sdStreams[XPrimeStream]; return _maxShuffles;
} }
/** /**
* Get InverseMass stream * Get number of shuffles
* *
* @return inverse mass stream * @return number of shuffles
* *
*/ */
BrookFloatStreamInternal* BrookRandomNumberGenerator::getInverseMassStream( void ) const { int BrookRandomNumberGenerator::_getNumberOfShuffles( void ) const {
return _sdStreams[InverseMassStream]; return _numberOfShuffles;
} }
/** /**
...@@ -984,110 +806,19 @@ int BrookRandomNumberGenerator::_initializeStreamSizes( int numberOfAtoms, const ...@@ -984,110 +806,19 @@ int BrookRandomNumberGenerator::_initializeStreamSizes( int numberOfAtoms, const
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_sdAtomStreamSize = getAtomStreamSize( platform ); setNumberOfAtoms( numberOfAtoms );
_sdAtomStreamWidth = getAtomStreamWidth( platform );
_sdAtomStreamHeight = getAtomStreamHeight( platform );
return DefaultReturnValue;
}
/**
* Initialize stream dimensions
*
* @param numberOfAtoms number of atoms
* @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
//std::string BrookRandomNumberGenerator::_getDerivedParametersString( BrookRandomNumberGenerator::DerivedParameters derivedParametersIndex ) const {
std::string BrookRandomNumberGenerator::_getDerivedParametersString( int derivedParametersIndex ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookRandomNumberGenerator::_getDerivedParametersString";
// ---------------------------------------------------------------------------------------
std::string returnString;
switch( derivedParametersIndex ){
case GDT:
returnString = "GDT";
break;
case EPH:
returnString = "EPH";
break;
case EMH:
returnString = "EMH";
break;
case EP:
returnString = "EP";
break;
case EM:
returnString = "EM";
break;
case B:
returnString = "B";
break;
case C:
returnString = "C";
break;
case D:
returnString = "D";
break;
case V:
returnString = "V";
break;
case X:
returnString = "X";
break;
case Yv:
returnString = "Yv";
break;
case Yx:
returnString = "Yx";
break;
case Sd1pc1:
returnString = "Sd1pc1";
break;
case Sd1pc2:
returnString = "Sd1pc2";
break;
case Sd1pc3:
returnString = "Sd1pc3";
break;
case Sd2pc1: // get randomNumber stream dimensions
returnString = "Sd2pc1";
break;
case Sd2pc2: const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
returnString = "Sd2pc2"; const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (platform.getDefaultStreamFactory() );
break; _randomNumberStreamWidth = brookStreamFactory.getDefaultRandomNumberStreamWidth();
_randomNumberStreamSize = brookStreamFactory.getDefaultRandomNumberStreamSize();
_randomNumberStreamHeight = (int) ( ((float) _randomNumberStreamSize)/( (float) _randomNumberStreamWidth) + 0.001);
default: _randomNumberStreamSize = _randomNumberStreamWidth*_randomNumberStreamHeight;
returnString = "Unknown";
break;
}
return returnString; return DefaultReturnValue;
} }
/** /**
...@@ -1105,138 +836,35 @@ int BrookRandomNumberGenerator::_initializeStreams( const Platform& platform ){ ...@@ -1105,138 +836,35 @@ int BrookRandomNumberGenerator::_initializeStreams( const Platform& platform ){
//static const std::string methodName = "BrookRandomNumberGenerator::_initializeStreams"; //static const std::string methodName = "BrookRandomNumberGenerator::_initializeStreams";
BrookOpenMMFloat dangleValue = (BrookOpenMMFloat) 0.0; float dangleValue = 0.0f;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
int sdAtomStreamSize = getStochasticDynamicsAtomStreamSize(); int randomNumberStreamSize = getRandomNumberStreamSize();
int sdAtomStreamWidth = getStochasticDynamicsAtomStreamWidth(); int randomNumberStreamWidth = getRandomNumberStreamWidth();
_sdStreams[SDPC1Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC1Stream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float2, dangleValue );
_sdStreams[SDPC2Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC2Stream, _auxiliaryStreams[ShuffleStream] = new BrookFloatStreamInternal( BrookCommon::ShuffleStream,
sdAtomStreamSize, sdAtomStreamWidth, randomNumberStreamSize, randomNumberStreamWidth,
BrookStreamInternal::Float2, dangleValue );
_sdStreams[SD2XStream] = new BrookFloatStreamInternal( BrookCommon::SD2XStream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[SD1VStream] = new BrookFloatStreamInternal( BrookCommon::SD1VStream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[VPrimeStream] = new BrookFloatStreamInternal( BrookCommon::VPrimeStream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[XPrimeStream] = new BrookFloatStreamInternal( BrookCommon::XPrimeStream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[InverseMassStream] = new BrookFloatStreamInternal( BrookCommon::InverseMassStream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float, dangleValue ); BrookStreamInternal::Float, dangleValue );
return DefaultReturnValue;
}
/**
* Update sd streams -- called after parameters change
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int BrookRandomNumberGenerator::_updateSdStreams( void ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookRandomNumberGenerator::_updateSdStreams";
// ---------------------------------------------------------------------------------------
int sdAtomStreamSize = getStochasticDynamicsAtomStreamSize();
BrookOpenMMFloat* sdpc[2];
for( int ii = 0; ii < 2; ii++ ){
sdpc[ii] = new BrookOpenMMFloat[2*sdAtomStreamSize];
memset( sdpc[ii], 0, 2*sdAtomStreamSize*sizeof( BrookOpenMMFloat ) );
}
const BrookOpenMMFloat* derivedParameters = getDerivedParameters( );
int numberOfAtoms = getNumberOfAtoms();
int index = 0;
for( int ii = 0; ii < numberOfAtoms; ii++, index += 2 ){
sdpc[0][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yv]) ); // insure number of random number streams is > 0
sdpc[0][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[V]) ); // delete if already allocated and then initialize
sdpc[1][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yx]) );
sdpc[1][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[X]) );
if( _numberOfRandomNumberStreams < 1 ){
_numberOfRandomNumberStreams = 1;
} }
_sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] ); if( _randomNumberGeneratorStreams ){
_sdStreams[SDPC2Stream]->loadFromArray( sdpc[1] ); delete[] _randomNumberGeneratorStreams;
for( int ii = 0; ii < 2; ii++ ){
delete[] sdpc[ii];
} }
// initialize SD2X _randomNumberGeneratorStreams = new BrookFloatStreamInternal*[_numberOfRandomNumberStreams];
BrookOpenMMFloat* sd2x = new BrookOpenMMFloat[3*sdAtomStreamSize];
//SimTKOpenMMUtilities::setRandomNumberSeed( (uint32_t) getRandomNumberSeed() );
memset( sd2x, 0, 3*sdAtomStreamSize*sizeof( BrookOpenMMFloat ) );
index = 0;
for( int ii = 0; ii < numberOfAtoms; ii++, index += 3 ){
sd2x[index] = _inverseSqrtMasses[ii]*derivedParameters[X]*( static_cast<BrookOpenMMFloat> (SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()) );
sd2x[index+1] = _inverseSqrtMasses[ii]*derivedParameters[X]*( static_cast<BrookOpenMMFloat> (SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()) );
sd2x[index+2] = _inverseSqrtMasses[ii]*derivedParameters[X]*( static_cast<BrookOpenMMFloat> (SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()) );
}
_sdStreams[SD2XStream]->loadFromArray( sd2x );
delete[] sd2x;
return DefaultReturnValue;
}
/**
* Set masses
*
* @param masses atomic masses
*
*/
int BrookRandomNumberGenerator::_setInverseSqrtMasses( const std::vector<double>& masses ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookRandomNumberGenerator::_setInverseSqrtMasses";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0;
// --------------------------------------------------------------------------------------- for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
_randomNumberGeneratorStreams[ii] = new BrookFloatStreamInternal( BrookCommon::ShuffleStream,
// setup inverse sqrt masses randomNumberStreamSize, randomNumberStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_inverseSqrtMasses = new BrookOpenMMFloat[masses.size()];
int index = 0;
for( std::vector<double>::const_iterator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
if( *ii != 0.0 ){
BrookOpenMMFloat value = static_cast<BrookOpenMMFloat>(*ii);
_inverseSqrtMasses[index] = ( SQRT( one/value ) );
} else {
_inverseSqrtMasses[index] = zero;
}
} }
return DefaultReturnValue; return DefaultReturnValue;
...@@ -1245,31 +873,26 @@ int BrookRandomNumberGenerator::_setInverseSqrtMasses( const std::vector<double> ...@@ -1245,31 +873,26 @@ int BrookRandomNumberGenerator::_setInverseSqrtMasses( const std::vector<double>
/* /*
* Setup of StochasticDynamics parameters * Setup of StochasticDynamics parameters
* *
* @param masses masses * @param numberOfAtoms number of atoms
* @param platform Brook platform * @param platform Brook platform
* *
* @return nonzero value if error * @return nonzero value if error
* *
* */ * */
int BrookRandomNumberGenerator::setup( const std::vector<double>& masses, const Platform& platform ){ int BrookRandomNumberGenerator::setup( int numberOfAtoms, const Platform& platform ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookRandomNumberGenerator::setup"; //static const std::string methodName = "BrookRandomNumberGenerator::setup";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
int numberOfAtoms = (int) masses.size();
setNumberOfAtoms( numberOfAtoms );
// set stream sizes and then create streams // set stream sizes and then create streams
_initializeStreamSizes( numberOfAtoms, platform ); _initializeStreamSizes( numberOfAtoms, platform );
_initializeStreams( platform ); _initializeStreams( platform );
_setInverseSqrtMasses( masses );
return DefaultReturnValue; return DefaultReturnValue;
} }
...@@ -1304,47 +927,44 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const { ...@@ -1304,47 +927,44 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) ); #define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );
#endif #endif
(void) LOCAL_SPRINTF( value, "%d", getNumberOfAtoms() ); (void) LOCAL_SPRINTF( value, "%s", (UseOriginalRng ? "Original Rng" : "Kiss Rng") );
message << _getLine( tab, "Number of atoms:", value ); message << _getLine( tab, "Random number generator:", value );
(void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() ); (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamWidth() );
message << _getLine( tab, "Atom stream width:", value ); message << _getLine( tab, "RandomNumber stream width:", value );
(void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() ); (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamWidth() );
message << _getLine( tab, "Atom stream height:", value ); message << _getLine( tab, "RandomNumber stream width:", value );
(void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() ); (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamHeight() );
message << _getLine( tab, "Atom stream size:", value ); message << _getLine( tab, "RandomNumber stream height:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getTau() ); (void) LOCAL_SPRINTF( value, "%d", getNumberOfRandomNumberStreams() );
message << _getLine( tab, "Tau:", value ); message << _getLine( tab, "Number of RandomNumber streams:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getTemperature() ); (void) LOCAL_SPRINTF( value, "%ld", getRandomNumberSeed() );
message << _getLine( tab, "Temperature:", value ); message << _getLine( tab, "Seed:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getStepSize() ); (void) LOCAL_SPRINTF( value, "%d", getRvStreamIndex() );
message << _getLine( tab, "Step size:", value ); message << _getLine( tab, "Rv stream index:", value );
const BrookOpenMMFloat* derivedParameters = getDerivedParameters(); (void) LOCAL_SPRINTF( value, "%d", getRvStreamOffset() );
for( int ii = 0; ii < MaxDerivedParameters; ii++ ){ message << _getLine( tab, "Rv stream offset:", value );
(void) LOCAL_SPRINTF( value, "%.5e", derivedParameters[ii] );
message << _getLine( tab, _getDerivedParametersString( ii ), value ); (void) LOCAL_SPRINTF( value, "%d", _getNumberOfShuffles() );
} message << _getLine( tab, "Number of rv stream shuffles:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getTemperature() ); (void) LOCAL_SPRINTF( value, "%d", getMaxShuffles() );
message << _getLine( tab, "Temperature:", value ); message << _getLine( tab, "Max number of rv stream shuffles:", value );
message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) ); message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) );
message << _getLine( tab, "SDPC1:", (getSDPC1Stream() ? Set : NotSet) ); message << _getLine( tab, "Shuffle:", (_getShuffleStream() ? Set : NotSet) );
message << _getLine( tab, "SDPC2:", (getSDPC2Stream() ? Set : NotSet) );
message << _getLine( tab, "SD2X:", (getSD2XStream() ? Set : NotSet) );
message << _getLine( tab, "SD1V:", (getSD1VStream() ? Set : NotSet) );
for( int ii = 0; ii < LastStreamIndex; ii++ ){ for( int ii = 0; ii < LastStreamIndex; ii++ ){
message << std::endl; message << std::endl;
if( _sdStreams[ii] ){ if( _auxiliaryStreams[ii] ){
message << _sdStreams[ii]->getContentsString( ); message << _auxiliaryStreams[ii]->getContentsString( );
} }
} }
......
...@@ -32,12 +32,6 @@ ...@@ -32,12 +32,6 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
#include <vector>
#include <set>
#include "BrookFloatStreamInternal.h"
#include "BrookShakeAlgorithm.h"
#include "BrookPlatform.h"
#include "BrookCommon.h" #include "BrookCommon.h"
namespace OpenMM { namespace OpenMM {
...@@ -52,6 +46,10 @@ class BrookRandomNumberGenerator : public BrookCommon { ...@@ -52,6 +46,10 @@ class BrookRandomNumberGenerator : public BrookCommon {
public: public:
// toggle between original rng & Kiss (Nvidia) code
static const int UseOriginalRng = 1;
/** /**
* Constructor * Constructor
* *
...@@ -67,98 +65,38 @@ class BrookRandomNumberGenerator : public BrookCommon { ...@@ -67,98 +65,38 @@ class BrookRandomNumberGenerator : public BrookCommon {
~BrookRandomNumberGenerator(); ~BrookRandomNumberGenerator();
/** /**
* Get tau * Get number of random number streams
*
* @return tau
*/
BrookOpenMMFloat getTau( void ) const;
/**
* Get friction
*
* @return friction
*/
BrookOpenMMFloat getFriction( void ) const;
/**
* Get temperature
*
* @return temperature
*/
BrookOpenMMFloat getTemperature( void ) const;
/**
* Get step size
*
* @return step size
*/
BrookOpenMMFloat getStepSize( void ) const;
/**
*
* Get array of derived parameters indexed by 'DerivedParameters' enums
*
* @return array
* *
*/ * @return number of random number streams
const BrookOpenMMFloat* getDerivedParameters( void ) const;
/**
* Get StochasticDynamics atom stream width
* *
* @return atom stream width
*/ */
int getStochasticDynamicsAtomStreamWidth( void ) const; int getNumberOfRandomNumberStreams( void ) const;
/** /**
* Get StochasticDynamics atom stream height * Get stream width
* *
* @return atom stream height * @return stream width
*/ */
int getStochasticDynamicsAtomStreamHeight( void ) const; int getRandomNumberStreamWidth( void ) const;
/** /**
* Get StochasticDynamics atom stream size * Get stream height
* *
* @return atom stream size * @return stream height
*/ */
int getStochasticDynamicsAtomStreamSize( void ) const; int getRandomNumberStreamHeight( void ) const;
/** /**
* Update parameters * Get stream size
*
* @param temperature temperature
* @param friction friction
* @param step size step size
*
* @return DefaultReturnValue
* *
* @return stream size
*/ */
int updateParameters( double temperature, double friction, double stepSize ); int getRandomNumberStreamSize( void ) const;
/**
* Update
*
* @param positions atom positions
* @param velocities atom velocities
* @param forces atom forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
*
* @return DefaultReturnValue
*
*/
int update( Stream& positions, Stream& velocities,
const Stream& forces, BrookShakeAlgorithm& brookShakeAlgorithm );
/** /**
* Get array of StochasticDynamics streams * Get array of StochasticDynamics streams
* *
...@@ -169,16 +107,16 @@ class BrookRandomNumberGenerator : public BrookCommon { ...@@ -169,16 +107,16 @@ class BrookRandomNumberGenerator : public BrookCommon {
BrookFloatStreamInternal** getStreams( void ); BrookFloatStreamInternal** getStreams( void );
/* /*
* Setup of StochasticDynamics parameters * Setup of RNG parameters
* *
* @param masses atom masses * @param numberOfAtoms number of atoms
* @param platform Brook platform * @param platform Brook platform
* *
* @return ErrorReturnValue value if error, else DefaultReturnValue * @return ErrorReturnValue value if error, else DefaultReturnValue
* *
* */ * */
int setup( const std::vector<double>& masses, const Platform& platform ); int setup( int numberOfAtoms, const Platform& platform );
/* /*
* Get contents of object * Get contents of object
...@@ -192,104 +130,93 @@ class BrookRandomNumberGenerator : public BrookCommon { ...@@ -192,104 +130,93 @@ class BrookRandomNumberGenerator : public BrookCommon {
std::string getContentsString( int level = 0 ) const; std::string getContentsString( int level = 0 ) const;
/** /**
* Get SDPC1 stream * Get random number stream
* *
* @return SDPC1 stream * @param index random number stream index
*
* @return random number stream
* *
*/ */
BrookFloatStreamInternal* getSDPC1Stream( void ) const; BrookFloatStreamInternal* getRandomNumberStream( int index ) const;
/** /**
* Get SDPC2 stream * Get random number seed
*
* @return SDPC2 stream
* *
* @return random number seed
*/ */
BrookFloatStreamInternal* getSDPC2Stream( void ) const; unsigned long int getRandomNumberSeed( void ) const;
/** /**
* Get shuffle stream * Increment random number seed
* *
* @return Shuffle stream * @param increment amount to increment random number seed; default = 1
* *
* @return updated random number seed
*/ */
BrookFloatStreamInternal* getShuffleStream( void ) const; unsigned long int incrementRandomNumberSeed( unsigned long int increment = 1 );
/** /**
* Generate a random number using algorithm in Gromacs * Set random number seed
*
* @param ig seed
* *
* @return random number * @param new random number seed; default = 1
* *
* @return random number seed
*/ */
BrookOpenMMFloat generateGromacsRandomNumber( int* ig ); unsigned long int setRandomNumberSeed( unsigned long int seed = 1 );
/** /**
* Generate a random number using algorithm in Nvidia code * Get index of rv texture
* http://www.helsbreth.org/random/rng_kiss.html
*
* @param randomV1 output random value
* @param randomV2 output random value
* @param randomV3 output random value
* @param state state
* *
* @return index of rv texture
*/ */
void generateRandomsAlaNvidia( float* randomV1, float* randomV2, float* randomV3, int getRvStreamIndex( void ) const;
unsigned int state[4] );
/** /**
* Load random number streams using Nvidia algorithm * Get max shuffles
* *
* @return max shuffles
* *
* @return DefaultReturnValue;
*/ */
int loadRandomNumberStreamsNvidia( void ); int getMaxShuffles( void ) const;
/** /**
* Get random number seed * Advance random values stream index
* *
* @return random number seed * @param numberOfEntriesToAdvance number of entries consumed in previous iteration
*/
unsigned long int getRandomNumberSeed( void ) const;
/**
* Increment random number seed
* *
* @param increment amount to increment random number seed; default = 1 * @return DefaultReturnValue
* *
* @return updated random number seed
*/ */
unsigned long int incrementRandomNumberSeed( unsigned long int increment = 1 ); int advanceGVCursor( int numberOfEntriesToAdvance );
/** /**
* Set random number seed * Get random value stream offset
* *
* @param new random number seed; default = 1 * @return random value stream offset
* *
* @return random number seed
*/ */
unsigned long int setRandomNumberSeed( unsigned long int seed = 1 ); int getRvStreamOffset( void ) const;
private: private:
// streams indices // streams indices
enum BrookRandomNumberGeneratorStreams { enum BrookRandomNumberGeneratorStreams {
RandomNumberStream,
ShuffleStream, ShuffleStream,
LastStreamIndex LastStreamIndex
}; };
BrookFloatStreamInternal* _auxiliaryStreams[LastStreamIndex];
BrookFloatStreamInternal** _randomNumberGeneratorStreams;
// randomNumberSeed // randomNumberSeed
unsigned long int _randomNumberSeed; unsigned long int _randomNumberSeed;
...@@ -304,6 +231,16 @@ class BrookRandomNumberGenerator : public BrookCommon { ...@@ -304,6 +231,16 @@ class BrookRandomNumberGenerator : public BrookCommon {
int _randomNumberStreamHeight; int _randomNumberStreamHeight;
int _randomNumberStreamSize; int _randomNumberStreamSize;
// control variables
int _rvStreamIndex;
int _rvStreamOffset;
int _numberOfShuffles;
int _maxShuffles;
float* _loadBuffer;
int* _shuffleIndices;
/* /*
* Setup of stream dimensions * Setup of stream dimensions
* *
...@@ -340,15 +277,125 @@ class BrookRandomNumberGenerator : public BrookCommon { ...@@ -340,15 +277,125 @@ class BrookRandomNumberGenerator : public BrookCommon {
int _initializeStreams( const Platform& platform ); int _initializeStreams( const Platform& platform );
/** /**
* Set masses * Increment random number offset
*
* @param increment increment for offset
*
* @return random number offset
*/
int _incrementRvOffset( int increment );
/**
* Get shuffle stream
*
* @return Shuffle stream
*
*/
BrookFloatStreamInternal* _getShuffleStream( void ) const;
/**
* Generate a random number using algorithm in Gromacs
*
* @param ig seed
*
* @return random number
*
*/
BrookOpenMMFloat _generateGromacsRandomNumber( unsigned long int* ig );
/**
* Generate a random number using Kiss (algorithm in Kiss code)
* http://www.helsbreth.org/random/rng_kiss.html
*
* @param randomV1 output random value
* @param randomV2 output random value
* @param randomV3 output random value
* @param state state
*
*/
void _generateRandomsKiss( float* randomV1, float* randomV2, float* randomV3,
unsigned int state[4] );
/**
* Load random number streams using Kiss algorithm
*
*
* @return DefaultReturnValue;
*/
int _loadRandomNumberStreamsKiss( void );
/**
* Load random number streams using original gpu algorithm
*
*
* @return DefaultReturnValue;
*/
int _loadGVStreamsOriginal( void );
/**
* Loads a permutation of indices from 0 to gvSize-1 in
* sdp->strShuffle. To make sure that the order of the
* permutation is atleast NGVSHUFFLE, we create the
* permutation by introducing a random number of p-cycles
* where p is randomly determined from 2,3,5,7 and 11.
* The LCM of these numbers is 2310.
* Ofcourse the p-cycles are not necessarily disjoint
* the way it's done here, but there's a good chance
* there will enough disjoint cycles to make the
* order of the permutation larger than NGVSHUFFLE
*
*
* This function is only called once at startup
*
* @return DefaultReturnValue;
**/
int _loadGVShuffle( void );
/**
* Get number of shuffles
* *
* @param masses atomic masses * @return number of shuffles
* *
*/ */
int _setInverseSqrtMasses( const std::vector<double>& masses ); int _getNumberOfShuffles( void ) const;
/**
* Load buffer
*
* @return ptr to load buffer
*
* @throw OpenMMException if rv stream size is < 1
*
**/
float* _getLoadBuffer( void );
/**
* Get ptr to shuffle indices
*
* @return ptr to shuffle indices
*
* @throw OpenMMException if size is < 1
*
**/
int* _getShuffleIndices( int size );
/**
* Shuffle streams
*
* @return DefaultReturnValue;
*/
int _shuffleGVStreams( void );
}; };
......
...@@ -34,6 +34,8 @@ ...@@ -34,6 +34,8 @@
#include "BrookPlatform.h" #include "BrookPlatform.h"
#include "OpenMMException.h" #include "OpenMMException.h"
#include "BrookStreamImpl.h" #include "BrookStreamImpl.h"
#include "kshakeh.h"
#include "kupdatesd.h"
// use random number generator // use random number generator
...@@ -330,6 +332,7 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti ...@@ -330,6 +332,7 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti
* @param velocities atom velocities * @param velocities atom velocities
* @param forces atom forces * @param forces atom forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference * @param brookShakeAlgorithm BrookShakeAlgorithm reference
* @param brookRandomNumberGenerator BrookRandomNumberGenerator reference
* *
* @return DefaultReturnValue * @return DefaultReturnValue
* *
...@@ -337,7 +340,8 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti ...@@ -337,7 +340,8 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti
int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
const Stream& forces, const Stream& forces,
BrookShakeAlgorithm& brookShakeAlgorithm ){ BrookShakeAlgorithm& brookShakeAlgorithm,
BrookRandomNumberGenerator& brookRandomNumberGenerator ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -349,26 +353,25 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -349,26 +353,25 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
BrookOpenMMFloat* derivedParameters = getDerivedParameters(); const BrookOpenMMFloat* derivedParameters = getDerivedParameters();
BrookStreamImpl& positionStream = dynamic_cast<BrookStreamImpl&> (positions.getImpl()); BrookStreamImpl& positionStream = dynamic_cast<BrookStreamImpl&> (positions.getImpl());
BrookStreamImpl& velocityStream = dynamic_cast<BrookStreamImpl&> (velocities.getImpl()); BrookStreamImpl& velocityStream = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
BrookStreamImpl& forceStreamC = dynamic_cast<BrookStreamImpl&> (forces.getImpl()); const BrookStreamImpl& forceStreamC = dynamic_cast<const BrookStreamImpl&> (forces.getImpl());
const BrookStreamImpl& forceStream = dynamic_cast<const BrookStreamImpl&> (forceStream.getImpl()); BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
// first integration step // first integration step
kupdate_sd1_fix1( kupdate_sd1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
(float) sdp->gvWidth, (float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
(float) sdp->gvOffset, (float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM], derivedParameters[EM],
derivedParameters[Sd1pc1], derivedParameters[Sd1pc1],
derivedParameters[Sd1pc2], derivedParameters[Sd1pc2],
derivedParameters[Sd1pc3], derivedParameters[Sd1pc3],
getSDPC1Stream()->getBrookStream(), getSDPC1Stream()->getBrookStream(),
sdp->strVGauss[ sdp->gvCurStream ], brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD2XStream()->getBrookStream(), getSD2XStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream(),
...@@ -379,100 +382,104 @@ sdp->strVGauss[ sdp->gvCurStream ], ...@@ -379,100 +382,104 @@ sdp->strVGauss[ sdp->gvCurStream ],
getXPrimeStream()->getBrookStream() getXPrimeStream()->getBrookStream()
); );
AdvanceGVCursor( sdp, gpu, step ); // advance random number cursor
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() );
// first Shake step // first Shake step
kshakeh_fix1( kshakeh_fix1(
10.0f, 10.0f,
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm->getInverseHydrogenMass(), brookShakeAlgorithm.getInverseHydrogenMass(),
omega, omega,
brookShakeAlgorithm->getShakeAtomIndicesStream()->getBrookStream(), brookShakeAlgorithm.getShakeAtomIndicesStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeAtomParameterStream()->getBrookStream(), brookShakeAlgorithm.getShakeAtomParameterStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream() ); brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// first Shake gather // first Shake gather
kshakeh_update1_fix1( kshakeh_update1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc1],
brookShakeAlgorithm->getShakeInverseMapStream()->getBrookStream(), brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(), getVPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() ); getXPrimeStream()->getBrookStream() );
// second integration step // second integration step
kupdate_sd2_fix1( kupdate_sd2_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
(float) sdp->gvWidth, (float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
(float) sdp->gvOffset, (float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc1],
derivedParameters[Sd2pc2], derivedParameters[Sd2pc2],
getSDPC2Stream()->getBrookStream(), getSDPC2Stream()->getBrookStream(),
sdp->strVGauss[ sdp->gvCurStream ], brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD1VStream()->getBrookStream(), getSD1VStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(), getVPrimeStream()->getBrookStream(),
getS2XStream()->getBrookStream(), getSD2XStream()->getBrookStream(),
velocityStream.getBrookStream(), velocityStream.getBrookStream(),
getXPrimeStream()->getBrookStream() getXPrimeStream()->getBrookStream()
); );
AdvanceGVCursor( sdp, gpu, step ); // advance random number cursor
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() );
// second Shake step // second Shake step
kshakeh_fix1( kshakeh_fix1(
10.0f, 10.0f,
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm->getInverseHydrogenMass(), brookShakeAlgorithm.getInverseHydrogenMass(),
omega, omega,
brookShakeAlgorithm->getShakeAtomIndicesStream()->getBrookStream(), brookShakeAlgorithm.getShakeAtomIndicesStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeAtomParameterStream()->getBrookStream(), brookShakeAlgorithm.getShakeAtomParameterStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream() ); brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// second Shake gather // second Shake gather
kshakeh_update1_fix1( kshakeh_update1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc1],
brookShakeAlgorithm->getShakeInverseMapStream()->getBrookStream(), brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(), getVPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() ); getXPrimeStream()->getBrookStream() );
kshakeh_update2_fix1( kshakeh_update2_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm->getShakeInverseMapStream()->getBrookStream(), brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons0Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm->getShakeXCons3Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
positionStream.getBrookStream() ); positionStream.getBrookStream() );
return DefaultReturnValue; return DefaultReturnValue;
......
...@@ -37,6 +37,7 @@ ...@@ -37,6 +37,7 @@
#include "BrookFloatStreamInternal.h" #include "BrookFloatStreamInternal.h"
#include "BrookShakeAlgorithm.h" #include "BrookShakeAlgorithm.h"
#include "BrookRandomNumberGenerator.h"
#include "BrookPlatform.h" #include "BrookPlatform.h"
#include "BrookCommon.h" #include "BrookCommon.h"
...@@ -152,13 +153,15 @@ class BrookStochasticDynamics : public BrookCommon { ...@@ -152,13 +153,15 @@ class BrookStochasticDynamics : public BrookCommon {
* @param velocities atom velocities * @param velocities atom velocities
* @param forces atom forces * @param forces atom forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference * @param brookShakeAlgorithm BrookShakeAlgorithm reference
* @param brookRandomNumberGenerator BrookRandomNumberGenerator reference
* *
* @return DefaultReturnValue * @return DefaultReturnValue
* *
*/ */
int update( Stream& positions, Stream& velocities, int update( Stream& positions, Stream& velocities,
const Stream& forces, BrookShakeAlgorithm& brookShakeAlgorithm ); const Stream& forces, BrookShakeAlgorithm& brookShakeAlgorithm,
BrookRandomNumberGenerator& brookRandomNumberGenerator );
/** /**
* Get array of StochasticDynamics streams * Get array of StochasticDynamics streams
* *
...@@ -227,6 +230,33 @@ class BrookStochasticDynamics : public BrookCommon { ...@@ -227,6 +230,33 @@ class BrookStochasticDynamics : public BrookCommon {
BrookFloatStreamInternal* getSD1VStream( void ) const; BrookFloatStreamInternal* getSD1VStream( void ) const;
/**
* Get V-prime stream
*
* @return V-prime stream
*
*/
BrookFloatStreamInternal* getVPrimeStream( void ) const;
/**
* Get X-prime stream
*
* @return X-prime stream
*
*/
BrookFloatStreamInternal* getXPrimeStream( void ) const;
/**
* Get inverse sqrt masses
*
* @return inverse sqrt masses stream
*
*/
BrookFloatStreamInternal* getInverseMassStream( void ) 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,
...@@ -385,8 +415,6 @@ class BrookStochasticDynamics : public BrookCommon { ...@@ -385,8 +415,6 @@ class BrookStochasticDynamics : public BrookCommon {
int _setInverseSqrtMasses( const std::vector<double>& masses ); int _setInverseSqrtMasses( const std::vector<double>& masses );
}; };
} // namespace OpenMM } // namespace OpenMM
......
...@@ -56,7 +56,9 @@ BrookStreamFactory::BrookStreamFactory( void ){ ...@@ -56,7 +56,9 @@ BrookStreamFactory::BrookStreamFactory( void ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_defaultDangleValue = 1.0e+38; _defaultDangleValue = 1.0e+38;
_defaultAtomStreamWidth = 32; _defaultAtomStreamWidth = DefaultStreamAtomWidth;
_defaultStreamRandomNumberWidth = DefaultStreamRandomNumberWidth;
_defaultStreamRandomNumberSize = DefaultStreamRandomNumberSize;
} }
...@@ -121,7 +123,111 @@ int BrookStreamFactory::setDefaultAtomStreamWidth( int atomStreamWidth ){ ...@@ -121,7 +123,111 @@ int BrookStreamFactory::setDefaultAtomStreamWidth( int atomStreamWidth ){
} }
/** /**
* get default dangle value * Get randomNumber stream width
*
* @return randomNumberStreamWidth
*
*/
int BrookStreamFactory::getDefaultRandomNumberStreamWidth( void ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookStreamFactory::getDefaultRandomNumberStreamWidth";
// ---------------------------------------------------------------------------------------
return _defaultStreamRandomNumberWidth;
}
/**
* Set randomNumber stream width
*
* @param randomNumberStreamWidth randomNumber stream width
*
* @return DefaultReturnValue
*
* @throw OpenMMException if randomNumberStreamWidth < 1
*
*/
int BrookStreamFactory::setDefaultRandomNumberStreamWidth( int randomNumberStreamWidth ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookStreamFactory::setDefaultRandomNumberStreamWidth";
// ---------------------------------------------------------------------------------------
// validate randomNumber stream width
if( randomNumberStreamWidth < 1 ){
std::stringstream message;
message << methodName << " randomNumberStreamWidth=" << randomNumberStreamWidth << " is less than 1.";
throw OpenMMException( message.str() );
return ErrorReturnValue;
}
_defaultStreamRandomNumberWidth = randomNumberStreamWidth;
return DefaultReturnValue;
}
/*
* Get randomNumber stream size
*
* @return randomNumberStreamSize
*
*/
int BrookStreamFactory::getDefaultRandomNumberStreamSize( void ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookStreamFactory::getDefaultRandomNumberStreamSize";
// ---------------------------------------------------------------------------------------
return _defaultStreamRandomNumberSize;
}
/**
* Set randomNumber stream size
*
* @param randomNumberStreamSize randomNumber stream size
*
* @return DefaultReturnValue
*
* @throw OpenMMException if randomNumberStreamSize < 1
*
*/
int BrookStreamFactory::setDefaultRandomNumberStreamSize( int randomNumberStreamSize ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookStreamFactory::setDefaultRandomNumberStreamSize";
// ---------------------------------------------------------------------------------------
// validate randomNumber stream size
if( randomNumberStreamSize < 1 ){
std::stringstream message;
message << methodName << " randomNumberStreamSize=" << randomNumberStreamSize << " is less than 1.";
throw OpenMMException( message.str() );
return ErrorReturnValue;
}
_defaultStreamRandomNumberSize = randomNumberStreamSize;
return DefaultReturnValue;
}
/**
* Get default dangle value
* *
* @return default dangle value * @return default dangle value
* *
......
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