Commit 6a4a17f8 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods

parent d0ff1456
...@@ -1056,6 +1056,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP ...@@ -1056,6 +1056,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookBonded::loadInvMaps"; static const std::string methodName = "BrookBonded::loadInvMaps";
static const int PrintOn = 0;
double dangleValue = 0.0; double dangleValue = 0.0;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -1100,7 +1101,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP ...@@ -1100,7 +1101,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
} }
} }
if( getLog() ){ if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s force stream strms=%d nbondeds=%d max counts=[%d %d %d %d] strSz&Wd=%d %d\n", methodName.c_str(), getNumberOfForceStreams(), (void) fprintf( getLog(), "%s force stream strms=%d nbondeds=%d max counts=[%d %d %d %d] strSz&Wd=%d %d\n", methodName.c_str(), getNumberOfForceStreams(),
nbondeds, getMaxInverseMapStreamCount(0), getMaxInverseMapStreamCount(1), getMaxInverseMapStreamCount(2), getMaxInverseMapStreamCount(3), nbondeds, getMaxInverseMapStreamCount(0), getMaxInverseMapStreamCount(1), getMaxInverseMapStreamCount(2), getMaxInverseMapStreamCount(3),
atomStreamSize, atomStreamWidth ); atomStreamSize, atomStreamWidth );
...@@ -1110,21 +1111,24 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP ...@@ -1110,21 +1111,24 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
// load data // load data
for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){ for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
for( int jj = 0; jj < 4*getMaxInverseMapStreamCount()*atomStreamSize; jj++ ){
block[jj] = -1.0f;
}
gpuCalcInvMap( ii, 4, nbondeds, natoms, atoms, getInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) ); gpuCalcInvMap( ii, 4, nbondeds, natoms, atoms, getInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) );
//gpuPrintInvMaps( _inverseMapStreamCount[ii], natoms, counts, invmaps, getLog() ); //gpuPrintInvMaps( _inverseMapStreamCount[ii], natoms, counts, invmaps, getLog() );
validateInverseMapStreamCount( ii, _inverseMapStreamCount[ii] ); validateInverseMapStreamCount( ii, _inverseMapStreamCount[ii] );
for( int jj = 0; jj < _inverseMapStreamCount[ii]; jj++ ){ for( int jj = 0; jj < _inverseMapStreamCount[ii]; jj++ ){
_inverseStreamMaps[ii][jj]->loadFromArray( invmaps[jj] ); _inverseStreamMaps[ii][jj]->loadFromArray( invmaps[jj] );
(void) fprintf( getLog(), "%s inverseMap stream strms=%d count=%d index=%d %d InverseMapStreamCount[ii]=%d max=%d\n", if( PrintOn && getLog() ){
methodName.c_str(), getNumberOfForceStreams(), _inverseMapStreamCount[ii], ii, jj, (void) fprintf( getLog(), "%s inverseMap stream strms=%d count=%d index=%d %d InverseMapStreamCount[ii]=%d max=%d\n",
getInverseMapStreamCount( ii ), getMaxInverseMapStreamCount( ii ) ); methodName.c_str(), getNumberOfForceStreams(), _inverseMapStreamCount[ii], ii, jj,
getInverseMapStreamCount( ii ), getMaxInverseMapStreamCount( ii ) );
/* for( int kk = 0; kk < atomStreamSize; kk++ ){
for( int kk = 0; kk < atomStreamSize; kk++ ){ (void) fprintf( getLog(), "%8d [ %.1f %.1f %.1f %.1f]\n", kk, invmaps[jj][kk].x, invmaps[jj][kk].y, invmaps[jj][kk].z, invmaps[jj][kk].w );
(void) fprintf( getLog(), "%8d [ %.1f %.1f %.1f %.1f]\n", kk, invmaps[jj][kk].x, invmaps[jj][kk].y, invmaps[jj][kk].z, invmaps[jj][kk].w ); }
} }
*/
} }
...@@ -1144,7 +1148,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP ...@@ -1144,7 +1148,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
// diagnostics // diagnostics
if( getLog() ){ if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s done\n", methodName.c_str() ); (void) fprintf( getLog(), "%s done\n", methodName.c_str() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1199,6 +1203,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1199,6 +1203,7 @@ int BrookBonded::setup( int numberOfAtoms,
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookBonded::setup"; static const std::string methodName = "BrookBonded::setup";
static const int PrintOn = 0;
double dangleValue = 0.0; double dangleValue = 0.0;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -1213,7 +1218,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1213,7 +1218,7 @@ int BrookBonded::setup( int numberOfAtoms,
std::stringstream message; std::stringstream message;
message << methodName << " number of harmonic bond atom indices=" << bondIndices.size() << " does not equal number of harmonic bond parameter entries=" << bondParameters.size(); message << methodName << " number of harmonic bond atom indices=" << bondIndices.size() << " does not equal number of harmonic bond parameter entries=" << bondParameters.size();
throw OpenMMException( message.str() ); throw OpenMMException( message.str() );
} else if( getLog() ){ } else if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s harmonic bonds=%d\n", methodName.c_str(), bondIndices.size() ); (void) fprintf( getLog(), "%s harmonic bonds=%d\n", methodName.c_str(), bondIndices.size() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1222,7 +1227,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1222,7 +1227,7 @@ int BrookBonded::setup( int numberOfAtoms,
std::stringstream message; std::stringstream message;
message << methodName << " number of angle atom indices=" << angleIndices.size() << " does not equal number of angle parameter entries=" << angleParameters.size(); message << methodName << " number of angle atom indices=" << angleIndices.size() << " does not equal number of angle parameter entries=" << angleParameters.size();
throw OpenMMException( message.str() ); throw OpenMMException( message.str() );
} else if( getLog() ){ } else if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s angle bonds=%d\n", methodName.c_str(), angleIndices.size() ); (void) fprintf( getLog(), "%s angle bonds=%d\n", methodName.c_str(), angleIndices.size() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1231,7 +1236,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1231,7 +1236,7 @@ int BrookBonded::setup( int numberOfAtoms,
std::stringstream message; std::stringstream message;
message << methodName << " number of periodicTorsion atom indices=" << periodicTorsionIndices.size() << " does not equal number of periodicTorsion parameter entries=" << periodicTorsionParameters.size(); message << methodName << " number of periodicTorsion atom indices=" << periodicTorsionIndices.size() << " does not equal number of periodicTorsion parameter entries=" << periodicTorsionParameters.size();
throw OpenMMException( message.str() ); throw OpenMMException( message.str() );
} else if( getLog() ){ } else if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s periodicTorsion bonds=%d\n", methodName.c_str(), periodicTorsionIndices.size() ); (void) fprintf( getLog(), "%s periodicTorsion bonds=%d\n", methodName.c_str(), periodicTorsionIndices.size() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1240,7 +1245,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1240,7 +1245,7 @@ int BrookBonded::setup( int numberOfAtoms,
std::stringstream message; std::stringstream message;
message << methodName << " number of rbTorsion atom indices=" << rbTorsionIndices.size() << " does not equal number of rbTorsion parameter entries=" << rbTorsionParameters.size(); message << methodName << " number of rbTorsion atom indices=" << rbTorsionIndices.size() << " does not equal number of rbTorsion parameter entries=" << rbTorsionParameters.size();
throw OpenMMException( message.str() ); throw OpenMMException( message.str() );
} else if( getLog() ){ } else if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s rbTorsion bonds=%d\n", methodName.c_str(), rbTorsionIndices.size() ); (void) fprintf( getLog(), "%s rbTorsion bonds=%d\n", methodName.c_str(), rbTorsionIndices.size() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1249,7 +1254,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1249,7 +1254,7 @@ int BrookBonded::setup( int numberOfAtoms,
std::stringstream message; std::stringstream message;
message << methodName << " number atoms=" << numberOfAtoms << " does not equal number of nb parameter entries=" << nonbondedParameters.size(); message << methodName << " number atoms=" << numberOfAtoms << " does not equal number of nb parameter entries=" << nonbondedParameters.size();
throw OpenMMException( message.str() ); throw OpenMMException( message.str() );
} else if( getLog() ){ } else if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s LJ 14 ixns=%d\n", methodName.c_str(), bonded14Indices.size() ); (void) fprintf( getLog(), "%s LJ 14 ixns=%d\n", methodName.c_str(), bonded14Indices.size() );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1321,7 +1326,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1321,7 +1326,7 @@ int BrookBonded::setup( int numberOfAtoms,
std::stringstream message; std::stringstream message;
message << methodName << " number of bonds=" << nbondeds << " is greater than maxBonds=" << maxBonds; message << methodName << " number of bonds=" << nbondeds << " is greater than maxBonds=" << maxBonds;
throw OpenMMException( message.str() ); throw OpenMMException( message.str() );
} else if( getLog() ){ } else if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s atoms=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfAtoms, nbondeds, maxBonds ); (void) fprintf( getLog(), "%s atoms=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfAtoms, nbondeds, maxBonds );
(void) fflush( getLog() ); (void) fflush( getLog() );
} }
...@@ -1364,7 +1369,7 @@ int BrookBonded::setup( int numberOfAtoms, ...@@ -1364,7 +1369,7 @@ int BrookBonded::setup( int numberOfAtoms,
// debug stuff // debug stuff
if( 1 && getLog() ){ if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s nbondeds=%d strDim [%d %d ] sz=%d\n", methodName.c_str(), nbondeds, (void) fprintf( getLog(), "%s nbondeds=%d strDim [%d %d ] sz=%d\n", methodName.c_str(), nbondeds,
_atomIndicesStream->getStreamWidth(), _atomIndicesStream->getStreamWidth(),
......
...@@ -169,7 +169,8 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -169,7 +169,8 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookCalcGBSAOBCForceFieldKernel::executeForces"; static const std::string methodName = "BrookCalcGBSAOBCForceFieldKernel::executeForces";
static const int PrintOn = 1;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -184,7 +185,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -184,7 +185,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// calculate Born radii first time thru and initialize on board // calculate Born radii first time thru and initialize on board
if( _brookGbsa->haveBornRadiiBeenInitialized() ){ if( !_brookGbsa->haveBornRadiiBeenInitialized() ){
_brookGbsa->calculateBornRadii( positions ); _brookGbsa->calculateBornRadii( positions );
} }
...@@ -210,6 +211,40 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -210,6 +211,40 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
gbsaForceStreams[3]->getBrookStream() gbsaForceStreams[3]->getBrookStream()
); );
// ---------------------------------------------------------------------------------------
// diagnostics
if( 1 && PrintOn ){
(void) fprintf( getLog(), "\nPost kObcLoop1: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
_brookGbsa->getNumberOfAtoms(),
_brookGbsa->getAtomSizeCeiling(),
_brookGbsa->getDuplicationFactor(),
_brookGbsa->getAtomStreamWidth( ),
_brookGbsa->getPartialForceStreamWidth( ),
_brookGbsa->getSoluteDielectric(),
_brookGbsa->getSolventDielectric(), includeAce );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nBornR\n" );
_brookGbsa->getObcBornRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\nAtomR\n" );
_brookGbsa->getObcAtomicRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStreams output\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
}
// ---------------------------------------------------------------------------------------
// gather for first loop // gather for first loop
kPostObcLoop1_nobranch( kPostObcLoop1_nobranch(
...@@ -228,6 +263,45 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -228,6 +263,45 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
_brookGbsa->getObcIntermediateForce()->getBrookStream(), _brookGbsa->getObcIntermediateForce()->getBrookStream(),
_brookGbsa->getObcBornRadii2()->getBrookStream() ); _brookGbsa->getObcBornRadii2()->getBrookStream() );
// ---------------------------------------------------------------------------------------
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: dup=%d aStrW=%d pStrW=%d no.atms=%3d ceil=%3d Unroll=%1d\n",
_brookGbsa->getDuplicationFactor(),
_brookGbsa->getAtomStreamWidth( ),
_brookGbsa->getPartialForceStreamWidth( ),
_brookGbsa->getNumberOfAtoms(),
_brookGbsa->getAtomSizeCeiling(),
_brookGbsa->getInnerLoopUnroll() );
(void) fprintf( getLog(), "\nForceStreams output\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
(void) fprintf( getLog(), "\nObcChain\n" );
_brookGbsa->getObcChain()->printToFile( getLog() );
(void) fprintf( getLog(), "\nBornR\n" );
_brookGbsa->getObcBornRadii()->printToFile( getLog() );
// output
(void) fprintf( getLog(), "\nObcIntermediateForce output\n" );
_brookGbsa->getObcIntermediateForce()->printToFile( getLog() );
// output
(void) fprintf( getLog(), "\ngetObcBornRadii2 output\n" );
_brookGbsa->getObcBornRadii2()->printToFile( getLog() );
}
// ---------------------------------------------------------------------------------------
// second major loop // second major loop
kObcLoop2( (float) _brookGbsa->getNumberOfAtoms(), kObcLoop2( (float) _brookGbsa->getNumberOfAtoms(),
...@@ -236,7 +310,6 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -236,7 +310,6 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
(float) _brookGbsa->getAtomStreamWidth( ), (float) _brookGbsa->getAtomStreamWidth( ),
(float) _brookGbsa->getPartialForceStreamWidth( ), (float) _brookGbsa->getPartialForceStreamWidth( ),
positionStream.getBrookStream(), positionStream.getBrookStream(),
_brookGbsa->getObcAtomicRadiiWithDielectricOffset()->getBrookStream(),
_brookGbsa->getObcScaledAtomicRadii()->getBrookStream(), _brookGbsa->getObcScaledAtomicRadii()->getBrookStream(),
_brookGbsa->getObcBornRadii2()->getBrookStream(), _brookGbsa->getObcBornRadii2()->getBrookStream(),
gbsaForceStreams[0]->getBrookStream(), gbsaForceStreams[0]->getBrookStream(),
...@@ -245,6 +318,38 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -245,6 +318,38 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
gbsaForceStreams[3]->getBrookStream() gbsaForceStreams[3]->getBrookStream()
); );
// ---------------------------------------------------------------------------------------
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\nPost kObcLoop2: no.atms=%5d ceil=%3d dup=%3d strW=%3d pStrW=%3d\n",
_brookGbsa->getNumberOfAtoms(),
_brookGbsa->getAtomSizeCeiling(),
_brookGbsa->getDuplicationFactor(),
_brookGbsa->getAtomStreamWidth( ),
_brookGbsa->getPartialForceStreamWidth( ) );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nObcScaledAtomicRadii\n" );
_brookGbsa->getObcScaledAtomicRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\ngetObcBornRadii2\n" );
_brookGbsa->getObcBornRadii2()->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
}
// ---------------------------------------------------------------------------------------
// gather for second loop // gather for second loop
float mergeNonObcForces = 1.0f; float mergeNonObcForces = 1.0f;
...@@ -270,6 +375,37 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S ...@@ -270,6 +375,37 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
forceStream.getBrookStream() forceStream.getBrookStream()
); );
// ---------------------------------------------------------------------------------------
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\nPost kObcLoop2: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
_brookGbsa->getNumberOfAtoms(),
_brookGbsa->getAtomSizeCeiling(),
_brookGbsa->getDuplicationFactor(),
_brookGbsa->getAtomStreamWidth( ),
_brookGbsa->getPartialForceStreamWidth( ),
_brookGbsa->getSoluteDielectric(),
_brookGbsa->getSolventDielectric(), includeAce );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nForceStream\n" );
brookStreamInternalF->printToFile( getLog() );
(void) fprintf( getLog(), "\nChain\n" );
_brookGbsa->getObcChain()->printToFile( getLog() );
(void) fprintf( getLog(), "\nBornR\n" );
_brookGbsa->getObcBornRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\nPartialForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
}
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
} }
......
...@@ -231,7 +231,6 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions ...@@ -231,7 +231,6 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
const BrookStreamImpl& positionStreamC = dynamic_cast<const BrookStreamImpl&> (positions.getImpl()); const BrookStreamImpl& positionStreamC = dynamic_cast<const BrookStreamImpl&> (positions.getImpl());
BrookStreamImpl& positionStream = const_cast<BrookStreamImpl&> (positionStreamC); BrookStreamImpl& positionStream = const_cast<BrookStreamImpl&> (positionStreamC);
...@@ -278,7 +277,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -278,7 +277,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
// diagnostics // diagnostics
if( PrintOn ){ if( 1 && PrintOn ){
(void) fprintf( getLog(), "\nPost knbforce_CDLJ4: atoms=%6d ceiling=%3d dupFac=%3d", _brookNonBonded->getNumberOfAtoms(), (void) fprintf( getLog(), "\nPost knbforce_CDLJ4: atoms=%6d ceiling=%3d dupFac=%3d", _brookNonBonded->getNumberOfAtoms(),
_brookNonBonded->getAtomSizeCeiling(), _brookNonBonded->getAtomSizeCeiling(),
_brookNonBonded->getDuplicationFactor() ); _brookNonBonded->getDuplicationFactor() );
...@@ -309,6 +308,8 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -309,6 +308,8 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
} }
} }
// ---------------------------------------------------------------------------------------
// gather forces // gather forces
kMergeFloat3_4_nobranch( (float) _brookNonBonded->getDuplicationFactor(), kMergeFloat3_4_nobranch( (float) _brookNonBonded->getDuplicationFactor(),
...@@ -323,6 +324,18 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -323,6 +324,18 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
nonbondedForceStreams[3]->getBrookStream(), nonbondedForceStreams[3]->getBrookStream(),
forceStream.getBrookStream() ); forceStream.getBrookStream() );
// diagnostics
if( 0 && PrintOn ){
(void) fprintf( getLog(), "\nNB forces" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
}
// ---------------------------------------------------------------------------------------
// bonded // bonded
epsfac = (float) (_brookBonded->getLJ_14Scale()*_brookBonded->getCoulombFactor()); epsfac = (float) (_brookBonded->getLJ_14Scale()*_brookBonded->getCoulombFactor());
...@@ -377,7 +390,14 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -377,7 +390,14 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
(void) fprintf( getLog(), "\nForce stream %d\n", ii ); (void) fprintf( getLog(), "\nForce stream %d\n", ii );
bondedForceStreams[ii]->printToFile( getLog() ); bondedForceStreams[ii]->printToFile( getLog() );
} }
(void) fprintf( getLog(), "\nInverse map streams\n" );
/*
(void) fprintf( getLog(), "\nNB1 forces" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
*/
(void) fprintf( getLog(), "\nInverse map streams -- K_Stream cnt=%d\n", _brookBonded->getInverseMapStreamCount( K_Stream ) );
for( int ii = 0; ii < 4; ii++ ){ for( int ii = 0; ii < 4; ii++ ){
for( int jj = 0; jj < countPrintInvMap[ii]; jj++ ){ for( int jj = 0; jj < countPrintInvMap[ii]; jj++ ){
(void) fprintf( getLog(), "\n Inverse map streams index=%d %d\n", ii, jj ); (void) fprintf( getLog(), "\n Inverse map streams index=%d %d\n", ii, jj );
...@@ -389,7 +409,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -389,7 +409,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
// gather forces // gather forces
if( _brookBonded->getInverseMapStreamCount( K_Stream ) <= 4 ){ if( _brookBonded->getInverseMapStreamCount( K_Stream ) <= 4 ){
// kinvmap_gather3_4( (float) bp->width, bp->strInvMapi[0], bp->strInvMapi[1], bp->strInvMapi[2], bp->fi, bp->strInvMapk[0], bp->strInvMapk[1], bp->strInvMapk[2], bp->strInvMapk[3], bp->fk, gpu->strF, gpu->strF );
kinvmap_gather3_4( width, kinvmap_gather3_4( width,
inverseStreamMaps[I_Stream][0]->getBrookStream(), inverseStreamMaps[I_Stream][0]->getBrookStream(),
...@@ -406,7 +426,6 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -406,7 +426,6 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
forceStream.getBrookStream(), forceStream.getBrookStream() ); forceStream.getBrookStream(), forceStream.getBrookStream() );
} else if( _brookBonded->getInverseMapStreamCount( K_Stream ) == 5 ){ } else if( _brookBonded->getInverseMapStreamCount( K_Stream ) == 5 ){
// kinvmap_gather3_5( (float) bp->width, bp->strInvMapi[0], bp->strInvMapi[1], bp->strInvMapi[2], bp->fi, bp->strInvMapk[0], bp->strInvMapk[1], bp->strInvMapk[2], bp->strInvMapk[3], bp->strInvMapk[4], bp->fk, gpu->strF, gpu->strF );
kinvmap_gather3_5( width, kinvmap_gather3_5( width,
inverseStreamMaps[I_Stream][0]->getBrookStream(), inverseStreamMaps[I_Stream][0]->getBrookStream(),
...@@ -436,6 +455,16 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -436,6 +455,16 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
} }
// diagnostics
if( 0 && PrintOn ){
(void) fprintf( getLog(), "\nPost 3_4/3_5 && NB forces" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
}
//kinvmap_gather5_2( (float) bp->width, bp->strInvMapj[0], bp->strInvMapj[1], bp->strInvMapj[2], bp->strInvMapj[3], bp->strInvMapj[4], bp->fj, //kinvmap_gather5_2( (float) bp->width, bp->strInvMapj[0], bp->strInvMapj[1], bp->strInvMapj[2], bp->strInvMapj[3], bp->strInvMapj[4], bp->fj,
// bp->strInvMapl[0], bp->strInvMapl[1], bp->fl, gpu->strF, gpu->strF ); // bp->strInvMapl[0], bp->strInvMapl[1], bp->fl, gpu->strF, gpu->strF );
kinvmap_gather5_2( width, kinvmap_gather5_2( width,
...@@ -450,6 +479,24 @@ nonbondedForceStreams[3]->fillWithValue( &zerof ); ...@@ -450,6 +479,24 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
bondedForceStreams[L_Stream]->getBrookStream(), bondedForceStreams[L_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() ); forceStream.getBrookStream(), forceStream.getBrookStream() );
// diagnostics
if( 1 && PrintOn ){
(void) fprintf( getLog(), "\nFinal NB & bonded forces" );
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
brookStreamInternalF->printToFile( getLog() );
/*
void* dataV = brookStreamInternalF->getData(1);
float* data = (float*) dataV;
(void) fprintf( getLog(), "\nFinal NB & bonded forces RAW\n" );
for( int ii = 0; ii < _brookNonBonded->getNumberOfAtoms()*3; ii += 3 ){
(void) fprintf( getLog(), "%d [%.6e %.6e %.6e]\n", ii, data[ii], data[ii+1], data[ii+2] );
}
*/
}
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
} }
......
...@@ -522,17 +522,20 @@ int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log ){ ...@@ -522,17 +522,20 @@ int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
void* dataArrayV = getDataArray( ); assert( log );
saveToArray( dataArrayV );
int streamSize = getStreamSize(); void* dataArrayV = getData( 1 );
int width = getWidth(); //void* dataArrayV = getDataArray( );
int index = 0; //saveToArray( dataArrayV );
float* dataArray = (float*) dataArrayV;
int streamSize = getStreamSize();
int width = getWidth();
int index = 0;
const float* dataArray = (float*) dataArrayV;
for( int ii = 0; ii < streamSize; ii++ ){ for( int ii = 0; ii < streamSize; ii++ ){
std::stringstream message; std::stringstream message;
message.width( 15 ); message.width( 15 );
message.precision( 5 ); message.precision( 8 );
message << ii << " [ "; message << ii << " [ ";
for( int jj = 0; jj < width; jj++ ){ for( int jj = 0; jj < width; jj++ ){
message << dataArray[index++] << " "; message << dataArray[index++] << " ";
...@@ -541,7 +544,7 @@ int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log ){ ...@@ -541,7 +544,7 @@ int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log ){
(void) fprintf( log, "%s", message.str().c_str() ); (void) fprintf( log, "%s", message.str().c_str() );
} }
delete[] dataArrayV; //delete[] dataArrayV;
return DefaultReturnValue; return DefaultReturnValue;
......
...@@ -223,6 +223,7 @@ int BrookGbsa::getAtomSizeCeiling( void ) const { ...@@ -223,6 +223,7 @@ int BrookGbsa::getAtomSizeCeiling( void ) const {
if( localThis->_atomSizeCeiling ){ if( localThis->_atomSizeCeiling ){
localThis->_atomSizeCeiling = localThis->getOuterLoopUnroll() - localThis->_atomSizeCeiling; localThis->_atomSizeCeiling = localThis->getOuterLoopUnroll() - localThis->_atomSizeCeiling;
} }
localThis->_atomSizeCeiling += localThis->getNumberOfAtoms();
} }
return _atomSizeCeiling; return _atomSizeCeiling;
...@@ -442,6 +443,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){ ...@@ -442,6 +443,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::calculateBornRadii"; static const std::string methodName = "BrookGbsa::calculateBornRadii";
static const int PrintOn = 1;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -463,11 +465,15 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){ ...@@ -463,11 +465,15 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
RealOpenMM* bornRadii = new RealOpenMM[streamSize]; RealOpenMM* bornRadii = new RealOpenMM[streamSize];
memset( bornRadii, 0, sizeof( RealOpenMM )*streamSize ); memset( bornRadii, 0, sizeof( RealOpenMM )*streamSize );
int index = 0; RealOpenMM* obcChain = new RealOpenMM[streamSize];
memset( obcChain, 0, sizeof( RealOpenMM )*streamSize );
int index = 0;
RealOpenMM* atomCoordinatesBlkPtr = atomCoordinatesBlk;
for( int ii = 0; ii < numberOfAtoms; ii++ ){ for( int ii = 0; ii < numberOfAtoms; ii++ ){
atomCoordinates[ii] = atomCoordinatesBlk; atomCoordinates[ii] = atomCoordinatesBlkPtr;
atomCoordinatesBlk += 3; atomCoordinatesBlkPtr += 3;
atomCoordinates[ii][0] = coordinates[index++]; atomCoordinates[ii][0] = coordinates[index++];
atomCoordinates[ii][1] = coordinates[index++]; atomCoordinates[ii][1] = coordinates[index++];
...@@ -477,11 +483,23 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){ ...@@ -477,11 +483,23 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
// calculate Born radii // calculate Born radii
_cpuObc->computeBornRadii( atomCoordinates, bornRadii ); _cpuObc->computeBornRadii( atomCoordinates, bornRadii, obcChain );
// diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\n%s: atms=%d\n", methodName.c_str(), numberOfAtoms );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
(void) fprintf( getLog(), "%d coord=[%.5e %.5e %.5e] bR=%.5e obcChain=%.6e\n", ii,
atomCoordinates[ii][0], atomCoordinates[ii][1], atomCoordinates[ii][2], bornRadii[ii], obcChain[ii] );
}
}
// write radii to board and set flag to indicate radii calculated once // write radii to board and set flag to indicate radii calculated once
_gbsaStreams[ObcBornRadiiStream]->loadFromArray( bornRadii ); _gbsaStreams[ObcBornRadiiStream]->loadFromArray( bornRadii );
_gbsaStreams[ObcChainStream]->loadFromArray( obcChain );
_bornRadiiInitialized = 1; _bornRadiiInitialized = 1;
// free memory // free memory
...@@ -489,6 +507,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){ ...@@ -489,6 +507,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
delete[] atomCoordinatesBlk; delete[] atomCoordinatesBlk;
delete[] atomCoordinates; delete[] atomCoordinates;
delete[] bornRadii; delete[] bornRadii;
delete[] obcChain;
return DefaultReturnValue; return DefaultReturnValue;
} }
...@@ -573,9 +592,9 @@ int BrookGbsa::initializeStreams( const Platform& platform ){ ...@@ -573,9 +592,9 @@ int BrookGbsa::initializeStreams( const Platform& platform ){
_gbsaStreams[ObcIntermediateForceStream] = new BrookFloatStreamInternal( BrookCommon::ObcIntermediateForceStream, _gbsaStreams[ObcIntermediateForceStream] = new BrookFloatStreamInternal( BrookCommon::ObcIntermediateForceStream,
gbsaAtomStreamSize, gbsaAtomStreamWidth, gbsaAtomStreamSize, gbsaAtomStreamWidth,
BrookStreamInternal::Float, dangleValue ); BrookStreamInternal::Float4, dangleValue );
// Born2 radii // Obc chain
_gbsaStreams[ObcChainStream] = new BrookFloatStreamInternal( BrookCommon::ObcChainStream, _gbsaStreams[ObcChainStream] = new BrookFloatStreamInternal( BrookCommon::ObcChainStream,
gbsaAtomStreamSize, gbsaAtomStreamWidth, gbsaAtomStreamSize, gbsaAtomStreamWidth,
...@@ -613,7 +632,7 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfAtomParam ...@@ -613,7 +632,7 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfAtomParam
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
static const int atomParametersSize = 4; static const int atomParametersSize = 3;
static const int maxErrors = 20; static const int maxErrors = 20;
static const std::string methodName = "BrookGbsa::setup"; static const std::string methodName = "BrookGbsa::setup";
...@@ -822,13 +841,6 @@ std::string BrookGbsa::getContentsString( int level ) const { ...@@ -822,13 +841,6 @@ std::string BrookGbsa::getContentsString( int level ) const {
message << _getLine( tab, "Partial force stream size:", value ); message << _getLine( tab, "Partial force stream size:", value );
message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) ); message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) );
/*
message << _getLine( tab, "ExclusionStream:", (getExclusionStream() ? Set : NotSet) );
message << _getLine( tab, "VdwStream:", (getOuterVdwStream() ? Set : NotSet) );
message << _getLine( tab, "ChargeStream:", (getChargeStream() ? Set : NotSet) );
message << _getLine( tab, "SigmaStream:", (getInnerSigmaStream() ? Set : NotSet) );
message << _getLine( tab, "EpsilonStream:", (getInnerEpsilonStream() ? Set : NotSet) );
*/
for( int ii = 0; ii < LastStreamIndex; ii++ ){ for( int ii = 0; ii < LastStreamIndex; ii++ ){
message << std::endl; message << std::endl;
......
...@@ -80,6 +80,8 @@ void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses, ...@@ -80,6 +80,8 @@ void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses,
_brookShakeAlgorithm = new BrookShakeAlgorithm( ); _brookShakeAlgorithm = new BrookShakeAlgorithm( );
_brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() ); _brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() );
// assert( (_brookShakeAlgorithm->getNumberOfConstraints() > 0) );
_brookRandomNumberGenerator = new BrookRandomNumberGenerator( ); _brookRandomNumberGenerator = new BrookRandomNumberGenerator( );
_brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() ); _brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
} }
...@@ -118,19 +120,12 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc ...@@ -118,19 +120,12 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
differences[0] = temperature - (double) _brookStochasticDynamics->getTemperature(); differences[0] = temperature - (double) _brookStochasticDynamics->getTemperature();
differences[1] = friction - (double) _brookStochasticDynamics->getFriction(); differences[1] = friction - (double) _brookStochasticDynamics->getFriction();
differences[2] = stepSize - (double) _brookStochasticDynamics->getStepSize(); differences[2] = stepSize - (double) _brookStochasticDynamics->getStepSize();
if( fabs( differences[0] ) < epsilon || fabs( differences[1] ) < epsilon || fabs( differences[2] ) < epsilon ){ if( fabs( differences[0] ) > epsilon || fabs( differences[1] ) > epsilon || fabs( differences[2] ) > epsilon ){
//printf( "%s calling updateParameters\n", methodName.c_str() );
_brookStochasticDynamics->updateParameters( temperature, friction, stepSize ); _brookStochasticDynamics->updateParameters( temperature, friction, stepSize );
} } else {
//printf( "%s NOT calling updateParameters\n", methodName.c_str() );
assert( _brookShakeAlgorithm ); }
/*
if( _brookShakeAlgorithm == NULL ){
std::stringstream message;
message << methodName << " _brookShakeAlgorithm is not set -- case not handled.";
throw OpenMMException( message.str() );
return NULL;
}
*/
_brookStochasticDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm, *_brookRandomNumberGenerator ); _brookStochasticDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm, *_brookRandomNumberGenerator );
......
...@@ -47,6 +47,11 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel { ...@@ -47,6 +47,11 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public: public:
// return values
static const int DefaultReturnValue = 0;
static const int ErrorReturnValue = -1;
/** /**
* BrookIntegrateLangevinStepKernel constructor * BrookIntegrateLangevinStepKernel constructor
* *
......
...@@ -144,7 +144,7 @@ void BrookPlatform::_initializeKernelFactory( void ){ ...@@ -144,7 +144,7 @@ void BrookPlatform::_initializeKernelFactory( void ){
BrookKernelFactory* factory = new BrookKernelFactory(); BrookKernelFactory* factory = new BrookKernelFactory();
registerKernelFactory( CalcStandardMMForceFieldKernel::Name(), factory); registerKernelFactory( CalcStandardMMForceFieldKernel::Name(), factory);
// registerKernelFactory( CalcGBSAOBCForceFieldKernel::Name(), factory); registerKernelFactory( CalcGBSAOBCForceFieldKernel::Name(), factory);
//registerKernelFactory( IntegrateVerletStepKernel::Name(), factory); //registerKernelFactory( IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory( IntegrateLangevinStepKernel::Name(), factory); registerKernelFactory( IntegrateLangevinStepKernel::Name(), factory);
//registerKernelFactory( IntegrateBrownianStepKernel::Name(), factory); //registerKernelFactory( IntegrateBrownianStepKernel::Name(), factory);
......
...@@ -36,6 +36,7 @@ ...@@ -36,6 +36,7 @@
#include "BrookStreamImpl.h" #include "BrookStreamImpl.h"
#include "kshakeh.h" #include "kshakeh.h"
#include "kupdatesd.h" #include "kupdatesd.h"
#include "kcommon.h"
// use random number generator // use random number generator
...@@ -222,13 +223,15 @@ int BrookStochasticDynamics::_setStepSize( BrookOpenMMFloat stepSize ){ ...@@ -222,13 +223,15 @@ int BrookStochasticDynamics::_setStepSize( BrookOpenMMFloat stepSize ){
* *
* @return DefaultReturnValue * @return DefaultReturnValue
* *
* @throw if tau too small
*
*/ */
int BrookStochasticDynamics::_updateDerivedParameters( void ){ int BrookStochasticDynamics::_updateDerivedParameters( void ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookStochasticDynamics::_updateDerivedParameters"; static const char* methodName = "\nBrookStochasticDynamics::_updateDerivedParameters";
static const BrookOpenMMFloat zero = 0.0; static const BrookOpenMMFloat zero = 0.0;
static const BrookOpenMMFloat one = 1.0; static const BrookOpenMMFloat one = 1.0;
...@@ -237,12 +240,21 @@ int BrookStochasticDynamics::_updateDerivedParameters( void ){ ...@@ -237,12 +240,21 @@ int BrookStochasticDynamics::_updateDerivedParameters( void ){
static const BrookOpenMMFloat four = 4.0; static const BrookOpenMMFloat four = 4.0;
static const BrookOpenMMFloat half = 0.5; static const BrookOpenMMFloat half = 0.5;
float epsilon = 1.0e-08f;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
BrookOpenMMFloat tau = getTau(); BrookOpenMMFloat tau = getTau();
BrookOpenMMFloat temperature = getTemperature(); BrookOpenMMFloat temperature = getTemperature();
BrookOpenMMFloat stepSize = getStepSize(); BrookOpenMMFloat stepSize = getStepSize();
if( fabsf( float( tau ) ) < epsilon ){
std::stringstream message;
message << methodName << " tau=" << tau << " too small.";
throw OpenMMException( message.str() );
}
_derivedParameters[GDT] = stepSize/tau; _derivedParameters[GDT] = stepSize/tau;
_derivedParameters[EPH] = EXP( half*_derivedParameters[GDT] ); _derivedParameters[EPH] = EXP( half*_derivedParameters[GDT] );
...@@ -293,7 +305,7 @@ int BrookStochasticDynamics::_updateDerivedParameters( void ){ ...@@ -293,7 +305,7 @@ int BrookStochasticDynamics::_updateDerivedParameters( void ){
return DefaultReturnValue; return DefaultReturnValue;
}; }
/** /**
* Update parameters -- only way parameters can be set * Update parameters -- only way parameters can be set
...@@ -310,17 +322,29 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti ...@@ -310,17 +322,29 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookStochasticDynamics::updateParameters"; static int showUpdate = 1;
static int maxShowUpdate = 3;
static const char* methodName = "\nBrookStochasticDynamics::updateParameters";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
_setStepSize( (BrookOpenMMFloat) stepSize ); _setStepSize( (BrookOpenMMFloat) stepSize );
_setFriction( (BrookOpenMMFloat) friction ); _setFriction( (BrookOpenMMFloat) friction );
//_setTau( (BrookOpenMMFloat) friction );
_setTemperature( (BrookOpenMMFloat) temperature ); _setTemperature( (BrookOpenMMFloat) temperature );
_updateDerivedParameters( ); _updateDerivedParameters( );
_updateSdStreams( ); _updateSdStreams( );
// show update
if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
std::string contents = getContentsString( );
(void) fprintf( getLog(), "%s contents %s", methodName, contents.c_str() );
(void) fflush( getLog() );
}
return DefaultReturnValue; return DefaultReturnValue;
}; };
...@@ -343,51 +367,103 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -343,51 +367,103 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
BrookShakeAlgorithm& brookShakeAlgorithm, BrookShakeAlgorithm& brookShakeAlgorithm,
BrookRandomNumberGenerator& brookRandomNumberGenerator ){ BrookRandomNumberGenerator& brookRandomNumberGenerator ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// unused Shake parameter // unused Shake parameter
float omega = 1.0f; float omega = 1.0f;
// static const char* methodName = "\nBrookStochasticDynamics::update"; static const char* methodName = "\nBrookStochasticDynamics::update";
// --------------------------------------------------------------------------------------- static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
const BrookOpenMMFloat* derivedParameters = getDerivedParameters(); const BrookOpenMMFloat* derivedParameters = getDerivedParameters();
BrookStreamImpl& positionStream = dynamic_cast<BrookStreamImpl&> (positions.getImpl()); BrookStreamImpl& positionStream = dynamic_cast<BrookStreamImpl&> (positions.getImpl());
BrookStreamImpl& velocityStream = dynamic_cast<BrookStreamImpl&> (velocities.getImpl()); BrookStreamImpl& velocityStream = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
const BrookStreamImpl& forceStreamC = dynamic_cast<const BrookStreamImpl&> (forces.getImpl()); const BrookStreamImpl& forceStreamC = dynamic_cast<const BrookStreamImpl&> (forces.getImpl());
BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC); BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
// first integration step if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() );
kupdate_sd1_fix1( (void) fflush( getLog() );
(float) getStochasticDynamicsAtomStreamWidth(),
(float) brookRandomNumberGenerator.getRandomNumberStreamWidth(), }
(float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM], // first integration step
derivedParameters[Sd1pc1],
derivedParameters[Sd1pc2], kupdate_sd1_fix1(
derivedParameters[Sd1pc3], (float) getStochasticDynamicsAtomStreamWidth(),
getSDPC1Stream()->getBrookStream(), (float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(), (float) brookRandomNumberGenerator.getRvStreamOffset(),
getSD2XStream()->getBrookStream(), derivedParameters[EM],
positionStream.getBrookStream(), derivedParameters[Sd1pc1],
forceStream.getBrookStream(), derivedParameters[Sd1pc2],
velocityStream.getBrookStream(), derivedParameters[Sd1pc3],
getInverseMassStream()->getBrookStream(), getSDPC1Stream()->getBrookStream(),
getSD1VStream()->getBrookStream(), brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getVPrimeStream()->getBrookStream(), getSD2XStream()->getBrookStream(),
getXPrimeStream()->getBrookStream() positionStream.getBrookStream(),
); forceStream.getBrookStream(),
velocityStream.getBrookStream(),
getInverseMassStream()->getBrookStream(),
getSD1VStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
getXPrimeStream()->getBrookStream()
);
// advance random number cursor // diagnostics
if( PrintOn ){
(void) fprintf( getLog(), "\nPost kupdate_sd1_fix1: atomStrW=%3d rngStrW=%3d rngOff=%3d"
"EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]",
getStochasticDynamicsAtomStreamWidth(),
brookRandomNumberGenerator.getRandomNumberStreamWidth(),
brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] );
(void) fprintf( getLog(), "\nSDPC1Stream\n" );
getSDPC1Stream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nSD2XStream\n" );
getSD2XStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nInverseMassStream\n" );
getInverseMassStream()->printToFile( getLog() );
//StreamImpl& positionStreamImpl = positionStream.getImpl();
//const BrookStreamImpl brookPositions = dynamic_cast<BrookStreamImpl&> (positionStreamImpl);
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(), "\nSD1VStream\n" );
getSD1VStream()->printToFile( getLog() );
(void) fprintf( getLog(), "\nVPrimeStream\n" );
getVPrimeStream()->printToFile( getLog() );
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() ); (void) fprintf( getLog(), "\nXPrimeStream\n" );
getXPrimeStream()->printToFile( getLog() );
}
// first Shake step // advance random number cursor
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() );
// first Shake step
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
kshakeh_fix1( kshakeh_fix1(
10.0f, 10.0f,
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
...@@ -401,9 +477,9 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -401,9 +477,9 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() ); brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// first Shake gather // first Shake gather
kshakeh_update1_fix1( kshakeh_update1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc1],
...@@ -416,32 +492,34 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -416,32 +492,34 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() ); getXPrimeStream()->getBrookStream() );
}
// second integration step // second integration step
kupdate_sd2_fix1( kupdate_sd2_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
(float) brookRandomNumberGenerator.getRandomNumberStreamWidth(), (float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
(float) brookRandomNumberGenerator.getRvStreamOffset(), (float) brookRandomNumberGenerator.getRvStreamOffset(),
derivedParameters[Sd2pc1], derivedParameters[Sd2pc1],
derivedParameters[Sd2pc2], derivedParameters[Sd2pc2],
getSDPC2Stream()->getBrookStream(), getSDPC2Stream()->getBrookStream(),
brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(), brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->getBrookStream(),
getSD1VStream()->getBrookStream(), getSD1VStream()->getBrookStream(),
positionStream.getBrookStream(), positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(), getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(), getVPrimeStream()->getBrookStream(),
getSD2XStream()->getBrookStream(), getSD2XStream()->getBrookStream(),
velocityStream.getBrookStream(), velocityStream.getBrookStream(),
getXPrimeStream()->getBrookStream() getXPrimeStream()->getBrookStream()
); );
// advance random number cursor // advance random number cursor
brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() ); brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() );
// second Shake step // second Shake step
if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
kshakeh_fix1( kshakeh_fix1(
10.0f, 10.0f,
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
...@@ -455,22 +533,9 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -455,22 +533,9 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() ); brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
// second Shake gather // second Shake gather
kshakeh_update1_fix1(
(float) getStochasticDynamicsAtomStreamWidth(),
derivedParameters[Sd2pc1],
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
positionStream.getBrookStream(),
getXPrimeStream()->getBrookStream(),
getVPrimeStream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
getXPrimeStream()->getBrookStream() );
kshakeh_update2_fix1( kshakeh_update2_fix1(
(float) getStochasticDynamicsAtomStreamWidth(), (float) getStochasticDynamicsAtomStreamWidth(),
brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(), brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
...@@ -481,6 +546,11 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities, ...@@ -481,6 +546,11 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(), brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
positionStream.getBrookStream() ); positionStream.getBrookStream() );
} else {
//kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
}
return DefaultReturnValue; return DefaultReturnValue;
...@@ -813,6 +883,8 @@ int BrookStochasticDynamics::_updateSdStreams( void ){ ...@@ -813,6 +883,8 @@ int BrookStochasticDynamics::_updateSdStreams( void ){
sdpc[ii] = new BrookOpenMMFloat[2*sdAtomStreamSize]; sdpc[ii] = new BrookOpenMMFloat[2*sdAtomStreamSize];
memset( sdpc[ii], 0, 2*sdAtomStreamSize*sizeof( BrookOpenMMFloat ) ); memset( sdpc[ii], 0, 2*sdAtomStreamSize*sizeof( BrookOpenMMFloat ) );
} }
BrookOpenMMFloat* inverseMass = new BrookOpenMMFloat[sdAtomStreamSize];
memset( inverseMass, 0, sdAtomStreamSize*sizeof( BrookOpenMMFloat ) );
const BrookOpenMMFloat* derivedParameters = getDerivedParameters( ); const BrookOpenMMFloat* derivedParameters = getDerivedParameters( );
int numberOfAtoms = getNumberOfAtoms(); int numberOfAtoms = getNumberOfAtoms();
...@@ -825,14 +897,18 @@ int BrookStochasticDynamics::_updateSdStreams( void ){ ...@@ -825,14 +897,18 @@ int BrookStochasticDynamics::_updateSdStreams( void ){
sdpc[1][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yx]) ); sdpc[1][index] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yx]) );
sdpc[1][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[X]) ); sdpc[1][index+1] = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[X]) );
inverseMass[ii] = _inverseSqrtMasses[ii]*_inverseSqrtMasses[ii];
} }
_sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] ); _sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] );
_sdStreams[SDPC2Stream]->loadFromArray( sdpc[1] ); _sdStreams[SDPC2Stream]->loadFromArray( sdpc[1] );
_sdStreams[InverseMassStream]->loadFromArray( inverseMass );
for( int ii = 0; ii < 2; ii++ ){ for( int ii = 0; ii < 2; ii++ ){
delete[] sdpc[ii]; delete[] sdpc[ii];
} }
delete[] inverseMass;
// initialize SD2X // initialize SD2X
...@@ -908,6 +984,9 @@ int BrookStochasticDynamics::setup( const std::vector<double>& masses, const Pla ...@@ -908,6 +984,9 @@ int BrookStochasticDynamics::setup( const std::vector<double>& masses, const Pla
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
setLog( brookPlatform.getLog() );
int numberOfAtoms = (int) masses.size(); int numberOfAtoms = (int) masses.size();
setNumberOfAtoms( numberOfAtoms ); setNumberOfAtoms( numberOfAtoms );
......
...@@ -281,3 +281,12 @@ brook::stream& BrookStreamImpl::getBrookStream( void ){ ...@@ -281,3 +281,12 @@ brook::stream& BrookStreamImpl::getBrookStream( void ){
return _brookStreamInternal->getBrookStream( ); return _brookStreamInternal->getBrookStream( );
} }
/**
* Get Brook stream impl
*
* @return Brook stream impl
*/
BrookStreamInternal* BrookStreamImpl::getBrookStreamImpl( void ) const {
return _brookStreamInternal;
}
...@@ -166,6 +166,14 @@ class BrookStreamImpl : public StreamImpl { ...@@ -166,6 +166,14 @@ class BrookStreamImpl : public StreamImpl {
//const std::string getContentsString( int level = 0 ) const; //const std::string getContentsString( int level = 0 ) const;
/**
* Get Brook stream impl
*
* @return Brook stream impl
*/
BrookStreamInternal* getBrookStreamImpl( void ) const;
protected: protected:
BrookStreamInternal* _brookStreamInternal; BrookStreamInternal* _brookStreamInternal;
......
...@@ -24,3 +24,12 @@ kernel void kzerof4( out float4 outstr<> ) { ...@@ -24,3 +24,12 @@ kernel void kzerof4( out float4 outstr<> ) {
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<> ) {
outstr = instr;
}
kernel void kadd3( float3 val<>, out float3 outstr<> ) {
outstr += val;
}
...@@ -75,28 +75,6 @@ void kAddAndMergeFloat4_4( ...@@ -75,28 +75,6 @@ void kAddAndMergeFloat4_4(
::brook::stream outputStream ::brook::stream outputStream
); );
void kPostObcLoop2(
const float repfac,
const float atomStrWidth,
const float pforceStrWidth,
const float natoms,
const float roundNatoms,
const float iUnroll,
const float conversion,
const float mergeNonObcForces,
::brook::stream inObcForces,
::brook::stream nonObcForces,
::brook::stream stream1,
::brook::stream stream2,
::brook::stream stream3,
::brook::stream stream4,
::brook::stream atomicRadii,
::brook::stream bornRadii,
::brook::stream obcChain,
::brook::stream outputStream
);
void kPostObcLoop2_nobranch( void kPostObcLoop2_nobranch(
const float repfac, const float repfac,
const float atomStrWidth, const float atomStrWidth,
...@@ -472,7 +450,6 @@ void kObcLoop2( ...@@ -472,7 +450,6 @@ void kObcLoop2(
float fstreamWidth, float fstreamWidth,
::brook::stream posq, ::brook::stream posq,
::brook::stream atomicRadii, ::brook::stream atomicRadii,
::brook::stream scaledAtomicRadii,
::brook::stream bornForceFactor, ::brook::stream bornForceFactor,
::brook::stream bornForce1, ::brook::stream bornForce1,
::brook::stream bornForce2, ::brook::stream bornForce2,
......
...@@ -56,6 +56,59 @@ kernel float3 do_gather( float strwidth, float4 invmap<>, float3 forces[][] ) { ...@@ -56,6 +56,59 @@ kernel float3 do_gather( float strwidth, float4 invmap<>, float3 forces[][] ) {
return f; return f;
} }
//helper function to make the unrolling look better
kernel float3 do_gather_nobranch( float strwidth, float4 invmap<>, float3 forces[][] ) {
float3 f1, f2, f3, f4, f, z, t;
float4 quotient, remainder;
float2 idx;
float m1, m2, m3, m4;
z = float3( 0.0f, 0.0f, 0.0f );
f=z;
//Convert from linear to 2D index
// quotient = floor( invmap / strwidth );
quotient = round( ( invmap - fmod(invmap, strwidth))/strwidth );
remainder = invmap - quotient * strwidth;
m1 = invmap.x;
m2 = invmap.y;
m3 = invmap.z;
m4 = invmap.w;
//Add each force only if non-negative
idx.y = quotient.x;
idx.x = remainder.x;
f1 = forces[ idx ];
idx.y = quotient.y;
idx.x = remainder.y;
f2 = forces[ idx ];
idx.y = quotient.z;
idx.x = remainder.z;
f3 = forces[ idx ];
idx.y = quotient.w;
idx.x = remainder.w;
f4 = forces[ idx ];
f = (m1 >= 0.0f) ? f1 : z;
t = (m2 >= 0.0f) ? f2 : z;
f = f+t;
t = (m3 >= 0.0f) ? f3 : z;
f = f+t;
t = (m4 >= 0.0f) ? f4 : z;
f = f+t;
return f;
}
//Simple version, takes only one index stream //Simple version, takes only one index stream
kernel void kinvmap_gather( kernel void kinvmap_gather(
float strwidth, //stream width of the dihedral forces float strwidth, //stream width of the dihedral forces
...@@ -67,7 +120,7 @@ kernel void kinvmap_gather( ...@@ -67,7 +120,7 @@ kernel void kinvmap_gather(
{ {
outforce = inforce; outforce = inforce;
outforce += do_gather( strwidth, invmap, forces ); outforce += do_gather_nobranch( strwidth, invmap, forces );
} }
//Takes two inverse maps //Takes two inverse maps
...@@ -82,8 +135,8 @@ kernel void kinvmap_gather2( ...@@ -82,8 +135,8 @@ kernel void kinvmap_gather2(
{ {
outforce = inforce; outforce = inforce;
outforce += do_gather( strwidth, invmap1, forces ); outforce += do_gather_nobranch( strwidth, invmap1, forces );
outforce += do_gather( strwidth, invmap2, forces ); outforce += do_gather_nobranch( strwidth, invmap2, forces );
} }
//Takes three inverse maps //Takes three inverse maps
...@@ -99,9 +152,9 @@ kernel void kinvmap_gather3( ...@@ -99,9 +152,9 @@ kernel void kinvmap_gather3(
{ {
outforce = inforce; outforce = inforce;
outforce += do_gather( strwidth, invmap1, forces ); outforce += do_gather_nobranch( strwidth, invmap1, forces );
outforce += do_gather( strwidth, invmap2, forces ); outforce += do_gather_nobranch( strwidth, invmap2, forces );
outforce += do_gather( strwidth, invmap3, forces ); outforce += do_gather_nobranch( strwidth, invmap3, forces );
} }
//Takes four inverse maps //Takes four inverse maps
...@@ -118,10 +171,10 @@ kernel void kinvmap_gather4( ...@@ -118,10 +171,10 @@ kernel void kinvmap_gather4(
{ {
outforce = inforce; outforce = inforce;
outforce += do_gather( strwidth, invmap1, forces ); outforce += do_gather_nobranch( strwidth, invmap1, forces );
outforce += do_gather( strwidth, invmap2, forces ); outforce += do_gather_nobranch( strwidth, invmap2, forces );
outforce += do_gather( strwidth, invmap3, forces ); outforce += do_gather_nobranch( strwidth, invmap3, forces );
outforce += do_gather( strwidth, invmap4, forces ); outforce += do_gather_nobranch( strwidth, invmap4, forces );
} }
//Takes five inverse maps //Takes five inverse maps
...@@ -139,11 +192,11 @@ kernel void kinvmap_gather5( ...@@ -139,11 +192,11 @@ kernel void kinvmap_gather5(
{ {
outforce = inforce; outforce = inforce;
outforce += do_gather( strwidth, invmap1, forces ); outforce += do_gather_nobranch( strwidth, invmap1, forces );
outforce += do_gather( strwidth, invmap2, forces ); outforce += do_gather_nobranch( strwidth, invmap2, forces );
outforce += do_gather( strwidth, invmap3, forces ); outforce += do_gather_nobranch( strwidth, invmap3, forces );
outforce += do_gather( strwidth, invmap4, forces ); outforce += do_gather_nobranch( strwidth, invmap4, forces );
outforce += do_gather( strwidth, invmap5, forces ); outforce += do_gather_nobranch( strwidth, invmap5, forces );
} }
//Takes six inverse maps - this is the last one! //Takes six inverse maps - this is the last one!
...@@ -162,15 +215,131 @@ kernel void kinvmap_gather6( ...@@ -162,15 +215,131 @@ kernel void kinvmap_gather6(
{ {
outforce = inforce; outforce = inforce;
outforce += do_gather( strwidth, invmap1, forces ); outforce += do_gather_nobranch( strwidth, invmap1, forces );
outforce += do_gather( strwidth, invmap2, forces ); outforce += do_gather_nobranch( strwidth, invmap2, forces );
outforce += do_gather( strwidth, invmap3, forces ); outforce += do_gather_nobranch( strwidth, invmap3, forces );
outforce += do_gather( strwidth, invmap4, forces ); outforce += do_gather_nobranch( strwidth, invmap4, forces );
outforce += do_gather( strwidth, invmap5, forces ); outforce += do_gather_nobranch( strwidth, invmap5, forces );
outforce += do_gather( strwidth, invmap6, forces ); outforce += do_gather_nobranch( strwidth, invmap6, forces );
} }
//Takes three + four inverse maps
kernel void kinvmap_gather3_4(
float strwidth, //stream width of the dihedral forces
float4 invmap3_1<>, //indices into the dihedral forces
float4 invmap3_2<>, //indices into the dihedral forces
float4 invmap3_3<>,
float3 forces3[][], //dihedral forces
float4 invmap4_1<>, //indices into the dihedral forces
float4 invmap4_2<>, //indices into the dihedral forces
float4 invmap4_3<>,
float4 invmap4_4<>,
float3 forces4[][], //dihedral forces
float3 inforce<>, //particle forces before
out float3 outforce<> //particle forces after
)
{
outforce = inforce;
outforce += do_gather_nobranch( strwidth, invmap3_1, forces3 );
outforce += do_gather_nobranch( strwidth, invmap3_2, forces3 );
outforce += do_gather_nobranch( strwidth, invmap3_3, forces3 );
outforce += do_gather_nobranch( strwidth, invmap4_1, forces4 );
outforce += do_gather_nobranch( strwidth, invmap4_2, forces4 );
outforce += do_gather_nobranch( strwidth, invmap4_3, forces4 );
outforce += do_gather_nobranch( strwidth, invmap4_4, forces4 );
}
//Takes three + five inverse maps
kernel void kinvmap_gather3_5(
float strwidth, //stream width of the dihedral forces
float4 invmap3_1[][], //indices into the dihedral forces
float4 invmap3_2[][], //indices into the dihedral forces
float4 invmap3_3[][],
float3 forces3[][], //dihedral forces
float4 invmap5_1[][], //indices into the dihedral forces
float4 invmap5_2[][], //indices into the dihedral forces
float4 invmap5_3[][],
float4 invmap5_4[][],
float4 invmap5_5[][],
float3 forces5[][], //dihedral forces
float3 inforce<>, //particle forces before
out float3 outforce<> //particle forces after
)
{
float2 idx = indexof(outforce);
outforce = inforce;
outforce += do_gather_nobranch( strwidth, invmap3_1[idx], forces3 );
outforce += do_gather_nobranch( strwidth, invmap3_2[idx], forces3 );
outforce += do_gather_nobranch( strwidth, invmap3_3[idx], forces3 );
outforce += do_gather_nobranch( strwidth, invmap5_1[idx], forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_2[idx], forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_3[idx], forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_4[idx], forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_5[idx], forces5 );
}
//Takes five + two inverse maps
kernel void kinvmap_gather5_2(
float strwidth, //stream width of the dihedral forces
float4 invmap5_1<>, //indices into the dihedral forces
float4 invmap5_2<>, //indices into the dihedral forces
float4 invmap5_3<>,
float4 invmap5_4<>,
float4 invmap5_5<>,
float3 forces5[][], //dihedral forces
float4 invmap2_1<>, //indices into the dihedral forces
float4 invmap2_2<>, //indices into the dihedral forces
float3 forces2[][], //dihedral forces
float3 inforce<>, //particle forces before
out float3 outforce<> //particle forces after
)
{
outforce = inforce;
outforce += do_gather_nobranch( strwidth, invmap5_1, forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_2, forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_3, forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_4, forces5 );
outforce += do_gather_nobranch( strwidth, invmap5_5, forces5 );
outforce += do_gather_nobranch( strwidth, invmap2_1, forces2 );
outforce += do_gather_nobranch( strwidth, invmap2_2, forces2 );
}
//Takes five + two inverse maps
kernel void kinvmap_gather4_2(
float strwidth, //stream width of the dihedral forces
float4 invmap4_1<>, //indices into the dihedral forces
float4 invmap4_2<>, //indices into the dihedral forces
float4 invmap4_3<>,
float4 invmap4_4<>,
float3 forces4[][], //dihedral forces
float4 invmap2_1<>, //indices into the dihedral forces
float4 invmap2_2<>, //indices into the dihedral forces
float3 forces2[][], //dihedral forces
float3 inforce<>, //particle forces before
out float3 outforce<> //particle forces after
)
{
outforce = inforce;
outforce += do_gather_nobranch( strwidth, invmap4_1, forces4 );
outforce += do_gather_nobranch( strwidth, invmap4_2, forces4 );
outforce += do_gather_nobranch( strwidth, invmap4_3, forces4 );
outforce += do_gather_nobranch( strwidth, invmap4_4, forces4 );
outforce += do_gather_nobranch( strwidth, invmap2_1, forces2 );
outforce += do_gather_nobranch( strwidth, invmap2_2, forces2 );
}
kernel float3 etch_force( kernel float3 etch_force(
float fpos, float fpos,
......
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