Skip to content

GitLab

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

Mods

parent 0e675a19
Hide whitespace changes
Inline Side-by-side
Showing with 938 additions and 881 deletions
platforms/brook/src/BrookBonded.cpp
View file @ adb4ec5a
...@@ -1805,13 +1805,13 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -1805,13 +1805,13 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
static const std::string methodName = "BrookBonded::computeForces"; static const std::string methodName = "BrookBonded::computeForces";
static const int PrintOn = 0;
static const int I_Stream = 0; static const int I_Stream = 0;
static const int J_Stream = 1; static const int J_Stream = 1;
static const int K_Stream = 2; static const int K_Stream = 2;
static const int L_Stream = 3; static const int L_Stream = 3;
static const int PrintOn = 0;
static const int MaxErrorMessages = 2; static const int MaxErrorMessages = 2;
static int ErrorMessages = 0; static int ErrorMessages = 0;
......
platforms/brook/src/BrookCalcKineticEnergyKernel.cpp
View file @ adb4ec5a
...@@ -141,34 +141,48 @@ double BrookCalcKineticEnergyKernel::execute( OpenMMContextImpl& context ){ ...@@ -141,34 +141,48 @@ double BrookCalcKineticEnergyKernel::execute( OpenMMContextImpl& context ){
void* dataV = _openMMBrookInterface.getParticleVelocities()->getData( 1 ); void* dataV = _openMMBrookInterface.getParticleVelocities()->getData( 1 );
float* velocity = (float*) dataV; float* velocity = (float*) dataV;
double energy = 0.0;
int index = 0;
/* if( 0 && _numberOfParticles ){
printf( "BrookCalcKineticEnergyKernel:\n" ); printf( "BrookCalcKineticEnergyKernel:\n" );
float com[3] = { 0.0, 0.0, 0.0 }; float com[3] = { 0.0, 0.0, 0.0 };
for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){ float localEnergy = 0.0f;
com[0] += velocity[index]; int localIndex = 0;
com[1] += velocity[index+1]; float massSum = 0.0f;
com[2] += velocity[index+2]; for ( int ii = 0; ii < _numberOfParticles; ii++, localIndex += 3 ){
printf( " %d %.3f [%12.5e %12.5e %12.5e]\n", ii, _masses[ii], velocity[index], velocity[index+1], velocity[index+2] );
com[0] += _masses[ii]*velocity[localIndex];
com[1] += _masses[ii]*velocity[localIndex+1];
com[2] += _masses[ii]*velocity[localIndex+2];
localEnergy += _masses[ii]*(velocity[localIndex]*velocity[localIndex] + velocity[localIndex + 1]*velocity[localIndex + 1] + velocity[localIndex + 2]*velocity[localIndex + 2]);
massSum += _masses[ii];
printf( " %d %.3f [%12.5e %12.5e %12.5e]\n", ii, _masses[ii], velocity[localIndex], velocity[localIndex+1], velocity[localIndex+2] );
} }
printf( "Com [%12.5e %12.5e %12.5e]\n", com[0], com[1], com[2] ); float inverseTotalMass = 1.0f/massSum;
index = 0; com[0] *= inverseTotalMass;
com[1] *= inverseTotalMass;
com[2] *= inverseTotalMass;
printf( "KE raw=%.5e Com [%12.5e %12.5e %12.5e]\n", 0.5f*localEnergy, com[0], com[1], com[2] );
float newcom[3] = { 0.0, 0.0, 0.0 }; float newcom[3] = { 0.0, 0.0, 0.0 };
for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){ localIndex = 0;
velocity[index] -= com[0]; for ( int ii = 0; ii < _numberOfParticles; ii++, localIndex += 3 ){
velocity[index+1] -= com[1]; velocity[localIndex] -= com[0];
velocity[index+2] -= com[2]; velocity[localIndex+1] -= com[1];
newcom[0] += velocity[index]; velocity[localIndex+2] -= com[2];
newcom[1] += velocity[index+1]; newcom[0] += velocity[localIndex];
newcom[2] += velocity[index+2]; newcom[1] += velocity[localIndex+1];
newcom[2] += velocity[localIndex+2];
printf( " %d %.3f [%12.5e %12.5e %12.5e]\n", ii, _masses[ii], velocity[localIndex], velocity[localIndex+1], velocity[localIndex+2] );
} }
printf( "NewCom [%12.5e %12.5e %12.5e]\n", newcom[0], newcom[1], newcom[2] ); printf( "NewCom [%12.5e %12.5e %12.5e]\n", newcom[0], newcom[1], newcom[2] );
index = 0;
*/
}
int index = 0;
double energy = 0.0;
for( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){ for( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
energy += _masses[ii]*(velocity[index]*velocity[index] + velocity[index + 1]*velocity[index + 1] + velocity[index + 2]*velocity[index + 2]); energy += _masses[ii]*(velocity[index]*velocity[index] + velocity[index + 1]*velocity[index + 1] + velocity[index + 2]*velocity[index + 2]);
} }
......
platforms/brook/src/BrookIntegrateVerletStepKernel.cpp
View file @ adb4ec5a
...@@ -118,7 +118,7 @@ int BrookIntegrateVerletStepKernel::setLog( FILE* log ){ ...@@ -118,7 +118,7 @@ int BrookIntegrateVerletStepKernel::setLog( FILE* log ){
* *
*/ */
void BrookIntegrateVerletStepKernel::initialize( const System& system, const VerletIntegrator& integrator ){ void BrookIntegrateVerletStepKernel::initialize( const System& system, const VerletIntegrator& integrator ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -166,9 +166,9 @@ void BrookIntegrateVerletStepKernel::initialize( const System& system, const Ve ...@@ -166,9 +166,9 @@ void BrookIntegrateVerletStepKernel::initialize( const System& system, const Ve
_brookShakeAlgorithm->setLog( log ); _brookShakeAlgorithm->setLog( log );
_brookShakeAlgorithm->setup( masses, constraintIndicesVector, constraintLengths, getPlatform() ); _brookShakeAlgorithm->setup( masses, constraintIndicesVector, constraintLengths, getPlatform() );
BrookOpenMMFloat tolerance = static_cast<BrookOpenMMFloat>( integrator.getConstraintTolerance() ); BrookOpenMMFloat tolerance = static_cast<BrookOpenMMFloat>( integrator.getConstraintTolerance() );
_brookShakeAlgorithm->setShakeTolerance( tolerance ); _brookShakeAlgorithm->setShakeTolerance( tolerance );
_brookShakeAlgorithm->setMaxIterations( 30 ); _brookShakeAlgorithm->setMaxIterations( 40 );
if( printOn && log ){ if( printOn && log ){
(void) fprintf( log, "%s done w/ setup: particles=%d const=%d\n", methodName.c_str(), numberOfParticles, numberOfConstraints ); (void) fprintf( log, "%s done w/ setup: particles=%d const=%d\n", methodName.c_str(), numberOfParticles, numberOfConstraints );
......
platforms/brook/src/BrookLangevinDynamics.cpp
View file @ adb4ec5a
...@@ -299,8 +299,22 @@ int BrookLangevinDynamics::_updateDerivedParameters( void ){ ...@@ -299,8 +299,22 @@ int BrookLangevinDynamics::_updateDerivedParameters( void ){
_derivedParameters[Sd1pc1] = tau*( one - _derivedParameters[EM] ); _derivedParameters[Sd1pc1] = tau*( one - _derivedParameters[EM] );
_derivedParameters[Sd1pc2] = tau*( _derivedParameters[EPH] - _derivedParameters[EMH] ); _derivedParameters[Sd1pc2] = tau*( _derivedParameters[EPH] - _derivedParameters[EMH] );
if( fabsf( _derivedParameters[Sd1pc2] ) < 1.0e-06 ){
_derivedParameters[Sd1pc2] = tau*_derivedParameters[GDT];
}
_derivedParameters[Sd1pc3] = _derivedParameters[D]/(tau*_derivedParameters[C] ); _derivedParameters[Sd1pc3] = _derivedParameters[D]/(tau*_derivedParameters[C] );
_derivedParameters[Sd2pc1] = one/_derivedParameters[Sd1pc2];
// if tau was greater than 20000, then _derivedParameters[Sd1pc2] was zero
// using tau*_derivedParameters[GDT] as approximation to tau*( _derivedParameters[EPH] - _derivedParameters[EMH] ) for values in this
// range
if( fabsf( _derivedParameters[Sd1pc2] ) > 1.0e-10 ){
_derivedParameters[Sd2pc1] = one/_derivedParameters[Sd1pc2];
} else {
std::stringstream message;
message << methodName << " Sd1pc2=" << _derivedParameters[Sd1pc2] << " is too small (Sd2pc1 is 1/Sd1pc2)";
throw OpenMMException( message.str() );
}
_derivedParameters[Sd2pc2] = tau*_derivedParameters[D]/( _derivedParameters[EM] - one ); _derivedParameters[Sd2pc2] = tau*_derivedParameters[D]/( _derivedParameters[EM] - one );
return DefaultReturnValue; return DefaultReturnValue;
...@@ -340,7 +354,7 @@ int BrookLangevinDynamics::updateParameters( double temperature, double friction ...@@ -340,7 +354,7 @@ int BrookLangevinDynamics::updateParameters( double temperature, double friction
if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){ if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
std::string contents = getContentsString( ); std::string contents = getContentsString( );
(void) fprintf( getLog(), "%s contents\n%s", methodName, contents.c_str() ); (void) fprintf( getLog(), "%s contents\n%s", methodName.c_str(), contents.c_str() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -349,468 +363,159 @@ int BrookLangevinDynamics::updateParameters( double temperature, double friction ...@@ -349,468 +363,159 @@ int BrookLangevinDynamics::updateParameters( double temperature, double friction
}; };
/**
* Update /**
*
* @param positions particle positions
* @param velocities particle velocities
* @param forces particle forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
* @param brookRandomNumberGenerator BrookRandomNumberGenerator reference
* *
* @return DefaultReturnValue * Get array of derived parameters indexed by 'DerivedParameters' enums
*
* @return array
* *
*/ */
const BrookOpenMMFloat* BrookLangevinDynamics::getDerivedParameters( void ) const {
int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamImpl& velocityStream, // ---------------------------------------------------------------------------------------
BrookStreamImpl& forceStream,
BrookShakeAlgorithm& brookShakeAlgorithm,
BrookRandomNumberGenerator& brookRandomNumberGenerator ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nBrookLangevinDynamics::update"; // static const std::string methodName = "\nBrookLangevinDynamics::getDerivedParameters";
static const int PrintOn = 0;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
const BrookOpenMMFloat* derivedParameters = getDerivedParameters(); return _derivedParameters;
}
if( (1 || PrintOn) && getLog() ){ /**
* Get Particle stream size
*
* @return Particle stream size
*
*/
static int showAux = 1; int BrookLangevinDynamics::getLangevinDynamicsParticleStreamSize( void ) const {
return _sdParticleStreamSize;
}
if( PrintOn ){ /**
(void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() ); * Get particle stream width
(void) fflush( getLog() ); *
} * @return particle stream width
*
*/
// show update int BrookLangevinDynamics::getLangevinDynamicsParticleStreamWidth( void ) const {
return _sdParticleStreamWidth;
if( showAux ){ }
showAux = 0;
std::string contents = brookRandomNumberGenerator.getContentsString( ); /**
(void) fprintf( getLog(), "%s RNG contents\n%s", methodName, contents.c_str() ); * Get particle stream height
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); *
* @return particle stream height
*/
contents = brookShakeAlgorithm.getContentsString( ); int BrookLangevinDynamics::getLangevinDynamicsParticleStreamHeight( void ) const {
(void) fprintf( getLog(), "%s Shake contents\n%s", methodName, contents.c_str() ); return _sdParticleStreamHeight;
}
(void) fflush( getLog() ); /**
} * Get SDPC1 stream
} *
* @return SDPC1 stream
*
*/
// first integration step BrookFloatStreamInternal* BrookLangevinDynamics::getSDPC1Stream( void ) const {
return _sdStreams[SDPC1Stream];
}
kupdate_sd1_fix1( /**
(float) getLangevinDynamicsParticleStreamWidth(), * Get SDPC2 stream
(float) brookRandomNumberGenerator.getRandomNumberStreamWidth(), *
(float) brookRandomNumberGenerator.getRvStreamOffset(), * @return SDPC2 stream
derivedParameters[EM], *
derivedParameters[Sd1pc1], */
derivedParameters[Sd1pc2],
derivedParameters[Sd1pc3],
getSDPC1Stream()->getBrookStream(),
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD2XStream()->getBrookStream(),
positionStream.getBrookStream(),
forceStream.getBrookStream(),
velocityStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
getSD1VStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getXPrimeStream()->getBrookStream()
);
// diagnostics
if( PrintOn ){ BrookFloatStreamInternal* BrookLangevinDynamics::getSDPC2Stream( void ) const {
(void) fprintf( getLog(), "\n%s Post kupdate_sd1_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d " return _sdStreams[SDPC2Stream];
"EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]", methodName, }
getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] );
(void) fprintf( getLog(), "\nSDPC1Stream\n" ); /**
getSDPC1Stream()->printToFile( getLog() ); * Get SD2X stream
*
* @return SD2X stream
*
*/
(void) fprintf( getLog(), "\nSD2XStream\n" ); BrookFloatStreamInternal* BrookLangevinDynamics::getSD2XStream( void ) const {
getSD2XStream()->printToFile( getLog() ); return _sdStreams[SD2XStream];
}
(void) fprintf( getLog(), "\nInverseMassStream\n" ); /**
getInverseMassStream()->printToFile( getLog() ); * Get SD1V stream
*
* @return SD1V stream
*
*/
//StreamImpl& positionStreamImpl = positionStream.getImpl(); BrookFloatStreamInternal* BrookLangevinDynamics::getSD1VStream( void ) const {
//const BrookStreamImpl brookPositions = dynamic_cast<BrookStreamImpl&> (positionStreamImpl); return _sdStreams[SD1VStream];
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl(); }
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStream\n" ); /**
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl(); * Get VPrime stream
brookStreamInternalF->printToFile( getLog() ); *
* @return Vprime stream
*
*/
BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl(); BrookFloatStreamInternal* BrookLangevinDynamics::getVPrimeStream( void ) const {
(void) fprintf( getLog(), "\nVelocityStream\n" ); return _sdStreams[VPrimeStream];
brookStreamInternalV->printToFile( getLog() ); }
(void) fprintf( getLog(), "\nSD1VStream\n" ); /**
getSD1VStream()->printToFile( getLog() ); * Get XPrime stream
*
* @return Xprime stream
*
*/
(void) fprintf( getLog(), "\nVPrimeStream\n" ); BrookFloatStreamInternal* BrookLangevinDynamics::getXPrimeStream( void ) const {
getVPrimeStream()->printToFile( getLog() ); return _sdStreams[XPrimeStream];
}
(void) fprintf( getLog(), "\nXPrimeStream\n" ); /**
getXPrimeStream()->printToFile( getLog() ); * Get InverseMass stream
*
* @return inverse mass stream
*
*/
(void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() ); BrookFloatStreamInternal* BrookLangevinDynamics::getInverseMassStream( void ) const {
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); return _sdStreams[InverseMassStream];
} }
// advance random number cursor /**
* Initialize stream dimensions
*
* @param numberOfParticles number of particles
* @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfParticles() ); int BrookLangevinDynamics::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
// first Shake step // ---------------------------------------------------------------------------------------
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){ //static const std::string methodName = "BrookLangevinDynamics::_initializeStreamSizes";
kshakeh_fix1( // ---------------------------------------------------------------------------------------
(float) brookShakeAlgorithm.getMaxIterations(),
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeTolerance(),
brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// first Shake gather
kshakeh_update1_fix1(
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() );
if( ( 0 || PrintOn) && getLog() ){ _sdParticleStreamSize = getParticleStreamSize( platform );
_sdParticleStreamWidth = getParticleStreamWidth( platform );
_sdParticleStreamHeight = getParticleStreamHeight( platform );
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d", return DefaultReturnValue;
methodName, getLangevinDynamicsParticleStreamWidth() ); }
(void) fprintf( getLog(), "\nShakeInverseMapStream\n" );
brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() );
}
}
// second integration step
kupdate_sd2_fix1(
(float) getLangevinDynamicsParticleStreamWidth(),
(float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
(float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1],
derivedParameters[Sd2pc2],
getSDPC2Stream()->getBrookStream(),
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD1VStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getSD2XStream()->getBrookStream(),
velocityStream.getBrookStream(),
getXPrimeStream()->getBrookStream()
);
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\n%s Post kupdate_sd2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"Sd2pc[]=[%12.5e %12.5e]", methodName,
getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
(void) fprintf( getLog(), "\nSDPC2Stream\n" );
getSDPC2Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nSD2XStream\n" );
getSD1VStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nVPrimeStream\n" );
getVPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\ngetSD2XStream\n" );
getSD2XStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalVel = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalVel->printToFile( getLog() );
(void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
}
// advance random number cursor
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfParticles() );
// second Shake step
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
kshakeh_fix1(
(float) brookShakeAlgorithm.getMaxIterations(),
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeTolerance(),
brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// second Shake gather
kshakeh_update2_fix1(
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
positionStream.getBrookStream() );
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"Sd2pc[]=[%12.5e %12.5e]", methodName,
getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
brookStreamInternalPos = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
}
} else {
//kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
}
return DefaultReturnValue;
};
/**
*
* Get array of derived parameters indexed by 'DerivedParameters' enums
*
* @return array
*
*/
const BrookOpenMMFloat* BrookLangevinDynamics::getDerivedParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nBrookLangevinDynamics::getDerivedParameters";
// ---------------------------------------------------------------------------------------
return _derivedParameters;
}
/**
* Get Particle stream size
*
* @return Particle stream size
*
*/
int BrookLangevinDynamics::getLangevinDynamicsParticleStreamSize( void ) const {
return _sdParticleStreamSize;
}
/**
* Get particle stream width
*
* @return particle stream width
*
*/
int BrookLangevinDynamics::getLangevinDynamicsParticleStreamWidth( void ) const {
return _sdParticleStreamWidth;
}
/**
* Get particle stream height
*
* @return particle stream height
*/
int BrookLangevinDynamics::getLangevinDynamicsParticleStreamHeight( void ) const {
return _sdParticleStreamHeight;
}
/**
* Get SDPC1 stream
*
* @return SDPC1 stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getSDPC1Stream( void ) const {
return _sdStreams[SDPC1Stream];
}
/**
* Get SDPC2 stream
*
* @return SDPC2 stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getSDPC2Stream( void ) const {
return _sdStreams[SDPC2Stream];
}
/**
* Get SD2X stream
*
* @return SD2X stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getSD2XStream( void ) const {
return _sdStreams[SD2XStream];
}
/**
* Get SD1V stream
*
* @return SD1V stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getSD1VStream( void ) const {
return _sdStreams[SD1VStream];
}
/**
* Get VPrime stream
*
* @return Vprime stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getVPrimeStream( void ) const {
return _sdStreams[VPrimeStream];
}
/**
* Get XPrime stream
*
* @return Xprime stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getXPrimeStream( void ) const {
return _sdStreams[XPrimeStream];
}
/**
* Get InverseMass stream
*
* @return inverse mass stream
*
*/
BrookFloatStreamInternal* BrookLangevinDynamics::getInverseMassStream( void ) const {
return _sdStreams[InverseMassStream];
}
/**
* Initialize stream dimensions
*
* @param numberOfParticles number of particles
* @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int BrookLangevinDynamics::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookLangevinDynamics::_initializeStreamSizes";
// ---------------------------------------------------------------------------------------
_sdParticleStreamSize = getParticleStreamSize( platform );
_sdParticleStreamWidth = getParticleStreamWidth( platform );
_sdParticleStreamHeight = getParticleStreamHeight( platform );
return DefaultReturnValue;
}
/** /**
* Initialize stream dimensions * Initialize stream dimensions
...@@ -911,274 +616,610 @@ std::string BrookLangevinDynamics::_getDerivedParametersString( int derivedParam ...@@ -911,274 +616,610 @@ std::string BrookLangevinDynamics::_getDerivedParametersString( int derivedParam
return returnString; return returnString;
} }
/** /**
* Initialize streams * Initialize streams
* *
* @param platform platform * @param platform platform
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int BrookLangevinDynamics::_initializeStreams( const Platform& platform ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookLangevinDynamics::_initializeStreams";
BrookOpenMMFloat dangleValue = (BrookOpenMMFloat) 0.0;
// ---------------------------------------------------------------------------------------
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize();
int sdParticleStreamWidth = getLangevinDynamicsParticleStreamWidth();
_sdStreams[SDPC1Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC1Stream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float2, dangleValue );
_sdStreams[SDPC2Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC2Stream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float2, dangleValue );
_sdStreams[SD2XStream] = new BrookFloatStreamInternal( BrookCommon::SD2XStream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[SD1VStream] = new BrookFloatStreamInternal( BrookCommon::SD1VStream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[VPrimeStream] = new BrookFloatStreamInternal( BrookCommon::VPrimeStream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[XPrimeStream] = new BrookFloatStreamInternal( BrookCommon::XPrimeStream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float3, dangleValue );
_sdStreams[InverseMassStream] = new BrookFloatStreamInternal( BrookCommon::InverseMassStream,
sdParticleStreamSize, sdParticleStreamWidth,
BrookStreamInternal::Float, dangleValue );
return DefaultReturnValue;
}
/**
* Update sd streams -- called after parameters change
*
* @return ErrorReturnValue if error, else DefaultReturnValue
*
*/
int BrookLangevinDynamics::_updateSdStreams( void ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookLangevinDynamics::_updateSdStreams";
// ---------------------------------------------------------------------------------------
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize();
BrookOpenMMFloat* sdpc[2];
for( int ii = 0; ii < 2; ii++ ){
sdpc[ii] = new BrookOpenMMFloat[2*sdParticleStreamSize];
memset( sdpc[ii], 0, 2*sdParticleStreamSize*sizeof( BrookOpenMMFloat ) );
}
BrookOpenMMFloat* inverseMass = new BrookOpenMMFloat[sdParticleStreamSize];
memset( inverseMass, 0, sdParticleStreamSize*sizeof( BrookOpenMMFloat ) );
const BrookOpenMMFloat* derivedParameters = getDerivedParameters( );
int numberOfParticles = getNumberOfParticles();
int index = 0;
for( int ii = 0; ii < numberOfParticles; 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]) );
inverseMass[ii] = _inverseSqrtMasses[ii]*_inverseSqrtMasses[ii];
}
_sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] );
_sdStreams[SDPC2Stream]->loadFromArray( sdpc[1] );
_sdStreams[InverseMassStream]->loadFromArray( inverseMass );
for( int ii = 0; ii < 2; ii++ ){
delete[] sdpc[ii];
}
delete[] inverseMass;
// initialize SD2X
BrookOpenMMFloat* sd2x = new BrookOpenMMFloat[3*sdParticleStreamSize];
//SimTKOpenMMUtilities::setRandomNumberSeed( (uint32_t) getRandomNumberSeed() );
memset( sd2x, 0, 3*sdParticleStreamSize*sizeof( BrookOpenMMFloat ) );
index = 0;
for( int ii = 0; ii < numberOfParticles; 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 particle masses
*
*/
int BrookLangevinDynamics::_setInverseSqrtMasses( const std::vector<double>& masses ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookLangevinDynamics::_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 LangevinDynamics parameters
*
* @param masses masses
* @param platform Brook platform
*
* @return nonzero value if error
*
* */
int BrookLangevinDynamics::setup( const std::vector<double>& masses, const Platform& platform ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookLangevinDynamics::setup";
// ---------------------------------------------------------------------------------------
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
setLog( brookPlatform.getLog() );
int numberOfParticles = (int) masses.size();
setNumberOfParticles( numberOfParticles );
// set stream sizes and then create streams
_initializeStreamSizes( numberOfParticles, platform );
_initializeStreams( platform );
_setInverseSqrtMasses( masses );
return DefaultReturnValue;
}
/*
* Get contents of object
* *
* @return ErrorReturnValue if error, else DefaultReturnValue * @param level level of dump
* *
*/ * @return string containing contents
*
* */
int BrookLangevinDynamics::_initializeStreams( const Platform& platform ){ std::string BrookLangevinDynamics::getContentsString( int level ) const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookLangevinDynamics::_initializeStreams"; static const std::string methodName = "BrookLangevinDynamics::getContentsString";
BrookOpenMMFloat dangleValue = (BrookOpenMMFloat) 0.0; static const unsigned int MAX_LINE_CHARS = 256;
char value[MAX_LINE_CHARS];
static const char* Set = "Set";
static const char* NotSet = "Not set";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize(); std::stringstream message;
int sdParticleStreamWidth = getLangevinDynamicsParticleStreamWidth(); std::string tab = " ";
_sdStreams[SDPC1Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC1Stream, #ifdef WIN32
sdParticleStreamSize, sdParticleStreamWidth, #define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );
BrookStreamInternal::Float2, dangleValue ); #else
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );
#endif
_sdStreams[SDPC2Stream] = new BrookFloatStreamInternal( BrookCommon::SDPC2Stream, (void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() );
sdParticleStreamSize, sdParticleStreamWidth, message << _getLine( tab, "Number of particles:", value );
BrookStreamInternal::Float2, dangleValue );
_sdStreams[SD2XStream] = new BrookFloatStreamInternal( BrookCommon::SD2XStream, (void) LOCAL_SPRINTF( value, "%d", getParticleStreamWidth() );
sdParticleStreamSize, sdParticleStreamWidth, message << _getLine( tab, "Particle stream width:", value );
BrookStreamInternal::Float3, dangleValue );
_sdStreams[SD1VStream] = new BrookFloatStreamInternal( BrookCommon::SD1VStream, (void) LOCAL_SPRINTF( value, "%d", getParticleStreamHeight() );
sdParticleStreamSize, sdParticleStreamWidth, message << _getLine( tab, "Particle stream height:", value );
BrookStreamInternal::Float3, dangleValue );
_sdStreams[VPrimeStream] = new BrookFloatStreamInternal( BrookCommon::VPrimeStream, (void) LOCAL_SPRINTF( value, "%d", getParticleStreamSize() );
sdParticleStreamSize, sdParticleStreamWidth, message << _getLine( tab, "Particle stream size:", value );
BrookStreamInternal::Float3, dangleValue );
_sdStreams[XPrimeStream] = new BrookFloatStreamInternal( BrookCommon::XPrimeStream, (void) LOCAL_SPRINTF( value, "%.5f", getTau() );
sdParticleStreamSize, sdParticleStreamWidth, message << _getLine( tab, "Tau:", value );
BrookStreamInternal::Float3, dangleValue );
_sdStreams[InverseMassStream] = new BrookFloatStreamInternal( BrookCommon::InverseMassStream, (void) LOCAL_SPRINTF( value, "%.5f", getTemperature() );
sdParticleStreamSize, sdParticleStreamWidth, message << _getLine( tab, "Temperature:", value );
BrookStreamInternal::Float, dangleValue );
return DefaultReturnValue; (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();
}
/** /**
* Update sd streams -- called after parameters change * Update
* *
* @return ErrorReturnValue if error, else DefaultReturnValue * @param positions particle positions
* @param velocities particle velocities
* @param forces particle forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
* @param brookRandomNumberGenerator BrookRandomNumberGenerator reference
*
* @return DefaultReturnValue
* *
*/ */
int BrookLangevinDynamics::_updateSdStreams( void ){ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamImpl& velocityStream,
BrookStreamImpl& forceStream,
BrookShakeAlgorithm& brookShakeAlgorithm,
BrookRandomNumberGenerator& brookRandomNumberGenerator ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookLangevinDynamics::_updateSdStreams"; static const std::string methodName = "\nBrookLangevinDynamics::update";
static const int PrintOn = 0;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize(); const BrookOpenMMFloat* derivedParameters = getDerivedParameters();
BrookOpenMMFloat* sdpc[2]; if( (0 || PrintOn) && getLog() ){
for( int ii = 0; ii < 2; ii++ ){
sdpc[ii] = new BrookOpenMMFloat[2*sdParticleStreamSize];
memset( sdpc[ii], 0, 2*sdParticleStreamSize*sizeof( BrookOpenMMFloat ) );
}
BrookOpenMMFloat* inverseMass = new BrookOpenMMFloat[sdParticleStreamSize];
memset( inverseMass, 0, sdParticleStreamSize*sizeof( BrookOpenMMFloat ) );
const BrookOpenMMFloat* derivedParameters = getDerivedParameters( ); static int showAux = 1;
int numberOfParticles = getNumberOfParticles();
int index = 0;
for( int ii = 0; ii < numberOfParticles; ii++, index += 2 ){
sdpc[0][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yv]) ); if( PrintOn ){
sdpc[0][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[V]) ); (void) fprintf( getLog(), "%s shake=%d\n", methodName.c_str(), brookShakeAlgorithm.getNumberOfConstraints() );
(void) fflush( getLog() );
}
sdpc[1][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yx]) ); // show update
sdpc[1][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[X]) );
if( showAux ){
showAux = 0;
inverseMass[ii] = _inverseSqrtMasses[ii]*_inverseSqrtMasses[ii]; std::string contents = brookRandomNumberGenerator.getContentsString( );
(void) fprintf( getLog(), "%s RNG contents\n%s", methodName.c_str(), contents.c_str() );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
contents = brookShakeAlgorithm.getContentsString( );
(void) fprintf( getLog(), "%s Shake contents\n%s", methodName.c_str(), contents.c_str() );
(void) fflush( getLog() );
}
} }
_sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] ); // first integration step
_sdStreams[SDPC2Stream]->loadFromArray( sdpc[1] );
_sdStreams[InverseMassStream]->loadFromArray( inverseMass );
for( int ii = 0; ii < 2; ii++ ){ kupdate_sd1_fix1(
delete[] sdpc[ii]; (float) getLangevinDynamicsParticleStreamWidth(),
} (float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
delete[] inverseMass; (float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM],
derivedParameters[Sd1pc1],
derivedParameters[Sd1pc2],
derivedParameters[Sd1pc3],
getSDPC1Stream()->getBrookStream(),
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD2XStream()->getBrookStream(),
positionStream.getBrookStream(),
forceStream.getBrookStream(),
velocityStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
getSD1VStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getXPrimeStream()->getBrookStream()
);
// diagnostics
// initialize SD2X if( PrintOn ){
(void) fprintf( getLog(), "\n%s Post kupdate_sd1_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]", methodName.c_str(),
getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] );
BrookOpenMMFloat* sd2x = new BrookOpenMMFloat[3*sdParticleStreamSize]; (void) fprintf( getLog(), "\nSDPC1Stream\n" );
//SimTKOpenMMUtilities::setRandomNumberSeed( (uint32_t) getRandomNumberSeed() ); getSDPC1Stream()->printToFile( getLog() );
memset( sd2x, 0, 3*sdParticleStreamSize*sizeof( BrookOpenMMFloat ) ); (void) fprintf( getLog(), "\nSD2XStream\n" );
getSD2XStream()->printToFile( getLog() );
index = 0; (void) fprintf( getLog(), "\nInverseMassStream\n" );
for( int ii = 0; ii < numberOfParticles; ii++, index += 3 ){ getInverseMassStream()->printToFile( getLog() );
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; //StreamImpl& positionStreamImpl = positionStream.getImpl();
//const BrookStreamImpl brookPositions = dynamic_cast<BrookStreamImpl&> (positionStreamImpl);
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
return DefaultReturnValue; (void) fprintf( getLog(), "\nForceStream\n" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalV->printToFile( getLog() );
(void) fprintf( getLog(), "\nSD1VStream\n" );
getSD1VStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nVPrimeStream\n" );
getVPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
}
// advance random number cursor
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfParticles() );
// first Shake step
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
kshakeh_fix1(
(float) brookShakeAlgorithm.getMaxIterations(),
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeTolerance(),
brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// first Shake gather
kshakeh_update1_fix1(
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() );
} if( ( 0 || PrintOn) && getLog() ){
/** (void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d",
* Set masses methodName.c_str(), getLangevinDynamicsParticleStreamWidth() );
*
* @param masses particle masses (void) fprintf( getLog(), "\nShakeInverseMapStream\n" );
* brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() );
*/
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
int BrookLangevinDynamics::_setInverseSqrtMasses( const std::vector<double>& masses ){ (void) fprintf( getLog(), "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() );
// --------------------------------------------------------------------------------------- (void) fprintf( getLog(), "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() );
//static const std::string methodName = "BrookLangevinDynamics::_setInverseSqrtMasses"; (void) fprintf( getLog(), "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() );
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0; (void) fprintf( getLog(), "\nShakeXCons3\n" );
BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0; brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() );
// --------------------------------------------------------------------------------------- }
// 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; // second integration step
}
/*
* Setup of LangevinDynamics parameters
*
* @param masses masses
* @param platform Brook platform
*
* @return nonzero value if error
*
* */
int BrookLangevinDynamics::setup( const std::vector<double>& masses, const Platform& platform ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookLangevinDynamics::setup";
// ---------------------------------------------------------------------------------------
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
setLog( brookPlatform.getLog() );
int numberOfParticles = (int) masses.size();
setNumberOfParticles( numberOfParticles );
// set stream sizes and then create streams
_initializeStreamSizes( numberOfParticles, platform ); kupdate_sd2_fix1(
_initializeStreams( platform ); (float) getLangevinDynamicsParticleStreamWidth(),
(float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
(float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1],
derivedParameters[Sd2pc2],
getSDPC2Stream()->getBrookStream(),
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD1VStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getSD2XStream()->getBrookStream(),
velocityStream.getBrookStream(),
getXPrimeStream()->getBrookStream()
);
_setInverseSqrtMasses( masses ); // diagnostics
return DefaultReturnValue; if( 0 || PrintOn ){
} (void) fprintf( getLog(), "\n%s Post kupdate_sd2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"Sd2pc[]=[%12.5e %12.5e]", methodName.c_str(),
getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
/* (void) fprintf( getLog(), "\nSDPC2Stream\n" );
* Get contents of object getSDPC2Stream()->printToFile( getLog() );
*
* @param level level of dump
*
* @return string containing contents
*
* */
std::string BrookLangevinDynamics::getContentsString( int level ) const { (void) fprintf( getLog(), "\nSD2XStream\n" );
getSD1VStream()->printToFile( getLog() );
// --------------------------------------------------------------------------------------- BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
static const std::string methodName = "BrookLangevinDynamics::getContentsString"; (void) fprintf( getLog(), "\nVPrimeStream\n" );
getVPrimeStream()->printToFile( getLog() );
static const unsigned int MAX_LINE_CHARS = 256; (void) fprintf( getLog(), "\nXPrimeStream\n" );
char value[MAX_LINE_CHARS]; getXPrimeStream()->printToFile( getLog() );
static const char* Set = "Set";
static const char* NotSet = "Not set";
// --------------------------------------------------------------------------------------- (void) fprintf( getLog(), "\ngetSD2XStream\n" );
getSD2XStream()->printToFile( getLog() );
std::stringstream message; BrookStreamInternal* brookStreamInternalVel = velocityStream.getBrookStreamImpl();
std::string tab = " "; (void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalVel->printToFile( getLog() );
#ifdef WIN32 (void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() );
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) ); // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
#else }
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );
#endif
(void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() ); // advance random number cursor
message << _getLine( tab, "Number of particles:", value );
(void) LOCAL_SPRINTF( value, "%d", getParticleStreamWidth() ); brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfParticles() );
message << _getLine( tab, "Particle stream width:", value );
(void) LOCAL_SPRINTF( value, "%d", getParticleStreamHeight() ); // second Shake step
message << _getLine( tab, "Particle stream height:", value );
(void) LOCAL_SPRINTF( value, "%d", getParticleStreamSize() ); if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
message << _getLine( tab, "Particle stream size:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getTau() ); kshakeh_fix1(
message << _getLine( tab, "Tau:", value ); (float) brookShakeAlgorithm.getMaxIterations(),
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeTolerance(),
brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// second Shake gather
kshakeh_update2_fix1(
(float) getLangevinDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
positionStream.getBrookStream() );
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
"Sd2pc[]=[%12.5e %12.5e]", methodName.c_str(),
getLangevinDynamicsParticleStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
brookStreamInternalPos = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
// brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
}
(void) LOCAL_SPRINTF( value, "%.5f", getTemperature() ); } else {
message << _getLine( tab, "Temperature:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getStepSize() ); // no constraints
message << _getLine( tab, "Step size:", value );
const BrookOpenMMFloat* derivedParameters = getDerivedParameters(); if( PrintOn ){
for( int ii = 0; ii < MaxDerivedParameters; ii++ ){ (void) fprintf( getLog(), "\n%s Pre ksetStr3 (no constraints)", methodName.c_str() );
(void) LOCAL_SPRINTF( value, "%.5e", derivedParameters[ii] ); // (void) fprintf( getLog(), "\nXPrimeStream\n" );
message << _getLine( tab, _getDerivedParametersString( ii ), value ); // getXPrimeStream()->printToFile( getLog() );
}
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
}
(void) LOCAL_SPRINTF( value, "%.5f", getTemperature() ); kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream(), positionStream.getBrookStream() );
message << _getLine( tab, "Temperature:", value ); //ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) ); // diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\n%s Post ksetStr3 (no constraints)", methodName.c_str() );
message << _getLine( tab, "SDPC1:", (getSDPC1Stream() ? Set : NotSet) ); (void) fprintf( getLog(), "\nXPrimeStream\n" );
message << _getLine( tab, "SDPC2:", (getSDPC2Stream() ? Set : NotSet) ); getXPrimeStream()->printToFile( getLog() );
message << _getLine( tab, "SD2X:", (getSD2XStream() ? Set : NotSet) );
message << _getLine( tab, "SD1V:", (getSD1VStream() ? Set : NotSet) ); BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
for( int ii = 0; ii < LastStreamIndex; ii++ ){ brookStreamInternalPos->printToFile( getLog() );
message << std::endl;
if( _sdStreams[ii] ){ // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() );
message << _sdStreams[ii]->getContentsString( ); }
}
} }
#undef LOCAL_SPRINTF return DefaultReturnValue;
return message.str();
} }
platforms/brook/src/BrookRemoveCMMotionKernel.cpp
View file @ adb4ec5a
...@@ -96,7 +96,9 @@ void BrookRemoveCMMotionKernel::initialize( const System& system, const CMMotion ...@@ -96,7 +96,9 @@ void BrookRemoveCMMotionKernel::initialize( const System& system, const CMMotion
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const std::string methodName = "BrookRemoveCMMotionKernel::initialize"; static const std::string methodName = "BrookRemoveCMMotionKernel::initialize";
static const int PrintOn = 0;
FILE* log = getLog();
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -110,6 +112,11 @@ void BrookRemoveCMMotionKernel::initialize( const System& system, const CMMotion ...@@ -110,6 +112,11 @@ void BrookRemoveCMMotionKernel::initialize( const System& system, const CMMotion
_brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval(); _brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval();
_brookVelocityCenterOfMassRemoval->setup( masses, getPlatform() ); _brookVelocityCenterOfMassRemoval->setup( masses, getPlatform() );
if( PrintOn && log ){
(void) fprintf( log, "%s\n", methodName.c_str() );
(void) fflush( log );
}
return; return;
} }
...@@ -162,7 +169,7 @@ void BrookRemoveCMMotionKernel::execute( OpenMMContextImpl& context ){ ...@@ -162,7 +169,7 @@ void BrookRemoveCMMotionKernel::execute( OpenMMContextImpl& context ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const std::string methodName = "BrookRemoveCMMotionKernel::execute"; static const std::string methodName = "BrookRemoveCMMotionKernel::execute";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
platforms/brook/src/BrookVerletDynamics.cpp
View file @ adb4ec5a
...@@ -140,7 +140,7 @@ int BrookVerletDynamics::updateParameters( double stepSize ){ ...@@ -140,7 +140,7 @@ int BrookVerletDynamics::updateParameters( double stepSize ){
static int showUpdate = 1; static int showUpdate = 1;
static int maxShowUpdate = 3; static int maxShowUpdate = 3;
static const std::string methodName = "\nBrookVerletDynamics::updateParameters"; static const std::string methodName = "BrookVerletDynamics::updateParameters";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -151,7 +151,7 @@ int BrookVerletDynamics::updateParameters( double stepSize ){ ...@@ -151,7 +151,7 @@ int BrookVerletDynamics::updateParameters( double stepSize ){
if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){ if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
std::string contents = getContentsString( ); std::string contents = getContentsString( );
(void) fprintf( getLog(), "%s contents\n%s", methodName, contents.c_str() ); (void) fprintf( getLog(), "%s contents\n%s", methodName.c_str(), contents.c_str() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -160,215 +160,6 @@ int BrookVerletDynamics::updateParameters( double stepSize ){ ...@@ -160,215 +160,6 @@ int BrookVerletDynamics::updateParameters( double stepSize ){
} }
/**
* Update
*
* @param positions particle positions
* @param velocities particle velocities
* @param forces particle forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
*
* @return DefaultReturnValue
*
*/
int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImpl& velocityStream,
const BrookStreamImpl& forceStreamC,
BrookShakeAlgorithm& brookShakeAlgorithm ){
// ---------------------------------------------------------------------------------------
static std::string methodName = "\nBrookVerletDynamics::update";
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
if( (1 || PrintOn) && getLog() ){
static int showAux = 1;
if( showAux ){
showAux = 0;
/*
std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
(void) fprintf( getLog(), "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
*/
(void) fprintf( getLog(), "%s Shake contents\n%s", methodName, brookShakeAlgorithm.getContentsString().c_str() );
(void) fflush( getLog() );
}
}
// To Shake or not to Shake
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
// integration step
kupdate_md1( (float) getStepSize(),
positionStream.getBrookStream(),
velocityStream.getBrookStream(),
forceStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
getXPrimeStream()->getBrookStream()
);
// diagnostics
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\n%s Post kupdate_md_verlet: particleStrW=%3d step=%.5f",
methodName.c_str(), getVerletDynamicsParticleStreamWidth(), getStepSize() );
(void) fprintf( getLog(), "\nInverseMassStream\n" );
getInverseMassStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStream\n" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalV->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
}
kshakeh_fix1(
(float) brookShakeAlgorithm.getMaxIterations(),
(float) getVerletDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeTolerance(),
brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
if( (1|| PrintOn) && getLog() ){
(void) fprintf( getLog(), "\n%s Post kshakeh_fix1: particleStrW=%3d", methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeParticleIndicesStream\n" );
brookShakeAlgorithm.getShakeParticleIndicesStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeParticleParameterStream\n" );
brookShakeAlgorithm.getShakeParticleParameterStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() );
}
// second Shake gather
kshakeh_update2_fix1(
(float) getVerletDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
positionStream.getBrookStream() );
if( ( 1 || PrintOn) && getLog() ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d",
methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeInverseMapStream\n" );
brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() );
}
// second integration step
float inverseStepSize = 1.0f/getStepSize();
kupdate_md2( inverseStepSize,
getXPrimeStream()->getBrookStream(),
positionStream.getBrookStream(),
velocityStream.getBrookStream(),
positionStream.getBrookStream() );
if( ( 1 || PrintOn) && getLog() ){
(void) fprintf( getLog(), "\n%s Post kupdate_md2: inverseStepSize=%3e",
methodName.c_str(), inverseStepSize );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
brookStreamInternalPos = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalPos->printToFile( getLog() );
}
} else {
kupdateMdNoShake( getStepSize(),
positionStream.getBrookStream(),
velocityStream.getBrookStream(),
forceStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
velocityStream.getBrookStream(),
positionStream.getBrookStream() );
}
//_brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
return DefaultReturnValue;
};
/** /**
* Get Particle stream size * Get Particle stream size
* *
...@@ -656,3 +447,213 @@ std::string BrookVerletDynamics::getContentsString( int level ) const { ...@@ -656,3 +447,213 @@ std::string BrookVerletDynamics::getContentsString( int level ) const {
return message.str(); return message.str();
} }
/**
* Update
*
* @param positions particle positions
* @param velocities particle velocities
* @param forces particle forces
* @param brookShakeAlgorithm BrookShakeAlgorithm reference
*
* @return DefaultReturnValue
*
*/
int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImpl& velocityStream,
const BrookStreamImpl& forceStreamC,
BrookShakeAlgorithm& brookShakeAlgorithm ){
// ---------------------------------------------------------------------------------------
static std::string methodName = "\nBrookVerletDynamics::update";
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
if( (1 || PrintOn) && getLog() ){
static int showAux = 1;
if( showAux ){
showAux = 0;
/*
std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
(void) fprintf( getLog(), "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
*/
(void) fprintf( getLog(), "%s Shake contents\n%s", methodName.c_str(), brookShakeAlgorithm.getContentsString().c_str() );
(void) fflush( getLog() );
}
}
// To Shake or not to Shake
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
// integration step
kupdate_md1( (float) getStepSize(),
positionStream.getBrookStream(),
velocityStream.getBrookStream(),
forceStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
getXPrimeStream()->getBrookStream()
);
// diagnostics
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\n%s Post kupdate_md1: particleStrW=%3d step=%.5f",
methodName.c_str(), getVerletDynamicsParticleStreamWidth(), getStepSize() );
(void) fprintf( getLog(), "\nInverseMassStream\n" );
getInverseMassStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStream\n" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalV = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalV->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
}
// Shake
kshakeh_fix1(
(float) brookShakeAlgorithm.getMaxIterations(),
(float) getVerletDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeTolerance(),
brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
if( (0|| PrintOn) && getLog() ){
(void) fprintf( getLog(), "\n%s Post kshakeh_fix1: particleStrW=%3d", methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeParticleIndicesStream\n" );
brookShakeAlgorithm.getShakeParticleIndicesStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeParticleParameterStream\n" );
brookShakeAlgorithm.getShakeParticleParameterStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() );
}
// Shake gather
kshakeh_update1_fix1(
(float) getVerletDynamicsParticleStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
getXPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() );
//positionStream.getBrookStream() );
if( ( 0 || PrintOn) && getLog() ){
(void) fprintf( getLog(), "\n%s Post kshakeh_update2_fix1: particleStrW=%3d",
methodName.c_str(), getVerletDynamicsParticleStreamWidth() );
(void) fprintf( getLog(), "\nShakeInverseMapStream\n" );
brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons0\n" );
brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons1\n" );
brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons2\n" );
brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nShakeXCons3\n" );
brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() );
}
// second integration step
float inverseStepSize = 1.0f/getStepSize();
kupdate_md2( inverseStepSize,
getXPrimeStream()->getBrookStream(),
positionStream.getBrookStream(),
velocityStream.getBrookStream(),
positionStream.getBrookStream() );
if( ( 0 || PrintOn) && getLog() ){
(void) fprintf( getLog(), "\n%s Post kupdate_md2: inverseStepSize=%3e",
methodName.c_str(), inverseStepSize );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
brookStreamInternalPos = velocityStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nVelocityStream\n" );
brookStreamInternalPos->printToFile( getLog() );
}
} else {
kupdateMdNoShake( getStepSize(),
positionStream.getBrookStream(),
velocityStream.getBrookStream(),
forceStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
velocityStream.getBrookStream(),
positionStream.getBrookStream() );
}
//_brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
return DefaultReturnValue;
}
platforms/brook/src/gpu/kcommon.br
View file @ adb4ec5a
...@@ -30,31 +30,29 @@ ...@@ -30,31 +30,29 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
//Inverse of above kernel void kgetxyz( float4 instr<>, out float3 outstr<> ){
kernel void kgetxyz( float4 instr<>, out float3 outstr<> ) {
outstr = instr.xyz; outstr = instr.xyz;
} }
//Zeroes out a stream //Zeroes out a stream
kernel void kzerof3( out float3 outstr<> ) { kernel void kzerof3( out float3 outstr<> ){
outstr = float3( 0.0f, 0.0f, 0.0f ); outstr = float3( 0.0f, 0.0f, 0.0f );
} }
//Zeros out a stream //Zeros out a stream
kernel void kzerof4( out float4 outstr<> ) { kernel void kzerof4( out float4 outstr<> ){
outstr = float4( 0.0f, 0.0f, 0.0f, 0.0f ); outstr = float4( 0.0f, 0.0f, 0.0f, 0.0f );
} }
kernel void ksetf4( float4 val, out float4 outstr<> ) { kernel void ksetf4( float4 val, out float4 outstr<> ){
outstr = val; outstr = val;
} }
kernel void ksetStr3( float3 instr<>, out float3 outstr<> ) { kernel void ksetStr3( float3 instr<>, out float3 outstr<> ){
outstr = instr; outstr = instr;
} }
kernel void kadd3( float3 val<>, out float3 outstr<> ) { kernel void kadd3( float3 in1<>, float3 in2<>, out float3 outstr<> ){
outstr += val; outstr = in1 + in2;
} }
platforms/brook/src/gpu/kcommon.h
View file @ adb4ec5a
...@@ -33,15 +33,14 @@ ...@@ -33,15 +33,14 @@
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
void kgetxyz (::brook::stream instr, void kgetxyz(::brook::stream instr, ::brook::stream outstr);
::brook::stream outstr);
void kzerof3 (::brook::stream outstr); void kzerof3(::brook::stream outstr);
void kzerof4 (::brook::stream outstr); void kzerof4(::brook::stream outstr);
void kzerof4 (::brook::stream outstr); void kzerof4(::brook::stream outstr);
void ksetf4 (const float4 val, ::brook::stream outstr); void ksetf4(const float4 val, ::brook::stream outstr);
void kadd3( ::brook::stream instr, ::brook::stream outstr ); void kadd3( ::brook::stream instr1, ::brook::stream instr2, ::brook::stream outstr );
void ksetStr3( ::brook::stream instr, ::brook::stream outstr ); void ksetStr3( ::brook::stream instr, ::brook::stream outstr );
#endif // __KCOMMON_H__ #endif // __KCOMMON_H__
platforms/brook/src/gpu/kshakeh.br
View file @ adb4ec5a
...@@ -145,7 +145,7 @@ kshakeh_fix1( ...@@ -145,7 +145,7 @@ kshakeh_fix1(
acor = ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ; acor = ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ;
diff = abs( ld1 - 2 * rrpr - rpsqij ) / (params.z * tolerance ); diff = abs( ld1 - 2 * rrpr - rpsqij ) / (params.z * tolerance );
acor = (diff < 1.0f) ? 0.0f : acor; acor = (diff < 1.0f) ? 0.0f : acor;
converged = acor; converged = abs( acor );
dr = rij1 * acor; dr = rij1 * acor;
xpi += dr * params.x; xpi += dr * params.x;
...@@ -160,7 +160,7 @@ kshakeh_fix1( ...@@ -160,7 +160,7 @@ kshakeh_fix1(
diff = abs( ld2 - 2.0f * rrpr - rpsqij ) / (params.z * tolerance ); diff = abs( ld2 - 2.0f * rrpr - rpsqij ) / (params.z * tolerance );
acor = mask2 * ( ld2 - 2.0f * rrpr - rpsqij ) * params.y / ( rrpr + rij2sq ) ; acor = mask2 * ( ld2 - 2.0f * rrpr - rpsqij ) * params.y / ( rrpr + rij2sq ) ;
acor = (diff < 1.0f) ? 0.0f : acor; acor = (diff < 1.0f) ? 0.0f : acor;
converged += acor; converged += abs( acor );
dr = rij2 * acor; dr = rij2 * acor;
xpi += dr * params.x; xpi += dr * params.x;
...@@ -174,7 +174,7 @@ kshakeh_fix1( ...@@ -174,7 +174,7 @@ kshakeh_fix1(
diff = abs( ld3 - 2.0f * rrpr - rpsqij ) / (params.z * tolerance ); diff = abs( ld3 - 2.0f * rrpr - rpsqij ) / (params.z * tolerance );
acor = mask3 * ( ld3 - 2.0f * rrpr - rpsqij ) * params.y / ( rrpr + rij3sq ) ; acor = mask3 * ( ld3 - 2.0f * rrpr - rpsqij ) * params.y / ( rrpr + rij3sq ) ;
acor = (diff < 1.0f) ? 0.0f : acor; acor = (diff < 1.0f) ? 0.0f : acor;
converged += acor; converged += abs( acor );
dr = rij3 * acor; dr = rij3 * acor;
xpi += dr * params.x; xpi += dr * params.x;
...@@ -195,7 +195,6 @@ kshakeh_fix1( ...@@ -195,7 +195,6 @@ kshakeh_fix1(
kernel void kshakeh_update1_fix1( kernel void kshakeh_update1_fix1(
float strwidth, //width of cposq streams float strwidth, //width of cposq streams
float2 invmap<>, //shakeh inverse map float2 invmap<>, //shakeh inverse map
float3 posq<>, //old positions
float3 posqp<>, //deltas from sd2 float3 posqp<>, //deltas from sd2
float3 cposq0[][], //constrained delta for heavy atom float3 cposq0[][], //constrained delta for heavy atom
float3 cposq1[][], //ditto for h1 float3 cposq1[][], //ditto for h1
......
platforms/brook/src/gpu/kshakeh.h
View file @ adb4ec5a
...@@ -45,7 +45,6 @@ void kshakeh_fix1 ( ...@@ -45,7 +45,6 @@ void kshakeh_fix1 (
void kshakeh_update1_fix1 ( void kshakeh_update1_fix1 (
const float strwidth, const float strwidth,
::brook::stream invmap, ::brook::stream invmap,
::brook::stream posq,
::brook::stream posqp, ::brook::stream posqp,
::brook::stream cposq0, ::brook::stream cposq0,
::brook::stream cposq1, ::brook::stream cposq1,
......
platforms/brook/src/gpu/kupdatemd.br
View file @ adb4ec5a
...@@ -48,8 +48,8 @@ kernel void kupdate_md2( ...@@ -48,8 +48,8 @@ kernel void kupdate_md2(
out float3 vnew<>, //Corrected velocities out float3 vnew<>, //Corrected velocities
out float3 posqnew<> //equal to posqp, avoids an extra call to copy out float3 posqnew<> //equal to posqp, avoids an extra call to copy
){ ){
posqnew = posq + posqp;
vnew = posqp * dtinv; vnew = posqp * dtinv;
posqnew = posq + posqp;
} }
kernel void kupdateMdNoShake( kernel void kupdateMdNoShake(
...@@ -64,4 +64,3 @@ kernel void kupdateMdNoShake( ...@@ -64,4 +64,3 @@ kernel void kupdateMdNoShake(
outv = v + dt*invmass*f; outv = v + dt*invmass*f;
posqp += dt*outv; posqp += dt*outv;
} }
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