Skip to content

GitLab

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

Daily commit

parent 4ed00931
Hide whitespace changes
Inline Side-by-side
Showing with 1942 additions and 0 deletions
platforms/brook/src/BrookCommon.cpp
View file @ a475f0a6
...@@ -75,6 +75,9 @@ const std::string BrookCommon::SDPC1Stream ...@@ -75,6 +75,9 @@ const std::string BrookCommon::SDPC1Stream
const std::string BrookCommon::SDPC2Stream = "SDPC2Stream"; const std::string BrookCommon::SDPC2Stream = "SDPC2Stream";
const std::string BrookCommon::SD2XStream = "SD2XStream"; const std::string BrookCommon::SD2XStream = "SD2XStream";
const std::string BrookCommon::SD1VStream = "SD1VStream"; const std::string BrookCommon::SD1VStream = "SD1VStream";
const std::string BrookCommon::VPrimeStream = "VPrimeStream";
const std::string BrookCommon::XPrimeStream = "XPrimeStream";
const std::string BrookCommon::InverseMassStream = "InverseMassStream";
// Shake streams // Shake streams
...@@ -108,6 +111,7 @@ BrookCommon::BrookCommon( ){ ...@@ -108,6 +111,7 @@ BrookCommon::BrookCommon( ){
_atomStreamSize = -1; _atomStreamSize = -1;
_log = NULL; _log = NULL;
_verbosity = 0;
} }
...@@ -333,6 +337,31 @@ int BrookCommon::setLog( FILE* log ){ ...@@ -333,6 +337,31 @@ int BrookCommon::setLog( FILE* log ){
return BrookCommon::DefaultReturnValue; return BrookCommon::DefaultReturnValue;
} }
/**
* Get verbosity
*
* @return verbosity
*
*/
int BrookCommon::getVerbosity( void ) const {
return _verbosity;
}
/**
* Set verbosity
*
* @param verbosity
*
* @return DefaultReturnValue
*
*/
int BrookCommon::setVerbosity( int verbosity ){
_verbosity = verbosity;
return BrookCommon::DefaultReturnValue;
}
/* /*
* Get contents of object * Get contents of object
* *
......
platforms/brook/src/BrookCommon.h
View file @ a475f0a6
...@@ -90,6 +90,9 @@ class BrookCommon { ...@@ -90,6 +90,9 @@ class BrookCommon {
static const std::string SDPC2Stream; static const std::string SDPC2Stream;
static const std::string SD2XStream; static const std::string SD2XStream;
static const std::string SD1VStream; static const std::string SD1VStream;
static const std::string VPrimeStream;
static const std::string XPrimeStream;
static const std::string InverseMassStream;
// Shake streams // Shake streams
...@@ -220,6 +223,24 @@ class BrookCommon { ...@@ -220,6 +223,24 @@ class BrookCommon {
FILE* getLog( void ) const; FILE* getLog( void ) const;
/**
* Be verbose flag
*
* @return verbosity flag
*
*/
int getVerbosity( void ) const;
/**
* Set verbosity flag
*
* @param verbosity flag
*
*/
int setVerbosity( int verbositym);
/* /*
* Given number of stream elements and width, returns the appropriate * Given number of stream elements and width, returns the appropriate
* height of the stream * height of the stream
...@@ -266,6 +287,10 @@ class BrookCommon { ...@@ -266,6 +287,10 @@ class BrookCommon {
FILE* _log; FILE* _log;
// verbosity
int _verbosity;
/** /**
* Set number of atoms * Set number of atoms
* *
......
platforms/brook/src/BrookRandomNumberGenerator.cpp 0 → 100644
View file @ a475f0a6
working on shuffleGVStreams
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include <sstream>
#include "BrookRandomNumberGenerator.h"
#include "BrookPlatform.h"
#include "OpenMMException.h"
#include "BrookStreamImpl.h"
// use random number generator
#include "SimTKOpenMMUtilities.h"
using namespace OpenMM;
using namespace std;
/**
*
* Constructor
*
*/
BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookRandomNumberGenerator::BrookRandomNumberGenerator";
// ---------------------------------------------------------------------------------------
// fixed for now
_numberOfRandomNumberStreams = 2;
// mark stream dimension variables as unset
_randomNumberStreamWidth = -1;
_randomNumberStreamHeight = -1;
_randomNumberStreamSize = -1;
for( int ii = 0; ii < LastStreamIndex; ii++ ){
_randomNumberStreams[ii] = NULL;
}
// set randomNumber seed
_randomNumberSeed = 1393;
//_randomNumberSeed = randomNumberSeed ? randomNumberSeed : 1393;
//SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
}
/**
* Destructor
*
*/
BrookRandomNumberGenerator::~BrookRandomNumberGenerator( ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookRandomNumberGenerator::~BrookRandomNumberGenerator";
// ---------------------------------------------------------------------------------------
for( int ii = 0; ii < LastStreamIndex; ii++ ){
delete _sdStreams[ii];
}
}
/**
* Get number of random number streams
*
* @return number of random number streams
*
*/
int BrookRandomNumberGenerator::getNumberOfRandomNumberStreams( void ) const {
return _numberOfRandomNumberStreams;
}
/**
* Get number of random number streams
*
* @return number of random number streams
*
*/
int BrookRandomNumberGenerator::getNumberOfRandomNumberStreams( void ) const {
return _numberOfRandomNumberStreams;
}
/**
* Get random number seed
*
* @return random number seed
*/
unsigned long int BrookRandomNumberGenerator::getRandomNumberSeed( void ) const {
return _randomNumberSeed;
}
/**
* Increment random number seed
*
* @param increment amount to increment random number seed; default = 1
*
* @return updated random number seed
*/
unsigned long int BrookRandomNumberGenerator::incrementRandomNumberSeed( unsigned long int increment ){
_randomNumberSeed += increment;
return _randomNumberSeed;
}
/**
* Set random number seed
*
* @param new random number seed; default = 1
*
* @return random number seed
*/
unsigned long int setRandomNumberSeed( unsigned long int seed = 1 );
_randomNumberSeed = seed;
return _randomNumberSeed;
}
/**
* Generate a random number using algorithm in Gromacs
*
* @param ig seed
*
* @return random number
*
*/
BrookOpenMMFloat BrookRandomNumberGenerator::generateGromacsRandomNumber( int* ig ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::generateGromacsRandomNumber";
int irand;
int m = 100000000;
float rm = 100000000.0; /* same number as m, but real format */
int m1 = 10000;
int mult = 31415821;
BrookOpenMMFloat r;
int irandh,irandl,multh,multl;
// ---------------------------------------------------------------------------------------
irand = abs(*ig) % m;
/* multiply irand by mult, but take into account that overflow
* must be discarded, and do not generate an error.
*/
irandh = irand / m1;
irandl = irand % m1;
multh = mult / m1;
multl = mult % m1;
irand = ((irandh*multl+irandl*multh) % m1) * m1 + irandl*multl;
irand = (irand + 1) % m;
/* convert irand to a real random number between 0 and 1. */
r = (BrookOpenMMFloat) (irand / 10);
r = r * 10 / rm;
if ((r <= 0) || (r > 1))
r = 0.0;
*ig = irand;
return r;
}
inline int MaxInt( unsigned int x, unsigned int y ){ return x > y ? x : y; }
/**
* Generate a random number using algorithm in Nvidia 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 BrookRandomNumberGenerator::generateRandomsAlaNvidia( float* randomV1, float* randomV2, float* randomV3,
unsigned int state[4] ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::generateRandomsAlaNvidia";
unsigned int carry = 0;
// ---------------------------------------------------------------------------------------
state[0] = state[0] * 69069 + 1;
state[1] ^= state[1] << 13;
state[1] ^= state[1] >> 17;
state[1] ^= state[1] << 5;
unsigned int k = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
unsigned int m = state[3] + state[3] + state[2] + carry;
state[2] = state[3];
state[3] = m;
carry = k >> 30;
unsigned int z1 = MaxInt(state[0] + state[1] + state[3], 0x00000001);
float x1 = (float) ( (double) z1 / (double)UINT_MAX );
/*
if( x1 < 0.0f ){
unsigned int z1 = MaxInt(state[0] + state[1] + state[3], 0x00000001);
double z2 = (double) z1/((double) UINT_MAX);
(void) fprintf( logFile, "x1=%.6e state[%u %u %u] sum %u %.6e den=%u z2=%.6e\n",
x1, state[0], state[1], state[3],
(state[0] + state[1] + state[3]), (float) (state[0] + state[1] + state[3]), z1, z2 );
}
*/
state[0] = state[0] * 69069 + 1;
state[1] ^= state[1] << 13;
state[1] ^= state[1] >> 17;
state[1] ^= state[1] << 5;
x1 = sqrt(-2.0f * log(x1));
k = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
m = state[3] + state[3] + state[2] + carry;
state[2] = state[3];
state[3] = m;
carry = k >> 30;
float x2 = (float)(state[0] + state[1] + state[3]) / (float)UINT_MAX;
state[0] = state[0] * 69069 + 1;
state[1] ^= state[1] << 13;
state[1] ^= state[1] >> 17;
state[1] ^= state[1] << 5;
*randomV1 = x1 * cos(2.0f * 3.14159265f * x2);
k = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
m = state[3] + state[3] + state[2] + carry;
state[2] = state[3];
state[3] = m;
carry = k >> 30;
unsigned int z3 = MaxInt(state[0] + state[1] + state[3], 0x00000001);
float x3 = (float) ( (double) z3 / (double)UINT_MAX );
//float x3 = (float)MaxInt(state[0] + state[1] + state[3], 0x00000001) / (float)UINT_MAX;
state[0] = state[0] * 69069 + 1;
state[1] ^= state[1] << 13;
state[1] ^= state[1] >> 17;
state[1] ^= state[1] << 5;
x3 = sqrt(-2.0f * log(x3));
k = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
m = state[3] + state[3] + state[2] + carry;
state[2] = state[3];
state[3] = m;
carry = k >> 30;
float x4 = (float)(state[0] + state[1] + state[3]) / (float)UINT_MAX;
state[0] = state[0] * 69069 + 1;
state[1] ^= state[1] << 13;
state[1] ^= state[1] >> 17;
state[1] ^= state[1] << 5;
*randomV2 = x3 * cos(2.0f * 3.14159265f * x4);
k = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
m = state[3] + state[3] + state[2] + carry;
state[2] = state[3];
state[3] = m;
carry = k >> 30;
//float x5 = (float)MaxInt(state[0] + state[1] + state[3], 0x00000001) / (float)UINT_MAX;
unsigned int z5 = MaxInt(state[0] + state[1] + state[3], 0x00000001);
float x5 = (float) ( (double) z5 / (double)UINT_MAX );
state[0] = state[0] * 69069 + 1;
state[1] ^= state[1] << 13;
state[1] ^= state[1] >> 17;
state[1] ^= state[1] << 5;
x5 = sqrt(-2.0f * log(x5));
k = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
m = state[3] + state[3] + state[2] + carry;
state[2] = state[3];
state[3] = m;
carry = k >> 30;
float x6 = (float)(state[0] + state[1] + state[3]) / (float)UINT_MAX;
*randomV3 = x5 * cos(2.0f * 3.14159265f * x6);
//(void) fprintf( logFile, "rv=%.6e %.6e %.6e\n", *randomV1, *randomV2, *randomV3 );
//exit(0);
return;
}
/**
* Load random number streams using Nvidia algorithm
*
*
* @return DefaultReturnValue;
*/
int BrookRandomNumberGenerator::loadRandomNumberStreamsNvidia( void ){
// ---------------------------------------------------------------------------------------
static float *buf = NULL;
static unsigned int state[4];
static int stateInitialized = 0;
static const int reseed = 5;
// static const char* methodName = "\nBrookRandomNumberGenerator::loadGVStreamsNvidia";
// ---------------------------------------------------------------------------------------
// periodically reset seeds
if( !stateInitialized || !(stateInitialized % reseed) ){
state[0] = rand();
state[1] = rand();
state[2] = rand();
state[3] = rand();
if( getVerbosity() && getLog() ){
(void) fprintf( getLog(), "LoadGVStreamsNvidia: reset state seeds stateInitialized=%d reseed=%d\n",
stateInitialized, reseed );
(void) fflush( getLog() );
}
/*
state[0] = 9578;
state[1] = 29245;
state[2] = 16266;
state[3] = 27587;
*/
}
stateInitialized++;
// allocate memory once for download of random nos.
if( buf == NULL ){
buf = (float*) malloc( sizeof(float) * 3 * getRandomNumberStreamSize() );
}
if( getVerbosity() && getLog() ){
static float count = 0.0f;
float block = (float) (3*sdp->gvSize);
count += 1.0f;
(void) fprintf( getLog(), "LoadGVStreamsNvidia: count=%.1f ttl=%.3e no./count=%.1f %d %d\n",
count, block*count, block, sdp->gvSize, NGVSTREAMS );
(void) fflush( getLog() );
}
for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){
float v1,v2,v3;
generateRandomsAlaNvidia( &v1, &v2, &v3, state, NULL );
buf[ii] = v1;
buf[ii+1] = v2;
buf[ii+2] = v3;
}
getRandomNumberStream( jj )->loadFromArray( buf );
}
return DefaultReturnValue;
}
/**
* Load random number streams using original gpu algorithm
*
*
* @return DefaultReturnValue;
*/
int BrookRandomNumberGenerator::loadGVStreamsOriginal( void ){
// ---------------------------------------------------------------------------------------
int i, j;
static float *buf = NULL;
unsigned long int jran;
// static const char* methodName = "\nBrookRandomNumberGenerator::LoadGVStreamsOriginal";
// ---------------------------------------------------------------------------------------
if( buf == NULL ){
buf = (float*) malloc( sizeof(float) * 3 * sdp->gvSize );
}
jran = getRandomNumberSeed();
for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){
//buf[i] = sdp->fgauss( &jran );
buf[i] = generateGromacsRandomNumber( &jran );
}
getRandomNumberStream( jj )->loadFromArray( buf );
}
incrementRandomNumberSeed( 1 );
return DefaultReturnValue;
}
/**
* 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 BrookRandomNumberGenerator::loadGVShuffle( void ){
// ---------------------------------------------------------------------------------------
const int p[] = { 2, 3, 5, 7, 11 };
const int np = sizeof(p) / sizeof(p[0]);
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";
// ---------------------------------------------------------------------------------------
if( buf == NULL ){
buf = (float*) malloc( sizeof(float) * sdp->gvSize );
indices = (int*) malloc( sizeof(int) * pmax );
}
int rvSize = getRandomNumberStreamSize();
for ( i = 0; i < rvSize; i++ ) {
buf[i] = (float) i;
}
//How to come up with this number here?
unsigned long int seed = getRandomNumberSeed();
for ( iter = 0; iter < 1000000; iter++ ) {
//for each p
for ( i = 0; i < np; i++ ){
//Come up with p random indices
//Not checking that they are distinct
//because that should be fairly rare
for ( j = 0; j < p[i]; j++ ) {
//indices[j] = (int) ( gmx_rando( &sdp->seed ) * sdp->gvSize );
indices[j] = (int) ( generateGromacsRandomNumber( &seed )*rvSize );
}
//do a p-cycle
tmp = buf[ indices[0] ];
for ( j = 0; j < p[i]-1; j++ ) {
buf[ indices[j] ] = buf[ indices[j+1] ];
}
buf[ indices[ p[i]-1 ] ] = tmp;
}
}
getShuffleStream()->loadFromArray( buf );
return DefaultReturnValue;
}
/**
* Shuffle streams
*
* @return DefaultReturnValue;
*/
int BrookRandomNumberGenerator::shuffleGVStreams( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::shuffleGVStreams";
// ---------------------------------------------------------------------------------------
int numberOfRvStreams = getNumberOfRandomNumberStreams();
for( int ii = 0; ii < numberOfRvStreams - 1; ii++ ){
kpermute_vectors( (float) , sdp->strShuffle, sdp->strVGauss[i+1],
sdp->strVGauss[i] );
}
for( int jj = 0; jj < numberOfRvStreams - 1; jj++ ){
kpermute_vectors( (float) sdp->gvWidth, sdp->strShuffle, sdp->strVGauss[0],
sdp->strVGauss[ NGVSTREAMS-1 ] );
}
return DefaultReturnValue;
}
static int
PrintGVStreams( gpuContext gpu )
{
int ii;
gpuSDParams sdp = gpu->sdparams;
char testName[128];
char fileName[128];
/* dump inputs */
(void) sprintf( testName, "ShuffleGVStreams" );
(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++ ){
(void) sprintf( testName, "ShuffleGVStreams%d", ii );
(void) sprintf( fileName, "%s.strVGaussIn.in", testName );
SaveStream( fileName, sdp->strVGauss[ii] );
}
return 1;
}
/* Advances the current position by 2*gpu->natoms
* If we run out of rand numbers, we shuffle and
* reuse a few times before reloading from the cpu
* */
static int
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 ){
time = clock();
}
sdp->gvOffset += 2*gpu->natoms;
//Check if we've used up this texture
if ( sdp->gvOffset > sdp->gvSize - 2*gpu->natoms ) {
//Next one if available
sdp->gvOffset = 0;
if ( sdp->gvCurStream < NGVSTREAMS - 1 ) {
sdp->gvCurStream++;
//printf( "Using stream %d\n", sdp->gvCurStream );
//fflush( stdout );
} else {
//No more textures, need to shuffle
if ( sdp->gvNumShuffles < NGVSHUFFLE ) {
//printf( "Shuffling streams\n" );
//fflush( stdout );
clock_t localTime;
if( trackTime ){
localTime = clock();
}
ShuffleGVStreams( sdp, gpu );
if( trackTime ){
shuffleTime += clock() - localTime;
}
sdp->gvNumShuffles++;
} else { //Need to refresh random numbers from cpu
static int reducePrint = 0;
fflush( gpu->log );
if( UseOriginalRng ){
LoadGVStreamsOriginal( sdp );
} else {
LoadGVStreamsNvidia( sdp, gpu->log );
}
sdp->gvNumShuffles = 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;
}
}
/*
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;
};
/**
* Update
*
* @param positions atom positions
* @param velocities atom velocities
* @param forces atom forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
*
* @return DefaultReturnValue
*
*/
int BrookRandomNumberGenerator::update( Stream& positions, Stream& velocities,
const Stream& forces,
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 array of derived parameters indexed by 'DerivedParameters' enums
*
* @return array
*
*/
const BrookOpenMMFloat* BrookRandomNumberGenerator::getDerivedParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::getDerivedParameters";
// ---------------------------------------------------------------------------------------
return _derivedParameters;
}
/**
* Get Atom stream size
*
* @return Atom stream size
*
*/
int BrookRandomNumberGenerator::getStochasticDynamicsAtomStreamSize( void ) const {
return _sdAtomStreamSize;
}
/**
* Get atom stream width
*
* @return atom stream width
*
*/
int BrookRandomNumberGenerator::getStochasticDynamicsAtomStreamWidth( void ) const {
return _sdAtomStreamWidth;
}
/**
* Get atom stream height
*
* @return atom stream height
*/
int BrookRandomNumberGenerator::getStochasticDynamicsAtomStreamHeight( void ) const {
return _sdAtomStreamHeight;
}
/**
* Get SDPC1 stream
*
* @return SDPC1 stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSDPC1Stream( void ) const {
return _sdStreams[SDPC1Stream];
}
/**
* Get SDPC2 stream
*
* @return SDPC2 stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSDPC2Stream( void ) const {
return _sdStreams[SDPC2Stream];
}
/**
* Get SD2X stream
*
* @return SD2X stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSD2XStream( void ) const {
return _sdStreams[SD2XStream];
}
/**
* Get SD1V stream
*
* @return SD1V stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getSD1VStream( void ) const {
return _sdStreams[SD1VStream];
}
/**
* Get VPrime stream
*
* @return Vprime stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getVPrimeStream( void ) const {
return _sdStreams[VPrimeStream];
}
/**
* Get XPrime stream
*
* @return Xprime stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getXPrimeStream( void ) const {
return _sdStreams[XPrimeStream];
}
/**
* Get InverseMass stream
*
* @return inverse mass stream
*
*/
BrookFloatStreamInternal* BrookRandomNumberGenerator::getInverseMassStream( void ) const {
return _sdStreams[InverseMassStream];
}
/**
* Initialize stream dimensions
*
* @param numberOfAtoms number of atoms
* @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int BrookRandomNumberGenerator::_initializeStreamSizes( int numberOfAtoms, const Platform& platform ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookRandomNumberGenerator::_initializeStreamSizes";
// ---------------------------------------------------------------------------------------
_sdAtomStreamSize = getAtomStreamSize( platform );
_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:
returnString = "Sd2pc1";
break;
case Sd2pc2:
returnString = "Sd2pc2";
break;
default:
returnString = "Unknown";
break;
}
return returnString;
}
/**
* Initialize streams
*
* @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int BrookRandomNumberGenerator::_initializeStreams( const Platform& platform ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookRandomNumberGenerator::_initializeStreams";
BrookOpenMMFloat dangleValue = (BrookOpenMMFloat) 0.0;
// ---------------------------------------------------------------------------------------
int sdAtomStreamSize = getStochasticDynamicsAtomStreamSize();
int sdAtomStreamWidth = getStochasticDynamicsAtomStreamWidth();
_sdStreams[SDPC1Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC1Stream,
sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float2, dangleValue );
_sdStreams[SDPC2Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC2Stream,
sdAtomStreamSize, sdAtomStreamWidth,
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 );
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]) );
sdpc[0][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[V]) );
sdpc[1][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yx]) );
sdpc[1][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[X]) );
}
_sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] );
_sdStreams[SDPC2Stream]->loadFromArray( sdpc[1] );
for( int ii = 0; ii < 2; ii++ ){
delete[] sdpc[ii];
}
// initialize SD2X
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;
// ---------------------------------------------------------------------------------------
// setup inverse sqrt masses
_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;
}
/*
* Setup of StochasticDynamics parameters
*
* @param masses masses
* @param platform Brook platform
*
* @return nonzero value if error
*
* */
int BrookRandomNumberGenerator::setup( const std::vector<double>& masses, const Platform& platform ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookRandomNumberGenerator::setup";
// ---------------------------------------------------------------------------------------
int numberOfAtoms = (int) masses.size();
setNumberOfAtoms( numberOfAtoms );
// set stream sizes and then create streams
_initializeStreamSizes( numberOfAtoms, platform );
_initializeStreams( platform );
_setInverseSqrtMasses( masses );
return DefaultReturnValue;
}
/*
* Get contents of object
*
* @param level level of dump
*
* @return string containing contents
*
* */
std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookRandomNumberGenerator::getContentsString";
static const unsigned int MAX_LINE_CHARS = 256;
char value[MAX_LINE_CHARS];
static const char* Set = "Set";
static const char* NotSet = "Not set";
// ---------------------------------------------------------------------------------------
std::stringstream message;
std::string tab = " ";
#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
(void) LOCAL_SPRINTF( value, "%d", getNumberOfAtoms() );
message << _getLine( tab, "Number of atoms:", value );
(void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() );
message << _getLine( tab, "Atom stream width:", value );
(void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() );
message << _getLine( tab, "Atom stream height:", value );
(void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
message << _getLine( tab, "Atom stream size:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getTau() );
message << _getLine( tab, "Tau:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getTemperature() );
message << _getLine( tab, "Temperature:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getStepSize() );
message << _getLine( tab, "Step size:", value );
const BrookOpenMMFloat* derivedParameters = getDerivedParameters();
for( int ii = 0; ii < MaxDerivedParameters; ii++ ){
(void) LOCAL_SPRINTF( value, "%.5e", derivedParameters[ii] );
message << _getLine( tab, _getDerivedParametersString( ii ), value );
}
(void) LOCAL_SPRINTF( value, "%.5f", getTemperature() );
message << _getLine( tab, "Temperature:", value );
message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) );
message << _getLine( tab, "SDPC1:", (getSDPC1Stream() ? 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++ ){
message << std::endl;
if( _sdStreams[ii] ){
message << _sdStreams[ii]->getContentsString( );
}
}
#undef LOCAL_SPRINTF
return message.str();
}
platforms/brook/src/BrookRandomNumberGenerator.h 0 → 100644
View file @ a475f0a6
#ifndef OPENMM_BROOK_RANDOM_NUMBER_GENERATOR_H_
#define OPENMM_BROOK_RANDOM_NUMBER_GENERATOR_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include <vector>
#include <set>
#include "BrookFloatStreamInternal.h"
#include "BrookShakeAlgorithm.h"
#include "BrookPlatform.h"
#include "BrookCommon.h"
namespace OpenMM {
/**
*
* Encapsulates stochastic dynamics algorithm
*
*/
class BrookRandomNumberGenerator : public BrookCommon {
public:
/**
* Constructor
*
*/
BrookRandomNumberGenerator( );
/**
* Destructor
*
*/
~BrookRandomNumberGenerator();
/**
* Get tau
*
* @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
*
*/
const BrookOpenMMFloat* getDerivedParameters( void ) const;
/**
* Get StochasticDynamics atom stream width
*
* @return atom stream width
*/
int getStochasticDynamicsAtomStreamWidth( void ) const;
/**
* Get StochasticDynamics atom stream height
*
* @return atom stream height
*/
int getStochasticDynamicsAtomStreamHeight( void ) const;
/**
* Get StochasticDynamics atom stream size
*
* @return atom stream size
*/
int getStochasticDynamicsAtomStreamSize( void ) const;
/**
* Update parameters
*
* @param temperature temperature
* @param friction friction
* @param step size step size
*
* @return DefaultReturnValue
*
*/
int updateParameters( double temperature, double friction, double stepSize );
/**
* 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
*
* @return array ofstreams
*
*/
BrookFloatStreamInternal** getStreams( void );
/*
* Setup of StochasticDynamics parameters
*
* @param masses atom masses
* @param platform Brook platform
*
* @return ErrorReturnValue value if error, else DefaultReturnValue
*
* */
int setup( const std::vector<double>& masses, const Platform& platform );
/*
* Get contents of object
*
* @param level of dump
*
* @return string containing contents
*
* */
std::string getContentsString( int level = 0 ) const;
/**
* Get SDPC1 stream
*
* @return SDPC1 stream
*
*/
BrookFloatStreamInternal* getSDPC1Stream( void ) const;
/**
* Get SDPC2 stream
*
* @return SDPC2 stream
*
*/
BrookFloatStreamInternal* getSDPC2Stream( void ) const;
/**
* 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( int* ig );
/**
* Generate a random number using algorithm in Nvidia 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 generateRandomsAlaNvidia( float* randomV1, float* randomV2, float* randomV3,
unsigned int state[4] );
/**
* Load random number streams using Nvidia algorithm
*
*
* @return DefaultReturnValue;
*/
int loadRandomNumberStreamsNvidia( void );
/**
* Get random number seed
*
* @return random number seed
*/
unsigned long int getRandomNumberSeed( void ) const;
/**
* Increment random number seed
*
* @param increment amount to increment random number seed; default = 1
*
* @return updated random number seed
*/
unsigned long int incrementRandomNumberSeed( unsigned long int increment = 1 );
/**
* Set random number seed
*
* @param new random number seed; default = 1
*
* @return random number seed
*/
unsigned long int setRandomNumberSeed( unsigned long int seed = 1 );
private:
// streams indices
enum BrookRandomNumberGeneratorStreams {
RandomNumberStream,
ShuffleStream,
LastStreamIndex
};
// randomNumberSeed
unsigned long int _randomNumberSeed;
// number of random number streams
int _numberOfRandomNumberStreams;
// random number stream dimensions
int _randomNumberStreamWidth;
int _randomNumberStreamHeight;
int _randomNumberStreamSize;
/*
* Setup of stream dimensions
*
* @param atomStreamSize atom stream size
* @param atomStreamWidth atom stream width
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
* */
int _initializeStreamSizes( int atomStreamSize, int atomStreamWidth );
/**
* Initialize stream dimensions
*
* @param numberOfAtoms number of atoms
* @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int _initializeStreamSizes( int numberOfAtoms, const Platform& platform );
/**
* Initialize stream dimensions and streams
*
* @param platform platform
*
* @return nonzero value if error
*
*/
int _initializeStreams( const Platform& platform );
/**
* Set masses
*
* @param masses atomic masses
*
*/
int _setInverseSqrtMasses( const std::vector<double>& masses );
};
} // namespace OpenMM
#endif /* OPENMM_BROOK_RANDOM_NUMBER_GENERATOR_H_ */
platforms/brook/src/BrookStochasticDynamics.cpp
View file @ a475f0a6
...@@ -341,10 +341,139 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -341,10 +341,139 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// unused Shake parameter
float omega = 1.0f;
// static const char* methodName = "\nBrookStochasticDynamics::update"; // static const char* methodName = "\nBrookStochasticDynamics::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; return DefaultReturnValue;
...@@ -445,6 +574,39 @@ BrookFloatStreamInternal* BrookStochasticDynamics::getSD1VStream( void ) const { ...@@ -445,6 +574,39 @@ BrookFloatStreamInternal* BrookStochasticDynamics::getSD1VStream( void ) const {
return _sdStreams[SD1VStream]; return _sdStreams[SD1VStream];
} }
/**
* Get VPrime stream
*
* @return Vprime stream
*
*/
BrookFloatStreamInternal* BrookStochasticDynamics::getVPrimeStream( void ) const {
return _sdStreams[VPrimeStream];
}
/**
* Get XPrime stream
*
* @return Xprime stream
*
*/
BrookFloatStreamInternal* BrookStochasticDynamics::getXPrimeStream( void ) const {
return _sdStreams[XPrimeStream];
}
/**
* Get InverseMass stream
*
* @return inverse mass stream
*
*/
BrookFloatStreamInternal* BrookStochasticDynamics::getInverseMassStream( void ) const {
return _sdStreams[InverseMassStream];
}
/** /**
* Initialize stream dimensions * Initialize stream dimensions
* *
...@@ -607,6 +769,18 @@ int BrookStochasticDynamics::_initializeStreams( const Platform& platform ){ ...@@ -607,6 +769,18 @@ int BrookStochasticDynamics::_initializeStreams( const Platform& platform ){
sdAtomStreamSize, sdAtomStreamWidth, sdAtomStreamSize, sdAtomStreamWidth,
BrookStreamInternal::Float3, dangleValue ); 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 );
return DefaultReturnValue; return DefaultReturnValue;
} }
......
platforms/brook/src/BrookStochasticDynamics.h
View file @ a475f0a6
...@@ -241,6 +241,9 @@ class BrookStochasticDynamics : public BrookCommon { ...@@ -241,6 +241,9 @@ class BrookStochasticDynamics : public BrookCommon {
SDPC2Stream, SDPC2Stream,
SD2XStream, SD2XStream,
SD1VStream, SD1VStream,
VPrimeStream,
XPrimeStream,
InverseMassStream,
LastStreamIndex LastStreamIndex
}; };
......
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