Skip to content

GitLab

Commit 75acdad0 authored by Peter Eastman's avatar Peter Eastman
Browse files

Converted reference implementation of Langevin integrator to use a different algorithm

parent 7f486d6b
Hide whitespace changes
Inline Side-by-side
Showing with 51 additions and 729 deletions
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp
View file @ 75acdad0
......@@ -69,7 +69,6 @@ ReferenceStochasticDynamics::ReferenceStochasticDynamics( int numberOfAtoms,
_tau = one;
}
_setFixedParameters( );
allocate2DArrays( numberOfAtoms, 3, Max2DArrays );
allocate1DArrays( numberOfAtoms, Max1DArrays );
......@@ -92,76 +91,6 @@ ReferenceStochasticDynamics::~ReferenceStochasticDynamics( ){
}
/**---------------------------------------------------------------------------------------
Set fixed parameters
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceStochasticDynamics::_setFixedParameters( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceStochasticDynamics::_setFixedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM four = 4.0;
static const RealOpenMM half = 0.5;
// ---------------------------------------------------------------------------------------
_fixedParameters[GDT] = getDeltaT()/getTau();
_fixedParameters[EPH] = EXP( half*_fixedParameters[GDT] );
_fixedParameters[EMH] = EXP( -half*_fixedParameters[GDT] );
_fixedParameters[EM] = EXP( -_fixedParameters[GDT] );
_fixedParameters[EP] = EXP( _fixedParameters[GDT] );
if( _fixedParameters[GDT] >= (RealOpenMM) 0.1 ){
RealOpenMM term1 = _fixedParameters[EPH] - one;
term1 *= term1;
_fixedParameters[B] = _fixedParameters[GDT]*(_fixedParameters[EP] - one) - four*term1;
_fixedParameters[C] = _fixedParameters[GDT] - three + four*_fixedParameters[EMH] - _fixedParameters[EM];
_fixedParameters[D] = two - _fixedParameters[EPH] - _fixedParameters[EMH];
} else {
// this has not been debugged
RealOpenMM term1 = half*_fixedParameters[GDT];
RealOpenMM term2 = term1*term1;
RealOpenMM term4 = term2*term2;
RealOpenMM third = (RealOpenMM) ( 1.0/3.0 );
RealOpenMM o7_9 = (RealOpenMM) ( 7.0/9.0 );
RealOpenMM o1_12 = (RealOpenMM) ( 1.0/12.0 );
RealOpenMM o17_90 = (RealOpenMM) ( 17.0/90.0 );
RealOpenMM o7_30 = (RealOpenMM) ( 7.0/30.0 );
RealOpenMM o31_1260 = (RealOpenMM) ( 31.0/1260.0 );
RealOpenMM o_360 = (RealOpenMM) ( 1.0/360.0 );
_fixedParameters[B] = term4*( third + term1*( third + term1*( o17_90 + term1*o7_9 )));
_fixedParameters[C] = term2*term1*( two*third + term1*( -half + term1*( o7_30 + term1*(-o1_12 + term1*o31_1260 ))));
_fixedParameters[D] = term2*( -one + term2*(-o1_12 - term2*o_360));
}
RealOpenMM kT = ((RealOpenMM) BOLTZ)*getTemperature();
_fixedParameters[V] = SQRT( kT*( one - _fixedParameters[EM]) );
_fixedParameters[X] = getTau()*SQRT( kT*_fixedParameters[C] );
_fixedParameters[Yv] = SQRT( kT*_fixedParameters[B]/_fixedParameters[C] );
_fixedParameters[Yx] = getTau()*SQRT( kT*_fixedParameters[B]/(one - _fixedParameters[EM]) );
return ReferenceDynamics::DefaultReturn;
};
/**---------------------------------------------------------------------------------------
Get tau
......@@ -181,25 +110,6 @@ RealOpenMM ReferenceStochasticDynamics::getTau( void ) const {
return _tau;
}
/**---------------------------------------------------------------------------------------
Get array of fixed parameters indexed by 'FixedParameters' enums
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* ReferenceStochasticDynamics::getFixedParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceStochasticDynamics::getFixedParameters";
// ---------------------------------------------------------------------------------------
return _fixedParameters;
}
/**---------------------------------------------------------------------------------------
Print parameters
......@@ -214,9 +124,7 @@ int ReferenceStochasticDynamics::printParameters( std::stringstream& message ) c
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceStochasticDynamics::printParameters";
static const char* parameterNames[MaxFixedParameters] = { "gdt", "ep", "eph", "emh", "em", "B", "C", "D",
"V", "X", "Yv", "Yx" };
//static const char* methodName = "\nReferenceStochasticDynamics::printParameters";
// ---------------------------------------------------------------------------------------
......@@ -225,14 +133,6 @@ int ReferenceStochasticDynamics::printParameters( std::stringstream& message ) c
ReferenceDynamics::printParameters( message );
message << " tau=" << getTau();
message << " T=" << getTemperature();
int cut = 3;
for( int ii = 0; ii < MaxFixedParameters; ii++ ){
message << " " << parameterNames[ii] << "=" << _fixedParameters[ii];
if( cut++ > 5 ){
cut = 0;
message << std::endl;
}
}
return ReferenceDynamics::DefaultReturn;
......@@ -248,9 +148,6 @@ int ReferenceStochasticDynamics::printParameters( std::stringstream& message ) c
@param forces forces
@param inverseMasses inverse atom masses
@param xPrime xPrime
@param oldVelocities previous velocities
@param xVector xVector
@param vVector vVector
@return ReferenceDynamics::DefaultReturn
......@@ -259,76 +156,29 @@ int ReferenceStochasticDynamics::printParameters( std::stringstream& message ) c
int ReferenceStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector ){
RealOpenMM** xPrime ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceStochasticDynamics::updatePart1";
static const RealOpenMM one = 1.0;
static int debug = 0;
//static const char* methodName = "\nReferenceStochasticDynamics::updatePart1";
// ---------------------------------------------------------------------------------------
// perform first update
const RealOpenMM* fixedParameters = getFixedParameters();
RealOpenMM tau = getTau();
RealOpenMM fix1 = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
for( int ii = 0; ii < numberOfAtoms; ii++ ){
RealOpenMM tau = getTau();
const RealOpenMM vscale = EXP(-getDeltaT()/tau);
const RealOpenMM fscale = (1-vscale)*tau;
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
for( int jj = 0; jj < 3; jj++ ){
oldVelocities[ii][jj] = velocities[ii][jj];
RealOpenMM Vmh = xVector[ii][jj]*fixedParameters[D]/(tau*fixedParameters[C]) +
sqrtInvMass*fixedParameters[Yv]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
vVector[ii][jj] = sqrtInvMass*fixedParameters[V]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
RealOpenMM vPrime = oldVelocities[ii][jj]*fixedParameters[EM] +
inverseMasses[ii]*forces[ii][jj]*tau*( one - fixedParameters[EM]) +
vVector[ii][jj] - fixedParameters[EM]*Vmh;
xPrime[ii][jj] = atomCoordinates[ii][jj] + vPrime*fix1;
for (int ii = 0; ii < numberOfAtoms; ii++) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
// diagnostics
if( debug ){
int maxPrint = 5;
std::stringstream message;
message << methodName << " Post SD1 atoms=" << numberOfAtoms << "\n";
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
for( int ii = 0; ii < maxPrint; ii++ ){
message << " mI=" << inverseMasses[ii];
RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
message << " sdpc[" << sqrtInvMass*fixedParameters[Yv] << " " << sqrtInvMass*fixedParameters[V];
message << " x[";
SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
message << "] xp[";
SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
message << "] v[";
SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
message << "] vV[";
SimTKOpenMMUtilities::formatRealStringStream( message, vVector[ii], 3, one );
message << "] ov[";
SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
message << "] xV[";
message << "] f[";
SimTKOpenMMUtilities::formatRealStringStream( message, forces[ii], 3, one );
message << "] xV[";
SimTKOpenMMUtilities::formatRealStringStream( message, xVector[ii], 3, one );
message << "]\n";
}
SimTKOpenMMLog::printMessage( message );
}
return ReferenceDynamics::DefaultReturn;
}
......@@ -350,68 +200,20 @@ int ReferenceStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** at
int ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector ){
RealOpenMM** xPrime ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceStochasticDynamics::updatePart2";
static const RealOpenMM one = 1.0;
static int debug = 0;
//static const char* methodName = "\nReferenceStochasticDynamics::updatePart2";
// ---------------------------------------------------------------------------------------
// perform second update
const RealOpenMM* fixedParameters = getFixedParameters();
RealOpenMM tau = getTau();
RealOpenMM fix1 = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
fix1 = one/fix1;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
for( int jj = 0; jj < 3; jj++ ){
velocities[ii][jj] = (xPrime[ii][jj] - atomCoordinates[ii][jj])*fix1;
RealOpenMM Xmh = vVector[ii][jj]*tau*fixedParameters[D]/(fixedParameters[EM] - one) +
sqrtInvMass*fixedParameters[Yx]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
xVector[ii][jj] = sqrtInvMass*fixedParameters[X]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
xPrime[ii][jj] += xVector[ii][jj] - Xmh;
}
}
// diagnostics
if( debug ){
int maxPrint = 5;
std::stringstream message;
message << methodName << " Post SD2 atoms=" << numberOfAtoms << "\n";
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
for( int ii = 0; ii < maxPrint; ii++ ){
message << " mI=" << inverseMasses[ii];
RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
message << " sdpc[" << sqrtInvMass*fixedParameters[Yx] << " " << sqrtInvMass*fixedParameters[X];
message << " x[";
message << " x[";
SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
message << "] xp[";
SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
message << "] v[";
SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
message << "] vV[";
SimTKOpenMMUtilities::formatRealStringStream( message, vVector[ii], 3, one );
message << "] oV[";
SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
message << "] xV[";
SimTKOpenMMUtilities::formatRealStringStream( message, xVector[ii], 3, one );
message << "]\n";
for (int ii = 0; ii < numberOfAtoms; ii++) {
for (int jj = 0; jj < 3; jj++) {
xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
}
SimTKOpenMMLog::printMessage( message );
}
return ReferenceDynamics::DefaultReturn;
......@@ -434,8 +236,7 @@ int ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, RealOpenMM** at
--------------------------------------------------------------------------------------- */
int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses ){
RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses ){
// ---------------------------------------------------------------------------------------
......@@ -444,17 +245,11 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static int debug = 0;
// ---------------------------------------------------------------------------------------
// get work arrays
RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D );
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
RealOpenMM** xVector = get2DArrayAtIndex( X2D );
RealOpenMM** vVector = get2DArrayAtIndex( V2D );
RealOpenMM* inverseMasses = get1DArrayAtIndex( InverseMasses );
// first-time-through initialization
......@@ -476,16 +271,6 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
}
}
// set xVector
const RealOpenMM* fixedParameters = getFixedParameters();
for( int ii = 0; ii < numberOfAtoms; ii++ ){
RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
for( int jj = 0; jj < 3; jj++ ){
xVector[ii][jj] = sqrtInverseMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
// exit if errors
if( errors ){
......@@ -495,80 +280,15 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
// 1st update
updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses,
xPrime, oldVelocities, xVector, vVector );
//writeState( numberOfAtoms, atomCoordinates, velocities, forces, masses, -1 , "Sd1" );
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
/*
std::stringstream message;
message << methodName;
message << " calling constrain1\n";
SimTKOpenMMLog::printMessage( message );
*/
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
inverseMasses );
}
updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime );
// 2nd update
if( debug ){
int maxPrint = 5;
std::stringstream message;
message << methodName << " Pre SD2 atoms=" << numberOfAtoms << std::endl;
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
for( int ii = 0; ii < maxPrint; ii++ ){
message << " x[";
SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
message << "] xp[";
SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
message << "] v[";
SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
message << "] ov[";
SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
message << "]\n";
}
SimTKOpenMMLog::printMessage( message );
}
updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses,
xPrime, oldVelocities, xVector, vVector );
if( debug ){
int maxPrint = 5;
std::stringstream message;
message << methodName << " Post SD2 atoms=" << numberOfAtoms << "\n";
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
for( int ii = 0; ii < maxPrint; ii++ ){
message << " x[";
SimTKOpenMMUtilities::formatRealStringStream( message, atomCoordinates[ii], 3, one );
message << "] xp[";
SimTKOpenMMUtilities::formatRealStringStream( message, xPrime[ii], 3, one );
message << "] v[";
SimTKOpenMMUtilities::formatRealStringStream( message, velocities[ii], 3, one );
message << "] ov[";
SimTKOpenMMUtilities::formatRealStringStream( message, oldVelocities[ii], 3, one );
message << "]\n";
}
SimTKOpenMMLog::printMessage( message );
}
updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime );
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
/*
std::stringstream message;
message << methodName;
message << " calling constrain2\n";
SimTKOpenMMLog::printMessage( message );
*/
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
inverseMasses );
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
}
// copy xPrime -> atomCoordinates
......@@ -584,150 +304,3 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
return ReferenceDynamics::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Write state
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param state 0 if initial state; otherwise nonzero
@param baseFileName base file name
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceStochasticDynamics::writeState( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses,
int state, const std::string& baseFileName ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceStochasticDynamics::writeState";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const int threeI = 3;
// ---------------------------------------------------------------------------------------
std::stringstream stateFileName;
stateFileName << baseFileName;
stateFileName << "_Step" << getTimeStep();
// stateFileName << "_State" << state;
stateFileName << ".txt";
// ---------------------------------------------------------------------------------------
// open file -- return if unsuccessful
FILE* stateFile = NULL;
#ifdef _MSC_VER
fopen_s( &stateFile, stateFileName.str().c_str(), "w" );
#else
stateFile = fopen( stateFileName.str().c_str(), "w" );
#endif
// ---------------------------------------------------------------------------------------
// diagnostics
if( stateFile != NULL ){
std::stringstream message;
message << methodName;
message << " Opened file=<" << stateFileName.str() << ">.\n";
SimTKOpenMMLog::printMessage( message );
} else {
std::stringstream message;
message << methodName;
message << " could not open file=<" << stateFileName.str() << "> -- abort output.\n";
SimTKOpenMMLog::printMessage( message );
return ReferenceDynamics::ErrorReturn;
}
// ---------------------------------------------------------------------------------------
StringVector scalarNameI;
IntVector scalarI;
StringVector scalarNameR;
RealOpenMMVector scalarR;
StringVector scalarNameR1;
RealOpenMMPtrVector scalarR1;
StringVector scalarNameR2;
RealOpenMMPtrPtrVector scalarR2;
scalarI.push_back( getNumberOfAtoms() );
scalarNameI.push_back( "Atoms" );
scalarI.push_back( getTimeStep() );
scalarNameI.push_back( "Timestep" );
if( state == 0 || state == -1 ){
scalarR.push_back( getDeltaT() );
scalarNameR.push_back( "delta_t" );
scalarR.push_back( getTemperature() );
scalarNameR.push_back( "T" );
scalarR.push_back( getTau() );
scalarNameR.push_back( "Tau" );
scalarR1.push_back( masses );
scalarNameR1.push_back( "mass" );
scalarR2.push_back( atomCoordinates );
scalarNameR2.push_back( "coord" );
scalarR2.push_back( velocities );
scalarNameR2.push_back( "velocities" );
scalarR2.push_back( forces );
scalarNameR2.push_back( "forces" );
if( state == -1 ){
RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D );
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
RealOpenMM** xVector = get2DArrayAtIndex( X2D );
RealOpenMM** vVector = get2DArrayAtIndex( V2D );
scalarR2.push_back( xPrime );
scalarNameR2.push_back( "xPrime" );
scalarR2.push_back( oldVelocities);
scalarNameR2.push_back( "vold" );
scalarR2.push_back( xVector );
scalarNameR2.push_back( "xVector" );
scalarR2.push_back( vVector );
scalarNameR2.push_back( "vVector" );
}
} else {
scalarR2.push_back( atomCoordinates );
scalarNameR2.push_back( "coord" );
scalarR2.push_back( velocities );
scalarNameR2.push_back( "velocities" );
}
writeStateToFile( stateFile, scalarNameI, scalarI, scalarNameR, scalarR, getNumberOfAtoms(), scalarNameR1, scalarR1, threeI, scalarNameR2, scalarR2 );
(void) fclose( stateFile );
return ReferenceDynamics::DefaultReturn;
}
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.h
View file @ 75acdad0
......@@ -33,22 +33,10 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
private:
enum FixedParameters { GDT, EPH, EMH, EP, EM, B, C, D, V, X, Yv, Yx, MaxFixedParameters };
enum TwoDArrayIndicies { X2D, V2D, OldV, xPrime2D, vPrime2D, Max2DArrays };
enum TwoDArrayIndicies { xPrime2D, Max2DArrays };
enum OneDArrayIndicies { InverseMasses, Max1DArrays };
RealOpenMM _tau;
RealOpenMM _fixedParameters[MaxFixedParameters];
/**---------------------------------------------------------------------------------------
Set fixed values
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int _setFixedParameters( void );
public:
......@@ -83,16 +71,6 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
RealOpenMM getTau( void ) const;
/**---------------------------------------------------------------------------------------
Get array of fixed parameters indexed by 'FixedParameters' enums
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* getFixedParameters( void ) const;
/**---------------------------------------------------------------------------------------
Print parameters
......@@ -132,18 +110,13 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
@param forces forces
@param inverseMasses inverse atom masses
@param xPrime xPrime
@param oldVelocities previous velocities
@param xVector xVector
@param vVector vVector
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector );
RealOpenMM** forces, RealOpenMM* inverseMasses, RealOpenMM** xPrime );
/**---------------------------------------------------------------------------------------
......@@ -160,29 +133,7 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
int updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector );
/**---------------------------------------------------------------------------------------
Write state
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param state 0 if initial state; otherwise nonzero
@param baseFileName base file name
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int writeState( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses,
int state, const std::string& baseFileName ) const;
RealOpenMM** forces, RealOpenMM* inverseMasses, RealOpenMM** xPrime );
};
......
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp
View file @ 75acdad0
......@@ -114,77 +114,6 @@ void ReferenceVariableStochasticDynamics::setAccuracy( RealOpenMM accuracy ) {
_accuracy = accuracy;
}
/**---------------------------------------------------------------------------------------
Set fixed parameters
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceVariableStochasticDynamics::_setFixedParameters( RealOpenMM timeStep, RealOpenMM prevTimeStep ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceVariableStochasticDynamics::_setFixedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM four = 4.0;
static const RealOpenMM half = 0.5;
// ---------------------------------------------------------------------------------------
_fixedParameters[GDT] = timeStep/getTau();
_fixedParameters[EPH] = EXP( half*_fixedParameters[GDT] );
_fixedParameters[EMH] = EXP( -half*_fixedParameters[GDT] );
_fixedParameters[EM] = EXP( -_fixedParameters[GDT] );
_fixedParameters[EM_V] = EXP( -half*(timeStep+prevTimeStep)/getTau() );
_fixedParameters[EP] = EXP( _fixedParameters[GDT] );
if( _fixedParameters[GDT] >= (RealOpenMM) 0.1 ){
RealOpenMM term1 = _fixedParameters[EPH] - one;
term1 *= term1;
_fixedParameters[B] = _fixedParameters[GDT]*(_fixedParameters[EP] - one) - four*term1;
_fixedParameters[C] = _fixedParameters[GDT] - three + four*_fixedParameters[EMH] - _fixedParameters[EM];
_fixedParameters[D] = two - _fixedParameters[EPH] - _fixedParameters[EMH];
} else {
// this has not been debugged
RealOpenMM term1 = half*_fixedParameters[GDT];
RealOpenMM term2 = term1*term1;
RealOpenMM term4 = term2*term2;
RealOpenMM third = (RealOpenMM) ( 1.0/3.0 );
RealOpenMM o7_9 = (RealOpenMM) ( 7.0/9.0 );
RealOpenMM o1_12 = (RealOpenMM) ( 1.0/12.0 );
RealOpenMM o17_90 = (RealOpenMM) ( 17.0/90.0 );
RealOpenMM o7_30 = (RealOpenMM) ( 7.0/30.0 );
RealOpenMM o31_1260 = (RealOpenMM) ( 31.0/1260.0 );
RealOpenMM o_360 = (RealOpenMM) ( 1.0/360.0 );
_fixedParameters[B] = term4*( third + term1*( third + term1*( o17_90 + term1*o7_9 )));
_fixedParameters[C] = term2*term1*( two*third + term1*( -half + term1*( o7_30 + term1*(-o1_12 + term1*o31_1260 ))));
_fixedParameters[D] = term2*( -one + term2*(-o1_12 - term2*o_360));
}
RealOpenMM kT = ((RealOpenMM) BOLTZ)*getTemperature();
_fixedParameters[V] = SQRT( kT*( one - _fixedParameters[EM]) );
_fixedParameters[X] = getTau()*SQRT( kT*_fixedParameters[C] );
_fixedParameters[Yv] = SQRT( kT*_fixedParameters[B]/_fixedParameters[C] );
_fixedParameters[Yx] = getTau()*SQRT( kT*_fixedParameters[B]/(one - _fixedParameters[EM]) );
return ReferenceDynamics::DefaultReturn;
};
/**---------------------------------------------------------------------------------------
Get tau
......@@ -204,25 +133,6 @@ RealOpenMM ReferenceVariableStochasticDynamics::getTau( void ) const {
return _tau;
}
/**---------------------------------------------------------------------------------------
Get array of fixed parameters indexed by 'FixedParameters' enums
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* ReferenceVariableStochasticDynamics::getFixedParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceVariableStochasticDynamics::getFixedParameters";
// ---------------------------------------------------------------------------------------
return _fixedParameters;
}
/**---------------------------------------------------------------------------------------
Print parameters
......@@ -237,9 +147,7 @@ int ReferenceVariableStochasticDynamics::printParameters( std::stringstream& mes
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableStochasticDynamics::printParameters";
static const char* parameterNames[MaxFixedParameters] = { "gdt", "ep", "eph", "emh", "em", "B", "C", "D",
"V", "X", "Yv", "Yx" };
//static const char* methodName = "\nReferenceVariableStochasticDynamics::printParameters";
// ---------------------------------------------------------------------------------------
......@@ -248,14 +156,6 @@ int ReferenceVariableStochasticDynamics::printParameters( std::stringstream& mes
ReferenceDynamics::printParameters( message );
message << " tau=" << getTau();
message << " T=" << getTemperature();
int cut = 3;
for( int ii = 0; ii < MaxFixedParameters; ii++ ){
message << " " << parameterNames[ii] << "=" << _fixedParameters[ii];
if( cut++ > 5 ){
cut = 0;
message << std::endl;
}
}
return ReferenceDynamics::DefaultReturn;
......@@ -272,9 +172,6 @@ int ReferenceVariableStochasticDynamics::printParameters( std::stringstream& mes
@param masses atom masses
@param inverseMasses inverse atom masses
@param xPrime xPrime
@param oldVelocities previous velocities
@param xVector xVector
@param vVector vVector
@param maxStepSize maximum time step
@return ReferenceDynamics::DefaultReturn
......@@ -284,18 +181,12 @@ int ReferenceVariableStochasticDynamics::printParameters( std::stringstream& mes
int ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector,
RealOpenMM maxStepSize ){
RealOpenMM** xPrime, RealOpenMM maxStepSize ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableStochasticDynamics::updatePart1";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static int debug = 0;
// ---------------------------------------------------------------------------------------
......@@ -310,11 +201,11 @@ int ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, RealOpe
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= zero ){
if( masses[ii] <= 0 ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = one/masses[ii];
inverseMasses[ii] = 1/masses[ii];
}
}
......@@ -326,7 +217,7 @@ int ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, RealOpe
}
// Select the step size to use
RealOpenMM error = zero;
RealOpenMM error = 0;
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
RealOpenMM xerror = inverseMasses[i]*forces[i][j];
......@@ -341,48 +232,20 @@ int ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, RealOpe
newStepSize = getDeltaT(); // Keeping dt constant between steps improves the behavior of the integrator.
if (newStepSize > maxStepSize)
newStepSize = maxStepSize;
_setFixedParameters(newStepSize, getDeltaT());
setDeltaT(newStepSize);
if( getTimeStep() == 0 ){
// Initialize xVector
const RealOpenMM* fixedParameters = getFixedParameters();
for( int ii = 0; ii < numberOfAtoms; ii++ ){
RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
for( int jj = 0; jj < 3; jj++ ){
xVector[ii][jj] = sqrtInverseMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
}
// perform first update
const RealOpenMM* fixedParameters = getFixedParameters();
RealOpenMM tau = getTau();
RealOpenMM fix1 = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
if (fix1 == zero)
fix1 = getDeltaT();
for( int ii = 0; ii < numberOfAtoms; ii++ ){
RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
for( int jj = 0; jj < 3; jj++ ){
oldVelocities[ii][jj] = velocities[ii][jj];
RealOpenMM Vmh = xVector[ii][jj]*fixedParameters[D]/(tau*fixedParameters[C]) +
sqrtInvMass*fixedParameters[Yv]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
vVector[ii][jj] = sqrtInvMass*fixedParameters[V]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
RealOpenMM tau = getTau();
const RealOpenMM vscale = EXP(-getDeltaT()/tau);
const RealOpenMM fscale = (1-vscale)*tau;
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
RealOpenMM vPrime = oldVelocities[ii][jj]*fixedParameters[EM_V] +
inverseMasses[ii]*forces[ii][jj]*tau*(one - fixedParameters[EM_V]) +
vVector[ii][jj] - fixedParameters[EM]*Vmh;
xPrime[ii][jj] = atomCoordinates[ii][jj] + vPrime*fix1;
for (int ii = 0; ii < numberOfAtoms; ii++) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
......@@ -407,41 +270,19 @@ int ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, RealOpe
int ReferenceVariableStochasticDynamics::updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates,
RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector ){
RealOpenMM** xPrime ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableStochasticDynamics::updatePart2";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static int debug = 0;
//static const char* methodName = "\nReferenceVariableStochasticDynamics::updatePart2";
// ---------------------------------------------------------------------------------------
// perform second update
const RealOpenMM* fixedParameters = getFixedParameters();
RealOpenMM tau = getTau();
RealOpenMM fix1 = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
if (fix1 == zero)
fix1 = getDeltaT();
fix1 = one/fix1;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
for( int jj = 0; jj < 3; jj++ ){
velocities[ii][jj] = (xPrime[ii][jj] - atomCoordinates[ii][jj])*fix1;
RealOpenMM Xmh = vVector[ii][jj]*tau*fixedParameters[D]/(fixedParameters[EM] - one) +
sqrtInvMass*fixedParameters[Yx]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
xVector[ii][jj] = sqrtInvMass*fixedParameters[X]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
xPrime[ii][jj] += xVector[ii][jj] - Xmh;
for (int ii = 0; ii < numberOfAtoms; ii++) {
for (int jj = 0; jj < 3; jj++) {
xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
}
}
......@@ -470,40 +311,24 @@ int ReferenceVariableStochasticDynamics::update( int numberOfAtoms, RealOpenMM**
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableStochasticDynamics::update";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static int debug = 0;
//static const char* methodName = "\nReferenceVariableStochasticDynamics::update";
// ---------------------------------------------------------------------------------------
// get work arrays
RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D );
RealOpenMM** oldVelocities = get2DArrayAtIndex( OldV );
RealOpenMM** xVector = get2DArrayAtIndex( X2D );
RealOpenMM** vVector = get2DArrayAtIndex( V2D );
RealOpenMM* inverseMasses = get1DArrayAtIndex( InverseMasses );
// 1st update
updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, masses, inverseMasses,
xPrime, oldVelocities, xVector, vVector, maxStepSize );
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
inverseMasses );
}
updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, masses, inverseMasses, xPrime, maxStepSize );
// 2nd update
updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses,
xPrime, oldVelocities, xVector, vVector );
updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime );
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
inverseMasses );
......
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.h
View file @ 75acdad0
......@@ -33,22 +33,10 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
private:
enum FixedParameters { GDT, EPH, EMH, EP, EM, EM_V, B, C, D, V, X, Yv, Yx, MaxFixedParameters };
enum TwoDArrayIndicies { X2D, V2D, OldV, xPrime2D, vPrime2D, Max2DArrays };
enum TwoDArrayIndicies { xPrime2D, Max2DArrays };
enum OneDArrayIndicies { InverseMasses, Max1DArrays };
RealOpenMM _tau, _accuracy;
RealOpenMM _fixedParameters[MaxFixedParameters];
/**---------------------------------------------------------------------------------------
Set fixed values
@return ReferenceDynamics::DefaultReturn
--------------------------------------------------------------------------------------- */
int _setFixedParameters( RealOpenMM timeStep, RealOpenMM prevTimeStep );
public:
......@@ -101,16 +89,6 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
void setAccuracy( RealOpenMM accuracy );
/**---------------------------------------------------------------------------------------
Get array of fixed parameters indexed by 'FixedParameters' enums
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* getFixedParameters( void ) const;
/**---------------------------------------------------------------------------------------
Print parameters
......@@ -152,9 +130,6 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
@param masses atom masses
@param inverseMasses inverse atom masses
@param xPrime xPrime
@param oldVelocities previous velocities
@param xVector xVector
@param vVector vVector
@param maxStepSize maximum time step
@return ReferenceDynamics::DefaultReturn
......@@ -163,8 +138,7 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
int updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector, RealOpenMM maxStepSize );
RealOpenMM** xPrime, RealOpenMM maxStepSize );
/**---------------------------------------------------------------------------------------
......@@ -182,8 +156,7 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
int updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
RealOpenMM** forces, RealOpenMM* inverseMasses,
RealOpenMM** xPrime, RealOpenMM** oldVelocities,
RealOpenMM** xVector, RealOpenMM** vVector );
RealOpenMM** xPrime );
};
......
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