Skip to content

GitLab

"platforms/vscode:/vscode.git/clone" did not exist on "a2f4ab8b47b87a5e7bbb7400b75d90ea600c58e6"
Commit 22ecd657 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

BrookCalcStandardMMForceFieldKernel: added _refForceField to calculate energy 

                                     overrides of bonded gathers for testing
BrookGbsa:                           OBC parameters fixed; added energy calculation
BrookCommon:                         added utility methods for allocating arrays
parent 18c9a78a
Hide whitespace changes
Inline Side-by-side
Showing with 499 additions and 51 deletions
platforms/brook/src/BrookCalcGBSAOBCForceFieldKernel.cpp
View file @ 22ecd657
......@@ -64,6 +64,7 @@ BrookCalcGBSAOBCForceFieldKernel::BrookCalcGBSAOBCForceFieldKernel( std::string
_numberOfAtoms = 0;
_brookGbsa = NULL;
_log = NULL;
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
if( brookPlatform.getLog() != NULL ){
......@@ -170,7 +171,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookCalcGBSAOBCForceFieldKernel::executeForces";
static const int PrintOn = 1;
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -215,7 +216,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// diagnostics
if( 1 && PrintOn ){
if( 1 && PrintOn && getLog() ){
(void) fprintf( getLog(), "\nPost kObcLoop1: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
_brookGbsa->getNumberOfAtoms(),
......@@ -267,7 +268,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// diagnostics
if( PrintOn ){
if( PrintOn && getLog()){
(void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: dup=%d aStrW=%d pStrW=%d no.atms=%3d ceil=%3d Unroll=%1d\n",
_brookGbsa->getDuplicationFactor(),
......@@ -277,7 +278,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
_brookGbsa->getAtomSizeCeiling(),
_brookGbsa->getInnerLoopUnroll() );
(void) fprintf( getLog(), "\nForceStreams output\n" );
(void) fprintf( getLog(), "\nForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
......@@ -295,7 +296,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// output
(void) fprintf( getLog(), "\ngetObcBornRadii2 output\n" );
(void) fprintf( getLog(), "\nObcBornRadii2 output\n" );
_brookGbsa->getObcBornRadii2()->printToFile( getLog() );
}
......@@ -322,7 +323,7 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// diagnostics
if( PrintOn ){
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\nPost kObcLoop2: no.atms=%5d ceil=%3d dup=%3d strW=%3d pStrW=%3d\n",
_brookGbsa->getNumberOfAtoms(),
......@@ -379,9 +380,9 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
// diagnostics
if( PrintOn ){
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\nPost kObcLoop2: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
(void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
_brookGbsa->getNumberOfAtoms(),
_brookGbsa->getAtomSizeCeiling(),
_brookGbsa->getDuplicationFactor(),
......@@ -390,6 +391,11 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
_brookGbsa->getSoluteDielectric(),
_brookGbsa->getSolventDielectric(), includeAce );
(void) fprintf( getLog(), "\nPartialForceStreams\n" );
for( int ii = 0; ii < 4; ii++ ){
gbsaForceStreams[ii]->printToFile( getLog() );
}
BrookStreamInternal* brookStreamInternalF = forceStream.getBrookStreamImpl();
(void) fprintf( getLog(), "\nForceStream\n" );
brookStreamInternalF->printToFile( getLog() );
......@@ -400,22 +406,17 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
(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() );
}
}
// ---------------------------------------------------------------------------------------
}
/**
* Execute the kernel to calculate the energy.
* Execute the kernel to calculate the OBC energy
*
* @param positions atom positions
*
* @return potential energy due to the StandardMMForceField
* Currently always return 0.0 since energies not calculated on gpu
* @return potential energy due to the OBC forces
*
*/
......@@ -427,5 +428,6 @@ double BrookCalcGBSAOBCForceFieldKernel::executeEnergy( const Stream& positions
// ---------------------------------------------------------------------------------------
return 0.0;
return (double) _brookGbsa->getEnergy( positions );
}
platforms/brook/src/BrookCalcStandardMMForceFieldKernel.cpp
View file @ 22ecd657
......@@ -38,6 +38,10 @@
#include "BrookCalcStandardMMForceFieldKernel.h"
#include "kforce.h"
#include "kinvmap_gather.h"
#include "ReferencePlatform.h"
#include "ReferenceFloatStreamImpl.h"
#include "VerletIntegrator.h"
#include "StandardMMForceField.h"
using namespace OpenMM;
using namespace std;
......@@ -59,9 +63,10 @@ BrookCalcStandardMMForceFieldKernel::BrookCalcStandardMMForceFieldKernel( std::s
// ---------------------------------------------------------------------------------------
_numberOfAtoms = 0;
_brookBonded = NULL;
_brookNonBonded = NULL;
_numberOfAtoms = 0;
_brookBonded = NULL;
_brookNonBonded = NULL;
_refForceField = NULL;
const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
if( brookPlatform.getLog() != NULL ){
......@@ -85,6 +90,10 @@ BrookCalcStandardMMForceFieldKernel::~BrookCalcStandardMMForceFieldKernel( ){
delete _brookBonded;
delete _brookNonBonded;
// deleted w/ kernel delete? If activated, program crashes
//delete _refForceField;
}
/**
......@@ -201,7 +210,94 @@ void BrookCalcStandardMMForceFieldKernel::initialize(
}
// ---------------------------------------------------------------------------------------
// used for calculating energy
/*
_referenceCalcStandardMMForceFieldKernel = new ReferenceCalcStandardMMForceFieldKernel( getName(), getPlatform() );
_referenceCalcStandardMMForceFieldKernel->initialize( bondIndices, bondParameters,
angleIndices, angleParameters,
periodicTorsionIndices, periodicTorsionParameters,
rbTorsionIndices, rbTorsionParameters,
bonded14Indices, lj14Scale, coulomb14Scale,
exclusions, nonbondedParameters,
nonbondedMethod, nonbondedCutoff, periodicBoxSize );
*/
_refForceField = new StandardMMForceField( _numberOfAtoms, bondIndices.size(), angleIndices.size(),
periodicTorsionIndices.size(), rbTorsionIndices.size() );
// bonds
vector<vector<int> >::const_iterator bondI_it = bondIndices.begin();
vector<vector<double> >::const_iterator bondP_it = bondParameters.begin();
for( unsigned int ii = 0; ii < bondIndices.size(); ii++, bondI_it++, bondP_it++ ){
const vector<int> bondInd = *bondI_it;
const vector<double> bondPrm = *bondP_it;
_refForceField->setBondParameters( ii, bondInd[0], bondInd[1], bondPrm[0], bondPrm[1] );
}
// angles
vector<vector<int> >::const_iterator angleI_it = angleIndices.begin();
vector<vector<double> >::const_iterator angleP_it = angleParameters.begin();
for( unsigned int ii = 0; ii < angleIndices.size(); ii++, angleI_it++, angleP_it++ ){
const vector<int> angleInd = *angleI_it;
const vector<double> anglePrm = *angleP_it;
_refForceField->setAngleParameters( ii, angleInd[0], angleInd[1], angleInd[2], anglePrm[0], anglePrm[1] );
}
// periodicTorsions
vector<vector<int> >::const_iterator periodicTorsionI_it = periodicTorsionIndices.begin();
vector<vector<double> >::const_iterator periodicTorsionP_it = periodicTorsionParameters.begin();
for( unsigned int ii = 0; ii < periodicTorsionIndices.size(); ii++, periodicTorsionI_it++, periodicTorsionP_it++ ){
const vector<int> periodicTorsionInd = *periodicTorsionI_it;
const vector<double> periodicTorsionPrm = *periodicTorsionP_it;
_refForceField->setPeriodicTorsionParameters( ii, periodicTorsionInd[0], periodicTorsionInd[1], periodicTorsionInd[2], periodicTorsionInd[3],
periodicTorsionPrm[2], periodicTorsionPrm[1], periodicTorsionPrm[0] );
/*
printf( "PeriodicTor: [%d %d %d %d] [%.5e %.5e %.5e]\n", periodicTorsionInd[0], periodicTorsionInd[1], periodicTorsionInd[2], periodicTorsionInd[3],
periodicTorsionPrm[2], periodicTorsionPrm[1], periodicTorsionPrm[0] ); fflush( stdout );
*/
}
// rbTorsions
vector<vector<int> >::const_iterator rbTorsionI_it = rbTorsionIndices.begin();
vector<vector<double> >::const_iterator rbTorsionP_it = rbTorsionParameters.begin();
for( unsigned int ii = 0; ii < rbTorsionIndices.size(); ii++, rbTorsionI_it++, rbTorsionP_it++ ){
const vector<int> rbTorsionInd = *rbTorsionI_it;
const vector<double> rbTorsionPrm = *rbTorsionP_it;
_refForceField->setRBTorsionParameters( ii, rbTorsionInd[0], rbTorsionInd[1], rbTorsionInd[2], rbTorsionInd[3],
rbTorsionPrm[0], rbTorsionPrm[1], rbTorsionPrm[2], rbTorsionPrm[3],
rbTorsionPrm[4], rbTorsionPrm[5] );
/*
printf( "RbTor: [%d %d %d %d] [%.5e %.5e %.5e %.5e %.5e %.5e]\n", rbTorsionInd[0], rbTorsionInd[1], rbTorsionInd[2], rbTorsionInd[3],
rbTorsionPrm[0], rbTorsionPrm[1], rbTorsionPrm[2], rbTorsionPrm[3], rbTorsionPrm[4], rbTorsionPrm[5] ); fflush( stdout );
*/
}
// nonbonded
for( unsigned int ii = 0; ii < nonbondedParameters.size(); ii++ ){
vector<double> nonbondedParameterVector
= nonbondedParameters[ii];
double c6 = nonbondedParameterVector[1];
double c12 = nonbondedParameterVector[2];
double charge = nonbondedParameterVector[0];
// int index, double charge, double radius, double depth
_refForceField->setAtomParameters( ii, charge, c6, c12 );
}
// ---------------------------------------------------------------------------------------
return;
}
/**
......@@ -408,7 +504,23 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
// gather forces
if( _brookBonded->getInverseMapStreamCount( K_Stream ) <= 4 ){
if( _brookBonded->getInverseMapStreamCount( I_Stream ) == 3 && _brookBonded->getInverseMapStreamCount( K_Stream ) == 3 ){
kinvmap_gather3_3( width,
inverseStreamMaps[I_Stream][0]->getBrookStream(),
inverseStreamMaps[I_Stream][1]->getBrookStream(),
inverseStreamMaps[I_Stream][2]->getBrookStream(),
bondedForceStreams[I_Stream]->getBrookStream(),
inverseStreamMaps[K_Stream][0]->getBrookStream(),
inverseStreamMaps[K_Stream][1]->getBrookStream(),
inverseStreamMaps[K_Stream][2]->getBrookStream(),
bondedForceStreams[K_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() );
} else if( _brookBonded->getInverseMapStreamCount( I_Stream ) == 3 && _brookBonded->getInverseMapStreamCount( K_Stream ) == 4 ){
kinvmap_gather3_4( width,
......@@ -425,7 +537,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
forceStream.getBrookStream(), forceStream.getBrookStream() );
} else if( _brookBonded->getInverseMapStreamCount( K_Stream ) == 5 ){
} else if( _brookBonded->getInverseMapStreamCount( I_Stream ) == 3 && _brookBonded->getInverseMapStreamCount( K_Stream ) == 5 ){
kinvmap_gather3_5( width,
inverseStreamMaps[I_Stream][0]->getBrookStream(),
......@@ -445,13 +557,31 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
// case not handled -- throw an exception
if( _brookBonded->getLog() ){
(void) fprintf( _brookBonded->getLog(), "%s nkmaps=%d -- not handled.", methodName.c_str(), _brookBonded->getInverseMapStreamCount( K_Stream ) );
(void) fprintf( _brookBonded->getLog(), "%s case: I-map=%d K-map=%d -- not handled.\n",
methodName.c_str(), _brookBonded->getInverseMapStreamCount( I_Stream ),
_brookBonded->getInverseMapStreamCount( K_Stream ) );
(void) fflush( _brookBonded->getLog() );
}
kinvmap_gather3_3( width,
inverseStreamMaps[I_Stream][0]->getBrookStream(),
inverseStreamMaps[I_Stream][1]->getBrookStream(),
inverseStreamMaps[I_Stream][2]->getBrookStream(),
bondedForceStreams[I_Stream]->getBrookStream(),
inverseStreamMaps[K_Stream][0]->getBrookStream(),
inverseStreamMaps[K_Stream][1]->getBrookStream(),
inverseStreamMaps[K_Stream][2]->getBrookStream(),
bondedForceStreams[K_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() );
/*
std::stringstream message;
message << methodName << "K-maps=" << _brookBonded->getInverseMapStreamCount( K_Stream ) << " not handled.";
message << methodName << "I-maps=" << _brookBonded->getInverseMapStreamCount( I_Stream ) << " and " <<
"K-maps=" << _brookBonded->getInverseMapStreamCount( K_Stream ) << " not handled.";
throw OpenMMException( message.str() );
*/
}
......@@ -465,19 +595,49 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
}
//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 );
kinvmap_gather5_2( width,
inverseStreamMaps[J_Stream][0]->getBrookStream(),
inverseStreamMaps[J_Stream][1]->getBrookStream(),
inverseStreamMaps[J_Stream][2]->getBrookStream(),
inverseStreamMaps[J_Stream][3]->getBrookStream(),
inverseStreamMaps[J_Stream][4]->getBrookStream(),
bondedForceStreams[J_Stream]->getBrookStream(),
inverseStreamMaps[L_Stream][0]->getBrookStream(),
inverseStreamMaps[L_Stream][1]->getBrookStream(),
bondedForceStreams[L_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() );
if( _brookBonded->getInverseMapStreamCount( J_Stream ) == 5 && _brookBonded->getInverseMapStreamCount( L_Stream ) == 2 ){
kinvmap_gather5_2( width,
inverseStreamMaps[J_Stream][0]->getBrookStream(),
inverseStreamMaps[J_Stream][1]->getBrookStream(),
inverseStreamMaps[J_Stream][2]->getBrookStream(),
inverseStreamMaps[J_Stream][3]->getBrookStream(),
inverseStreamMaps[J_Stream][4]->getBrookStream(),
bondedForceStreams[J_Stream]->getBrookStream(),
inverseStreamMaps[L_Stream][0]->getBrookStream(),
inverseStreamMaps[L_Stream][1]->getBrookStream(),
bondedForceStreams[L_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() );
} else {
// case not handled -- throw an exception
if( _brookBonded->getLog() ){
(void) fprintf( _brookBonded->getLog(), "%s case: J-map=%d L-map=%d -- not handled.\n",
methodName.c_str(), _brookBonded->getInverseMapStreamCount( J_Stream ),
_brookBonded->getInverseMapStreamCount( L_Stream ) );
(void) fflush( _brookBonded->getLog() );
}
kinvmap_gather5_2( width,
inverseStreamMaps[J_Stream][0]->getBrookStream(),
inverseStreamMaps[J_Stream][1]->getBrookStream(),
inverseStreamMaps[J_Stream][2]->getBrookStream(),
inverseStreamMaps[J_Stream][3]->getBrookStream(),
inverseStreamMaps[J_Stream][4]->getBrookStream(),
bondedForceStreams[J_Stream]->getBrookStream(),
inverseStreamMaps[L_Stream][0]->getBrookStream(),
inverseStreamMaps[L_Stream][1]->getBrookStream(),
bondedForceStreams[L_Stream]->getBrookStream(),
forceStream.getBrookStream(), forceStream.getBrookStream() );
/*
std::stringstream message;
message << methodName << "J-maps=" << _brookBonded->getInverseMapStreamCount( J_Stream ) << " and " <<
"L-maps=" << _brookBonded->getInverseMapStreamCount( L_Stream ) << " not handled.";
throw OpenMMException( message.str() );
*/
}
// diagnostics
......@@ -510,8 +670,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
*
*/
double BrookCalcStandardMMForceFieldKernel::executeEnergy( const Stream& positions ){
double BrookCalcStandardMMForceFieldKernel::executeEnergy( const Stream& atomPositions ){
// ---------------------------------------------------------------------------------------
......@@ -519,5 +678,32 @@ double BrookCalcStandardMMForceFieldKernel::executeEnergy( const Stream& positio
// ---------------------------------------------------------------------------------------
return 0.0;
const BrookStreamImpl& positionStreamC = dynamic_cast<const BrookStreamImpl&> (atomPositions.getImpl());
BrookStreamImpl& positionStream = const_cast<BrookStreamImpl&> (positionStreamC);
BrookOpenMMFloat* positionsF = (BrookOpenMMFloat*) positionStream.getData();
ReferencePlatform refPlatform;
System system( atomPositions.getSize(), 0);
VerletIntegrator integrator( 0.01 );
system.addForce( _refForceField );
OpenMMContext context(system, integrator, refPlatform);
vector<Vec3> positions( positionStream.getSize() );
int index = 0;
for( int ii = 0; ii < positionStream.getSize(); ii++, index += 3 ){
positions[ii] = Vec3( positionsF[index], positionsF[index+1], positionsF[index+2] );
}
context.setPositions(positions);
State state = context.getState( State::Energy );
double energy = state.getPotentialEnergy();
//delete forceField;
// (void) fprintf( stdout, "BrookCalcStandardMMForceFieldKernel::executeEnergy E=%.5e\n", energy ); fflush( stdout );
return energy;
}
platforms/brook/src/BrookCalcStandardMMForceFieldKernel.h
View file @ 22ecd657
......@@ -36,6 +36,7 @@
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include "BrookBonded.h"
#include "BrookNonBonded.h"
#include "StandardMMForceField.h"
namespace OpenMM {
......@@ -151,6 +152,10 @@ class BrookCalcStandardMMForceFieldKernel : public CalcStandardMMForceFieldKerne
BrookBonded* _brookBonded;
BrookNonBonded* _brookNonBonded;
// used to calculate energy
StandardMMForceField* _refForceField;
};
} // namespace OpenMM
......
platforms/brook/src/BrookCommon.cpp
View file @ 22ecd657
......@@ -486,4 +486,83 @@ void BrookCommon::getStreamDimensions( int streamSize, int *streamWidth, int *st
return;
}
/*
* Allocate array
*
* @param length length of array
* @param width width of array
*
* @return ptr to array
*
*/
RealOpenMM** BrookCommon::allocateRealArray( int length, int width ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookCommon::allocateRealArray";
// ---------------------------------------------------------------------------------------
RealOpenMM** array = new RealOpenMM*[length];
RealOpenMM* buffer = new RealOpenMM[length*width];
for( int ii = 0; ii < length; ii++ ){
array[ii] = buffer;
buffer += width;
}
return array;
}
/*
* Free array
*
* @param array array to be freed (assumed allocated using BrookCommon::allocateRealArray
*
* @return DefaultReturnValue
*
*/
int BrookCommon::disposeRealArray( RealOpenMM** array ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookCommon::disposeRealArray";
// ---------------------------------------------------------------------------------------
delete[] array[0];
delete[] array;
return DefaultReturnValue;
}
/*
* Copy 1D BrookOpenMMFloat* array to 2D array of RealOpenMM
*
* @param length length of array
* @param width width of array
* @param array1D array to copy
*
* @return ptr to array
*
*/
RealOpenMM** BrookCommon::copy1DArrayTo2DArray( int length, int width, BrookOpenMMFloat* array1D ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "BrookCommon::copy1DArrayTo2DArray";
// ---------------------------------------------------------------------------------------
RealOpenMM** array = allocateRealArray( length, width );
int index = 0;
for( int ii = 0; ii < length; ii++ ){
for( int jj = 0; jj < width; jj++ ){
array[ii][jj] = static_cast<RealOpenMM> (array1D[index++]);
}
}
return array;
}
platforms/brook/src/BrookCommon.h
View file @ 22ecd657
......@@ -270,6 +270,42 @@ class BrookCommon {
*/
static void getStreamDimensions( int streamSize, int *streamWidth, int *streamHeight );
/*
* Allocate array
*
* @param length length of array
* @param width width of array
*
* @return ptr to array
*
*/
RealOpenMM** allocateRealArray( int length, int width ) const;
/*
* Free array
*
* @param array array to be freed (assumed allocated using BrookCommon::allocateRealArray
*
* @return DefaultReturnValue
*
*/
int disposeRealArray( RealOpenMM** array ) const;
/*
* Copy 1D BrookOpenMMFloat* array to 2D array of RealOpenMM
*
* @param length length of array
* @param width width of array
* @param array1D array to copy
*
* @return ptr to array
*
*/
RealOpenMM** copy1DArrayTo2DArray( int length, int width, BrookOpenMMFloat* array1D ) const;
protected:
......
platforms/brook/src/BrookFloatStreamInternal.cpp
View file @ 22ecd657
......@@ -173,8 +173,6 @@ BrookFloatStreamInternal::~BrookFloatStreamInternal( ){
// ---------------------------------------------------------------------------------------
//delete _aStream;
//printf( "%s %s data=%u stream=%d [%d %d] width=%d\n", methodName.c_str(), getName().c_str(), (unsigned int) _data, getStreamSize(), _streamHeight, _streamWidth, _width );
//fflush( stdout );
......
platforms/brook/src/BrookGbsa.cpp
View file @ 22ecd657
......@@ -68,7 +68,7 @@ BrookGbsa::BrookGbsa( ){
_includeAce = 1;
_solventDielectric = 78.3;
_soluteDielectric = 1.0;
_dielectricOffset = 0.09;
_dielectricOffset = 0.009;
for( int ii = 0; ii < LastStreamIndex; ii++ ){
_gbsaStreams[ii] = NULL;
......@@ -80,6 +80,8 @@ BrookGbsa::BrookGbsa( ){
_bornRadiiInitialized = 0;
_cpuObc = NULL;
_charges = NULL;
}
/**
......@@ -105,6 +107,8 @@ BrookGbsa::~BrookGbsa( ){
delete _cpuObc;
delete[] _charges;
}
/**
......@@ -443,7 +447,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::calculateBornRadii";
static const int PrintOn = 1;
static const int PrintOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -487,7 +491,7 @@ int BrookGbsa::calculateBornRadii( const Stream& positions ){
// diagnostics
if( PrintOn ){
if( PrintOn && getLog() ){
(void) fprintf( getLog(), "\n%s: atms=%d\n", methodName.c_str(), numberOfAtoms );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
......@@ -534,6 +538,26 @@ int BrookGbsa::initializeStreamSizes( int numberOfAtoms, const Platform& platfor
_gbsaAtomStreamWidth = getAtomStreamWidth( platform );
_gbsaAtomStreamHeight = getAtomStreamHeight( platform );
int innerUnroll = getInnerLoopUnroll();
if( innerUnroll < 1 ){
std::stringstream message;
message << methodName << " innerUnrolls=" << innerUnroll << " is less than 1.";
throw OpenMMException( message.str() );
return ErrorReturnValue;
}
if( _partialForceStreamWidth < 1 ){
std::stringstream message;
message << methodName << " partial force stream width=" << _partialForceStreamWidth << " is less than 1.";
throw OpenMMException( message.str() );
return ErrorReturnValue;
}
_partialForceStreamSize = _gbsaAtomStreamSize*getDuplicationFactor()/innerUnroll;
_partialForceStreamHeight = _partialForceStreamSize/_partialForceStreamWidth;
_partialForceStreamHeight += ( (_partialForceStreamSize % _partialForceStreamWidth) ? 1 : 0);
_partialForceStreamSize = _partialForceStreamHeight*_partialForceStreamWidth;
return DefaultReturnValue;
}
......@@ -607,7 +631,7 @@ int BrookGbsa::initializeStreams( const Platform& platform ){
std::stringstream name;
name << partialForceStream << ii;
_gbsaForceStreams[ii] = new BrookFloatStreamInternal( name.str(), getPartialForceStreamSize(),
getPartialForceStreamWidth(), BrookStreamInternal::Float4, dangleValue );
getPartialForceStreamWidth(), BrookStreamInternal::Float4, dangleValue );
}
return DefaultReturnValue;
......@@ -649,8 +673,13 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfAtomParam
initializeStreamSizes( numberOfAtoms, platform );
initializeStreams( platform );
BrookOpenMMFloat* radiiAndCharge = new BrookOpenMMFloat[getNumberOfAtoms()*2];
BrookOpenMMFloat* scaledRadiiAndOffset = new BrookOpenMMFloat[getNumberOfAtoms()*2];
int atomStreamSize = getGbsaAtomStreamSize();
BrookOpenMMFloat* radiiAndCharge = new BrookOpenMMFloat[atomStreamSize*2];
BrookOpenMMFloat* scaledRadiiAndOffset = new BrookOpenMMFloat[atomStreamSize*2];
memset( radiiAndCharge, 0, atomStreamSize*2*sizeof( BrookOpenMMFloat ) );
memset( scaledRadiiAndOffset, 0, atomStreamSize*2*sizeof( BrookOpenMMFloat ) );
_charges = new RealOpenMM[atomStreamSize];
// used by CpuObc to calculate initial Born radii
......@@ -688,12 +717,16 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfAtomParam
atomicRadii[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);
scaledRadiiAndOffset[streamIndex] = static_cast<BrookOpenMMFloat> (radius*scalingFactor);
scaledRadiiAndOffset[streamIndex+1] = static_cast<BrookOpenMMFloat> (radius - dielectricOffset);
scaledRadiiAndOffset[streamIndex] = static_cast<BrookOpenMMFloat> (scaledRadiiAndOffset[streamIndex+1]*scalingFactor);
// scaledRadiiAndOffset[streamIndex] = static_cast<BrookOpenMMFloat> (radius - dielectricOffset);
// scaledRadiiAndOffset[streamIndex+1] = static_cast<BrookOpenMMFloat> (scaledRadiiAndOffset[streamIndex]*scalingFactor);
}
......@@ -714,10 +747,10 @@ int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfAtomParam
delete[] radiiAndCharge;
delete[] scaledRadiiAndOffset;
// setup for Born radii
// setup for Born radii calculation
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
obcParameters->setAtomicRadii( atomicRadii, SimTKOpenMMCommon::MdUnits);
obcParameters->setAtomicRadii( atomicRadii);
obcParameters->setScaledRadiusFactors( scaleFactors );
obcParameters->setSolventDielectric( static_cast<RealOpenMM>(solventDielectric) );
......@@ -868,3 +901,53 @@ std::string BrookGbsa::getContentsString( int level ) const {
return message.str();
}
/*
* Calculate OBC energy
*
* @param atomPositions atom positions
* @return energy
*
* @throw OpenMMException if _cpuObc or charges are not set
*
* */
double BrookGbsa::getEnergy( const Stream& atomPositions ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "BrookGbsa::getEnergy";
// ---------------------------------------------------------------------------------------
// validate initialization
if( _cpuObc == NULL ){
std::stringstream message;
message << methodName << " _cpuObc not set.";
throw OpenMMException( message.str() );
return ErrorReturnValue;
}
if( _charges == NULL ){
std::stringstream message;
message << methodName << " charges not set.";
throw OpenMMException( message.str() );
return ErrorReturnValue;
}
const BrookStreamImpl& positionStreamC = dynamic_cast<const BrookStreamImpl&> (atomPositions.getImpl());
BrookStreamImpl& positionStream = const_cast<BrookStreamImpl&> (positionStreamC);
BrookOpenMMFloat* positionsF = (BrookOpenMMFloat*) positionStream.getData();
RealOpenMM** positions = copy1DArrayTo2DArray( positionStream.getSize(), 3, positionsF );
RealOpenMM** forces = allocateRealArray( positionStream.getSize(), 3 );
_cpuObc->computeImplicitSolventForces( positions, _charges, forces, 1 );
disposeRealArray( forces );
disposeRealArray( positions );
return _cpuObc->getEnergy();
}
platforms/brook/src/BrookGbsa.h
View file @ 22ecd657
......@@ -337,6 +337,19 @@ class BrookGbsa : public BrookCommon {
std::string getContentsString( int level = 0 ) const;
/*
* Calculate energy
*
* @param atomPositions atom positions
* @return energy
*
* @throw OpenMMException if _cpuObc or charges are not set
*
* */
double getEnergy( const Stream& atomPositions );
private:
// fixed number of force streams
......@@ -401,6 +414,10 @@ class BrookGbsa : public BrookCommon {
int _bornRadiiInitialized;
// atom charges
RealOpenMM* _charges;
// CpuObc reference
CpuObc* _cpuObc;
......
platforms/brook/src/BrookStreamInternal.cpp
View file @ 22ecd657
......@@ -60,6 +60,9 @@ BrookStreamInternal::~BrookStreamInternal( ){
//static const std::string methodName = "BrookStreamInternal::~BrookStreamInternal";
// ---------------------------------------------------------------------------------------
// _aStream is not ptr
// delete _aStream;
}
......
platforms/brook/src/gpu/kinvmap_gather.br
View file @ 22ecd657
......@@ -224,6 +224,33 @@ kernel void kinvmap_gather6(
}
//Takes three + four inverse maps
kernel void kinvmap_gather3_3(
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<>,
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 );
}
//Takes three + four inverse maps
kernel void kinvmap_gather3_4(
float strwidth, //stream width of the dihedral forces
......
platforms/brook/src/gpu/kinvmap_gather.h
View file @ 22ecd657
......@@ -50,6 +50,18 @@ void kinvmap_gather6 (const float strwidth,
::brook::stream outforce);
void kinvmap_gather3_3 (const float strwidth,
::brook::stream invmap3_1,
::brook::stream invmap3_2,
::brook::stream invmap3_3,
::brook::stream forces3,
::brook::stream invmap4_1,
::brook::stream invmap4_2,
::brook::stream invmap4_3,
::brook::stream forces4,
::brook::stream inforce,
::brook::stream outforce);
void kinvmap_gather3_4 (const float strwidth,
::brook::stream invmap3_1,
::brook::stream invmap3_2,
......
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