Skip to content

GitLab

Commit 1ce7d813 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods

parent 495a43c6
Hide whitespace changes
Inline Side-by-side
Showing with 485 additions and 350 deletions
platforms/brook/src/BrookBonded.cpp
View file @ 1ce7d813
......@@ -736,7 +736,10 @@ int BrookBonded::_addRBTorsions( int *nbondeds, int *particles, float *params[],
(*nbondeds)++;
}
index = 0;
// note -- we are starting w/ index c1 not c0 -- c0 is not used
// to calculate the forces!
index = 1;
PARAMS( ibonded, 0, 0 ) = (BrookOpenMMFloat) rbParameters[index++];
PARAMS( ibonded, 0, 1 ) = (BrookOpenMMFloat) rbParameters[index++];
PARAMS( ibonded, 0, 2 ) = (BrookOpenMMFloat) rbParameters[index++];
......@@ -746,7 +749,7 @@ int BrookBonded::_addRBTorsions( int *nbondeds, int *particles, float *params[],
if( debug && getLog() ){
(void) fprintf( getLog(), " %d [%d %d %d %d] %.3e %.3e %.3e %.3e\n", ibonded, i, j, k, l,
rbParameters[0], rbParameters[1],
rbParameters[2], rbParameters[3], rbParameters[4] );
rbParameters[2], rbParameters[3], rbParameters[4], rbParameters[5] );
}
}
......@@ -1200,7 +1203,7 @@ int BrookBonded::setup( int numberOfParticles,
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookBonded::setup";
static const int PrintOn = 1;
static const int PrintOn = 0;
double dangleValue = 0.0;
// ---------------------------------------------------------------------------------------
......@@ -1991,7 +1994,17 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
}
if( getInverseMapStreamCount( J_Stream ) == 5 && getInverseMapStreamCount( L_Stream ) == 2 ){
if( getInverseMapStreamCount( J_Stream ) == 1 && getInverseMapStreamCount( L_Stream ) == 1 ){
kinvmap_gather1_1( width,
inverseStreamMaps[J_Stream][0]->getBrookStream(),
bondedForceStreams[J_Stream]->getBrookStream(),
inverseStreamMaps[L_Stream][0]->getBrookStream(),
bondedForceStreams[L_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() );
} else if( getInverseMapStreamCount( J_Stream ) == 5 && getInverseMapStreamCount( L_Stream ) == 2 ){
kinvmap_gather5_2( width,
inverseStreamMaps[J_Stream][0]->getBrookStream(),
inverseStreamMaps[J_Stream][1]->getBrookStream(),
......
platforms/brook/src/BrookBrownianDynamics.cpp
View file @ 1ce7d813
......@@ -349,7 +349,7 @@ int BrookBrownianDynamics::updateParameters( double temperature, double friction
static int showUpdate = 1;
static int maxShowUpdate = 3;
static const char* methodName = "\nBrookBrownianDynamics::updateParameters";
static const std::string methodName = "\nBrookBrownianDynamics::updateParameters";
// ---------------------------------------------------------------------------------------
......@@ -393,14 +393,12 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
// ---------------------------------------------------------------------------------------
// unused Shake parameter
float omega = 1.0f;
float one = 1.0f;
float omega = 1.0f;
float one = 1.0f;
static const std::string methodName = "\nBrookBrownianDynamics::update";
static const char* methodName = "\nBrookBrownianDynamics::update";
static const int PrintOn = 0;
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -794,7 +792,7 @@ int BrookBrownianDynamics::setup( const std::vector<double>& masses, const Platf
// ---------------------------------------------------------------------------------------
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
setLog( brookPlatform.getLog() );
int numberOfParticles = (int) masses.size();
......
platforms/brook/src/BrookCalcHarmonicBondForceKernel.cpp
View file @ 1ce7d813
......@@ -126,6 +126,7 @@ void BrookCalcHarmonicBondForceKernel::initialize( const System& system, const H
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookCalcHarmonicBondForceKernel::initialize";
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -133,8 +134,6 @@ void BrookCalcHarmonicBondForceKernel::initialize( const System& system, const H
// ---------------------------------------------------------------------------------------
//(void) fprintf( getLog(), "%s\n", methodName.c_str() ); (void) fflush( getLog() );
// create _brookBondParameters object containing particle indices/parameters
int numberOfBonds = force.getNumBonds();
......@@ -165,7 +164,7 @@ void BrookCalcHarmonicBondForceKernel::initialize( const System& system, const H
_openMMBrookInterface.setTriggerForceKernel( this );
_openMMBrookInterface.setTriggerEnergyKernel( this );
if( log ){
if( PrintOn && log ){
std::string contents = _brookBondParameters->getContentsString( );
(void) fprintf( log, "%s contents\n%s", methodName.c_str(), contents.c_str() );
(void) fflush( log );
......
platforms/brook/src/BrookCalcKineticEnergyKernel.cpp
View file @ 1ce7d813
......@@ -173,6 +173,6 @@ index = 0;
energy += _masses[ii]*(velocity[index]*velocity[index] + velocity[index + 1]*velocity[index + 1] + velocity[index + 2]*velocity[index + 2]);
}
printf( " KQ=%12.5e\n", 0.5*energy );
//printf( " KQ=%12.5e\n", 0.5*energy );
return 0.5*energy;
}
platforms/brook/src/BrookCalcRBTorsionForceKernel.cpp
View file @ 1ce7d813
......@@ -128,7 +128,7 @@ void BrookCalcRBTorsionForceKernel::initialize( const System& system, const RBTo
// ---------------------------------------------------------------------------------------
FILE* log = getLog();
FILE* log = getLog();
// ---------------------------------------------------------------------------------------
......@@ -150,6 +150,7 @@ void BrookCalcRBTorsionForceKernel::initialize( const System& system, const RBTo
double parameters[NumberOfParametersInBond];
force.getTorsionParameters( ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
......@@ -170,7 +171,7 @@ void BrookCalcRBTorsionForceKernel::initialize( const System& system, const RBTo
if( log ){
std::string contents = _brookBondParameters->getContentsString( );
(void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
(void) fprintf( log, "%s contents\n%s", methodName.c_str(), contents.c_str() );
(void) fflush( log );
}
......
platforms/brook/src/BrookGbsa.cpp
View file @ 1ce7d813
......@@ -45,7 +45,7 @@ using namespace std;
*
*/
BrookGbsa::BrookGbsa( ){
BrookGbsa::BrookGbsa( ){
// ---------------------------------------------------------------------------------------
......@@ -81,8 +81,7 @@ BrookGbsa::BrookGbsa( ){
}
_bornRadiiInitialized = 0;
_cpuObc = NULL;
_charges = NULL;
//_cpuObc = NULL;
}
......@@ -107,9 +106,7 @@ BrookGbsa::~BrookGbsa( ){
delete _gbsaForceStreams[ii];
}
delete _cpuObc;
delete[] _charges;
//delete _cpuObc;
}
......@@ -443,7 +440,7 @@ int BrookGbsa::haveBornRadiiBeenInitialized( void ) const {
* @return calculate Born radii
*
*/
/*
int BrookGbsa::calculateBornRadii( const Stream& positions ){
// ---------------------------------------------------------------------------------------
......@@ -517,7 +514,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
return DefaultReturnValue;
}
*/
/**
* Initialize stream dimensions
*
......@@ -528,11 +525,11 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
*
*/
int BrookGbsa::initializeStreamSizes( int numberOfParticles, const Platform& platform ){
int BrookGbsa::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::initializeStreamSizes";
static const std::string methodName = "BrookGbsa::_initializeStreamSizes";
// ---------------------------------------------------------------------------------------
......@@ -572,11 +569,11 @@ int BrookGbsa::initializeStreamSizes( int numberOfParticles, const Platform& pla
*
*/
int BrookGbsa::initializeStreams( const Platform& platform ){
int BrookGbsa::_initializeStreams( const Platform& platform ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::initializeStreams";
static const std::string methodName = "BrookGbsa::_initializeStreams";
static const double dangleValue = 0.0;
// ---------------------------------------------------------------------------------------
......@@ -672,8 +669,8 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfParticleP
// initialize stream sizes and then Brook streams
initializeStreamSizes( numberOfParticles, platform );
initializeStreams( platform );
_initializeStreamSizes( numberOfParticles, platform );
_initializeStreams( platform );
int particleStreamSize = getGbsaParticleStreamSize();
BrookOpenMMFloat* radiiAndCharge = new BrookOpenMMFloat[particleStreamSize*2];
......@@ -681,8 +678,6 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfParticleP
memset( radiiAndCharge, 0, particleStreamSize*2*sizeof( BrookOpenMMFloat ) );
memset( scaledRadiiAndOffset, 0, particleStreamSize*2*sizeof( BrookOpenMMFloat ) );
_charges = new RealOpenMM[particleStreamSize];
// used by CpuObc to calculate initial Born radii
vector<RealOpenMM> particleRadii(numberOfParticles);
......@@ -719,7 +714,6 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfParticleP
particleRadii[vectorIndex] = static_cast<RealOpenMM> (radius);
scaleFactors[vectorIndex] = static_cast<RealOpenMM> (scalingFactor);
_charges[vectorIndex] = static_cast<RealOpenMM> (charge);
radiiAndCharge[streamIndex] = static_cast<BrookOpenMMFloat> (radius);
radiiAndCharge[streamIndex+1] = static_cast<BrookOpenMMFloat> (charge);
......@@ -758,8 +752,8 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfParticleP
obcParameters->setSolventDielectric( static_cast<RealOpenMM>(solventDielectric) );
obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(soluteDielectric) );
_cpuObc = new CpuObc(obcParameters);
_cpuObc->setIncludeAceApproximation( true );
//_cpuObc = new CpuObc(obcParameters);
//_cpuObc->setIncludeAceApproximation( true );
return DefaultReturnValue;
}
......@@ -774,11 +768,11 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfParticleP
*
* */
int BrookGbsa::initializePartialForceStreamSize( int particleStreamSize, int particleStreamWidth ){
int BrookGbsa::_initializePartialForceStreamSize( int particleStreamSize, int particleStreamWidth ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::initializePartialForceStreamSize";
static const std::string methodName = "BrookGbsa::_initializePartialForceStreamSize";
//static const int debug = 1;
// ---------------------------------------------------------------------------------------
......@@ -842,6 +836,12 @@ std::string BrookGbsa::getContentsString( int level ) const {
(void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() );
message << _getLine( tab, "Number of particles:", value );
(void) LOCAL_SPRINTF( value, "%d", includeAce() );
message << _getLine( tab, "ACE included:", value );
(void) LOCAL_SPRINTF( value, "%.5f", getDielectricOffset() );
message << _getLine( tab, "Dielectric offset:", value );
(void) LOCAL_SPRINTF( value, "%d", getNumberOfForceStreams() );
message << _getLine( tab, "Number of force streams:", value );
......@@ -866,6 +866,15 @@ std::string BrookGbsa::getContentsString( int level ) const {
(void) LOCAL_SPRINTF( value, "%d", getParticleStreamSize() );
message << _getLine( tab, "Particle stream size:", value );
(void) LOCAL_SPRINTF( value, "%d", getGbsaParticleStreamWidth() );
message << _getLine( tab, "Gbsa stream width:", value );
(void) LOCAL_SPRINTF( value, "%d", getGbsaParticleStreamHeight() );
message << _getLine( tab, "Gbsa stream height:", value );
(void) LOCAL_SPRINTF( value, "%d", getGbsaParticleStreamSize() );
message << _getLine( tab, "Gbsa stream size:", value );
(void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamWidth() );
message << _getLine( tab, "Partial force stream width:", value );
......@@ -875,6 +884,7 @@ std::string BrookGbsa::getContentsString( int level ) const {
(void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamSize() );
message << _getLine( tab, "Partial force stream size:", value );
message << _getLine( tab, "Log:", (getLog() ? Set : NotSet) );
for( int ii = 0; ii < LastStreamIndex; ii++ ){
......@@ -914,7 +924,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::executeForces";
static const int PrintOn = 1;
static const int PrintOn = 0;
float mergeNonObcForces = 1.0f;
float kcalMolTokJNM = -0.4184f;
......@@ -925,8 +935,6 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
// calculate Born radii
(void) fprintf( getLog(), "\nPost kCalculateBornRadii: obcParticleRadiiWithDielectricOffset & obcScaledParticleRadii not set correctly!!!!!\n" );
kCalculateBornRadii( (float) getNumberOfParticles(),
(float) getParticleSizeCeiling(),
(float) getDuplicationFactor(),
......@@ -940,7 +948,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
// diagnostics
if( 0 && PrintOn && getLog() ){
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\n%s Post kCalculateBornRadii: atms=%d ceil=%d dup=%d particleStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
methodName.c_str(), getNumberOfParticles(),
......@@ -1010,7 +1018,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
// ---------------------------------------------------------------------------------------
// seecond major loop
// first major OBC loop
kObcLoop1( (float) getNumberOfParticles(),
(float) getParticleSizeCeiling(),
......@@ -1020,12 +1028,9 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
getSoluteDielectric(),
getSolventDielectric(),
includeAceTerm,
positionStream.getBrookStream(),
getObcBornRadii()->getBrookStream(),
getObcParticleRadii()->getBrookStream(),
gbsaForceStreams[0]->getBrookStream(),
gbsaForceStreams[1]->getBrookStream(),
gbsaForceStreams[2]->getBrookStream(),
......@@ -1036,7 +1041,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
// diagnostics
if( 1 && PrintOn && getLog() ){
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\nPost kObcLoop1: atms=%d ceil=%d dup=%d particleStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
getNumberOfParticles(),
......@@ -1048,17 +1053,18 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
getSolventDielectric(), includeAceTerm );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
(void) fprintf( getLog(), "\nPost kObcLoop1 PositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nBornR\n" );
(void) fprintf( getLog(), "\nPost kObcLoop1 BornR\n" );
getObcBornRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\nParticleR\n" );
(void) fprintf( getLog(), "\nPost kObcLoop1 ParticleR\n" );
getObcParticleRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStreams output\n" );
(void) fprintf( getLog(), "\nPost kObcLoop1 ForceStreams output\n" );
for( int ii = 0; ii < 4; ii++ ){
(void) fprintf( getLog(), "\nPost kObcLoop1 ForceStream %d output\n", ii );
gbsaForceStreams[ii]->printToFile( getLog() );
}
......@@ -1098,34 +1104,33 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
getParticleSizeCeiling(),
getInnerLoopUnroll() );
(void) fprintf( getLog(), "\nForceStreams\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: ForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: %d ForceStreams\n", ii );
gbsaForceStreams[ii]->printToFile( getLog() );
}
(void) fprintf( getLog(), "\nObcChain\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: ObcChain\n" );
getObcChain()->printToFile( getLog() );
(void) fprintf( getLog(), "\nBornR\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: BornR\n" );
getObcBornRadii()->printToFile( getLog() );
// output
(void) fprintf( getLog(), "\nObcIntermediateForce output\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: ObcIntermediateForce output\n" );
getObcIntermediateForce()->printToFile( getLog() );
// output
(void) fprintf( getLog(), "\nObcBornRadii2 output\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: ObcBornRadii2 output\n" );
getObcBornRadii2()->printToFile( getLog() );
}
// ---------------------------------------------------------------------------------------
// second major loop
(void) fprintf( getLog(), "\nkObcLoop2 messed up see /home/friedrim/src/openmmWork/trunk/OpenMM/platforms/brook/src/gpu/kObcBaseD2.br\n" );
// second major OBC loop
kObcLoop2( (float) getNumberOfParticles(),
(float) getParticleSizeCeiling(),
......@@ -1135,7 +1140,6 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
positionStream.getBrookStream(),
getObcScaledParticleRadii()->getBrookStream(),
getObcBornRadii2()->getBrookStream(),
getObcBornRadii2()->getBrookStream(),
gbsaForceStreams[0]->getBrookStream(),
gbsaForceStreams[1]->getBrookStream(),
gbsaForceStreams[2]->getBrookStream(),
......@@ -1156,16 +1160,16 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
getPartialForceStreamWidth( ) );
BrookStreamInternal* brookStreamInternalPos = positionStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nPositionStream\n" );
(void) fprintf( getLog(), "\nPost kObcLoop2: PositionStream\n" );
brookStreamInternalPos->printToFile( getLog() );
(void) fprintf( getLog(), "\nObcScaledParticleRadii\n" );
(void) fprintf( getLog(), "\nPost kObcLoop2: ObcScaledParticleRadii\n" );
getObcScaledParticleRadii()->printToFile( getLog() );
(void) fprintf( getLog(), "\ngetObcBornRadii2\n" );
(void) fprintf( getLog(), "\nPost kObcLoop2: ObcBornRadii2\n" );
getObcBornRadii2()->printToFile( getLog() );
(void) fprintf( getLog(), "\nForceStreams\n" );
(void) fprintf( getLog(), "\nPost kObcLoop2: ForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
......@@ -1212,19 +1216,20 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
getSoluteDielectric(),
getSolventDielectric(), includeAceTerm );
(void) fprintf( getLog(), "\nPartialForceStreams\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: PartialForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
(void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: PartialForceStreams %d\n", ii );
gbsaForceStreams[ii]->printToFile( getLog() );
}
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nForceStream\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: ForceStream\n" );
brookStreamInternalF->printToFile( getLog() );
(void) fprintf( getLog(), "\nChain\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: Chain\n" );
getObcChain()->printToFile( getLog() );
(void) fprintf( getLog(), "\nBornR\n" );
(void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: BornR\n" );
getObcBornRadii()->printToFile( getLog() );
}
......
platforms/brook/src/BrookGbsa.h
View file @ 1ce7d813
......@@ -33,7 +33,6 @@
* -------------------------------------------------------------------------- */
#include <vector>
#include <set>
#include "BrookStreamImpl.h"
#include "BrookPlatform.h"
......@@ -307,12 +306,12 @@ class BrookGbsa : public BrookCommon {
*
*/
int calculateBornRadii( const Stream& positions );
// int calculateBornRadii( const Stream& positions );
/*
* Setup of Gbsa parameters
*
* @param particleParameters vector of OBC parameters [particleI][0=charge]
* @param particleParameters vector of OBC parameters [particleI][0=charge]
* [particleI][1=radius]
* [particleI][2=scaling factor]
* @param solventDielectric solvent dielectric
......@@ -337,19 +336,6 @@ class BrookGbsa : public BrookCommon {
std::string getContentsString( int level = 0 ) const;
/*
* Calculate energy
*
* @param particlePositions particle positions
* @return energy
*
* @throw OpenMMException if _cpuObc or charges are not set
*
* */
double getEnergy( const Stream& particlePositions );
/**
* Compute forces
*
......@@ -421,13 +407,10 @@ class BrookGbsa : public BrookCommon {
int _bornRadiiInitialized;
// particle charges
RealOpenMM* _charges;
// CpuObc reference
// CpuObc reference -- was used to calculate initial Born radii
// no longer used -- to be removed?
CpuObc* _cpuObc;
// CpuObc* _cpuObc;
/*
* Setup of stream dimensions
......@@ -439,7 +422,7 @@ class BrookGbsa : public BrookCommon {
*
* */
int initializeStreamSizes( int particleStreamSize, int particleStreamWidth );
int _initializeStreamSizes( int particleStreamSize, int particleStreamWidth );
/**
* Initialize stream dimensions
......@@ -451,7 +434,7 @@ class BrookGbsa : public BrookCommon {
*
*/
int initializeStreamSizes( int numberOfParticles, const Platform& platform );
int _initializeStreamSizes( int numberOfParticles, const Platform& platform );
/**
* Initialize stream dimensions and streams
......@@ -462,7 +445,7 @@ class BrookGbsa : public BrookCommon {
*
*/
int initializeStreams( const Platform& platform );
int _initializeStreams( const Platform& platform );
/*
* Setup of stream dimensions for partial force streams
......@@ -474,7 +457,7 @@ class BrookGbsa : public BrookCommon {
*
* */
int initializePartialForceStreamSize( int particleStreamSize, int particleStreamWidth );
int _initializePartialForceStreamSize( int particleStreamSize, int particleStreamWidth );
};
......
platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp
View file @ 1ce7d813
......@@ -123,7 +123,7 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
// ---------------------------------------------------------------------------------------
static const int printOn = 1;
static const int printOn = 0;
static const std::string methodName = "BrookIntegrateLangevinStepKernel::initialize";
// ---------------------------------------------------------------------------------------
......
platforms/brook/src/BrookLangevinDynamics.cpp
View file @ 1ce7d813
......@@ -231,7 +231,7 @@ int BrookLangevinDynamics::_updateDerivedParameters( void ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookLangevinDynamics::_updateDerivedParameters";
static const std::string methodName = "\nBrookLangevinDynamics::_updateDerivedParameters";
static const BrookOpenMMFloat zero = 0.0;
static const BrookOpenMMFloat one = 1.0;
......@@ -322,9 +322,9 @@ int BrookLangevinDynamics::updateParameters( double temperature, double friction
// ---------------------------------------------------------------------------------------
static int showUpdate = 1;
static int maxShowUpdate = 3;
static const char* methodName = "\nBrookLangevinDynamics::updateParameters";
static int showUpdate = 1;
static int maxShowUpdate = 3;
static const std::string methodName = "\nBrookLangevinDynamics::updateParameters";
// ---------------------------------------------------------------------------------------
......@@ -369,8 +369,8 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookLangevinDynamics::update";
static const int PrintOn = 0;
static const std::string methodName = "\nBrookLangevinDynamics::update";
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -671,7 +671,7 @@ const BrookOpenMMFloat* BrookLangevinDynamics::getDerivedParameters( void ) cons
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookLangevinDynamics::getDerivedParameters";
// static const std::string methodName = "\nBrookLangevinDynamics::getDerivedParameters";
// ---------------------------------------------------------------------------------------
......@@ -979,7 +979,7 @@ int BrookLangevinDynamics::_updateSdStreams( void ){
// ---------------------------------------------------------------------------------------
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize();
int sdParticleStreamSize = getLangevinDynamicsParticleStreamSize();
BrookOpenMMFloat* sdpc[2];
for( int ii = 0; ii < 2; ii++ ){
......@@ -1087,7 +1087,7 @@ int BrookLangevinDynamics::setup( const std::vector<double>& masses, const Platf
// ---------------------------------------------------------------------------------------
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
setLog( brookPlatform.getLog() );
int numberOfParticles = (int) masses.size();
......
platforms/brook/src/BrookRandomNumberGenerator.cpp
View file @ 1ce7d813
......@@ -66,8 +66,8 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
_rvStreamIndex = 0;
_rvStreamOffset = 0;
_numberOfShuffles = 0;
_maxShuffles = 0;
//_maxShuffles = 100;
//_maxShuffles = 0;
_maxShuffles = 100;
_loadBuffer = NULL;
_shuffleIndices = NULL;
......@@ -79,10 +79,9 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
// set randomNumber seed
_randomNumberSeed = 1393;
_randomNumberGenerator = Mersenne;
//_randomNumberGenerator = Mersenne;
_randomNumberGenerator = Kiss;
//_randomNumberSeed = randomNumberSeed ? randomNumberSeed : 1393;
//SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
}
......@@ -191,7 +190,7 @@ BrookOpenMMFloat BrookRandomNumberGenerator::_generateGromacsRandomNumber( unsig
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::_generateGromacsRandomNumber";
// static const std::string methodName = "\nBrookRandomNumberGenerator::_generateGromacsRandomNumber";
int irand;
......@@ -248,7 +247,7 @@ void BrookRandomNumberGenerator::_generateRandomsKiss( float* randomV1, float* r
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::_generateRandomsKiss";
// static const std::string methodName = "\nBrookRandomNumberGenerator::_generateRandomsKiss";
unsigned int carry = 0;
......@@ -363,9 +362,10 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
static unsigned int state[4];
static int stateInitialized = 0;
static int PrintOn = 0;
static const int reseed = 10000;
// static const char* methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsKiss";
static std::string methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsKiss";
// ---------------------------------------------------------------------------------------
......@@ -378,9 +378,9 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
state[2] = rand();
state[3] = rand();
if( getLog() ){
(void) fprintf( getLog(), "LoadGVStreamsKiss: reset state seeds stateInitialized=%d reseed=%d [%u %u %u %u]\n",
stateInitialized, reseed, state[0], state[1], state[2], state[3] );
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s reset state seeds stateInitialized=%d reseed=%d [%u %u %u %u]\n",
methodName.c_str(), stateInitialized, reseed, state[0], state[1], state[2], state[3] );
(void) fflush( getLog() );
}
......@@ -398,11 +398,11 @@ state[3] = 27587;
float* loadBuffer = _getLoadBuffer();
if( getLog() ){
if( PrintOn && getLog() ){
static float count = 0.0f;
float block = (float) (3*getRandomNumberStreamSize() );
count += 1.0f;
(void) fprintf( getLog(), "LoadGVStreamsKiss: count=%.1f ttl=%.3e no./count=%.1f %d %d\n",
(void) fprintf( getLog(), "%s: count=%.1f ttl=%.3e no./count=%.1f %d %d\n", methodName.c_str(),
count, block*count, block, getRandomNumberStreamSize(), getNumberOfRandomNumberStreams() );
(void) fflush( getLog() );
}
......@@ -418,8 +418,8 @@ state[3] = 27587;
getRandomNumberStream( jj )->loadFromArray( loadBuffer );
}
if( getLog() ){
(void) fprintf( getLog(), "LoadGVStreamsKiss: stats\n%s\n", getStatisticsString().c_str() );
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
(void) fflush( getLog() );
}
......@@ -438,7 +438,8 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne";
static const std::string methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne";
static int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -455,8 +456,8 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){
getRandomNumberStream( jj )->loadFromArray( loadBuffer );
}
if( getLog() ){
(void) fprintf( getLog(), "_loadRandomNumberStreamsMersenne: stats\n%s\n", getStatisticsString().c_str() );
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
(void) fflush( getLog() );
}
......@@ -476,17 +477,15 @@ int BrookRandomNumberGenerator::_loadGVStreamsOriginal( void ){
unsigned long int jran;
// static const char* methodName = "\nBrookRandomNumberGenerator::_loadGVStreamsOriginal";
// static const std::string methodName = "\nBrookRandomNumberGenerator::_loadGVStreamsOriginal";
// ---------------------------------------------------------------------------------------
float* loadBuffer = _getLoadBuffer();
jran = getRandomNumberSeed();
jran = getRandomNumberSeed();
for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii++ ){
//loadBuffer[i] = sdp->fgauss( &jran );
loadBuffer[ii] = _generateGromacsRandomNumber( &jran );
}
getRandomNumberStream( jj )->loadFromArray( loadBuffer );
......@@ -510,7 +509,7 @@ float* BrookRandomNumberGenerator::_getLoadBuffer( void ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookRandomNumberGenerator::_getLoadBuffer";
static const std::string methodName = "\nBrookRandomNumberGenerator::_getLoadBuffer";
// ---------------------------------------------------------------------------------------
......@@ -542,7 +541,7 @@ int* BrookRandomNumberGenerator::_getShuffleIndices( int size ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookRandomNumberGenerator::_getShuffleIndices";
static const std::string methodName = "\nBrookRandomNumberGenerator::_getShuffleIndices";
// ---------------------------------------------------------------------------------------
......@@ -587,7 +586,7 @@ int BrookRandomNumberGenerator::_loadGVShuffle( void ){
const int np = sizeof(p) / sizeof(p[0]);
const int pmax = p[np-1];
// static const char* methodName = "\nBrookRandomNumberGenerator::loadGVShuffle";
// static const std::string methodName = "\nBrookRandomNumberGenerator::loadGVShuffle";
// ---------------------------------------------------------------------------------------
......@@ -641,7 +640,7 @@ int BrookRandomNumberGenerator::_shuffleGVStreams( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookRandomNumberGenerator::_shuffleGVStreams";
// static const std::string methodName = "\nBrookRandomNumberGenerator::_shuffleGVStreams";
// ---------------------------------------------------------------------------------------
......@@ -676,7 +675,7 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
// ---------------------------------------------------------------------------------------
static const char* methodName = "BrookRandomNumberGenerator::advanceGVCursor";
static const std::string methodName = "BrookRandomNumberGenerator::advanceGVCursor";
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -773,7 +772,7 @@ BrookFloatStreamInternal* BrookRandomNumberGenerator::getRandomNumberStream( int
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nBrookRandomNumberGenerator::getRandomNumberStream";
static const std::string methodName = "\nBrookRandomNumberGenerator::getRandomNumberStream";
// ---------------------------------------------------------------------------------------
......@@ -897,7 +896,6 @@ int BrookRandomNumberGenerator::_initializeStreams( const Platform& platform ){
randomNumberStreamSize, randomNumberStreamWidth,
BrookStreamInternal::Float, dangleValue );
// insure number of random number streams is > 0
// delete if already allocated and then initialize
......
platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp
View file @ 1ce7d813
......@@ -147,8 +147,8 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
// ---------------------------------------------------------------------------------------
static const char* methodName = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
static const int debug = 0;
static const std::string methodName = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
static const int debug = 0;
// ---------------------------------------------------------------------------------------
......@@ -229,7 +229,7 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nBrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass";
// static const std::string methodName = "\nBrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0;
......
platforms/brook/src/BrookVerletDynamics.cpp
View file @ 1ce7d813
......@@ -138,9 +138,9 @@ int BrookVerletDynamics::updateParameters( double stepSize ){
// ---------------------------------------------------------------------------------------
static int showUpdate = 1;
static int maxShowUpdate = 3;
static const char* methodName = "\nBrookVerletDynamics::updateParameters";
static int showUpdate = 1;
static int maxShowUpdate = 3;
static const std::string methodName = "\nBrookVerletDynamics::updateParameters";
// ---------------------------------------------------------------------------------------
......
platforms/brook/src/gpu/kforce_CDLJ.br
View file @ 1ce7d813
......@@ -578,207 +578,356 @@ kernel void kforce14_CDLJ(
fj = -fi;
}
kernel void kbonded_CDLJ (
float epsfac,
float xstrwidth,
float4 nbparams,
float3 posq[][],
float charge<>,
float4 atoms<>,
float4 parm0<>,
float4 parm1<>,
float4 parm2<>,
float4 parm3<>,
float4 parm4<>,
out float3 fi<>,
out float3 fj<>,
out float3 fk<>,
out float3 fl<>
)
{
float3 rij, rkj, rkl, ril;
float2 ai, aj, ak, al;
float3 m, n;
float sgnphi;
float cosfac;
float ddphi, ddphi_rb;
float3 u, v, s;
float invlkj, invlij, invlkl, msq, nsq, cos_phi, sin_phi;
float3 uij, ukj, ukl;
float phi, mdphi;
float rij_d_ukj, ukj_d_rkl;
float normfac;
float qq;
float3 fi_ang, fj_ang, fk_ang;
float3 fi_bond;
float cii, ckk, cik;
float fs, st, sth;
float3 fi_pair;
float r2 ;
float torsion_numerator;
float theta1, costheta1, invsintheta1;
float theta2, costheta2, invsintheta2;
float3 posqi, posqj, posqk, posql;
float sig, eps;
ai.y = round( (atoms.x - fmod( atoms.x, xstrwidth ))/xstrwidth );
ai.x = atoms.x - ai.y * xstrwidth;
posqi = posq[ ai ];
if ( atoms.y > -0.5f ) {
aj.y = round( (atoms.y - fmod( atoms.y, xstrwidth ))/xstrwidth );
aj.x = atoms.y - aj.y * xstrwidth;
posqj = posq[ aj ];
}
else
posqj = float4( 0.0f, 0.0f, 1.0f, 0.0f );
if ( atoms.z > -0.5f ) {
ak.y = round( (atoms.z - fmod( atoms.z, xstrwidth ))/xstrwidth );
ak.x = atoms.z - ak.y * xstrwidth;
posqk = posq[ ak ];
}
else
posqk = float4( 0.0f, 1.0f, 1.0f, 0.0f );
if ( atoms.w > -0.5f ) {
al.y = round( (atoms.w - fmod( atoms.w, xstrwidth ))/xstrwidth );
al.x = atoms.w - al.y * xstrwidth;
posql = posq[ al ];
} else {
posql = float4( 1.0f, 1.0f, 1.0f, 0.0f );
}
rij = posqi.xyz - posqj.xyz;
rkj = posqk.xyz - posqj.xyz;
rkl = posqk.xyz - posql.xyz;
ril = posqi.xyz - posql.xyz;
invlij = rsqrt( dot( rij, rij ) );
invlkj = rsqrt( dot( rkj, rkj ) );
invlkl = rsqrt( dot( rkl, rkl ) );
uij = rij * invlij;
ukj = rkj * invlkj;
ukl = rkl * invlkl;
rij_d_ukj = dot( rij, ukj );
ukj_d_rkl = dot( ukj, rkl );
m = cross( uij, ukj );
n = cross( ukj, ukl );
costheta1 = clamp( rij_d_ukj * invlij, -1.0f, 1.0f );
theta1 = acos( costheta1 );
invsintheta1 = rsqrt( 1.0f - costheta1 * costheta1 );
costheta2 = clamp( ukj_d_rkl * invlkl, -1.0f, 1.0f );
theta2 = acos( costheta2 );
invsintheta2 = rsqrt( 1.0f - costheta2 * costheta2 );
cos_phi = clamp( dot(m,n) * invsintheta1 * invsintheta2, -1.0f, 1.0f );
sgnphi = sign( dot(rij, n) );
phi = sgnphi * acos( cos_phi );
mdphi = parm1.w * phi - parm1.z;
ddphi = -parm1.y * parm1.w * sin( mdphi );
cos_phi = -cos_phi;
sin_phi = -sgnphi*sqrt( clamp( 1.0f - cos_phi * cos_phi, 0.0f, 1.0f) );
ddphi_rb = 5.0f * parm1.x;
ddphi_rb = 4.0f * parm0.w + ddphi_rb * cos_phi;
ddphi_rb = 3.0f * parm0.z + ddphi_rb * cos_phi;
ddphi_rb = 2.0f * parm0.y + ddphi_rb * cos_phi;
ddphi_rb = parm0.x + ddphi_rb * cos_phi;
ddphi_rb = -ddphi_rb * sin_phi;
ddphi += ddphi_rb;
fi = (-ddphi * invlij * invsintheta1 * invsintheta1 ) * m;
fl = ( ddphi * invlkl * invsintheta2 * invsintheta2 ) * n;
u = rij_d_ukj * invlkj * fi;
v = ukj_d_rkl * invlkj * fl;
s = u - v;
fj = s - fi;
fk = -(s + fl);
kernel void kAcosV( float x<>, out float acosVal<> ){
float x2, x3, x5, x7, x9, x11, x13, x15, x17, x19;
float pi2 = 1.57079632679489661923f;
float coefficient3 = 1.0f/6.0f;
float coefficient5 = 3.0f/40.0f;
float coefficient7 = 5.0f/112.0f;
float coefficient9 = 35.0f/1152.0f;
float coefficient11 = 63.0f/2816.0f;
float coefficient13 = 231.0f/13312.0f;
float coefficient15 = 143.0f/10240.0f;
float coefficient17 = 6435.0f/557056.0f;
float coefficient19 = 12155.0f/1245184.0f;
// ---------------------------------------------------------------------------------------
x2 = x*x;
x3 = x*x2;
x5 = x3*x2;
x7 = x5*x2;
x9 = x7*x2;
x11 = x9*x2;
x13 = x11*x2;
x15 = x13*x2;
x17 = x15*x2;
x19 = x17*x2;
acosVal = x19*coefficient19;
acosVal += x17*coefficient17;
acosVal += x15*coefficient15;
acosVal += x13*coefficient13;
acosVal += x11*coefficient11;
acosVal += x9*coefficient9;
acosVal += x7*coefficient7;
acosVal += x5*coefficient5;
acosVal += x3*coefficient3 + x;
acosVal = pi2 - acosVal;
}
fs = -parm2.y * ( theta1 - parm2.x );
st = fs * invsintheta1;
st = clamp( st, -100000.0f, 100000.0f );
kernel void kAcosTanV( float xIn<>, out float acosVal<> ){
const float c1 = 48.70107004404898384f;
const float c2 = 49.5326263772254345f;
const float c3 = 9.40604244231624f;
const float c4 = 48.70107004404996166f;
const float c5 = 65.7663163908956299f;
const float c6 = 21.587934067020262f;
const float tantwelfthpi = 0.26794919243112270647255365849413f;
const float tansixthpi = 0.57735026918962576450914878050196f;
const float sixthpi = 0.52359877559829887307710723054658f;
const float halfpi = 1.5707963267948966192313216916398f;
float x2, y, x;
float sign, complement, region;
// ---------------------------------------------------------------------------------------
complement = 0.0f;
region = 0.0f;
x = xIn;
y = rsqrt( (1.0f + x)*(1.0f - x) );
x *= y;
if( x < 0.0f ){
x = -x;
sign = -1.0f;
} else {
sign = 1.0f;
}
if( x > 1.0f ){
x = 1.0f/x;
complement = 1.0f;
} else {
complement = 0.0f;
}
if( x > tantwelfthpi ){
x = (x-tansixthpi)/(1.0f + tansixthpi*x);
region = sixthpi;
} else {
region = 0.0f;
}
x2 = x*x;
y = (x*(c1 + x2*(c2 + x2*c3))/(c4 + x2*(c5 + x2*(c6 + x2))));
y += region;
y = (complement > 0.0f) ? (halfpi- y) : y;
acosVal = halfpi - sign*y;
sth = st * costheta1;
}
cik = st * invlkj * invlij;
cii = sth * invlij;
ckk = sth * invlkj;
kernel void kAcos( float xIn, out float acosVal<> ){
fi_ang = -( cik * rkj - cii * uij );
fk_ang = -( cik * rij - ckk * ukj );
fj_ang = -fi_ang - fk_ang;
const float c1 = 48.70107004404898384f;
const float c2 = 49.5326263772254345f;
const float c3 = 9.40604244231624f;
const float c4 = 48.70107004404996166f;
const float c5 = 65.7663163908956299f;
const float c6 = 21.587934067020262f;
const float tantwelfthpi = 0.26794919243112270647255365849413f;
const float tansixthpi = 0.57735026918962576450914878050196f;
const float sixthpi = 0.52359877559829887307710723054658f;
const float halfpi = 1.5707963267948966192313216916398f;
fi += fi_ang;
fj += fj_ang;
fk += fk_ang;
fs = -parm2.w * ( theta2 - parm2.z );
st = fs * invsintheta2;
st = clamp( st, -100000.0f, 100000.0f );
float x2, y, x;
float sign, complement, region;
sth = st * costheta2;
// ---------------------------------------------------------------------------------------
cik = st * invlkj * invlkl;
cii = sth * invlkj;
ckk = sth * invlkl;
// acos(x) = atan(x/sqrt(1-x**2))
fi_ang = ( cik * rkl - cii * ukj );
fk_ang = ( cik * rkj - ckk * ukl );
fj_ang = -fi_ang - fk_ang;
x = xIn;
y = rsqrt( (1.0f + x)*(1.0f - x) );
x *= y;
fj += fi_ang;
fk += fj_ang;
fl += fk_ang;
if( x < 0.0f ){
x = -x;
sign = -1.0f;
} else {
sign = 1.0f;
}
if( x > 1.0f ){
x = 1.0f/x;
complement = 1.0f;
} else {
complement = 0.0f;
}
fi_bond = -parm3.y * ( 1.0f - parm3.x * invlij ) * rij;
fi += fi_bond;
fj += -fi_bond;
if( x > tantwelfthpi ){
x = (x - tansixthpi)/(1.0f + tansixthpi*x);
region = sixthpi;
} else {
region = 0.0f;
}
x2 = x*x;
y = (x*(c1 + x2*(c2 + x2*c3))/(c4 + x2*(c5 + x2*(c6 + x2))));
fi_bond = parm3.w * ( 1.0f - parm3.z * invlkj ) * rkj;
fj += fi_bond;
fk += -fi_bond;
y += region;
y = (complement > 0.0f) ? (halfpi- y) : y;
acosVal = halfpi - sign*y;
fi_bond = -parm4.y * ( 1.0f - parm4.x * invlkl ) * rkl;
fk += fi_bond;
fl += -fi_bond;
}
if ( parm4.z > -0.5f ) {
r2 = dot( ril, ril );
fs = scalar_force_single_CDLJ( charge, epsfac, parm4.z, parm4.w, r2, nbparams );
fi_pair = fs * ril;
fi += fi_pair;
fl -= fi_pair;
}
kernel void kbonded_CDLJ (
float epsfac,
float xstrwidth,
float4 nbparams,
float3 posq[][],
float charge<>,
float4 atoms<>,
float4 parm0<>,
float4 parm1<>,
float4 parm2<>,
float4 parm3<>,
float4 parm4<>,
out float3 fi<>,
out float3 fj<>,
out float3 fk<>,
out float3 fl<>
){
float3 rij, rkj, rkl, ril;
float2 ai, aj, ak, al;
float3 m, n;
float sgnphi;
float cosfac;
float ddphi, ddphi_rb;
float3 u, v, s;
float invlkj, invlij, invlkl, msq, nsq, cos_phi, sin_phi;
float3 uij, ukj, ukl;
float phi, mdphi;
float rij_d_ukj, ukj_d_rkl;
float normfac;
float qq;
float3 fi_ang, fj_ang, fk_ang;
float3 fi_bond;
float cii, ckk, cik;
float fs, st, sth;
float3 fi_pair;
float r2 ;
float torsion_numerator;
float theta1, costheta1, invsintheta1;
float theta2, costheta2, invsintheta2;
float3 posqi, posqj, posqk, posql;
float sig, eps;
// i-coordinates
ai.y = round( (atoms.x - fmod( atoms.x, xstrwidth ))/xstrwidth );
ai.x = atoms.x - ai.y * xstrwidth;
posqi = posq[ ai ];
// j-coordinates
if ( atoms.y > -0.5f ) {
aj.y = round( (atoms.y - fmod( atoms.y, xstrwidth ))/xstrwidth );
aj.x = atoms.y - aj.y * xstrwidth;
posqj = posq[ aj ];
} else {
posqj = float3( 0.0f, 1.0f, 0.0f );
}
// k-coordinates
if ( atoms.z > -0.5f ) {
ak.y = round( (atoms.z - fmod( atoms.z, xstrwidth ))/xstrwidth );
ak.x = atoms.z - ak.y * xstrwidth;
posqk = posq[ ak ];
} else {
posqk = float3( 1.0f, 1.0f, 0.0f );
}
// l-coordinates
if ( atoms.w > -0.5f ) {
al.y = round( (atoms.w - fmod( atoms.w, xstrwidth ))/xstrwidth );
al.x = atoms.w - al.y * xstrwidth;
posql = posq[ al ];
} else {
posql = float3( 1.0f, 1.0f, 1.0f );
}
rij = posqi - posqj;
rkj = posqk - posqj;
rkl = posqk - posql;
ril = posqi - posql;
invlij = rsqrt( dot( rij, rij ) );
invlkj = rsqrt( dot( rkj, rkj ) );
invlkl = rsqrt( dot( rkl, rkl ) );
uij = rij * invlij;
ukj = rkj * invlkj;
ukl = rkl * invlkl;
rij_d_ukj = dot( rij, ukj );
ukj_d_rkl = dot( ukj, rkl );
m = cross( uij, ukj );
n = cross( ukj, ukl );
costheta1 = clamp( rij_d_ukj * invlij, -1.0f, 1.0f );
//theta1 = acos( costheta1 );
kAcos( costheta1, theta1 );
invsintheta1 = rsqrt( (1.0f - costheta1) * (1.0f + costheta1) );
costheta2 = clamp( ukj_d_rkl * invlkl, -1.0f, 1.0f );
//theta2 = acos( costheta2 );
kAcos( costheta2, theta2 );
invsintheta2 = rsqrt( (1.0f - costheta2) * (1.0f + costheta2) );
cos_phi = clamp( dot(m,n) * invsintheta1 * invsintheta2, -1.0f, 1.0f );
sgnphi = sign( dot(rij, n) );
phi = sgnphi * acos( cos_phi );
mdphi = parm1.w * phi - parm1.z;
ddphi = -parm1.y * parm1.w * sin( mdphi );
cos_phi = -cos_phi;
sin_phi = -sgnphi*sqrt( (1.0f - cos_phi) * (1.0f + cos_phi) );
ddphi_rb = 5.0f * parm1.x;
ddphi_rb = 4.0f * parm0.w + ddphi_rb * cos_phi;
ddphi_rb = 3.0f * parm0.z + ddphi_rb * cos_phi;
ddphi_rb = 2.0f * parm0.y + ddphi_rb * cos_phi;
ddphi_rb = parm0.x + ddphi_rb * cos_phi;
ddphi_rb = -ddphi_rb * sin_phi;
ddphi += ddphi_rb;
/*
fi = float3( cos_phi, sgnphi, ddphi_rb );
fj = float3( sin_phi, parm1.x, parm0.w );
fk = float3( parm0.z, parm0.y, parm0.x );
fl = float3( ddphi, 0.0f, 0.0f );
*/
fi = (-ddphi * invlij * invsintheta1 * invsintheta1 ) * m;
fl = ( ddphi * invlkl * invsintheta2 * invsintheta2 ) * n;
u = rij_d_ukj * invlkj * fi;
v = ukj_d_rkl * invlkj * fl;
s = u - v;
fj = s - fi;
fk = -(s + fl);
fs = -parm2.y * ( theta1 - parm2.x );
st = fs * invsintheta1;
st = clamp( st, -100000.0f, 100000.0f );
sth = st * costheta1;
cik = st * invlkj * invlij;
cii = sth * invlij;
ckk = sth * invlkj;
fi_ang = -( cik * rkj - cii * uij );
fk_ang = -( cik * rij - ckk * ukj );
fj_ang = -fi_ang - fk_ang;
fi += fi_ang;
fj += fj_ang;
fk += fk_ang;
fs = -parm2.w * ( theta2 - parm2.z );
st = fs * invsintheta2;
st = clamp( st, -100000.0f, 100000.0f );
sth = st * costheta2;
cik = st * invlkj * invlkl;
cii = sth * invlkj;
ckk = sth * invlkl;
fi_ang = ( cik * rkl - cii * ukj );
fk_ang = ( cik * rkj - ckk * ukl );
fj_ang = -fi_ang - fk_ang;
fj += fi_ang;
fk += fj_ang;
fl += fk_ang;
fi_bond = -parm3.y * ( 1.0f - parm3.x * invlij ) * rij;
fi += fi_bond;
fj += -fi_bond;
fi_bond = parm3.w * ( 1.0f - parm3.z * invlkj ) * rkj;
fj += fi_bond;
fk += -fi_bond;
fi_bond = -parm4.y * ( 1.0f - parm4.x * invlkl ) * rkl;
fk += fi_bond;
fl += -fi_bond;
if( parm4.z > -0.5f ){
r2 = dot( ril, ril );
fs = scalar_force_single_CDLJ( charge, epsfac, parm4.z, parm4.w, r2, nbparams );
fi_pair = fs * ril;
fi += fi_pair;
fl -= fi_pair;
}
}
platforms/brook/src/gpu/kgbsa.br
View file @ 1ce7d813
......@@ -824,33 +824,28 @@ kernel void kPostObcLoop1_nobranch(
atomIndex = pindex.x + pindex.y*atomStreamWidth;
forceIndex = atomIndex;
outstream = float4( 0.0f, 0.0f, 0.0f, 0.0f );
outstream = float4( 0.0f, 0.0f, 0.0f, 0.0f );
for( i = 0.0f; i < repfac; i += 1.0f ){
// qIndex = floor( forceIndex/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qIndex = round( (forceIndex - fmod( forceIndex, iUnroll))/iUnroll );
qOff = forceIndex - iUnroll*qIndex;
qOff = forceIndex - iUnroll*qIndex;
// pindex.y = floor( qIndex/ pStreamWidth );
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
// pindex.x = qIndex - pindex.y*pStreamWidth + qOff;
pindex.x = qIndex - pindex.y*pStreamWidth;
pindex.y = round( (qIndex - fmod( qIndex, pStreamWidth ))/pStreamWidth );
pindex.x = qIndex - pindex.y*pStreamWidth;
o1 = pstream1[ pindex ];
o2 = pstream2[ pindex ];
o3 = pstream3[ pindex ];
o4 = pstream4[ pindex ];
o1 = pstream1[ pindex ];
o2 = pstream2[ pindex ];
o3 = pstream3[ pindex ];
o4 = pstream4[ pindex ];
tmp = qOff < 0.5f ? o1 : o2;
tmp = qOff < 1.5f ? tmp : o3;
tmp = qOff < 2.5f ? tmp : o4;
tmp = qOff < 0.5f ? o1 : o2;
tmp = qOff < 1.5f ? tmp : o3;
tmp = qOff < 2.5f ? tmp : o4;
outstream += tmp;
outstream += tmp;
forceIndex += roundNatoms;
forceIndex += roundNatoms;
}
bornRadii2Force = obcChain*bornRadii*bornRadii*outstream.w;
......@@ -1140,7 +1135,7 @@ kernel void loop2InternalNoSum( float3 d1, float3 d2, float3 d3, float3 d4, floa
//1376 FLOPS
kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicationFactor,
float streamWidth, float fstreamWidth, float3 posq[][],
float scaledAtomicRadii[][], float bornForceFactor[][],
float2 scaledAtomicRadii[][], float bornForceFactor[][],
out float4 bornForce1<>, out float4 bornForce2<>,
out float4 bornForce3<>, out float4 bornForce4<> ){
......@@ -1221,7 +1216,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
i1pos = posq[ iAtom ];
i1BornForceFactor = bornForceFactor[ iAtom ];
//i1AtomicRadii = atomicRadii[ iAtom ];
i1AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i1ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce1 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
......@@ -1231,7 +1225,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
iAtom.x += 1;
i2pos = posq[ iAtom ];
i2BornForceFactor = bornForceFactor[ iAtom ];
//i2AtomicRadii = atomicRadii[ iAtom ];
i2AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i2ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce2 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
......@@ -1241,7 +1234,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
iAtom.x += 1;
i3pos = posq[ iAtom ];
i3BornForceFactor = bornForceFactor[ iAtom ];
//i3AtomicRadii = atomicRadii[ iAtom ];
i3AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i3ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce3 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
......@@ -1251,24 +1243,22 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
iAtom.x += 1;
i4pos = posq[ iAtom ];
i4BornForceFactor = bornForceFactor[ iAtom ];
//i4AtomicRadii = atomicRadii[ iAtom ];
i4AtomicRadii = scaledAtomicRadii[ iAtom ].y;
i4ScaledAtomicRadii = scaledAtomicRadii[ iAtom ].x;
bornForce4 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
// create float4 for main vars
//iBornForceFactor = float4( i1BornForceFactor, i2BornForceFactor, i3BornForceFactor, i4BornForceFactor );
iAtomicRadii = float4( i1AtomicRadii, i2AtomicRadii, i3AtomicRadii, i4AtomicRadii );
iScaledAtomicRadii = float4( i1ScaledAtomicRadii, i2ScaledAtomicRadii, i3ScaledAtomicRadii, i4ScaledAtomicRadii );
// inner loop indices setup
//changed the following instruction for rounding issues on some ASICs
//whichRep = floor( forceIndex / roundedUpAtoms );
tmp = fmod(forceIndex, roundedUpAtoms);
whichRep = round((forceIndex - tmp)/roundedUpAtoms);
//jBlock = 1 + floor( numberOfAtoms/(duplicationFactor*streamWidth ) );
jBlock = 1 + round( (numberOfAtoms-fmod( numberOfAtoms, (duplicationFactor*streamWidth) )) / (duplicationFactor*streamWidth) );
jStart = whichRep*jBlock;
......@@ -1314,7 +1304,6 @@ kernel void kObcLoop2( float numberOfAtoms, float roundedUpAtoms, float duplicat
jAtom.x += 1.0f;
//jBornForceFactor = float4( j1BornForceFactor, j2BornForceFactor, j3BornForceFactor, j4BornForceFactor );
jAtomicRadii = float4( j1AtomicRadii, j2AtomicRadii, j3AtomicRadii, j4AtomicRadii );
jScaledAtomicRadii = float4( j1ScaledAtomicRadii, j2ScaledAtomicRadii, j3ScaledAtomicRadii, j4ScaledAtomicRadii );
......
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