Skip to content

GitLab

Commit 326627a8 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Cleaned up code

parent a587c652
Hide whitespace changes
Inline Side-by-side
Showing with 1444 additions and 5224 deletions
platforms/reference/src/ReferenceKernels.cpp
View file @ 326627a8
......@@ -847,12 +847,13 @@ double ReferenceCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool inclu
vector<RealVec>& forceData = extractForces(context);
if (isPeriodic)
obc->getObcParameters()->setPeriodic(extractBoxSize(context));
obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
return obc->getEnergy();
return obc->computeBornEnergyForces(posData, charges, forceData);
}
ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
if (gbvi) {
GBVIParameters * gBVIParameters = gbvi->getGBVIParameters();
delete gBVIParameters;
delete gbvi;
}
}
......@@ -887,18 +888,21 @@ void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIFo
}
double ReferenceCalcGBVIForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealOpenMM> bornRadii(context.getSystem().getNumParticles());
if (isPeriodic)
gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
gbvi->computeBornRadii(posData, bornRadii );
RealOpenMM energy;
if (includeForces) {
vector<RealVec>& forceData = extractForces(context);
gbvi->computeBornForces(bornRadii, posData, &charges[0], forceData);
gbvi->computeBornForces(posData, charges, forceData);
energy = 0.0;
}
if( includeEnergy ){
energy = gbvi->computeBornEnergy(posData, charges);
}
RealOpenMM energy = 0.0;
if (includeEnergy)
energy = gbvi->computeBornEnergy(bornRadii ,posData, &charges[0]);
return static_cast<double>(energy);
}
......
platforms/reference/src/gbsa/CpuGBVI.cpp
View file @ 326627a8
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
* Contributors: Pande Group
*
......@@ -23,899 +22,577 @@
*/
#include <string.h>
#include <math.h>
#include <sstream>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "CpuGBVI.h"
#include "../SimTKReference/ReferenceForce.h"
#include <math.h>
#include "CpuGBVI.h"
using namespace std;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
CpuGBVI constructor
gbviParameters gbviParameters object
--------------------------------------------------------------------------------------- */
CpuGBVI constructor
CpuGBVI::CpuGBVI( ImplicitSolventParameters* gbviParameters ) : CpuImplicitSolvent( gbviParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVI::CpuGBVI";
// ---------------------------------------------------------------------------------------
_initializeGBVIDataMembers( );
_gbviParameters = static_cast<GBVIParameters*> (gbviParameters);
gbviParameters gbviParameters object
--------------------------------------------------------------------------------------- */
CpuGBVI::CpuGBVI( GBVIParameters* gbviParameters ) : _gbviParameters(gbviParameters) {
_switchDeriviative.resize( gbviParameters->getNumberOfAtoms() );
}
/**---------------------------------------------------------------------------------------
CpuGBVI destructor
CpuGBVI destructor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
CpuGBVI::~CpuGBVI( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVI::~CpuGBVI";
// ---------------------------------------------------------------------------------------
//if( _gbviParameters != NULL ){
// delete _gbviParameters;
//}
}
/**---------------------------------------------------------------------------------------
Initialize data members
--------------------------------------------------------------------------------------- */
void CpuGBVI::_initializeGBVIDataMembers( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVI::initializeDataMembers";
// ---------------------------------------------------------------------------------------
_gbviParameters = NULL;
}
/**---------------------------------------------------------------------------------------
Get GBVIParameters reference
Get GBVIParameters reference
@return GBVIParameters reference
@return GBVIParameters reference
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
GBVIParameters* CpuGBVI::getGBVIParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVI::getGBVIParameters";
// ---------------------------------------------------------------------------------------
return _gbviParameters;
return _gbviParameters;
}
/**---------------------------------------------------------------------------------------
Set GBVIParameters reference
Set GBVIParameters reference
@param GBVIParameters reference
@param GBVIParameters reference
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void CpuGBVI::setGBVIParameters( GBVIParameters* gbviParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVI::setGBVIParameters";
// ---------------------------------------------------------------------------------------
_gbviParameters = gbviParameters;
if( _switchDeriviative.size() <= 0 ){
_switchDeriviative.resize(_gbviParameters->getNumberOfAtoms());
}
void CpuGBVI::setGBVIParameters( GBVIParameters* gbviParameters ){
_gbviParameters = gbviParameters;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
std::vector<RealOpenMM>& CpuGBVI::getSwitchDeriviative( void ){
@return array
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
// static const char* methodName = "\nCpuGBVI::getSwitchDeriviative";
// ---------------------------------------------------------------------------------------
if( _switchDeriviative.size() <= 0 ){
_switchDeriviative.resize(_gbviParameters->getNumberOfAtoms());
}
return _switchDeriviative;
}
/**---------------------------------------------------------------------------------------
Return switching function derivative
@return array
--------------------------------------------------------------------------------------- */
/*
std::vector<RealOpenMM> CpuGBVI::getSwitchDeriviativeConst( void ) const {
return _switchDeriviative;
RealOpenMMVector& CpuGBVI::getSwitchDeriviative( void ){
return _switchDeriviative;
}
*/
/**---------------------------------------------------------------------------------------
Compute quintic spline value and associated derviative
@param x value to compute spline at
@param rl lower cutoff value
@param ru upper cutoff value
@param outValue value of spline at x
@param outDerivative value of derivative of spline at x
Compute quintic spline value and associated derviative
--------------------------------------------------------------------------------------- */
@param x value to compute spline at
@param rl lower cutoff value
@param ru upper cutoff value
@param outValue value of spline at x
@param outDerivative value of derivative of spline at x
#define GBVIDebug 0
--------------------------------------------------------------------------------------- */
void CpuGBVI::quinticSpline( RealOpenMM x, RealOpenMM rl, RealOpenMM ru,
RealOpenMM* outValue, RealOpenMM* outDerivative ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM minusSix = (RealOpenMM) -6.0;
static const RealOpenMM minusTen = (RealOpenMM) -10.0;
static const RealOpenMM minusThirty = (RealOpenMM) -30.0;
static const RealOpenMM fifteen = (RealOpenMM) 15.0;
static const RealOpenMM sixty = (RealOpenMM) 60.0;
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::quinticSpline";
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM minusSix = static_cast<RealOpenMM>( -6.0 );
static const RealOpenMM minusTen = static_cast<RealOpenMM>( -10.0 );
static const RealOpenMM minusThirty = static_cast<RealOpenMM>( -30.0 );
static const RealOpenMM fifteen = static_cast<RealOpenMM>( 15.0 );
static const RealOpenMM sixty = static_cast<RealOpenMM>( 60.0 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
RealOpenMM numerator = x - rl;
RealOpenMM denominator = ru - rl;
RealOpenMM ratio = numerator/denominator;
RealOpenMM ratio2 = ratio*ratio;
RealOpenMM ratio3 = ratio2*ratio;
RealOpenMM numerator = x - rl;
RealOpenMM denominator = ru - rl;
RealOpenMM ratio = numerator/denominator;
RealOpenMM ratio2 = ratio*ratio;
RealOpenMM ratio3 = ratio2*ratio;
*outValue = one + ratio3*(minusTen + fifteen*ratio + minusSix*ratio2);
*outDerivative = ratio2*(minusThirty + sixty*ratio + minusThirty*ratio2)/denominator;
*outValue = one + ratio3*(minusTen + fifteen*ratio + minusSix*ratio2);
*outDerivative = ratio2*(minusThirty + sixty*ratio + minusThirty*ratio2)/denominator;
}
/**---------------------------------------------------------------------------------------
Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
and quintic splice switching function
Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
and quintic splice switching function
@param atomicRadius3 atomic radius cubed
@param bornSum Born sum (volume integral)
@param gbviParameters Gbvi parameters (parameters used in spline
QuinticLowerLimitFactor & QuinticUpperBornRadiusLimit)
@param bornRadius output Born radius
@param switchDeriviative output switching function deriviative
@param atomicRadius3 atomic radius cubed
@param bornSum Born sum (volume integral)
@param gbviParameters Gbvi parameters (parameters used in spline
QuinticLowerLimitFactor & QuinticUpperBornRadiusLimit)
@param bornRadius output Born radius
@param switchDeriviative output switching function deriviative
--------------------------------------------------------------------------------------- */
#define GBVIDebug 0
--------------------------------------------------------------------------------------- */
void CpuGBVI::computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadius3, RealOpenMM bornSum,
GBVIParameters* gbviParameters,
RealOpenMM* bornRadius, RealOpenMM* switchDeriviative ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM minusOne = (RealOpenMM) -1.0;
static const RealOpenMM minusThree = (RealOpenMM) -3.0;
static const RealOpenMM oneEighth = (RealOpenMM) 0.125;
static const RealOpenMM minusOneThird = (RealOpenMM) (-1.0/3.0);
static const RealOpenMM three = (RealOpenMM) 3.0;
// ---------------------------------------------------------------------------------------
static const char* methodName = "CpuGBVI::computeBornRadiiUsingQuinticSpline";
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM minusOne = static_cast<RealOpenMM>( -1.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
static const RealOpenMM oneEighth = static_cast<RealOpenMM>( 0.125 );
static const RealOpenMM minusOneThird = static_cast<RealOpenMM>( (-1.0/3.0) );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
#if( GBVIDebug == 1 )
FILE* logFile = stderr;
#endif
// R = [ S(V)*(A - V) ]**(-1/3)
// R = [ S(V)*(A - V) ]**(-1/3)
// S(V) = 1 V < L
// S(V) = qSpline + U/(A-V) L < V < A
// S(V) = U/(A-V) U < V
// S(V) = 1 V < L
// S(V) = qSpline + U/(A-V) L < V < A
// S(V) = U/(A-V) U < V
// dR/dr = (-1/3)*[ S(V)*(A - V) ]**(-4/3)*[ d{ S(V)*(A-V) }/dr
// dR/dr = (-1/3)*[ S(V)*(A - V) ]**(-4/3)*[ d{ S(V)*(A-V) }/dr
// d{ S(V)*(A-V) }/dr = (dV/dr)*[ (A-V)*dS/dV - S(V) ]
// d{ S(V)*(A-V) }/dr = (dV/dr)*[ (A-V)*dS/dV - S(V) ]
// (A - V)*dS/dV - S(V) = 0 - 1 V < L
// (A - V)*dS/dV - S(V) = 0 - 1 V < L
// (A - V)*dS/dV - S(V) = (A-V)*d(qSpline) + (A-V)*U/(A-V)**2 - qSpline - U/(A-V)
// (A - V)*dS/dV - S(V) = (A-V)*d(qSpline) + (A-V)*U/(A-V)**2 - qSpline - U/(A-V)
// = (A-V)*d(qSpline) - qSpline L < V < A**(-3)
// (A - V)*dS/dV - S(V) = (A-V)*U*/(A-V)**2 - U/(A-V) = 0 U < V
RealOpenMM splineL = gbviParameters->getQuinticLowerLimitFactor()*atomicRadius3;
RealOpenMM sum;
if( bornSum > splineL ){
if( bornSum < atomicRadius3 ){
RealOpenMM splineValue, splineDerivative;
quinticSpline( bornSum, splineL, atomicRadius3, &splineValue, &splineDerivative );
sum = (atomicRadius3 - bornSum)*splineValue + gbviParameters->getQuinticUpperBornRadiusLimit();
*switchDeriviative = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
#if( GBVIDebug == 1 )
(void) fprintf( logFile, " Qv=%14.6e splnDrvtv=%14.6e spline[%10.3e %10.3e] ", splineValue, splineDerivative,
splineL, gbviParameters->getQuinticUpperBornRadiusLimit() );
#endif
} else {
sum = gbviParameters->getQuinticUpperBornRadiusLimit();
*switchDeriviative = zero;
}
} else {
sum = atomicRadius3 - bornSum;
*switchDeriviative = one;
}
*bornRadius = POW( sum, minusOneThird );
// (A - V)*dS/dV - S(V) = (A-V)*U*/(A-V)**2 - U/(A-V) = 0 U < V
RealOpenMM splineL = gbviParameters->getQuinticLowerLimitFactor()*atomicRadius3;
RealOpenMM sum;
if( bornSum > splineL ){
if( bornSum < atomicRadius3 ){
RealOpenMM splineValue, splineDerivative;
quinticSpline( bornSum, splineL, atomicRadius3, &splineValue, &splineDerivative );
sum = (atomicRadius3 - bornSum)*splineValue + gbviParameters->getQuinticUpperBornRadiusLimit();
*switchDeriviative = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
} else {
sum = gbviParameters->getQuinticUpperBornRadiusLimit();
*switchDeriviative = zero;
}
} else {
sum = atomicRadius3 - bornSum;
*switchDeriviative = one;
}
*bornRadius = POW( sum, minusOneThird );
}
/**---------------------------------------------------------------------------------------
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param chain not used here
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param chain not used here
#define GBVIDebug 0
--------------------------------------------------------------------------------------- */
void CpuGBVI::computeBornRadii( vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii ){
void CpuGBVI::computeBornRadii( const vector<RealVec>& atomCoordinates, RealOpenMMVector& bornRadii ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM minusThree = (RealOpenMM) -3.0;
static const RealOpenMM oneEighth = (RealOpenMM) 0.125;
static const RealOpenMM minusOneThird = (RealOpenMM) (-1.0/3.0);
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
static const RealOpenMM oneEighth = static_cast<RealOpenMM>( 0.125 );
static const RealOpenMM minusOneThird = static_cast<RealOpenMM>( (-1.0/3.0) );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const char* methodName = "CpuGBVI::computeBornRadii";
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
GBVIParameters* gbviParameters = getGBVIParameters();
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& scaledRadii = gbviParameters->getScaledRadii();
GBVIParameters* gbviParameters = getGBVIParameters();
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMMVector& switchDeriviatives = getSwitchDeriviative();
std::vector<RealOpenMM>& switchDeriviatives = getSwitchDeriviative();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// calculate Born radii
#if( GBVIDebug == 1 )
FILE* logFile = stderr;
(void) fprintf( logFile, "\n%s scalingMethdd=%d\n", methodName, _gbviParameters->getBornRadiusScalingMethod() );
(void) fprintf( logFile, "Lower=%14.6e Upper=%10.3e\n", gbviParameters->getQuinticLowerLimitFactor(), gbviParameters->getQuinticUpperBornRadiusLimit() );
#endif
// calculate Born radii
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM sum = zero;
// sum over volumes
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
continue;
sum += CpuGBVI::getVolume( r, radiusI, scaledRadii[atomJ] );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM sum = zero;
// sum over volumes
#if( GBVIDebug == 1 )
if( atomI == 1568 || atomJ == 1568 ){
(void) fprintf( logFile, "%d addJ=%d scR=%14.6e %14.6e sum=%14.6e rI=%14.6e r=%14.6e S-R=%14.6e\n",
atomI, atomJ, scaledRadii[atomJ], getVolume( r, radiusI, scaledRadii[atomJ] ), sum,
radiusI, r, (scaledRadii[atomJ]-radiusI) );
}
#endif
}
}
#if( GBVIDebug == 1 )
(void) fprintf( logFile, "%d Born radius sum=%14.6e %14.6e %14.6e q=%d ",
atomI, sum, POW( radiusI, minusThree ), (POW( radiusI, minusThree ) - sum), GBVIParameters::QuinticSpline );
#endif
RealOpenMM atomicRadius3 = POW( radiusI, minusThree );
if( _gbviParameters->getBornRadiusScalingMethod() == GBVIParameters::QuinticSpline ){
sum = atomicRadius3 - sum;
bornRadii[atomI] = POW( sum, minusOneThird );
switchDeriviatives[atomI] = one;
} else {
RealOpenMM bornRadius, switchDeriviative;
computeBornRadiiUsingQuinticSpline( atomicRadius3, sum, gbviParameters,
&bornRadius, &switchDeriviative );
bornRadii[atomI] = bornRadius;
switchDeriviatives[atomI] = switchDeriviative;
}
#if( GBVIDebug == 1 )
(void) fprintf( logFile, " br=%14.6e swDeriv=%14.6e\n", atomI, bornRadii[atomI], switchDeriviatives[atomI] );
#endif
}
#if( GBVIDebug == 2 )
std::string fileName = "br.txt";
FILE* brFile = NULL;
#ifdef WIN32
fopen_s( &brFile, fileName.c_str(), "w" );
#else
brFile = fopen( fileName.c_str(), "w" );
#endif
(void) fprintf( brFile, "%5d\n", numberOfAtoms );
for( unsigned int ii = 0; ii < numberOfAtoms; ii++ ){
(void) fprintf( brFile, "%6u %14.6e %9.4f %9.4f %14.6e %14.6e %14.6e\n",
ii, bornRadii[ii], atomicRadii[ii], scaledRadii[ii],
atomCoordinates[ii][0], atomCoordinates[ii][1], atomCoordinates[ii][2] );
}
(void) fclose( brFile );
#endif
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
continue;
sum += CpuGBVI::getVolume( r, radiusI, scaledRadii[atomJ] );
}
}
RealOpenMM atomicRadius3 = POW( radiusI, minusThree );
if( _gbviParameters->getBornRadiusScalingMethod() == GBVIParameters::QuinticSpline ){
sum = atomicRadius3 - sum;
bornRadii[atomI] = POW( sum, minusOneThird );
switchDeriviatives[atomI] = one;
} else {
RealOpenMM bornRadius, switchDeriviative;
computeBornRadiiUsingQuinticSpline( atomicRadius3, sum, gbviParameters,
&bornRadius, &switchDeriviative );
bornRadii[atomI] = bornRadius;
switchDeriviatives[atomI] = switchDeriviative;
}
}
}
#undef GBVIDebug
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
@return volume
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVI::getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::getVolume";
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM minusThree = (RealOpenMM) -3.0;
// ---------------------------------------------------------------------------------------
RealOpenMM diff = (S - R);
if( FABS( diff ) < r ){
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
return (CpuGBVI::getL( r, (r + S), S ) -
CpuGBVI::getL( r, lowerBound, S ));
RealOpenMM diff = (S - R);
if( FABS( diff ) < r ){
} else if( r <= diff ){
RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
return (CpuGBVI::getL( r, (r + S), S ) -
CpuGBVI::getL( r, lowerBound, S ));
} else if( r <= diff ){
return CpuGBVI::getL( r, (r + S), S ) -
CpuGBVI::getL( r, (r - S), S ) +
POW( R, minusThree );
return CpuGBVI::getL( r, (r + S), S ) -
CpuGBVI::getL( r, (r - S), S ) +
POW( R, minusThree );
} else {
return zero;
}
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
Get L (used in analytical solution for volume integrals)
Get L (used in analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVI::getL( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::getL";
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM threeHalves = (RealOpenMM) 1.5;
static const RealOpenMM third = (RealOpenMM) (1.0/3.0);
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125;
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM threeHalves = static_cast<RealOpenMM>( 1.5 );
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0) );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
RealOpenMM diff2 = (r + S)*(r - S);
return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
}
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVI::dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVI::dL_dr";
// static const char* methodName = "\nCpuGBVI::dL_dr";
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM threeHalves = (RealOpenMM) 1.5;
static const RealOpenMM threeEights = (RealOpenMM) 0.375;
static const RealOpenMM third = (RealOpenMM) (1.0/3.0);
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125;
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM threeHalves = static_cast<RealOpenMM>( 1.5 );
static const RealOpenMM threeEights = static_cast<RealOpenMM>( 0.375 );
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0) );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM rInv2 = rInv*rInv;
RealOpenMM rInv = one/r;
RealOpenMM rInv2 = rInv*rInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
RealOpenMM diff2 = (r + S)*(r - S);
return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
}
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVI::dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::dL_dx";
// static const char* methodName = "CpuGBVI::dL_dx";
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM threeHalvesM = (RealOpenMM) -1.5;
static const RealOpenMM third = (RealOpenMM) (1.0/3.0);
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM threeHalvesM = static_cast<RealOpenMM>( -1.5 );
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0) );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff = (r + S)*(r - S);
RealOpenMM diff = (r + S)*(r - S);
return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
}
/**---------------------------------------------------------------------------------------
Sgb function
Sgb function
@param t r*r*G_i*G_j
@param t r*r*G_i*G_j
@return Sgb (top of p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
@return Sgb (top of p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVI::Sgb( RealOpenMM t ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::Sgb";
// static const char* methodName = "CpuGBVI::Sgb";
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
return ( (t != zero) ? one/SQRT( (one + (fourth*EXP( -t ))/t) ) : zero);
return ( (t != zero) ? one/SQRT( (one + (fourth*EXP( -t ))/t) ) : zero);
}
#define GBVIDebug 0
/**---------------------------------------------------------------------------------------
Get GB/VI energy
Get GB/VI energy
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@return energy
@return energy
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVI::computeBornEnergy( const vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges ){
RealOpenMM CpuGBVI::computeBornEnergy( const vector<RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
static const char* methodName = "CpuGBVI::computeBornEnergy";
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM two = (RealOpenMM) 2.0;
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM four = (RealOpenMM) 4.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
const GBVIParameters* gbviParameters = getGBVIParameters();
const RealOpenMM preFactor = gbviParameters->getElectricConstant();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& gammaParameters = gbviParameters->getGammaParameters();
const GBVIParameters* gbviParameters = getGBVIParameters();
const RealOpenMM preFactor = gbviParameters->getElectricConstant();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* gammaParameters = gbviParameters->getGammaParameters();
// compute Born radii
#if( GBVIDebug == 1 )
FILE* logFile = stderr;
(void) fprintf( logFile, "\n%s\n", methodName );
(void) fflush( logFile );
#endif
RealOpenMMVector bornRadii( numberOfAtoms );
computeBornRadii( atomCoordinates, bornRadii );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
// to minimze roundoff error sum cavityEnergy separately since in general much
// smaller than other contributions
// Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
// to minimze roundoff error sum cavityEnergy separately since in general much
// smaller than other contributions
RealOpenMM energy = zero;
RealOpenMM cavityEnergy = zero;
RealOpenMM energy = zero;
RealOpenMM cavityEnergy = zero;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = partialCharges[atomI];
// self-energy term
RealOpenMM atomIEnergy = half*partialChargeI/bornRadii[atomI];
// cavity term
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
cavityEnergy += gammaParameters[atomI]*ratio*ratio*ratio;
/*
RealOpenMM e1 = partialChargeI*partialCharges[atomI]/bornRadii[atomI];
RealOpenMM e2 = gammaParameters[atomI]*ratio*ratio*ratio;
(void) fprintf( stderr, "E %d self=%.4e gamma=%.4e e=%.4e\n", atomI, e1, e2, energy );
*/
for( int atomJ = atomI + 1; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM t = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);
atomIEnergy += partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
/*
RealOpenMM e3 = -partialChargeI*partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
(void) fprintf( stderr, "E %d %d e3=%.4e r2=%4e t=%.3e sgb=%.4e e=%.5e\n", atomI, atomJ, e3, r2, t, Sgb( t ), energy );
*/
}
energy += two*partialChargeI*atomIEnergy;
}
energy *= preFactor;
energy -= cavityEnergy;
#if( GBVIDebug == 1 )
(void) fprintf( logFile, "ElectricConstant=%.4e Tau=%.4e e=%.5e eOut=%.5e\n", preFactor, gbviParameters->getTau(), energy, gbviParameters->getTau()*energy );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "bR %d bR=%16.8e\n", atomI, bornRadii[atomI] );
}
(void) fflush( logFile );
#endif
RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
return (conversion*energy);
RealOpenMM partialChargeI = partialCharges[atomI];
}
#undef GBVIDebug
// self-energy term
RealOpenMM atomIEnergy = half*partialChargeI/bornRadii[atomI];
// cavity term
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
cavityEnergy += gammaParameters[atomI]*ratio*ratio*ratio;
for( int atomJ = atomI + 1; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM t = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);
atomIEnergy += partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
}
energy += two*partialChargeI*atomIEnergy;
}
energy *= preFactor;
energy -= cavityEnergy;
#define GBVIDebug 0
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
return (conversion*energy);
}
/**---------------------------------------------------------------------------------------
Get GB/VI forces
Get GB/VI forces
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void CpuGBVI::computeBornForces( const std::vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& inputForces){
void CpuGBVI::computeBornForces( std::vector<RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges,
std::vector<OpenMM::RealVec>& inputForces){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
static const char* methodName = "CpuGBVI::computeBornEnergyForces";
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM oneThird = static_cast<RealOpenMM>( (1.0/3.0) );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM two = (RealOpenMM) 2.0;
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM four = (RealOpenMM) 4.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM oneThird = (RealOpenMM) (1.0/3.0);
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
const GBVIParameters* gbviParameters = getGBVIParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& gammaParameters = gbviParameters->getGammaParameters();
#if( GBVIDebug == 1 )
FILE* logFile = stderr;
(void) fprintf( logFile, "\n%s\n", methodName );
(void) fflush( logFile );
#endif
// ---------------------------------------------------------------------------------------
const GBVIParameters* gbviParameters = getGBVIParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* gammaParameters = gbviParameters->getGammaParameters();
// constants
// ---------------------------------------------------------------------------------------
const RealOpenMM preFactor = two*gbviParameters->getElectricConstant();
// constants
// ---------------------------------------------------------------------------------------
const RealOpenMM preFactor = two*gbviParameters->getElectricConstant();
// compute Born radii
// ---------------------------------------------------------------------------------------
RealOpenMMVector bornRadii( numberOfAtoms );
computeBornRadii( atomCoordinates, bornRadii );
// set energy/forces to zero
// set energy/forces to zero
RealOpenMM** forces = new RealOpenMM*[numberOfAtoms];
RealOpenMM* block = new RealOpenMM[numberOfAtoms*3];
memset( block, 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
RealOpenMM* blockPtr = block;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forces[ii] = blockPtr;
blockPtr += 3;
}
std::vector<OpenMM::RealVec> forces( numberOfAtoms );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forces[ii][0] = zero;
forces[ii][1] = zero;
forces[ii][2] = zero;
}
vector<RealOpenMM>& bornForces = getBornForce();
bornForces.assign(numberOfAtoms, 0.0);
RealOpenMMVector bornForces( numberOfAtoms, 0.0);
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// first main loop
// first main loop
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// partial of polar term wrt Born radius
// and (dGpol/dr)(dr/dx)
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
if( atomI != atomJ ){
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
// 3 FLOP
#if 0
if( atomI == 0 ){
(void) fprintf( logFile, "bFCalc: %6d %6d %14.6e %14.6e %14.6e %14.6e\n", atomI, atomJ, dGpol_dalpha2_ij, bornRadii[atomJ], bornForces[atomI], bornRadii[atomI] );
}
#endif
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
#if( GBVIDebug == 1 )
{
double stupidFactor = three;
RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
int maxPrint = 10;
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
(void) fprintf( logFile, "Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, 1, gbviParameters->getTau() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI];
// partial of cavity term wrt Born radius
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
bornForces[atomI] += (stupidFactor*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
double xx = bornForces[atomI]*oneThird*b2*b2;
// xx*conversion should agree w/ values pulled out of kReduceGBVIBornForces_kernel in kForces.cu
(void) fprintf( logFile, "F1 %6d r/sclR[%14.6e %14.6e] bR=%14.6e bF=%14.6e %14.6e f[%14.6e %14.6e %14.6e](cnvrtd)"
" x[%14.6e %14.6e %14.6e]\n",
atomI, atomicRadii[atomI], scaledRadii[atomI], bornRadii[atomI], bornForces[atomI], xx*conversion,
// forces[atomI][0], forces[atomI][1], forces[atomI][2],
conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*forces[atomI][2],
atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
if( atomI == maxPrint ){
atomI = numberOfAtoms - maxPrint;
if( atomI < maxPrint )atomI = maxPrint;
}
}
(void) fflush( logFile );
}
#endif
// ---------------------------------------------------------------------------------------
// second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
// dGpol/dBornRadius) = bornForces[]
// dBornRadius/dr = (1/3)*(bR**4)*(dV/dr)
#if 0
(void) fprintf( logFile, "Clearing forces before loop2 periodic=%d cutoff=%d cutoffR=%14.7e\n",
_gbviParameters->getPeriodic(), _gbviParameters->getUseCutoff(), _gbviParameters->getCutoffDistance() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = zero;
forces[atomI][1] = zero;
forces[atomI][2] = zero;
}
(void) fflush( logFile );
#endif
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
std::vector<RealOpenMM>& switchDeriviative = getSwitchDeriviative();
RealOpenMM stupidFactor = three;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// partial of polar term wrt Born radius
// and (dGpol/dr)(dr/dx)
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM R = atomicRadii[atomI];
// partial of cavity term wrt Born radius
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
bornForces[atomI] += (stupidFactor*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
bornForces[atomI] *= switchDeriviative[atomI]*oneThird*b2*b2;
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaX = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
RealOpenMM deltaY = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
RealOpenMM deltaZ = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
......@@ -923,420 +600,288 @@ if( atomI == 0 ){
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
deltaX = deltaR[ReferenceForce::XIndex];
deltaY = deltaR[ReferenceForce::YIndex];
deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = SQRT( r2 );
RealOpenMM S = scaledRadii[atomJ];
RealOpenMM diff = (S - R);
RealOpenMM de = zero;
// find dRb/dr, where Rb is the Born radius
de = CpuGBVI::dL_dr( r, r+S, S ) + CpuGBVI::dL_dx( r, r+S, S );
if( FABS( diff ) < r ){
if( R > (r - S) ){
de -= CpuGBVI::dL_dr( r, R, S );
} else {
de -= ( CpuGBVI::dL_dr( r, (r-S), S ) + CpuGBVI::dL_dx( r, (r-S), S ) );
}
} else if( r < (S - R) ){
de -= ( CpuGBVI::dL_dr( r, r-S, S ) + CpuGBVI::dL_dx( r, r-S, S ) );
}
#if 0
RealOpenMM delta = (RealOpenMM) 1.0e-02;
(void) fprintf( stderr, "\n" );
for( int kk = 0; kk < 5; kk++ ){
RealOpenMM V1 = CpuGBVI::getVolume( r, R, S );
RealOpenMM V2 = CpuGBVI::getVolume( r+delta, R, S );
RealOpenMM df = (V2-V1)/delta;
(void) fprintf( stderr, "df %d %d [%14.6e %14.6e] V[%14.6e %14.6e] %.2e\n", atomI, atomJ, de, df, V2, V1, delta );
delta *= (RealOpenMM) 0.1;
}
double deltaD = 1.0e-02;
double ded = CpuGBVI::dL_drD( (double) r, r+S, S ) + CpuGBVI::dL_dxD( r, r+S, S ) - ( CpuGBVI::dL_drD( r, (r-S), S ) + CpuGBVI::dL_dxD( r, (r-S), S ) );
for( int kk = 0; kk < 5; kk++ ){
double V1 = CpuGBVI::getVolumeD( r, R, S );
double V2 = CpuGBVI::getVolumeD( r+deltaD, R, S );
double df = (V2-V1)/deltaD;
(void) fprintf( stderr, "df %d %d [%14.6e %14.6e] V[%14.6e %14.6e] %.2e\n", atomI, atomJ, ded, df, V2, V1, deltaD );
deltaD *= 0.1;
}
#endif
// de = (dG/dRb)(dRb/dr)
de *= bornForces[atomI]/r;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
if( atomI != atomJ ){
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
}
}
#if( GBVIDebug == 1 )
(void) fprintf( logFile, "Atom BornRadii BornForce Forces\n" );
double forceSum[3] = { 0.0, 0.0, 0.0 };
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forceSum[0] += forces[atomI][0];
forceSum[1] += forces[atomI][1];
forceSum[2] += forces[atomI][2];
(void) fprintf( logFile, "%4d %14.6e %14.6e [%14.6e %14.6e %14.6e]\n", atomI, bornRadii[atomI], bornForces[atomI], forces[atomI][0], forces[atomI][1], forces[atomI][2] );
}
(void) fprintf( logFile, "F sum=[%14.6e %14.6e %14.6e]\n", forceSum[0], forceSum[1], forceSum[2] );
(void) fflush( logFile );
#endif
// convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1];
inputForces[atomI][2] += conversion*forces[atomI][2];
}
#if( GBVIDebug == 1 )
{
(void) fprintf( logFile, "\nPost conversion\n" );
(void) fprintf( logFile, "Atom BornRadii BornForce Forces\n" );
int maxPrint = 10;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%4d %14.6e %14.6e [%14.6e %14.6e %14.6e]\n", atomI, bornRadii[atomI], conversion*bornForces[atomI],
inputForces[atomI][0], inputForces[atomI][1], inputForces[atomI][2] );
if( atomI == maxPrint ){
atomI = numberOfAtoms - maxPrint;
if( atomI < maxPrint )atomI = numberOfAtoms;
}
}
(void) fflush( logFile );
}
#endif
#undef GBVIDebug
delete[] forces;
delete[] block;
}
/**---------------------------------------------------------------------------------------
Get string w/ state
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
@param title title (optional)
@return string containing state
--------------------------------------------------------------------------------------- */
std::string CpuGBVI::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << CpuImplicitSolvent::getStateString( title );
return message.str();
}
/**---------------------------------------------------------------------------------------
Write Born energy and forces (Simbios)
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@param resultsFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
}
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
int CpuGBVI::writeBornEnergyForces( vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& forces,
const std::string& resultsFileName ) const {
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
static const char* methodName = "\nCpuGBVI::writeBornEnergyForces";
// dGpol/dBornRadius) = bornForces[]
// dBornRadius/dr = (1/3)*(bR**4)*(dV/dr)
// ---------------------------------------------------------------------------------------
/*
ImplicitSolventParameters* implicitSolventParameters = getImplicitSolventParameters();
const GBVIParameters* gbviParameters = static_cast<const GBVIParameters*>(implicitSolventParameters);
const RealOpenMMVector& scaledRadii = gbviParameters->getScaledRadii();
const RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI];
// partial of cavity term wrt Born radius
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
bornForces[atomI] += (three*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
bornForces[atomI] *= switchDeriviative[atomI]*oneThird*b2*b2;
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaX = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
RealOpenMM deltaY = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
RealOpenMM deltaZ = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
deltaX = deltaR[ReferenceForce::XIndex];
deltaY = deltaR[ReferenceForce::YIndex];
deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = SQRT( r2 );
RealOpenMM S = scaledRadii[atomJ];
RealOpenMM diff = (S - R);
RealOpenMM de = zero;
// find dRb/dr, where Rb is the Born radius
de = CpuGBVI::dL_dr( r, r+S, S ) + CpuGBVI::dL_dx( r, r+S, S );
if( FABS( diff ) < r ){
if( R > (r - S) ){
de -= CpuGBVI::dL_dr( r, R, S );
} else {
de -= ( CpuGBVI::dL_dr( r, (r-S), S ) + CpuGBVI::dL_dx( r, (r-S), S ) );
}
} else if( r < (S - R) ){
de -= ( CpuGBVI::dL_dr( r, r-S, S ) + CpuGBVI::dL_dx( r, r-S, S ) );
}
// de = (dG/dRb)(dRb/dr)
de *= bornForces[atomI]/r;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
}
}
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* bornRadii = getBornRadiiConst();
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadiusFactors();
const RealOpenMM* gbviChain = getObcChainConst();
const RealOpenMM energy = getEnergy();
// ---------------------------------------------------------------------------------------
// open file -- return if unsuccessful
FILE* implicitSolventResultsFile = NULL;
#ifdef WIN32
fopen_s( &implicitSolventResultsFile, resultsFileName.c_str(), "w" );
#else
implicitSolventResultsFile = fopen( resultsFileName.c_str(), "w" );
#endif
// diganostics
std::stringstream message;
message << methodName;
if( implicitSolventResultsFile != NULL ){
std::stringstream message;
message << methodName;
message << " Opened file=<" << resultsFileName << ">.";
SimTKOpenMMLog::printMessage( message );
} else {
std::stringstream message;
message << methodName;
message << " could not open file=<" << resultsFileName << "> -- abort output.";
SimTKOpenMMLog::printMessage( message );
return SimTKOpenMMCommon::ErrorReturn;
}
// header
(void) fprintf( implicitSolventResultsFile, "# %d atoms E=%.7e format: coords(3) bornRadii(input) q atomicRadii scaleFactors forces gbviChain\n",
numberOfAtoms, energy );
RealOpenMM forceConversion = (RealOpenMM) 1.0;
RealOpenMM lengthConversion = (RealOpenMM) 1.0;
// output
if( forces != NULL && atomCoordinates != NULL && partialCharges != NULL && atomicRadii != NULL ){
for( int ii = 0; ii < numberOfAtoms; ii++ ){
(void) fprintf( implicitSolventResultsFile, "%.7e %.7e %.7e %.7e %.5f %.5f %.5f %.7e %.7e %.7e %.7e\n",
lengthConversion*atomCoordinates[ii][0],
lengthConversion*atomCoordinates[ii][1],
lengthConversion*atomCoordinates[ii][2],
(bornRadii != NULL ? lengthConversion*bornRadii[ii] : 0.0),
partialCharges[ii], lengthConversion*atomicRadii[ii], scaledRadii[ii],
forceConversion*forces[ii][0],
forceConversion*forces[ii][1],
forceConversion*forces[ii][2],
forceConversion*gbviChain[ii]
);
}
}
(void) fclose( implicitSolventResultsFile );
*/
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get Obc Born energy and forces -- used debugging
@param bornRadii Born radii -- optional; if NULL, then GBVIParameters
entry is used
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
The array bornRadii is also updated and the obcEnergy
--------------------------------------------------------------------------------------- */
void CpuGBVI::computeBornEnergyForces( RealOpenMM* bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& forces ){
// ---------------------------------------------------------------------------------------
// convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
// static const char* methodName = "\nCpuGBVI::computeBornEnergyForcesPrint";
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1];
inputForces[atomI][2] += conversion*forces[atomI][2];
}
}
/**---------------------------------------------------------------------------------------
Use double precision
Use double precision
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
@return volume
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
double CpuGBVI::getVolumeD( double r, double R, double S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::getVolume";
// ---------------------------------------------------------------------------------------
static const double zero = 0.0;
static const double minusThree = -3.0;
static const double zero = 0.0;
static const double minusThree = -3.0;
double diff = (S - R);
if( fabs( diff ) < r ){
double diff = (S - R);
if( fabs( diff ) < r ){
double lowerBound = (R > (r - S)) ? R : (r - S);
double lowerBound = (R > (r - S)) ? R : (r - S);
return (CpuGBVI::getLD( r, (r + S), S ) -
CpuGBVI::getLD( r, lowerBound, S ));
return (CpuGBVI::getLD( r, (r + S), S ) -
CpuGBVI::getLD( r, lowerBound, S ));
} else if( r < diff ){
} else if( r < diff ){
return CpuGBVI::getLD( r, (r + S), S ) -
CpuGBVI::getLD( r, (r - S), S ) +
pow( R, minusThree );
return CpuGBVI::getLD( r, (r + S), S ) -
CpuGBVI::getLD( r, (r - S), S ) +
pow( R, minusThree );
} else {
return zero;
}
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
Use double precision
Use double precision
Get L (used in analytical solution for volume integrals)
Get L (used in analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
double CpuGBVI::getLD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::getL";
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double rInv = one/r;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff2 = (r + S)*(r - S);
double diff2 = (r + S)*(r - S);
return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
}
/**---------------------------------------------------------------------------------------
Use double precision
Use double precision
Get partial derivative of L wrt r
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
double CpuGBVI::dL_drD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::dL_dr";
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double threeEights = 0.375;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double threeEights = 0.375;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double rInv2 = rInv*rInv;
double rInv = one/r;
double rInv2 = rInv*rInv;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff2 = (r + S)*(r - S);
double diff2 = (r + S)*(r - S);
return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
}
/**---------------------------------------------------------------------------------------
Use double precision
Use double precision
Get partial derivative of L wrt x
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
double CpuGBVI::dL_dxD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVI::dL_dx";
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double half = 0.5;
static const double threeHalvesM = -1.5;
static const double third = 1.0/3.0;
static const double one = 1.0;
static const double half = 0.5;
static const double threeHalvesM = -1.5;
static const double third = 1.0/3.0;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double rInv = one/r;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff = (r + S)*(r - S);
double diff = (r + S)*(r - S);
return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
}
platforms/reference/src/gbsa/CpuGBVI.h
View file @ 326627a8
......@@ -25,26 +25,21 @@
#ifndef __CpuGBVI_H__
#define __CpuGBVI_H__
#include <vector>
#include "../SimTKUtilities/RealVec.h"
#include "GBVIParameters.h"
#include "CpuImplicitSolvent.h"
#include <vector>
// ---------------------------------------------------------------------------------------
class CpuGBVI : public CpuImplicitSolvent {
class CpuGBVI {
private:
// GB/VI parameters
GBVIParameters* _gbviParameters;
std::vector<RealOpenMM> _switchDeriviative;
// initialize data members (more than
// one constructor, so centralize intialization here)
void _initializeGBVIDataMembers( void );
RealOpenMMVector _switchDeriviative;
public:
......@@ -58,7 +53,7 @@ class CpuGBVI : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
CpuGBVI( ImplicitSolventParameters* gbviParameters );
CpuGBVI( GBVIParameters* gbviParameters );
/**---------------------------------------------------------------------------------------
......@@ -98,50 +93,7 @@ class CpuGBVI : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii );
/**---------------------------------------------------------------------------------------
Get Born energy and forces (not used)
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
void computeBornEnergyForces( RealOpenMM* bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& forces );
/**---------------------------------------------------------------------------------------
Get state
title title (optional)
@return state string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Write Born energy and forces (Simbios)
@param atomCoordinates atomic coordinates
@param partialCharges partial atom charges
@param forces force array
@param resultsFileName output file name
@return SimTKOpenMMCommon::DefaultReturn if file opened; else return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int writeBornEnergyForces( std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& forces,
const std::string& resultsFileName ) const;
void computeBornRadii( const std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMMVector& bornRadii );
/**---------------------------------------------------------------------------------------
......@@ -155,7 +107,7 @@ class CpuGBVI : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
static RealOpenMM getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S );
static RealOpenMM getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S );
/**---------------------------------------------------------------------------------------
......@@ -215,7 +167,6 @@ class CpuGBVI : public CpuImplicitSolvent {
Get GB/VI energy
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
......@@ -223,22 +174,20 @@ class CpuGBVI : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
RealOpenMM computeBornEnergy( const std::vector<RealOpenMM>& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges );
RealOpenMM computeBornEnergy( const std::vector<OpenMM::RealVec>& atomCoordinates, const RealOpenMMVector& partialCharges );
/**---------------------------------------------------------------------------------------
Get GB/VI forces
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces output forces
--------------------------------------------------------------------------------------- */
void computeBornForces( const std::vector<RealOpenMM>& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& inputForces );
void computeBornForces( std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges, std::vector<OpenMM::RealVec>& inputForces );
/**---------------------------------------------------------------------------------------
......@@ -305,7 +254,7 @@ class CpuGBVI : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getSwitchDeriviative( void );
RealOpenMMVector& getSwitchDeriviative( void );
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/gbsa/CpuImplicitSolvent.cpp deleted 100644 → 0
View file @ a587c652
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "CpuImplicitSolvent.h"
#include <cstdio>
using namespace OpenMM;
using namespace std;
//#define UseGromacsMalloc 1
// Replacement new/delete w/ Gromac's smalloc() and sfree()
// static data member created by call to cpuSetImplicitSolventParameters()
// stores parameter settings, ...
// used to calculate GBSA forces/energy
CpuImplicitSolvent* CpuImplicitSolvent::_cpuImplicitSolvent = NULL;
/**---------------------------------------------------------------------------------------
CpuImplicitSolvent constructor
@param implicitSolventParameters implicitSolventParameters object
--------------------------------------------------------------------------------------- */
CpuImplicitSolvent::CpuImplicitSolvent( ImplicitSolventParameters* implicitSolventParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::CpuImplicitSolvent(2)";
// ---------------------------------------------------------------------------------------
_initializeDataMembers( );
_implicitSolventParameters = implicitSolventParameters;
}
/**---------------------------------------------------------------------------------------
CpuImplicitSolvent destructor
--------------------------------------------------------------------------------------- */
CpuImplicitSolvent::~CpuImplicitSolvent( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::~CpuImplicitSolvent";
// ---------------------------------------------------------------------------------------
delete _implicitSolventParameters;
}
/**---------------------------------------------------------------------------------------
Initialize data members
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::_initializeDataMembers( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::_initializeDataMembers";
// ---------------------------------------------------------------------------------------
_includeAceApproximation = 0;
_forceConversionFactor = (RealOpenMM) 1.0;
_energyConversionFactor = (RealOpenMM) 1.0;
_forceCallIndex = 0;
_implicitSolventEnergy = (RealOpenMM) 0.0;
}
/**---------------------------------------------------------------------------------------
Return number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int CpuImplicitSolvent::getNumberOfAtoms( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getNumberOfAtoms";
// ---------------------------------------------------------------------------------------
return _implicitSolventParameters->getNumberOfAtoms();
}
/**---------------------------------------------------------------------------------------
Return energy
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM CpuImplicitSolvent::getEnergy( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getEnergy";
// ---------------------------------------------------------------------------------------
return _implicitSolventEnergy;
}
/**---------------------------------------------------------------------------------------
Delete static _cpuImplicitSolvent object if set
@return SimTKOpenMMCommon::DefaultReturn if _cpuImplicitSolvent was set;
otherwise return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int CpuImplicitSolvent::deleteCpuImplicitSolvent( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::deleteCpuImplicitSolvent";
// ---------------------------------------------------------------------------------------
if( _cpuImplicitSolvent != NULL ){
delete _cpuImplicitSolvent;
_cpuImplicitSolvent = NULL;
return SimTKOpenMMCommon::DefaultReturn;
} else {
return SimTKOpenMMCommon::ErrorReturn;
}
}
/**---------------------------------------------------------------------------------------
Set static member _cpuImplicitSolvent
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::setCpuImplicitSolvent( CpuImplicitSolvent* cpuImplicitSolvent ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::setCpuImplicitSolvent";
// ---------------------------------------------------------------------------------------
_cpuImplicitSolvent = cpuImplicitSolvent;
}
/**---------------------------------------------------------------------------------------
Get static member cpuImplicitSolvent
@return static member cpuImplicitSolvent
--------------------------------------------------------------------------------------- */
CpuImplicitSolvent* CpuImplicitSolvent::getCpuImplicitSolvent( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getCpuImplicitSolvent";
// ---------------------------------------------------------------------------------------
return _cpuImplicitSolvent;
}
/**---------------------------------------------------------------------------------------
Set energy
@param energy
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::setEnergy( RealOpenMM energy ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::setEnergy";
// ---------------------------------------------------------------------------------------
_implicitSolventEnergy = energy;
}
/**---------------------------------------------------------------------------------------
Return ImplicitSolventParameters
@return ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters* CpuImplicitSolvent::getImplicitSolventParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
return _implicitSolventParameters;
}
/**---------------------------------------------------------------------------------------
Set ImplicitSolventParameters
@param ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::setImplicitSolventParameters( ImplicitSolventParameters* implicitSolventParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::setImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
_implicitSolventParameters = implicitSolventParameters;
}
/**---------------------------------------------------------------------------------------
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return flag
--------------------------------------------------------------------------------------- */
int CpuImplicitSolvent::includeAceApproximation( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::includeAceApproximation";
// ---------------------------------------------------------------------------------------
return _includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether AceApproximation is to be included
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::setIncludeAceApproximation( int includeAceApproximation ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::setImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
_includeAceApproximation = includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Return ForceConversionFactor for units
@return ForceConversionFactor
--------------------------------------------------------------------------------------- */
RealOpenMM CpuImplicitSolvent::getForceConversionFactor( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getForceConversionFactor";
// ---------------------------------------------------------------------------------------
return _forceConversionFactor;
}
/**---------------------------------------------------------------------------------------
Set ForceConversionFactor
@param ForceConversionFactor (units conversion)
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::setForceConversionFactor( RealOpenMM forceConversionFactor ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::setForceConversionFactor";
// ---------------------------------------------------------------------------------------
_forceConversionFactor = forceConversionFactor;
}
/**---------------------------------------------------------------------------------------
Return EnergyConversionFactor for units
@return EnergyConversionFactor
--------------------------------------------------------------------------------------- */
RealOpenMM CpuImplicitSolvent::getEnergyConversionFactor( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getEnergyConversionFactor";
// ---------------------------------------------------------------------------------------
return _energyConversionFactor;
}
/**---------------------------------------------------------------------------------------
Set EnergyConversionFactor
@param EnergyConversionFactor (units conversion)
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::setEnergyConversionFactor( RealOpenMM energyConversionFactor ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::setEnergyConversionFactor";
// ---------------------------------------------------------------------------------------
_energyConversionFactor = energyConversionFactor;
}
/**---------------------------------------------------------------------------------------
Return ForceCallIndex -- number of times forces have been calculated
@return ForceCallIndex
--------------------------------------------------------------------------------------- */
int CpuImplicitSolvent::getForceCallIndex( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getForceCallIndex";
// ---------------------------------------------------------------------------------------
return _forceCallIndex;
}
/**---------------------------------------------------------------------------------------
Increment ForceCallIndex
@return incremented forceCallIndex
--------------------------------------------------------------------------------------- */
int CpuImplicitSolvent::incrementForceCallIndex( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::incrementForceCallIndex";
// ---------------------------------------------------------------------------------------
_forceCallIndex++;
return _forceCallIndex;
}
/**---------------------------------------------------------------------------------------
Return bornForce, a work array of size _implicitSolventParameters->getNumberOfAtoms()*sizeof( RealOpenMM )
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& CpuImplicitSolvent::getBornForce( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getImplicitSolventBornForce";
// ---------------------------------------------------------------------------------------
if( _bornForce.size() == 0 ){
_bornForce.resize(_implicitSolventParameters->getNumberOfAtoms());
}
return _bornForce;
}
/**---------------------------------------------------------------------------------------
Return Born radii: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if it is not set
@return array
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& CpuImplicitSolvent::getBornRadii( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getBornRadii";
// ---------------------------------------------------------------------------------------
if( _bornRadii.size() == 0 ){
_bornRadii.resize(_implicitSolventParameters->getNumberOfAtoms(), 0.0);
}
return _bornRadii;
}
/**---------------------------------------------------------------------------------------
Return Born radii: size = _implicitSolventParameters->getNumberOfAtoms()
@return array
--------------------------------------------------------------------------------------- */
const vector<RealOpenMM>& CpuImplicitSolvent::getBornRadiiConst( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getBornRadii";
// ---------------------------------------------------------------------------------------
return _bornRadii;
}
/**---------------------------------------------------------------------------------------
Return Born radii temp work array of size=_implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& CpuImplicitSolvent::getBornRadiiTemp( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getImplicitSolventBornRadiiTemp";
// ---------------------------------------------------------------------------------------
if( _tempBornRadii.size() == 0 ){
_tempBornRadii.resize(_implicitSolventParameters->getNumberOfAtoms(), 0.0);
}
return _tempBornRadii;
}
/**---------------------------------------------------------------------------------------
Compute Born radii
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param obcChain output array of Obc chain derivatives
--------------------------------------------------------------------------------------- */
/*
void CpuImplicitSolvent::computeBornRadii( vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nCpuImplicitSolvent::computeBornRadii";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << methodName;
message << " Error: calling from base class.";
SimTKOpenMMLog::printError( message );
} */
/**---------------------------------------------------------------------------------------
Get Born energy and forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces (output)
@param updateBornRadii if set, then Born radii are updated for current configuration;
otherwise radii correspond to configuration from previous iteration
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::computeImplicitSolventForces( vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges,
vector<RealVec>& forces, int updateBornRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nCpuImplicitSolvent::computeImplicitSolventForces";
// ---------------------------------------------------------------------------------------
int callId = incrementForceCallIndex();
// check parameters have been initialized
ImplicitSolventParameters* implicitSolventParameters = getImplicitSolventParameters();
if( !implicitSolventParameters ){
std::stringstream message;
message << methodName;
message << " implicitSolventParameters has not been initialized!";
SimTKOpenMMLog::printError( message );
}
// check to see if Born radii have been previously calculated
// if not, then calculate;
// logic here assumes that the radii are intitialized to zero
// and then once computed, always greater than zero.
// after first iteration Born radii are updated in force calculation (computeBornEnergyForces())
// unless updateBornRadii is set
vector<RealOpenMM>& bornRadii = getBornRadii();
if( updateBornRadii || bornRadii[0] < (RealOpenMM) 0.0001 || callId == 1 ){
computeBornRadii( atomCoordinates, bornRadii );
}
// compute forces
computeBornEnergyForces( getBornRadii(), atomCoordinates, partialCharges, forces );
}
/**---------------------------------------------------------------------------------------
Get Born energy and forces based on J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
@param bornRadii Born radii -- optional; if NULL, then ImplicitSolventParameters
entry is used
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::computeBornEnergyForces( vector<RealOpenMM>& bornRadii,
vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges,
vector<RealVec>& forces ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nCpuImplicitSolvent::computeBornEnergyForces";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << methodName;
message << " Error: calling from base class.";
SimTKOpenMMLog::printError( message );
}
/**---------------------------------------------------------------------------------------
Get nonpolar solvation force constribution via ACE approximation
@param implicitSolventParameters parameters
@param vdwRadii Vdw radii
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
--------------------------------------------------------------------------------------- */
void CpuImplicitSolvent::computeAceNonPolarForce( const ImplicitSolventParameters* implicitSolventParameters,
const vector<RealOpenMM>& bornRadii, RealOpenMM* energy,
vector<RealOpenMM>& forces ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::computeAceNonPolarForce";
static const RealOpenMM minusSix = -6.0;
// ---------------------------------------------------------------------------------------
// compute the nonpolar solvation via ACE approximation
const RealOpenMM probeRadius = implicitSolventParameters->getProbeRadius();
const RealOpenMM surfaceAreaFactor = implicitSolventParameters->getPi4Asolv();
const RealOpenMM* atomicRadii = implicitSolventParameters->getAtomicRadii();
int numberOfAtoms = implicitSolventParameters->getNumberOfAtoms();
// 1 + 1 + pow + 3 + 1 + 2 FLOP
// the original ACE equation is based on Eq.2 of
// M. Schaefer, C. Bartels and M. Karplus, "Solution Conformations
// and Thermodynamics of Structured Peptides: Molecular Dynamics
// Simulation with an Implicit Solvation Model", J. Mol. Biol.,
// 284, 835-848 (1998) (ACE Method)
// The original equation includes the factor (atomicRadii[atomI]/bornRadii[atomI]) to the first power,
// whereas here the ratio is raised to the sixth power: (atomicRadii[atomI]/bornRadii[atomI])**6
// This modification was made by Jay Ponder who observed it gave better correlations w/
// observed values. He did not think it was important enough to write up, so there is
// no paper to cite.
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
if( bornRadii[atomI] > 0.0 ){
RealOpenMM r = atomicRadii[atomI] + probeRadius;
RealOpenMM ratio6 = POW( atomicRadii[atomI]/bornRadii[atomI], (RealOpenMM) 6.0 );
RealOpenMM saTerm = surfaceAreaFactor*r*r*ratio6;
*energy += saTerm;
forces[atomI] += minusSix*saTerm/bornRadii[atomI];
}
}
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string containing state
--------------------------------------------------------------------------------------- */
std::string CpuImplicitSolvent::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
if( title ){
message << title;
}
message << "\nImplicitSolvent info:";
message << "\nForce conversion=" << getForceConversionFactor() << " Energy conversion=" << getEnergyConversionFactor();
message << "\nInclude ACE approximation=";
if( includeAceApproximation() ){
message << "Y";
} else {
message << "N";
}
// get parameters
message << getImplicitSolventParameters()->getStateString( NULL );
return message.str();
}
platforms/reference/src/gbsa/CpuImplicitSolvent.h deleted 100755 → 0
View file @ a587c652
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __CpuImplicitSolvent_H__
#define __CpuImplicitSolvent_H__
#include "ImplicitSolventParameters.h"
#include <vector>
// ---------------------------------------------------------------------------------------
class OPENMM_EXPORT CpuImplicitSolvent {
private:
// used for direct calls
static CpuImplicitSolvent* _cpuImplicitSolvent;
// parameters
ImplicitSolventParameters* _implicitSolventParameters;
// flag to signal whether ACE approximation
// is to be included
int _includeAceApproximation;
// force index call
int _forceCallIndex;
// work arrays
std::vector<RealOpenMM> _bornForce;
// Born radii and force
std::vector<RealOpenMM> _bornRadii;
std::vector<RealOpenMM> _tempBornRadii;
// convert units for energy/force
RealOpenMM _forceConversionFactor;
RealOpenMM _energyConversionFactor;
// Ed, 2007-04-27: Store the energy internally
RealOpenMM _implicitSolventEnergy;
/**---------------------------------------------------------------------------------------
Initialize data members -- potentially more than
one constructor, so centralize intialization here
--------------------------------------------------------------------------------------- */
void _initializeDataMembers( void );
protected:
/**---------------------------------------------------------------------------------------
Return implicitSolventBornForce, a work array of size _implicitSolventParameters->getNumberOfAtoms()*sizeof( RealOpenMM )
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getBornForce( void );
/**---------------------------------------------------------------------------------------
Return Born radii temp work array of size=_implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getBornRadiiTemp( void );
/**---------------------------------------------------------------------------------------
Set energy
@param energy new energy
--------------------------------------------------------------------------------------- */
void setEnergy( RealOpenMM energy );
public:
/**---------------------------------------------------------------------------------------
Constructor
@param implicitSolventParameters ImplicitSolventParameters reference
@return CpuImplicitSolvent object
--------------------------------------------------------------------------------------- */
CpuImplicitSolvent( ImplicitSolventParameters* implicitSolventParameters );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
virtual ~CpuImplicitSolvent( );
/**---------------------------------------------------------------------------------------
Delete static _cpuImplicitSolvent object if set
@return SimTKOpenMMCommon::DefaultReturn if _cpuImplicitSolvent was set;
otherwise return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
static int deleteCpuImplicitSolvent( void );
/**---------------------------------------------------------------------------------------
Set static member _cpuImplicitSolvent
--------------------------------------------------------------------------------------- */
static void setCpuImplicitSolvent( CpuImplicitSolvent* cpuImplicitSolvent );
/**---------------------------------------------------------------------------------------
Get static member cpuImplicitSolvent
@return static member cpuImplicitSolvent
--------------------------------------------------------------------------------------- */
static CpuImplicitSolvent* getCpuImplicitSolvent( void );
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int getNumberOfAtoms( void ) const;
/**---------------------------------------------------------------------------------------
Get energy
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM getEnergy( void ) const;
/**---------------------------------------------------------------------------------------
Return ImplicitSolventParameters
@return ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters* getImplicitSolventParameters( void ) const;
/**---------------------------------------------------------------------------------------
Set ImplicitSolventParameters
@param ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
void setImplicitSolventParameters( ImplicitSolventParameters* implicitSolventParameters );
/**---------------------------------------------------------------------------------------
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return flag
--------------------------------------------------------------------------------------- */
int includeAceApproximation( void ) const;
/**---------------------------------------------------------------------------------------
Set flag indicating whether AceApproximation is to be included
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void setIncludeAceApproximation( int includeAceApproximation );
/**---------------------------------------------------------------------------------------
Return ForceConversionFactor for units
@return ForceConversionFactor
--------------------------------------------------------------------------------------- */
RealOpenMM getForceConversionFactor( void ) const;
/**---------------------------------------------------------------------------------------
Set ForceConversionFactor
@param ForceConversionFactor (units conversion)
--------------------------------------------------------------------------------------- */
void setForceConversionFactor( RealOpenMM forceConversionFactor );
/**---------------------------------------------------------------------------------------
Return EnergyConversionFactor for units
@return EnergyConversionFactor
--------------------------------------------------------------------------------------- */
RealOpenMM getEnergyConversionFactor( void ) const;
/**---------------------------------------------------------------------------------------
Set EnergyConversionFactor
@param EnergyConversionFactor (units conversion)
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
void setEnergyConversionFactor( RealOpenMM energyConversionFactor );
/**---------------------------------------------------------------------------------------
Return ForceCallIndex -- number of times forces have been calculated
@return ForceCallIndex
--------------------------------------------------------------------------------------- */
int getForceCallIndex( void ) const;
/**---------------------------------------------------------------------------------------
Increment ForceCallIndex
@return incremented forceCallIndex
--------------------------------------------------------------------------------------- */
int incrementForceCallIndex( void );
/**---------------------------------------------------------------------------------------
Return Born radii: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getBornRadii( void );
/**---------------------------------------------------------------------------------------
Return Born radii: size = _implicitSolventParameters->getNumberOfAtoms()
@return array
--------------------------------------------------------------------------------------- */
const std::vector<RealOpenMM>& getBornRadiiConst( void ) const;
/**---------------------------------------------------------------------------------------
Get Born energy and forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces (output); if not set on input, then memory is allocated
@param updateBornRadii if set, then Born radii are updated for current configuration;
otherwise radii correspond to configuration from previous iteration
--------------------------------------------------------------------------------------- */
void computeImplicitSolventForces( std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges,
std::vector<OpenMM::RealVec>& forces, int updateBornRadii = 0 );
/**---------------------------------------------------------------------------------------
Get Born radii based on J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
@param atomCoordinates atomic coordinates dimension: [0-numberAtoms-1][0-2]
@param bornRadii output array of Born radii
@param obcChain output array of OBC chain derivative
@return array of Born radii
--------------------------------------------------------------------------------------- */
virtual void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii ) = 0;
/**---------------------------------------------------------------------------------------
Get Born energy and forces based on J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@return force array
--------------------------------------------------------------------------------------- */
virtual void computeBornEnergyForces( std::vector<RealOpenMM>& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& forces );
/**---------------------------------------------------------------------------------------
Get nonpolar solvation force constribution via ACE approximation
@param implicitSolventParameters parameters
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
--------------------------------------------------------------------------------------- */
void computeAceNonPolarForce( const ImplicitSolventParameters* implicitSolventParameters,
const std::vector<RealOpenMM>& bornRadii, RealOpenMM* energy,
std::vector<RealOpenMM>& forces ) const;
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string containing state
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
};
// ---------------------------------------------------------------------------------------
#endif // __CpuImplicitSolvent_H__
platforms/reference/src/gbsa/CpuObc.cpp
View file @ 326627a8
......@@ -25,624 +25,475 @@
#include <string.h>
#include <stdlib.h>
#include <sstream>
#include <cmath>
#include <cstdio>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "CpuObc.h"
#include "../SimTKReference/ReferenceForce.h"
#include <cmath>
#include <cstdio>
#include "CpuObc.h"
using namespace OpenMM;
using namespace std;
/**---------------------------------------------------------------------------------------
CpuObc constructor
CpuObc constructor
obcParameters obcParameters object
--------------------------------------------------------------------------------------- */
CpuObc::CpuObc( ImplicitSolventParameters* obcParameters ) : CpuImplicitSolvent( obcParameters ){
// ---------------------------------------------------------------------------------------
obcParameters obcParameters object
--------------------------------------------------------------------------------------- */
// static const char* methodName = "\nCpuObc::CpuObc";
CpuObc::CpuObc( ObcParameters* obcParameters ) : _obcParameters(obcParameters), _includeAceApproximation(1) {
_obcChain.resize(_obcParameters->getNumberOfAtoms());
}
// ---------------------------------------------------------------------------------------
/**---------------------------------------------------------------------------------------
_initializeObcDataMembers( );
CpuObc destructor
_obcParameters = static_cast<ObcParameters*> (obcParameters);
--------------------------------------------------------------------------------------- */
CpuObc::~CpuObc( ){
}
/**---------------------------------------------------------------------------------------
CpuObc destructor
Get ObcParameters reference
--------------------------------------------------------------------------------------- */
CpuObc::~CpuObc( ){
// ---------------------------------------------------------------------------------------
@return ObcParameters reference
// static const char* methodName = "\nCpuObc::~CpuObc";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
ObcParameters* CpuObc::getObcParameters( void ) const {
return _obcParameters;
}
/**---------------------------------------------------------------------------------------
Initialize data members
--------------------------------------------------------------------------------------- */
Set ObcParameters reference
void CpuObc::_initializeObcDataMembers( void ){
@param ObcParameters reference
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
// static const char* methodName = "\nCpuObc::initializeDataMembers";
// ---------------------------------------------------------------------------------------
_obcParameters = NULL;
void CpuObc::setObcParameters( ObcParameters* obcParameters ){
_obcParameters = obcParameters;
}
/**---------------------------------------------------------------------------------------
Get ObcParameters reference
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return ObcParameters reference
@return flag
--------------------------------------------------------------------------------------- */
ObcParameters* CpuObc::getObcParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::getObcParameters";
// ---------------------------------------------------------------------------------------
return _obcParameters;
int CpuObc::includeAceApproximation( void ) const {
return _includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Set ObcParameters reference
Set flag indicating whether AceApproximation is to be included
@param ObcParameters reference
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void CpuObc::setObcParameters( ObcParameters* obcParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setObcParameters";
// ---------------------------------------------------------------------------------------
_obcParameters = obcParameters;
void CpuObc::setIncludeAceApproximation( int includeAceApproximation ){
_includeAceApproximation = includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
@return array
@return array
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& CpuObc::getObcChain( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::getObcChain";
// ---------------------------------------------------------------------------------------
if( _obcChain.size() == 0 ){
_obcChain.resize(_obcParameters->getNumberOfAtoms());
}
return _obcChain;
return _obcChain;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
Get Born radii based on papers:
@return array
J. Phys. Chem. 1996 100, 19824-19839 (HCT paper)
Proteins: Structure, Function, and Bioinformatcis 55:383-394 (2004) (OBC paper)
--------------------------------------------------------------------------------------- */
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
const vector<RealOpenMM>& CpuObc::getObcChainConst( void ) const {
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
void CpuObc::computeBornRadii( const vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii ){
// static const char* methodName = "\nCpuObc::getObcChain";
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
return _obcChain;
}
// ---------------------------------------------------------------------------------------
/**---------------------------------------------------------------------------------------
ObcParameters* obcParameters = getObcParameters();
Return OBC chain temp work array of size=_obcParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
int numberOfAtoms = obcParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = obcParameters->getAtomicRadii();
const RealOpenMMVector& scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
RealOpenMMVector& obcChain = getObcChain();
@return array
RealOpenMM dielectricOffset = obcParameters->getDielectricOffset();
RealOpenMM alphaObc = obcParameters->getAlphaObc();
RealOpenMM betaObc = obcParameters->getBetaObc();
RealOpenMM gammaObc = obcParameters->getGammaObc();
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
vector<RealOpenMM>& CpuObc::getObcChainTemp( void ){
// calculate Born radii
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
// static const char* methodName = "\nCpuObc::getImplicitSolventObcChainTemp";
RealOpenMM radiusIInverse = one/offsetRadiusI;
RealOpenMM sum = zero;
// ---------------------------------------------------------------------------------------
// HCT code
return _obcChainTemp;
}
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
/**---------------------------------------------------------------------------------------
if( atomJ != atomI ){
Get Born radii based on papers:
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_obcParameters->getUseCutoff() && r > _obcParameters->getCutoffDistance())
continue;
J. Phys. Chem. 1996 100, 19824-19839 (HCT paper)
Proteins: Structure, Function, and Bioinformatcis 55:383-394 (2004) (OBC paper)
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM rInverse = one/r;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
--------------------------------------------------------------------------------------- */
RealOpenMM u_ij = one/rScaledRadiusJ;
void CpuObc::computeBornRadii( vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii ){
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM ratio = LN( (u_ij/l_ij) );
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + ( half*rInverse*ratio) + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
// ---------------------------------------------------------------------------------------
// this case (atom i completely inside atom j) is not considered in the original paper
// Jay Ponder and the authors of Tinker recognized this and
// worked out the details
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM two = (RealOpenMM) 2.0;
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM fourth = (RealOpenMM) 0.25;
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( radiusIInverse - l_ij);
}
sum += term;
static const char* methodName = "\nCpuObc::computeBornRadii";
}
}
}
// OBC-specific code (Eqs. 6-8 in paper)
sum *= half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
bornRadii[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
obcChain[atomI] = offsetRadiusI*( alphaObc - two*betaObc*sum + three*gammaObc*sum2 );
obcChain[atomI] = (one - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
// ---------------------------------------------------------------------------------------
}
}
ObcParameters* obcParameters = getObcParameters();
/**---------------------------------------------------------------------------------------
int numberOfAtoms = obcParameters->getNumberOfAtoms();
RealOpenMM* atomicRadii = obcParameters->getAtomicRadii();
const RealOpenMM* scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
vector<RealOpenMM>& obcChain = getObcChain();
Get nonpolar solvation force constribution via ACE approximation
RealOpenMM dielectricOffset = obcParameters->getDielectricOffset();
RealOpenMM alphaObc = obcParameters->getAlphaObc();
RealOpenMM betaObc = obcParameters->getBetaObc();
RealOpenMM gammaObc = obcParameters->getGammaObc();
@param obcParameters parameters
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
// calculate Born radii
void CpuObc::computeAceNonPolarForce( const ObcParameters* obcParameters,
const RealOpenMMVector& bornRadii,
RealOpenMM* energy,
RealOpenMMVector& forces ) const {
//FILE* logFile = SimTKOpenMMLog::getSimTKOpenMMLogFile( );
//FILE* logFile = NULL;
//FILE* logFile = fopen( "bR", "w" );
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM minusSix = -6.0;
static const RealOpenMM six = static_cast<RealOpenMM>( 6.0 );
RealOpenMM radiusIInverse = one/offsetRadiusI;
RealOpenMM sum = zero;
// ---------------------------------------------------------------------------------------
// HCT code
// compute the nonpolar solvation via ACE approximation
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
const RealOpenMM probeRadius = obcParameters->getProbeRadius();
const RealOpenMM surfaceAreaFactor = obcParameters->getPi4Asolv();
if( atomJ != atomI ){
const RealOpenMMVector& atomicRadii = obcParameters->getAtomicRadii();
int numberOfAtoms = obcParameters->getNumberOfAtoms();
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_obcParameters->getUseCutoff() && r > _obcParameters->getCutoffDistance())
continue;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM rInverse = one/r;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM ratio = LN( (u_ij/l_ij) );
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + ( half*rInverse*ratio) + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
// this case (atom i completely inside atom j) is not considered in the original paper
// Jay Ponder and the authors of Tinker recognized this and
// worked out the details
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( radiusIInverse - l_ij);
}
sum += term;
/*
if( logFile && atomI == 0 ){
(void) fprintf( logFile, "\nRR %d %d r=%.4f rads[%.6f %.6f] scl=[%.3f %.3f] sum=%12.6e %12.6e %12.6e %12.6e",
atomI, atomJ, r, offsetRadiusI, offsetRadiusJ, scaledRadiusFactor[atomI], scaledRadiusFactor[atomJ], 0.5f*sum,
l_ij, u_ij, term );
}
*/
// the original ACE equation is based on Eq.2 of
}
}
}
// OBC-specific code (Eqs. 6-8 in paper)
// M. Schaefer, C. Bartels and M. Karplus, "Solution Conformations
// and Thermodynamics of Structured Peptides: Molecular Dynamics
// Simulation with an Implicit Solvation Model", J. Mol. Biol.,
// 284, 835-848 (1998) (ACE Method)
sum *= half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
bornRadii[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
obcChain[atomI] = offsetRadiusI*( alphaObc - two*betaObc*sum + three*gammaObc*sum2 );
obcChain[atomI] = (one - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
// The original equation includes the factor (atomicRadii[atomI]/bornRadii[atomI]) to the first power,
// whereas here the ratio is raised to the sixth power: (atomicRadii[atomI]/bornRadii[atomI])**6
/*
if( logFile ){
(void) fprintf( logFile, "\nRRQ %d sum %12.6e tanhS %12.6e radI %.5f %.5f born %18.10e obc %12.6e",
atomI, sum, tanhSum, radiusI, offsetRadiusI, bornRadii[atomI], obcChain[atomI] );
}
*/
// This modification was made by Jay Ponder who observed it gave better correlations w/
// observed values. He did not think it was important enough to write up, so there is
// no paper to cite.
}
/*
if( logFile ){
(void) fclose( logFile );
} */
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
if( bornRadii[atomI] > zero ){
RealOpenMM r = atomicRadii[atomI] + probeRadius;
RealOpenMM ratio6 = POW( atomicRadii[atomI]/bornRadii[atomI], six );
RealOpenMM saTerm = surfaceAreaFactor*r*r*ratio6;
*energy += saTerm;
forces[atomI] += minusSix*saTerm/bornRadii[atomI];
}
}
}
/**---------------------------------------------------------------------------------------
Get Obc Born energy and forces
@param bornRadii Born radii -- optional; if NULL, then ObcParameters
entry is used
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
Get Obc Born energy and forces
The array bornRadii is also updated and the obcEnergy
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
The array bornRadii is also updated and the obcEnergy
void CpuObc::computeBornEnergyForces( vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges, vector<RealVec>& inputForces ){
// static const char* methodName = "\nCpuObc::computeBornEnergyForces";
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM two = (RealOpenMM) 2.0;
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM four = (RealOpenMM) 4.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125;
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
const ObcParameters* obcParameters = getObcParameters();
const int numberOfAtoms = obcParameters->getNumberOfAtoms();
const ObcParameters* obcParameters = getObcParameters();
const int numberOfAtoms = obcParameters->getNumberOfAtoms();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// constants
// constants
const RealOpenMM preFactor = obcParameters->getPreFactor();
const RealOpenMM dielectricOffset = obcParameters->getDielectricOffset();
RealOpenMM preFactor;
if( obcParameters->getSoluteDielectric() != zero && obcParameters->getSolventDielectric() != zero ){
preFactor = two*obcParameters->getElectricConstant()*( (one/obcParameters->getSoluteDielectric()) - (one/obcParameters->getSolventDielectric()) );
} else {
preFactor = zero;
}
// ---------------------------------------------------------------------------------------
const RealOpenMM dielectricOffset = obcParameters->getDielectricOffset();
#if 0
{
RealOpenMM* atomicRadii = obcParameters->getAtomicRadii();
const RealOpenMM* scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
RealOpenMM* obcChain = getObcChain();
FILE* logFile = fopen( "bornParameters", "w" );
(void) fprintf( logFile, "%5d dielOff=%.4e rad::hct::q::bR::Chain::coords\n", numberOfAtoms, dielectricOffset );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%5d %10.5f %10.5f %10.5f %14.7e %14.7e %14.7e %14.7e %14.7e\n", atomI,
atomicRadii[atomI], scaledRadiusFactor[atomI], partialCharges[atomI], bornRadii[atomI], obcChain[atomI],
atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
}
(void) fclose( logFile );
}
#endif
// ---------------------------------------------------------------------------------------
// set energy/forces to zero
// compute Born radii
RealOpenMM obcEnergy = zero;
RealOpenMMVector bornRadii( numberOfAtoms );
computeBornRadii( atomCoordinates, bornRadii );
RealOpenMM** forces = (RealOpenMM**) malloc( sizeof( RealOpenMM* )*numberOfAtoms );
RealOpenMM* block = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms*3 );
memset( block, 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
RealOpenMM* blockPtr = block;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forces[ii] = blockPtr;
blockPtr += 3;
}
// set energy/forces to zero
vector<RealOpenMM>& bornForces = getBornForce();
bornForces.assign(numberOfAtoms, 0.0);
RealOpenMM obcEnergy = zero;
RealOpenMMVector bornForces( numberOfAtoms, 0.0 );
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// N*( 8 + pow) ACE
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcParameters, bornRadii, &obcEnergy, bornForces );
}
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcParameters, bornRadii, &obcEnergy, bornForces );
}
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// first main loop
// first main loop
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
continue;
// 3 FLOP
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
// exp + 2 + sqrt FLOP
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
// 6 FLOP
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
if( atomI != atomJ ){
// 5 FLOP
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
// 11 FLOP
inputForces[atomI][0] += deltaX;
inputForces[atomI][1] += deltaY;
inputForces[atomI][2] += deltaZ;
if( atomI != atomJ ){
inputForces[atomJ][0] -= deltaX;
inputForces[atomJ][1] -= deltaY;
inputForces[atomJ][2] -= deltaZ;
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
} else {
Gpol *= half;
}
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
obcEnergy += Gpol;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
}
}
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
// ---------------------------------------------------------------------------------------
} else {
Gpol *= half;
}
// second main loop
// 3 FLOP
const RealOpenMMVector& obcChain = getObcChain();
const RealOpenMMVector& atomicRadii = obcParameters->getAtomicRadii();
obcEnergy += Gpol;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
//obcEnergy *= getEnergyConversionFactor();
// ---------------------------------------------------------------------------------------
// second main loop
// initialize Born radii & ObcChain temp arrays -- contain values
// used in next iteration
vector<RealOpenMM>& bornRadiiTemp = getBornRadiiTemp();
bornRadiiTemp.assign(numberOfAtoms, 0.0);
vector<RealOpenMM>& obcChainTemp = getObcChainTemp();
obcChainTemp.assign(numberOfAtoms, 0.0);
vector<RealOpenMM>& obcChain = getObcChain();
const RealOpenMM* atomicRadii = obcParameters->getAtomicRadii();
const RealOpenMM alphaObc = obcParameters->getAlphaObc();
const RealOpenMM betaObc = obcParameters->getBetaObc();
const RealOpenMM gammaObc = obcParameters->getGammaObc();
const RealOpenMM* scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
const RealOpenMM alphaObc = obcParameters->getAlphaObc();
const RealOpenMM betaObc = obcParameters->getBetaObc();
const RealOpenMM gammaObc = obcParameters->getGammaObc();
const RealOpenMMVector& scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
// compute factor that depends only on the outer loop index
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bornForces[atomI] *= bornRadii[atomI]*bornRadii[atomI]*obcChain[atomI];
}
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bornForces[atomI] *= bornRadii[atomI]*bornRadii[atomI]*obcChain[atomI];
}
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// radius w/ dielectric offset applied
// radius w/ dielectric offset applied
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
// used to compute Born radius for next iteration
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM bornSum = zero;
if( atomJ != atomI ){
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
// radius w/ dielectric offset applied
// radius w/ dielectric offset applied
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
if( offsetRadiusI < rScaledRadiusJ ){
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
RealOpenMM de = bornForces[atomI]*t3*rInverse;
RealOpenMM de = bornForces[atomI]*t3*rInverse;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
forces[atomI][0] -= deltaX;
forces[atomI][1] -= deltaY;
forces[atomI][2] -= deltaZ;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
inputForces[atomI][0] -= deltaX;
inputForces[atomI][1] -= deltaY;
inputForces[atomI][2] -= deltaZ;
forces[atomJ][0] += deltaX;
forces[atomJ][1] += deltaY;
forces[atomJ][2] += deltaZ;
inputForces[atomJ][0] += deltaX;
inputForces[atomJ][1] += deltaY;
inputForces[atomJ][2] += deltaZ;
// Born radius term
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + (half*rInverse)*LN(u_ij/l_ij) +
(fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2-u_ij2);
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( (one/offsetRadiusI) - l_ij);
}
bornSum += term;
}
}
}
// OBC-specific code (Eqs. 6-8 in paper)
bornSum *= half*offsetRadiusI;
RealOpenMM sum2 = bornSum*bornSum;
RealOpenMM sum3 = bornSum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*bornSum - betaObc*sum2 + gammaObc*sum3 );
bornRadiiTemp[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
obcChainTemp[atomI] = offsetRadiusI*( alphaObc - two*betaObc*bornSum + three*gammaObc*sum2 );
obcChainTemp[atomI] = (one - tanhSum*tanhSum)*obcChainTemp[atomI]/radiusI;
}
// cal to Joule conversion
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += forces[atomI][0];
inputForces[atomI][1] += forces[atomI][1];
inputForces[atomI][2] += forces[atomI][2];
}
setEnergy( obcEnergy );
// copy new Born radii and obcChain values into permanent array
bornRadii = bornRadiiTemp;
obcChain = obcChainTemp;
free( (char*) block );
free( (char*) forces );
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string containing state
--------------------------------------------------------------------------------------- */
std::string CpuObc::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getStateString";
// ---------------------------------------------------------------------------------------
}
}
}
std::stringstream message;
message << CpuImplicitSolvent::getStateString( title );
}
return message.str();
return obcEnergy;
}
platforms/reference/src/gbsa/CpuObc.h
View file @ 326627a8
......@@ -26,11 +26,10 @@
#define __CpuObc_H__
#include "ObcParameters.h"
#include "CpuImplicitSolvent.h"
// ---------------------------------------------------------------------------------------
class CpuObc : public CpuImplicitSolvent {
class CpuObc {
private:
......@@ -40,13 +39,13 @@ class CpuObc : public CpuImplicitSolvent {
// arrays containing OBC chain derivative
std::vector<RealOpenMM> _obcChain;
std::vector<RealOpenMM> _obcChainTemp;
RealOpenMMVector _obcChain;
// initialize data members (more than
// one constructor, so centralize intialization here)
// flag to signal whether ACE approximation
// is to be included
int _includeAceApproximation;
void _initializeObcDataMembers( void );
public:
......@@ -60,7 +59,7 @@ class CpuObc : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
CpuObc( ImplicitSolventParameters* obcParameters );
CpuObc( ObcParameters* obcParameters );
/**---------------------------------------------------------------------------------------
......@@ -92,26 +91,33 @@ class CpuObc : public CpuImplicitSolvent {
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return array
@return flag
--------------------------------------------------------------------------------------- */
int includeAceApproximation( void ) const;
/**---------------------------------------------------------------------------------------
Set flag indicating whether AceApproximation is to be included
std::vector<RealOpenMM>& getObcChain( void );
const std::vector<RealOpenMM>& getObcChainConst( void ) const;
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void setIncludeAceApproximation( int includeAceApproximation );
/**---------------------------------------------------------------------------------------
Return OBC chain temp work array of size=_implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
Return OBC chain derivative: size = _implicitSolventParameters->getNumberOfAtoms()
@return array
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getObcChainTemp( void );
RealOpenMMVector& getObcChain( void );
/**---------------------------------------------------------------------------------------
......@@ -119,38 +125,38 @@ class CpuObc : public CpuImplicitSolvent {
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param obcChain output array of OBC chain derivative factors; if NULL,
then ignored
--------------------------------------------------------------------------------------- */
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii );
void computeBornRadii( const std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMMVector& bornRadii );
/**---------------------------------------------------------------------------------------
Get nonpolar solvation force constribution via ACE approximation
@param obcParameters parameters
@param vdwRadii Vdw radii
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
--------------------------------------------------------------------------------------- */
void computeAceNonPolarForce( const ObcParameters* obcParameters, const RealOpenMMVector& bornRadii,
RealOpenMM* energy, RealOpenMMVector& forces ) const;
/**---------------------------------------------------------------------------------------
Get Born energy and forces based on OBC
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
void computeBornEnergyForces( std::vector<RealOpenMM>& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& forces );
/**---------------------------------------------------------------------------------------
Get state
title title (optional)
@return state string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
RealOpenMM computeBornEnergyForces( const std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges, std::vector<OpenMM::RealVec>& forces );
};
......
platforms/reference/src/gbsa/GBVIParameters.cpp
View file @ 326627a8
......@@ -26,590 +26,399 @@
#include <sstream>
#include <string.h>
#include "openmm/OpenMMException.h"
#include "GBVIParameters.h"
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
using std::vector;
using OpenMM::RealVec;
const std::string GBVIParameters::ParameterFileName = std::string( "params.agb" );
/**---------------------------------------------------------------------------------------
GBVIParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in Obc equation gPol
(4) list of atoms that should be excluded in calculating
force -- nonbonded atoms (1-2, and 1-3 atoms)
Implementation:
Slightly different sequence of calls when running on CPU vs GPU.
Difference arise because the CPU-side data arrays for the Brook
streams are allocated by the BrookStreamWrapper objects. These
arrays are then used by GBVIParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
GBVIParameters* gb_VIParameters = new GBVIParameters( numberOfAtoms, log );
gb_VIParameters->initializeParameters( top );
Gpu:
gb_VIParameters = new GBVIParameters( gpu->natoms, log );
// set arrays for cpu using stream data field;
// initializeParameters() only allocates space for arrays if they are not set (==NULL)
// also set flag so that GBVIParameters destructor does not free arrays
gb_VIParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuObc::gb_VIVdwRadii )->getData() );
gb_VIParameters->setVolume( getBrookStreamWrapperAtIndex( GpuObc::gb_VIVolume )->getData() );
gb_VIParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuObc::gb_VIGpolFixed )->getData() );
gb_VIParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuObc::gb_VIBornRadii )->getData() );
gb_VIParameters->setFreeArrays( false );
gb_VIParameters->initializeParameters( top );
Issues:
Tinker's atom radii are used.
The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
only tested for generic proteins
see mapGmxAtomNameToTinkerAtomNumber()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
GBVIParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
GBVIParameters constructor
GBVIParameters::GBVIParameters( int numberOfAtoms ) : ImplicitSolventParameters( numberOfAtoms ), cutoff(false), periodic(false) {
@param numberOfAtoms number of atoms
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
// static const char* methodName = "\nGBVIParameters::GBVIParameters";
// ---------------------------------------------------------------------------------------
GBVIParameters::GBVIParameters( int numberOfAtoms ) : _numberOfAtoms(numberOfAtoms),
_soluteDielectric(1.0),
_solventDielectric(78.3),
_electricConstant(-0.5*ONE_4PI_EPS0),
_cutoff(false),
_periodic(false),
_bornRadiusScalingMethod(0),
_quinticLowerLimitFactor(0.8),
_quinticUpperBornRadiusLimit(5.0) {
_ownScaledRadii = 0;
_scaledRadii = NULL;
_ownGammaParameters = 0;
_gammaParameters = NULL;
_bornRadiusScalingMethod = 0;
_quinticLowerLimitFactor = 0.8f;
_quinticUpperBornRadiusLimit = 5.0f;
_atomicRadii.resize( numberOfAtoms );
_scaledRadii.resize( numberOfAtoms );
_gammaParameters.resize( numberOfAtoms );
}
/**---------------------------------------------------------------------------------------
GBVIParameters destructor (Simbios)
GBVIParameters destructor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
GBVIParameters::~GBVIParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVIParameters::~GBVIParameters";
// ---------------------------------------------------------------------------------------
// in GPU runs, arrays may be 'owned' by BrookStreamWrapper -- hence they should not
// be freed here, i.e., _freeArrays should be 'false'
#ifdef UseGromacsMalloc
/*
if( _freeArrays ){
if( _vdwRadii != NULL ){
save_free( "_vdwRadii", __FILE__, __LINE__, _vdwRadii );
}
} */
#else
if( _ownScaledRadii ){
delete[] _scaledRadii;
}
delete[] _gammaParameters;
/*
if( getFreeArrays() ){
} */
#endif
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
Get number of atoms
@return array of atomic radii
@return number of atoms
--------------------------------------------------------------------------------------- */
RealOpenMM* GBVIParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
return atomicRadii;
int GBVIParameters::getNumberOfAtoms( void ) const {
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Get electric constant
@param atomicRadii array of atomic radii
@return electric constant
--------------------------------------------------------------------------------------- */
void GBVIParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVIParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
RealOpenMM GBVIParameters::getElectricConstant( void ) const {
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Get solvent dielectric
@param atomicRadii vector of atomic radii
@return solvent dielectric
--------------------------------------------------------------------------------------- */
void GBVIParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGBVIParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
RealOpenMM GBVIParameters::getSolventDielectric( void ) const {
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Return scaled radii
If not previously set, allocate space
Set solvent dielectric
@return array
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
const RealOpenMM* GBVIParameters::getScaledRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVIParameters::getScaledRadii";
// ---------------------------------------------------------------------------------------
if( _scaledRadii == NULL ){
GBVIParameters* localThis = const_cast<GBVIParameters* const>(this);
localThis->_scaledRadii = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownScaledRadii = true;
memset( _scaledRadii, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _scaledRadii;
void GBVIParameters::setSolventDielectric( RealOpenMM solventDielectric ){
_solventDielectric = solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
Get solute dielectric
@param ownScaledRadii flag indicating whether scale factors
array should be deleted
@return soluteDielectric
--------------------------------------------------------------------------------------- */
void GBVIParameters::setOwnScaledRadii( int ownScaledRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVIParameters::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownScaledRadii = ownScaledRadii;
RealOpenMM GBVIParameters::getSoluteDielectric( void ) const {
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set scaled radii
Set solute dielectric
@param scaledRadii scaledRadii
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void GBVIParameters::setScaledRadii( RealOpenMM* scaledRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadii";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadii && _scaledRadii != scaledRadii ){
delete[] _scaledRadii;
_ownScaledRadii = false;
}
_scaledRadii = scaledRadii;
void GBVIParameters::setSoluteDielectric( RealOpenMM soluteDielectric ){
_soluteDielectric = soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set scaled radii
@param scaledRadii scaledRadii
--------------------------------------------------------------------------------------- */
void GBVIParameters::setScaledRadii( const RealOpenMMVector& scaledRadii ){
// ---------------------------------------------------------------------------------------
Get AtomicRadii array
// static const char* methodName = "\nCpuObc::setScaledRadii";
@return array of atomic radii
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
if( _ownScaledRadii && _scaledRadii != NULL ){
delete[] _scaledRadii;
}
_ownScaledRadii = true;
_scaledRadii = new RealOpenMM[getNumberOfAtoms()];
for( int ii = 0; ii < (int) scaledRadii.size(); ii++ ){
_scaledRadii[ii] = scaledRadii[ii];
}
const RealOpenMMVector& GBVIParameters::getAtomicRadii( void ) const {
return _atomicRadii;
}
/**---------------------------------------------------------------------------------------
Return gamma parameters
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
Set AtomicRadii array
RealOpenMM* GBVIParameters::getGammaParameters( void ) const {
@param atomicRadii vector of atomic radii
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
// static const char* methodName = "\nGBVIParameters::getGammaParameters";
void GBVIParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
if( atomicRadii.size() == _atomicRadii.size() ){
for( unsigned int ii = 0; ii < atomicRadii.size(); ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
} else {
std::stringstream msg;
msg << "GBVIParameters: input size for atomic radii does not agree w/ current size: input=";
msg << atomicRadii.size();
msg << " current size=" << _atomicRadii.size();
throw OpenMM::OpenMMException(msg.str());
}
if( _gammaParameters == NULL ){
GBVIParameters* localThis = const_cast<GBVIParameters* const>(this);
localThis->_gammaParameters = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownGammaParameters = true;
memset( _gammaParameters, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _gammaParameters;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
@param ownGammaParameters flag indicating whether gamma parameter
array should be deleted
--------------------------------------------------------------------------------------- */
void GBVIParameters::setOwnGammaParameters( int ownGammaParameters ){
// ---------------------------------------------------------------------------------------
Return scaled radii
@return array
// static const char* methodName = "\nGBVIParameters::setOwnScaleFactors";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
_ownGammaParameters = ownGammaParameters;
const RealOpenMMVector& GBVIParameters::getScaledRadii( void ) const {
return _scaledRadii;
}
/**---------------------------------------------------------------------------------------
Set gamma parameters
@param gammas gamma parameters
--------------------------------------------------------------------------------------- */
void GBVIParameters::setGammaParameters( RealOpenMM* gammas ){
Set scaled radii
// ---------------------------------------------------------------------------------------
@param scaledRadii scaledRadii
// static const char* methodName = "\nCpuObc::setGammas";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
void GBVIParameters::setScaledRadii( const RealOpenMMVector& scaledRadii ){
if( _ownGammaParameters && _gammaParameters != gammas ){
delete[] _gammaParameters;
_ownGammaParameters = false;
}
if( scaledRadii.size() == _scaledRadii.size() ){
for( unsigned int ii = 0; ii < scaledRadii.size(); ii++ ){
_scaledRadii[ii] = scaledRadii[ii];
}
} else {
std::stringstream msg;
msg << "GBVIParameters: input size for scaled radii does not agree w/ current size: input=";
msg << scaledRadii.size();
msg << " current size=" << _scaledRadii.size();
throw OpenMM::OpenMMException(msg.str());
}
_gammaParameters = gammas;
}
/**---------------------------------------------------------------------------------------
Set gamma parameters
@param gammas gammas
--------------------------------------------------------------------------------------- */
void GBVIParameters::setGammaParameters( const RealOpenMMVector& gammas ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setGammas";
Return gamma parameters
If not previously set, allocate space
// ---------------------------------------------------------------------------------------
@return array
if( _ownGammaParameters && _gammaParameters != NULL ){
delete[] _gammaParameters;
}
_ownGammaParameters = true;
--------------------------------------------------------------------------------------- */
_gammaParameters = new RealOpenMM[getNumberOfAtoms()];
for( int ii = 0; ii < (int) gammas.size(); ii++ ){
_gammaParameters[ii] = gammas[ii];
}
const RealOpenMMVector& GBVIParameters::getGammaParameters( void ) const {
return _gammaParameters;
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string GBVIParameters::getStateString( const char* title ) const {
Set gamma parameters
// ---------------------------------------------------------------------------------------
@param gammas gammas
// static const char* methodName = "\nGBVIParameters::getStateString";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << ImplicitSolventParameters::getStateString( title );
std::string tab = getStringTab();
void GBVIParameters::setGammaParameters( const RealOpenMMVector& gammas ){
return message.str();
if( gammas.size() == _gammaParameters.size() ){
for( unsigned int ii = 0; ii < gammas.size(); ii++ ){
_gammaParameters[ii] = gammas[ii];
}
} else {
std::stringstream msg;
msg << "GBVIParameters: input size for gammas does not agree w/ current size: input=";
msg << gammas.size();
msg << " current size=" << _gammaParameters.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
Set the force to use a cutoff.
@param distance the cutoff distance
@param distance the cutoff distance
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void GBVIParameters::setUseCutoff( RealOpenMM distance ) {
cutoff = true;
cutoffDistance = distance;
_cutoff = true;
_cutoffDistance = distance;
}
/**---------------------------------------------------------------------------------------
Get whether to use a cutoff.
Get whether to use a cutoff.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
bool GBVIParameters::getUseCutoff() {
return cutoff;
return _cutoff;
}
/**---------------------------------------------------------------------------------------
Get the cutoff distance.
Get the cutoff distance.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getCutoffDistance() {
return cutoffDistance;
return _cutoffDistance;
}
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param boxSize the X, Y, and Z widths of the periodic box
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void GBVIParameters::setPeriodic( RealVec& boxSize ) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
assert(_cutoff);
assert(boxSize[0] >= 2.0*_cutoffDistance);
assert(boxSize[1] >= 2.0*_cutoffDistance);
assert(boxSize[2] >= 2.0*_cutoffDistance);
_periodic = true;
_periodicBoxSize[0] = boxSize[0];
_periodicBoxSize[1] = boxSize[1];
_periodicBoxSize[2] = boxSize[2];
}
/**---------------------------------------------------------------------------------------
Get whether to use periodic boundary conditions.
Get whether to use periodic boundary conditions.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
bool GBVIParameters::getPeriodic() {
return periodic;
return _periodic;
}
/**---------------------------------------------------------------------------------------
Get the periodic box dimension
Get the periodic box dimension
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
const RealOpenMM* GBVIParameters::getPeriodicBox() {
return periodicBoxSize;
return _periodicBoxSize;
}
/**---------------------------------------------------------------------------------------
Get tau prefactor
Get tau prefactor
@return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric
@return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getTau( void ) const {
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVIParameters::getTau:";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
RealOpenMM tau;
if( getSoluteDielectric() != zero && getSolventDielectric() != zero ){
tau = (one/getSoluteDielectric()) - (one/getSolventDielectric());
} else {
tau = zero;
}
RealOpenMM tau;
if( getSoluteDielectric() != zero && getSolventDielectric() != zero ){
tau = (one/getSoluteDielectric()) - (one/getSolventDielectric());
} else {
tau = zero;
}
return tau;
return tau;
}
/**---------------------------------------------------------------------------------------
Get bornRadiusScalingMethod
Get bornRadiusScalingMethod
@return bornRadiusScalingMethod
@return bornRadiusScalingMethod
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
int GBVIParameters::getBornRadiusScalingMethod( void ) const {
return _bornRadiusScalingMethod;
return _bornRadiusScalingMethod;
}
/**---------------------------------------------------------------------------------------
Set bornRadiusScalingMethod
Set bornRadiusScalingMethod
@param bornRadiusScalingMethod bornRadiusScalingMethod
@param bornRadiusScalingMethod bornRadiusScalingMethod
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void GBVIParameters::setBornRadiusScalingMethod( int bornRadiusScalingMethod ){
_bornRadiusScalingMethod = bornRadiusScalingMethod;
_bornRadiusScalingMethod = bornRadiusScalingMethod;
}
/**---------------------------------------------------------------------------------------
Get quinticLowerLimitFactor
Get quinticLowerLimitFactor
@return quinticLowerLimitFactor
@return quinticLowerLimitFactor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getQuinticLowerLimitFactor( void ) const {
return _quinticLowerLimitFactor;
return _quinticLowerLimitFactor;
}
/**---------------------------------------------------------------------------------------
Set quinticLowerLimitFactor
Set quinticLowerLimitFactor
@param quinticLowerLimitFactor quinticLowerLimitFactor
@param quinticLowerLimitFactor quinticLowerLimitFactor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void GBVIParameters::setQuinticLowerLimitFactor( RealOpenMM quinticLowerLimitFactor ){
_quinticLowerLimitFactor = quinticLowerLimitFactor;
_quinticLowerLimitFactor = quinticLowerLimitFactor;
}
/**---------------------------------------------------------------------------------------
Get quinticUpperBornRadiusLimit
Get quinticUpperBornRadiusLimit
@return quinticUpperBornRadiusLimit
@return quinticUpperBornRadiusLimit
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getQuinticUpperBornRadiusLimit( void ) const {
return _quinticUpperBornRadiusLimit;
return _quinticUpperBornRadiusLimit;
}
/**---------------------------------------------------------------------------------------
Set quinticUpperBornRadiusLimit
Set quinticUpperBornRadiusLimit
@param quinticUpperBornRadiusLimit quinticUpperBornRadiusLimit
@param quinticUpperBornRadiusLimit quinticUpperBornRadiusLimit
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void GBVIParameters::setQuinticUpperBornRadiusLimit( RealOpenMM quinticUpperBornRadiusLimit ){
_quinticUpperBornRadiusLimit = quinticUpperBornRadiusLimit;
_quinticUpperBornRadiusLimit = quinticUpperBornRadiusLimit;
}
platforms/reference/src/gbsa/GBVIParameters.h
View file @ 326627a8
......@@ -26,16 +26,13 @@
#define __GBVIParameters_H__
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "ImplicitSolventParameters.h"
// ---------------------------------------------------------------------------------------
class GBVIParameters : public ImplicitSolventParameters {
class GBVIParameters {
public:
static const std::string ParameterFileName;
/**
* This is an enumeration of the different methods that may be used for scaling of the Born radii.
*/
......@@ -44,10 +41,6 @@ class GBVIParameters : public ImplicitSolventParameters {
* No scaling method is applied.
*/
NoScaling = 0,
/**
* Use the method outlined in Proteins 55, 383-394 (2004), Eq. 6
*/
Tanh = 2,
/**
* Use quintic spline scaling function
*/
......@@ -56,21 +49,24 @@ class GBVIParameters : public ImplicitSolventParameters {
private:
// scaled radii
int _numberOfAtoms;
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _electricConstant;
int _ownScaledRadii;
RealOpenMM* _scaledRadii;
// parameter vectors
// gamma parameters
int _ownGammaParameters;
RealOpenMM* _gammaParameters;
RealOpenMMVector _atomicRadii;
RealOpenMMVector _scaledRadii;
RealOpenMMVector _gammaParameters;
// cutoff and periodic boundary conditions
bool cutoff;
bool periodic;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
bool _cutoff;
bool _periodic;
RealOpenMM _periodicBoxSize[3];
RealOpenMM _cutoffDistance;
int _bornRadiusScalingMethod;
RealOpenMM _quinticLowerLimitFactor;
......@@ -98,118 +94,123 @@ class GBVIParameters : public ImplicitSolventParameters {
/**---------------------------------------------------------------------------------------
Return scaled radii
Get number of atoms
@return array
@return number of atoms
--------------------------------------------------------------------------------------- */
const RealOpenMM* getScaledRadii( void ) const;
int getNumberOfAtoms( void ) const;
/**---------------------------------------------------------------------------------------
Return scaled radii
Get electric constant
@return array
@return electric constant
--------------------------------------------------------------------------------------- */
void setScaledRadii( RealOpenMM* scaledRadii );
void setScaledRadii( const RealOpenMMVector& scaledRadii );
RealOpenMM getElectricConstant( void ) const;
/**---------------------------------------------------------------------------------------
Set flag indicating whether scaled radii array should be deleted
Get solvent dielectric
@param ownScaledRadiusFactors flag indicating whether scaled radii
array should be deleted
@return solvent dielectric
--------------------------------------------------------------------------------------- */
void setOwnScaledRadii( int ownScaledRadii );
RealOpenMM getSolventDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
Set solvent dielectric
@return array of atom volumes
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM* getAtomicRadii( void ) const;
void setSolventDielectric( RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Get solute dielectric
@param atomicRadii array of atomic radii
@return soluteDielectric
--------------------------------------------------------------------------------------- */
void setAtomicRadii( RealOpenMM* atomicRadii );
RealOpenMM getSoluteDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Set solute dielectric
@param atomicRadii vector of atomic radii
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void setAtomicRadii( const RealOpenMMVector& atomicRadii );
void setSoluteDielectric( RealOpenMM soluteDielectric );
/**---------------------------------------------------------------------------------------
Set flag indicating whether gamma parameter array should be deleted
Return scaled radii
@param ownGammaParameters flag indicating whether gamma parameter
array should be deleted
@return array
--------------------------------------------------------------------------------------- */
void setOwnGammaParameters( int ownGammaParameters );
const RealOpenMMVector& getScaledRadii( void ) const;
/**---------------------------------------------------------------------------------------
Get GammaParameters array
Return scaled radii
@return array of gamma values
@return array
--------------------------------------------------------------------------------------- */
void setScaledRadii( const RealOpenMMVector& scaledRadii );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
@return array of atom volumes
--------------------------------------------------------------------------------------- */
RealOpenMM* getGammaParameters( void ) const;
const RealOpenMMVector& getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set GammaParameters array
Set AtomicRadii array
@param gammaParameters array of gamma parameters
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
void setGammaParameters( RealOpenMM* gammaParameters );
void setAtomicRadii( const RealOpenMMVector& atomicRadii );
/**---------------------------------------------------------------------------------------
Set GammaParameters array
Get GammaParameters array
@param gammaParameters array of gamma parameters
@return array of gamma values
--------------------------------------------------------------------------------------- */
void setGammaParameters( const RealOpenMMVector& gammaParameters );
const RealOpenMMVector& getGammaParameters( void ) const;
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
Set GammaParameters array
@param gammaParameters array of gamma parameters
--------------------------------------------------------------------------------------- */
void setGammaParameters( const RealOpenMMVector& gammaParameters );
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/gbsa/ImplicitSolventParameters.cpp deleted 100644 → 0
View file @ a587c652
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ImplicitSolventParameters.h"
#include <sstream>
#include <string.h>
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in ImplicitSolvent equation gPol
(4) list of atoms that should be excluded in calculating
force -- nonbonded atoms (1-2, and 1-3 atoms)
Implementation:
Slightly different sequence of calls when running on CPU vs GPU.
Difference arise because the CPU-side data arrays for the Brook
streams are allocated by the BrookStreamWrapper objects. These
arrays are then used by ImplicitSolventParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
ImplicitSolventParameters* implicitSolventParameters = new ImplicitSolventParameters( numberOfAtoms, log );
implicitSolventParameters->initializeParameters( top );
Gpu:
implicitSolventParameters = new ImplicitSolventParameters( gpu->natoms, log );
// set arrays for cpu using stream data field;
// initializeParameters() only allocates space for arrays if they are not set (==NULL)
// also set flag so that ImplicitSolventParameters destructor does not free arrays
implicitSolventParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventVdwRadii )->getData() );
implicitSolventParameters->setVolume( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventVolume )->getData() );
implicitSolventParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventGpolFixed )->getData() );
implicitSolventParameters->setAtomicRadii( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventAtomicRadii )->getData() );
implicitSolventParameters->setFreeArrays( false );
implicitSolventParameters->initializeParameters( top );
Issues:
Tinker's atom radii are used.
The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
only tested for generic proteins
see mapGmxAtomNameToTinkerAtomNumber()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters::ImplicitSolventParameters( int numberOfAtoms ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::ImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
_numberOfAtoms = numberOfAtoms;
_ownAtomicRadii = false;
_atomicRadii = NULL;
// see comments in ~ImplicitSolventParameters for explanation
_freeArrays = false;
_initializeImplicitSolventConstants();
}
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters::~ImplicitSolventParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::~ImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
if( _ownAtomicRadii ){
delete[] _atomicRadii;
}
}
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::getNumberOfAtoms( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getNumberOfAtoms:";
// ---------------------------------------------------------------------------------------
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getSolventDielectric( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getSolventDielectric:";
// ---------------------------------------------------------------------------------------
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setSolventDielectric( RealOpenMM solventDielectric ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setSolventDielectric:";
// ---------------------------------------------------------------------------------------
_solventDielectric = solventDielectric;
_resetPreFactor();
}
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getSoluteDielectric( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getSoluteDielectric:";
// ---------------------------------------------------------------------------------------
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setSoluteDielectric( RealOpenMM soluteDielectric ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setSoluteDielectric:";
// ---------------------------------------------------------------------------------------
_soluteDielectric = soluteDielectric;
_resetPreFactor();
}
/**---------------------------------------------------------------------------------------
Get electric constant (Simbios)
@return electricConstant
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getElectricConstant( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getElectricConstant:";
// ---------------------------------------------------------------------------------------
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Set electric constant (Simbios)
@param electricConstant electric constant
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setElectricConstant( RealOpenMM electricConstant ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setElectricConstant:";
// ---------------------------------------------------------------------------------------
_electricConstant = electricConstant;
_resetPreFactor();
}
/**---------------------------------------------------------------------------------------
Get probe radius (Simbios)
@return probeRadius
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getProbeRadius( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getProbeRadius:";
// ---------------------------------------------------------------------------------------
return _probeRadius;
}
/**---------------------------------------------------------------------------------------
Set probe radius (Simbios)
@param probeRadius probe radius
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setProbeRadius( RealOpenMM probeRadius ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setProbeRadius:";
// ---------------------------------------------------------------------------------------
_probeRadius = probeRadius;
}
/**---------------------------------------------------------------------------------------
Get pi*4*Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049*1000.0; (Still)
((RealOpenMM) M_PI)*4.0f*0.0054*1000.0; (OBC)
@return pi4Asolv
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getPi4Asolv( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getPi4Asolv:";
// ---------------------------------------------------------------------------------------
return _pi4Asolv;
}
/**---------------------------------------------------------------------------------------
Get prefactor
@return prefactor
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getPreFactor( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getPreFactor:";
// ---------------------------------------------------------------------------------------
return _preFactor;
}
/**---------------------------------------------------------------------------------------
Set pi4Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049*1000.0; (Still)
((RealOpenMM) M_PI)*4.0f*0.0054*1000.0; (OBC)
@param pi4Asolv probe radius
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setPi4Asolv( RealOpenMM pi4Asolv ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setPi4Asolv:";
// ---------------------------------------------------------------------------------------
_pi4Asolv = pi4Asolv;
}
/**---------------------------------------------------------------------------------------
Get conversion factor for kcal/A -> kJ/nm (Simbios)
@return kcalA_To_kJNm factor
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getKcalA_To_kJNm( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getKcalA_To_kJNm:";
// ---------------------------------------------------------------------------------------
return _kcalA_To_kJNm;
}
/**---------------------------------------------------------------------------------------
Set KcalA_To_kJNm
@param kcalA_To_kJNm probe radius
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setKcalA_To_kJNm( RealOpenMM kcalA_To_kJNm ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setKcalA_To_kJNm:";
// ---------------------------------------------------------------------------------------
_kcalA_To_kJNm = kcalA_To_kJNm;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atomic radii
--------------------------------------------------------------------------------------- */
RealOpenMM* ImplicitSolventParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
if( _atomicRadii == NULL ){
ImplicitSolventParameters* localThis = const_cast<ImplicitSolventParameters* const>(this);
localThis->_atomicRadii = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownAtomicRadii = true;
memset( localThis->_atomicRadii, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _atomicRadii;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
if( _ownAtomicRadii && atomicRadii != _atomicRadii ){
delete[] _atomicRadii;
_ownAtomicRadii = false;
}
_atomicRadii = atomicRadii;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nImplicitSolventParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
int numberOfAtoms = getNumberOfAtoms();
if( _ownAtomicRadii ){
delete[] _atomicRadii;
_atomicRadii = new RealOpenMM[numberOfAtoms];
} else if( _atomicRadii == NULL ){
_atomicRadii = new RealOpenMM[numberOfAtoms];
_ownAtomicRadii = true;
}
if( numberOfAtoms != (int) atomicRadii.size() ){
std::stringstream message;
message << methodName;
message << " number of object atoms=" << numberOfAtoms << " does not match number in input vector=" << atomicRadii.size();
SimTKOpenMMLog::printWarning( message );
numberOfAtoms = numberOfAtoms < (int) atomicRadii.size() ? numberOfAtoms : (int) atomicRadii.size();
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setOwnAtomicRadii( int ownAtomicRadii ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nImplicitSolventParameters::setOwnAtomicRadii:";
// ---------------------------------------------------------------------------------------
_ownAtomicRadii = ownAtomicRadii;
}
/**---------------------------------------------------------------------------------------
Get free array flag -- if set then work arrays are freed when destructor is called
@return freeArrays flag
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::getFreeArrays( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setFreeArrays:";
// ---------------------------------------------------------------------------------------
return _freeArrays;
}
/**---------------------------------------------------------------------------------------
Set free array flag -- if set then work arrays are freed when destructor is called
@param freeArrays flag
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::setFreeArrays( int freeArrays ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setFreeArrays:";
// ---------------------------------------------------------------------------------------
_freeArrays = freeArrays;
}
/**---------------------------------------------------------------------------------------
Initialize ImplicitSolvent Parameters (Simbios)
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::_initializeImplicitSolventConstants( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::_initializeImplicitSolventConstants:";
// ---------------------------------------------------------------------------------------
_soluteDielectric = (RealOpenMM) 1.0;
_solventDielectric = (RealOpenMM) 78.3;
_kcalA_To_kJNm = (RealOpenMM) 0.4184;
_probeRadius = (RealOpenMM) 0.14;
_electricConstant = (RealOpenMM) (-0.5*ONE_4PI_EPS0);
//_pi4Asolv = (RealOpenMM) PI_M*4.0*0.0049*1000.0;
//_pi4Asolv = (RealOpenMM) PI_M*19.6;
//_pi4Asolv = (RealOpenMM) PI_M*4.0*0.0054;
_pi4Asolv = (RealOpenMM) 28.3919551;
//_pi4Asolv = (RealOpenMM) -400.71504079;
//_pi4Asolv = (RealOpenMM) 0.0;
_resetPreFactor();
}
/**---------------------------------------------------------------------------------------
Reset prefactor (Simbios)
called when _electricConstant, _soluteDielectric, or _solventDielectric are modified
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::_resetPreFactor( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::_resetPreFactor:";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
if( getSoluteDielectric() != zero && getSolventDielectric() != zero ){
_preFactor = two*getElectricConstant()*( (one/getSoluteDielectric()) - (one/getSolventDielectric()) );
} else {
_preFactor = zero;
}
}
/**---------------------------------------------------------------------------------------
Get state (Simbios)
@param title title (optional)
--------------------------------------------------------------------------------------- */
std::string ImplicitSolventParameters::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
if( title ){
message << title;
}
std::string tab = getStringTab();
message << "\nImplicitSolventParameters :";
message << tab << "Solute dielectric: " << getSoluteDielectric();
message << tab << "Solvent dielectric: " << getSolventDielectric();
message << tab << "Electric constant: " << getElectricConstant();
message << tab << "Probe radius: " << getProbeRadius();
message << tab << "Number of atoms: " << getNumberOfAtoms();
message << tab << "Free arrays: " << getFreeArrays();
return message.str();
}
/**---------------------------------------------------------------------------------------
Get string tab -- centralized
@return tab string
--------------------------------------------------------------------------------------- */
std::string ImplicitSolventParameters::getStringTab( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getStringTab";
// ---------------------------------------------------------------------------------------
return std::string( "\n " );
}
platforms/reference/src/gbsa/ImplicitSolventParameters.h deleted 100755 → 0
View file @ a587c652
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ImplicitSolventParameters_H__
#define __ImplicitSolventParameters_H__
#include "../SimTKUtilities/SimTKOpenMMRealType.h"
#include "../../../../openmmapi/include/openmm/internal/windowsExport.h"
#include <string>
// ---------------------------------------------------------------------------------------
class OPENMM_EXPORT ImplicitSolventParameters {
protected:
// parameters common to implicit solvent models
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _kcalA_To_kJNm;
RealOpenMM _electricConstant;
RealOpenMM _probeRadius;
RealOpenMM _pi4Asolv;
RealOpenMM _preFactor;
// ---------------------------------------------------------------------------------------
// atom count
int _numberOfAtoms;
// flag signalling whether arrays are to be freed
int _freeArrays;
// atomic radii
int _ownAtomicRadii;
RealOpenMM* _atomicRadii;
/**---------------------------------------------------------------------------------------
Initialize ImplicitSolvent Parameters (Simbios)
--------------------------------------------------------------------------------------- */
void _initializeImplicitSolventConstants( void );
/**---------------------------------------------------------------------------------------
Set KcalA_To_kJNm
@param kcalA_To_kJNm probe radius
--------------------------------------------------------------------------------------- */
void setKcalA_To_kJNm( RealOpenMM kcalA_To_kJNm );
/**---------------------------------------------------------------------------------------
Set free array flag -- if set then work arrays are freed when destructor is called
@param freeArrays flag
--------------------------------------------------------------------------------------- */
void setFreeArrays( int freeArrays );
/**---------------------------------------------------------------------------------------
Reset prefactor (Simbios)
called when _electricConstant, _soluteDielectric, or _solventDielectric are modified
--------------------------------------------------------------------------------------- */
void _resetPreFactor( void );
public:
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters( int numberOfAtoms );
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
virtual ~ImplicitSolventParameters( );
// override of new/delete
//static void* operator new( size_t size );
//static void operator delete( void *p );
//static void* operator new[]( size_t size );
//static void operator delete[]( void *p );
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int getNumberOfAtoms( void ) const;
/**---------------------------------------------------------------------------------------
Get free array flag -- if set then work arrays are freed when destructor is called
@return freeArrays flag
--------------------------------------------------------------------------------------- */
int getFreeArrays( void ) const;
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSolventDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void setSolventDielectric( RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSoluteDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void setSoluteDielectric( RealOpenMM soluteDielectric );
/**---------------------------------------------------------------------------------------
Get conversion factor for kcal/A -> kJ/nm (Simbios)
@return kcalA_To_kJNm factor
--------------------------------------------------------------------------------------- */
RealOpenMM getKcalA_To_kJNm( void ) const;
/**---------------------------------------------------------------------------------------
Get electric constant (Simbios)
@return electricConstant
--------------------------------------------------------------------------------------- */
RealOpenMM getElectricConstant( void ) const;
/**---------------------------------------------------------------------------------------
Set electric constant (Simbios)
@param electricConstant electric constant
--------------------------------------------------------------------------------------- */
void setElectricConstant( RealOpenMM electricConstant );
/**---------------------------------------------------------------------------------------
Get probe radius (Simbios)
@return probeRadius
--------------------------------------------------------------------------------------- */
RealOpenMM getProbeRadius( void ) const;
/**---------------------------------------------------------------------------------------
Set probe radius (Simbios)
@param probeRadius probe radius
--------------------------------------------------------------------------------------- */
void setProbeRadius( RealOpenMM probeRadius );
/**---------------------------------------------------------------------------------------
Get pi4Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049f*1000.0f; (Simbios)
@return pi4Asolv
--------------------------------------------------------------------------------------- */
RealOpenMM getPi4Asolv( void ) const;
/**---------------------------------------------------------------------------------------
Get prefactor
@returnpreFactor
--------------------------------------------------------------------------------------- */
RealOpenMM getPreFactor( void ) const;
/**---------------------------------------------------------------------------------------
Set values used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049f*1000.0f; (Simbios)
@param pi4Asolv see above
--------------------------------------------------------------------------------------- */
void setPi4Asolv( RealOpenMM pi4Asolv );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atom volumes
--------------------------------------------------------------------------------------- */
virtual RealOpenMM* getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
--------------------------------------------------------------------------------------- */
virtual void setAtomicRadii( RealOpenMM* atomicRadii );
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
virtual void setAtomicRadii( const RealOpenMMVector& atomicRadii );
/**---------------------------------------------------------------------------------------
Set flag indicating whether AtomicRadii array is owned by
object; if flag is set, then when the object is deleted,
the array will be freed
@param ownAtomicRadii flag indicating whether array of atomic radii
should be freed
--------------------------------------------------------------------------------------- */
void setOwnAtomicRadii( int ownAtomicRadii );
/**---------------------------------------------------------------------------------------
Print state to log file (Simbios)
@param title title (optional)
--------------------------------------------------------------------------------------- */
virtual std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Get string tab -- centralized
@return tab string
--------------------------------------------------------------------------------------- */
std::string getStringTab( void ) const;
};
// ---------------------------------------------------------------------------------------
#endif // __ImplicitSolventParameters_H__
platforms/reference/src/gbsa/ObcParameters.cpp
View file @ 326627a8
......@@ -26,597 +26,385 @@
#include <sstream>
#include <string.h>
#include "openmm/OpenMMException.h"
#include "ObcParameters.h"
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
using std::vector;
using OpenMM::RealVec;
const std::string ObcParameters::ParameterFileName = std::string( "params.agb" );
/**---------------------------------------------------------------------------------------
ObcParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in Obc equation gPol
(4) list of atoms that should be excluded in calculating
force -- nonbonded atoms (1-2, and 1-3 atoms)
Implementation:
Slightly different sequence of calls when running on CPU vs GPU.
Difference arise because the CPU-side data arrays for the Brook
streams are allocated by the BrookStreamWrapper objects. These
arrays are then used by ObcParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, log );
obcParameters->initializeParameters( top );
Gpu:
obcParameters = new ObcParameters( gpu->natoms, log );
// set arrays for cpu using stream data field;
// initializeParameters() only allocates space for arrays if they are not set (==NULL)
// also set flag so that ObcParameters destructor does not free arrays
obcParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuObc::obcVdwRadii )->getData() );
obcParameters->setVolume( getBrookStreamWrapperAtIndex( GpuObc::obcVolume )->getData() );
obcParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuObc::obcGpolFixed )->getData() );
obcParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuObc::obcBornRadii )->getData() );
obcParameters->setFreeArrays( false );
obcParameters->initializeParameters( top );
Issues:
Tinker's atom radii are used.
The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
only tested for generic proteins
see mapGmxAtomNameToTinkerAtomNumber()
--------------------------------------------------------------------------------------- */
ObcParameters constructor
@param numberOfAtoms number of atoms
@param obcType OBC type (Eq. 7 or 8 in paper)
/**---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
ObcParameters constructor (Simbios)
@param numberOfAtoms number of atoms
@param obcType OBC type (Eq. 7 or 8 in paper)
--------------------------------------------------------------------------------------- */
ObcParameters::ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType ) :
_numberOfAtoms(numberOfAtoms),
_solventDielectric( 78.3 ),
_soluteDielectric( 1.0 ),
_electricConstant(-0.5*ONE_4PI_EPS0),
_probeRadius(0.14),
_pi4Asolv( 28.3919551),
_dielectricOffset( 0.009 ),
_obcType( obcType ),
_cutoff(false),
_periodic(false) {
ObcParameters::ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType ) : ImplicitSolventParameters( numberOfAtoms ), cutoff(false), periodic(false) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::ObcParameters";
// ---------------------------------------------------------------------------------------
_obcType = obcType;
_dielectricOffset = 0.009f;
_ownScaledRadiusFactors = 0;
_scaledRadiusFactors = NULL;
setObcTypeParameters( obcType );
_atomicRadii.resize( numberOfAtoms );
_scaledRadiusFactors.resize( numberOfAtoms );
setObcTypeParameters( obcType );
}
/**---------------------------------------------------------------------------------------
ObcParameters destructor (Simbios)
ObcParameters destructor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
ObcParameters::~ObcParameters( ){
}
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::~ObcParameters";
// ---------------------------------------------------------------------------------------
// in GPU runs, arrays may be 'owned' by BrookStreamWrapper -- hence they should not
// be freed here, i.e., _freeArrays should be 'false'
#ifdef UseGromacsMalloc
/*
if( _freeArrays ){
if( _vdwRadii != NULL ){
save_free( "_vdwRadii", __FILE__, __LINE__, _vdwRadii );
}
} */
#else
/**---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors ){
delete[] _scaledRadiusFactors;
}
/*
if( getFreeArrays() ){
Get number of atoms
} */
@return number of atoms
#endif
--------------------------------------------------------------------------------------- */
int ObcParameters::getNumberOfAtoms( void ) const {
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Get OBC type
Get OBC type
@return OBC type
@return OBC type
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
ObcParameters::ObcType ObcParameters::getObcType( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getObcType:";
// ---------------------------------------------------------------------------------------
return _obcType;
return _obcType;
}
/**---------------------------------------------------------------------------------------
Set OBC type specific parameters
Set OBC type specific parameters
@param obcType OBC type (ObcTypeI or ObcTypeII -- Eq. 7 or 8)
@param obcType OBC type (ObcTypeI or ObcTypeII -- Eq. 7 or 8)
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcParameters::setObcTypeParameters( ObcParameters::ObcType obcType ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::setObcTypeParameters:";
// ---------------------------------------------------------------------------------------
if( obcType == ObcTypeI ){
_alphaObc = 0.8f;
_betaObc = 0.0f;
_gammaObc = 2.91f;
} else {
_alphaObc = 1.0f;
_betaObc = 0.8f;
_gammaObc = 4.85f;
}
_obcType = obcType;
if( obcType == ObcTypeI ){
_alphaObc = 0.8f;
_betaObc = 0.0f;
_gammaObc = 2.91f;
} else {
_alphaObc = 1.0f;
_betaObc = 0.8f;
_gammaObc = 4.85f;
}
_obcType = obcType;
}
/**---------------------------------------------------------------------------------------
Get dielectricOffset
Get dielectricOffset
@return _dielectricOffset
@return _dielectricOffset
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getDielectricOffset( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getDielectricOffset:";
// ---------------------------------------------------------------------------------------
return _dielectricOffset;
return _dielectricOffset;
}
/**---------------------------------------------------------------------------------------
Get alpha OBC (Eqs. 6 & 7) in Proteins paper
Get alpha OBC (Eqs. 6 & 7) in Proteins paper
@return alphaObc
@return alphaObc
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getAlphaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getAlphaObc:";
// ---------------------------------------------------------------------------------------
return _alphaObc;
return _alphaObc;
}
/**---------------------------------------------------------------------------------------
Get beta OBC (Eqs. 6 & 7) in Proteins paper
Get beta OBC (Eqs. 6 & 7) in Proteins paper
@return betaObc
@return betaObc
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getBetaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getBetaObc:";
// ---------------------------------------------------------------------------------------
return _betaObc;
return _betaObc;
}
/**---------------------------------------------------------------------------------------
Get gamma OBC (Eqs. 6 & 7) in Proteins paper
Get gamma OBC (Eqs. 6 & 7) in Proteins paper
@return gammaObc
@return gammaObc
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getGammaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getGammaObc:";
// ---------------------------------------------------------------------------------------
return _gammaObc;
return _gammaObc;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
Get solvent dielectric
@return array of atomic radii
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM* ObcParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
// if dielectric offset applied, then unapply
return atomicRadii;
RealOpenMM ObcParameters::getSolventDielectric( void ) const {
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Set solvent dielectric
@param atomicRadii array of atomic radii
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void ObcParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
void ObcParameters::setSolventDielectric( RealOpenMM solventDielectric ){
_solventDielectric = solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Get solute dielectric
@param atomicRadii vector of atomic radii
@return soluteDielectric
--------------------------------------------------------------------------------------- */
void ObcParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nObcParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
RealOpenMM ObcParameters::getSoluteDielectric( void ) const {
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Return OBC scale factors
If not previously set, allocate space
Set solute dielectric
@return array
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcParameters::getScaledRadiusFactors( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::getScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _scaledRadiusFactors == NULL ){
ObcParameters* localThis = const_cast<ObcParameters* const>(this);
localThis->_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownScaledRadiusFactors = true;
memset( _scaledRadiusFactors, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _scaledRadiusFactors;
void ObcParameters::setSoluteDielectric( RealOpenMM soluteDielectric ){
_soluteDielectric = soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
Get electric constant
@param ownScaledRadiusFactors flag indicating whether scale factors
array should be deleted
@return electricConstant
--------------------------------------------------------------------------------------- */
void ObcParameters::setOwnScaleFactors( int ownScaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownScaledRadiusFactors = ownScaledRadiusFactors;
RealOpenMM ObcParameters::getElectricConstant( void ) const {
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Set OBC scale factors
Get probe radius
@param scaledRadiusFactors scaledRadiusFactors
@return probeRadius
--------------------------------------------------------------------------------------- */
void ObcParameters::setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != scaledRadiusFactors ){
delete[] _scaledRadiusFactors;
_ownScaledRadiusFactors = false;
}
_scaledRadiusFactors = scaledRadiusFactors;
RealOpenMM ObcParameters::getProbeRadius( void ) const {
return _probeRadius;
}
/**---------------------------------------------------------------------------------------
Set OBC scale factors
Set probe radius
@param scaledRadiusFactors scaledRadiusFactors
@param probeRadius probe radius
--------------------------------------------------------------------------------------- */
void ObcParameters::setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != NULL ){
delete[] _scaledRadiusFactors;
}
_ownScaledRadiusFactors = true;
_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
for( int ii = 0; ii < (int) scaledRadiusFactors.size(); ii++ ){
_scaledRadiusFactors[ii] = scaledRadiusFactors[ii];
}
void ObcParameters::setProbeRadius( RealOpenMM probeRadius ){
_probeRadius = probeRadius;
}
/**---------------------------------------------------------------------------------------
Map Gmx atom name to Tinker atom number (Simbios)
@param atomName atom name (CA, HA, ...); upper and lower case should both work
@param log if set, then print error messages to log file
Get pi*4*Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049*1000.0; (Still)
((RealOpenMM) M_PI)*4.0f*0.0054*1000.0; (OBC)
@return Tinker atom number if atom name is valid; else return -1
@return pi4Asolv
--------------------------------------------------------------------------------------- */
int ObcParameters::mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const {
RealOpenMM ObcParameters::getPi4Asolv( void ) const {
return _pi4Asolv;
}
// ---------------------------------------------------------------------------------------
static int mapCreated = 0;
static int atomNameMap[26];
// ---------------------------------------------------------------------------------------
/**---------------------------------------------------------------------------------------
// set up atomNameMap array on first call to this method
Get AtomicRadii array
// atomNameMap[ii] = Tinker atom number
// where ii = (the ASCII index - 65) of the first character in the
// input atom name; name may be lower case
@return array of atomic radii
if( !mapCreated ){
--------------------------------------------------------------------------------------- */
mapCreated = 1;
const RealOpenMMVector& ObcParameters::getAtomicRadii( void ) const {
return _atomicRadii;
}
for( int ii = 0; ii < 26; ii++ ){
atomNameMap[ii] = -1;
}
/**---------------------------------------------------------------------------------------
// H
atomNameMap[7] = 1;
Set AtomicRadii array
// C
atomNameMap[2] = 6;
@param atomicRadii vector of atomic radii
// N
atomNameMap[13] = 7;
--------------------------------------------------------------------------------------- */
// O
atomNameMap[14] = 8;
void ObcParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// S
atomNameMap[18] = 16;
}
if( atomicRadii.size() == _atomicRadii.size() ){
for( unsigned int ii = 0; ii < atomicRadii.size(); ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
} else {
std::stringstream msg;
msg << "ObcParameters: input size for atomic radii does not agree w/ current size: input=";
msg << atomicRadii.size();
msg << " current size=" << _atomicRadii.size();
throw OpenMM::OpenMMException(msg.str());
}
// map first letter in atom name to Tinker atom number
}
int firstAsciiValue = ((int) atomName[0]) - 65;
/**---------------------------------------------------------------------------------------
// check for lower case
Return OBC scale factors
if( firstAsciiValue > 25 ){
firstAsciiValue -= 32;
}
@return array
// validate
--------------------------------------------------------------------------------------- */
if( firstAsciiValue < 0 || firstAsciiValue > 25 ){
if( log != NULL ){
(void) fprintf( log, "Atom name=<%s> unrecognized.", atomName );
}
(void) fprintf( stderr, "Atom name=<%s> unrecognized.", atomName );
return -1;
}
return atomNameMap[firstAsciiValue];
const RealOpenMMVector& ObcParameters::getScaledRadiusFactors( void ) const {
return _scaledRadiusFactors;
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string ObcParameters::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
Set OBC scale factors
// static const char* methodName = "\nObcParameters::getStateString";
@param scaledRadiusFactors scaledRadiusFactors
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
std::stringstream message;
message << ImplicitSolventParameters::getStateString( title );
std::string tab = getStringTab();
if( getObcType() == ObcTypeI ){
message << tab << "OBC type: Type I";
} else {
message << tab << "OBC type: Type II";
}
message << tab << "Alpha: " << getAlphaObc();
message << tab << "Beta: " << getBetaObc();
message << tab << "Gamma: " << getGammaObc();
return message.str();
void ObcParameters::setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors ){
if( scaledRadiusFactors.size() == _scaledRadiusFactors.size() ){
for( unsigned int ii = 0; ii < scaledRadiusFactors.size(); ii++ ){
_scaledRadiusFactors[ii] = scaledRadiusFactors[ii];
}
} else {
std::stringstream msg;
msg << "ObcParameters: input size for scaled radius factors does not agree w/ current size: input=";
msg << scaledRadiusFactors.size();
msg << " current size=" << _scaledRadiusFactors.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
Set the force to use a cutoff.
@param distance the cutoff distance
@param distance the cutoff distance
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcParameters::setUseCutoff( RealOpenMM distance ) {
cutoff = true;
cutoffDistance = distance;
_cutoff = true;
_cutoffDistance = distance;
}
/**---------------------------------------------------------------------------------------
Get whether to use a cutoff.
Get whether to use a cutoff.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
bool ObcParameters::getUseCutoff() {
return cutoff;
return _cutoff;
}
/**---------------------------------------------------------------------------------------
Get the cutoff distance.
Get the cutoff distance.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getCutoffDistance() {
return cutoffDistance;
return _cutoffDistance;
}
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
--------------------------------------------------------------------------------------- */
@param boxSize the X, Y, and Z widths of the periodic box
void ObcParameters::setPeriodic( const OpenMM::RealVec& boxSize ) {
--------------------------------------------------------------------------------------- */
assert(_cutoff);
void ObcParameters::setPeriodic( RealVec& boxSize ) {
assert(boxSize[0] >= 2.0*_cutoffDistance);
assert(boxSize[1] >= 2.0*_cutoffDistance);
assert(boxSize[2] >= 2.0*_cutoffDistance);
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
_periodic = true;
_periodicBoxSize[0] = boxSize[0];
_periodicBoxSize[1] = boxSize[1];
_periodicBoxSize[2] = boxSize[2];
}
/**---------------------------------------------------------------------------------------
Get whether to use periodic boundary conditions.
Get whether to use periodic boundary conditions.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
bool ObcParameters::getPeriodic() {
return periodic;
return _periodic;
}
/**---------------------------------------------------------------------------------------
Get the periodic box dimension
Get the periodic box dimension
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcParameters::getPeriodicBox() {
return periodicBoxSize;
return _periodicBoxSize;
}
platforms/reference/src/gbsa/ObcParameters.h
View file @ 326627a8
......@@ -26,11 +26,10 @@
#define __ObcParameters_H__
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "ImplicitSolventParameters.h"
// ---------------------------------------------------------------------------------------
class ObcParameters : public ImplicitSolventParameters {
class ObcParameters {
public:
......@@ -38,12 +37,18 @@ class ObcParameters : public ImplicitSolventParameters {
enum ObcType { ObcTypeI, ObcTypeII };
static const std::string ParameterFileName;
private:
// OBC constants & parameters
int _numberOfAtoms;
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _electricConstant;
RealOpenMM _probeRadius;
RealOpenMM _pi4Asolv;
RealOpenMM _dielectricOffset;
RealOpenMM _alphaObc;
RealOpenMM _betaObc;
......@@ -52,20 +57,19 @@ class ObcParameters : public ImplicitSolventParameters {
// scaled radius factors (S_kk in HCT paper)
int _ownScaledRadiusFactors;
RealOpenMM* _scaledRadiusFactors;
RealOpenMMVector _atomicRadii;
RealOpenMMVector _scaledRadiusFactors;
// cutoff and periodic boundary conditions
bool cutoff;
bool periodic;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
bool _cutoff;
bool _periodic;
RealOpenMM _periodicBoxSize[3];
RealOpenMM _cutoffDistance;
/**---------------------------------------------------------------------------------------
Set solvent dielectric (Simbios)
Set solvent dielectric
@param dielectricOffset solvent dielectric
......@@ -77,7 +81,7 @@ class ObcParameters : public ImplicitSolventParameters {
/**---------------------------------------------------------------------------------------
ObcParameters constructor (Simbios)
ObcParameters constructor
@param numberOfAtoms number of atoms
......@@ -87,19 +91,102 @@ class ObcParameters : public ImplicitSolventParameters {
/**---------------------------------------------------------------------------------------
ObcParameters destructor (Simbios)
ObcParameters destructor
--------------------------------------------------------------------------------------- */
~ObcParameters( );
// override of new/delete
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int getNumberOfAtoms( void ) const;
/**---------------------------------------------------------------------------------------
Get electric constant
@return electric constant
--------------------------------------------------------------------------------------- */
RealOpenMM getElectricConstant( void ) const;
/**---------------------------------------------------------------------------------------
Get probe radius (Simbios)
@return probeRadius
--------------------------------------------------------------------------------------- */
RealOpenMM getProbeRadius( void ) const;
/**---------------------------------------------------------------------------------------
Set probe radius (Simbios)
@param probeRadius probe radius
--------------------------------------------------------------------------------------- */
void setProbeRadius( RealOpenMM probeRadius );
/**---------------------------------------------------------------------------------------
Get pi4Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049f*1000.0f; (Simbios)
@return pi4Asolv
--------------------------------------------------------------------------------------- */
RealOpenMM getPi4Asolv( void ) const;
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSolventDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void setSolventDielectric( RealOpenMM solventDielectric );
//static void* operator new( size_t size );
//static void operator delete( void *p );
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
//static void* operator new[]( size_t size );
//static void operator delete[]( void *p );
RealOpenMM getSoluteDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void setSoluteDielectric( RealOpenMM soluteDielectric );
/**---------------------------------------------------------------------------------------
......@@ -153,7 +240,7 @@ class ObcParameters : public ImplicitSolventParameters {
/**---------------------------------------------------------------------------------------
Get solvent dielectric (Simbios)
Get solvent dielectric
@return dielectricOffset dielectric offset
......@@ -169,7 +256,7 @@ class ObcParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
const RealOpenMM* getScaledRadiusFactors( void ) const;
const RealOpenMMVector& getScaledRadiusFactors( void ) const;
/**---------------------------------------------------------------------------------------
......@@ -179,41 +266,8 @@ class ObcParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
void setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors );
void setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors );
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors arra should be deleted
@param ownScaledRadiusFactors flag indicating whether scale factors
array should be deleted
--------------------------------------------------------------------------------------- */
void setOwnScaleFactors( int ownScaledRadiusFactors );
/**---------------------------------------------------------------------------------------
Assign standard radii for GB/SA methods other than ACE;
taken from Macromodel and OPLS-AA, except for hydrogens (Simbios)
Logic based on logic in Tinker's ksolv.f
Currently only works for standard amino acid atoms
If invalid atom name is encountered, a message is printed to log file and the
radius for that atom is set to 1.0f
@param numberOfAtoms number of atoms
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
--------------------------------------------------------------------------------------- */
void getMacroModelAtomicRadii( int numberOfAtoms,
char*** atomNames, RealOpenMM* radii, FILE* log );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
......@@ -222,17 +276,7 @@ class ObcParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
RealOpenMM* getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
--------------------------------------------------------------------------------------- */
void setAtomicRadii( RealOpenMM* atomicRadii );
const RealOpenMMVector& getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
......@@ -244,30 +288,6 @@ class ObcParameters : public ImplicitSolventParameters {
void setAtomicRadii( const RealOpenMMVector& atomicRadii );
/**---------------------------------------------------------------------------------------
Map Gmx atom name to Tinker atom number (Simbios)
@param atomName atom name (CA, HA, ...); upper and lower case should both work
@param log if set, then print error messages to log file
return Tinker atom number if atom name is valid; else return -1
--------------------------------------------------------------------------------------- */
int mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
......@@ -305,7 +325,7 @@ class ObcParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
void setPeriodic( OpenMM::RealVec& boxSize );
void setPeriodic( const OpenMM::RealVec& boxSize );
/**---------------------------------------------------------------------------------------
......@@ -325,19 +345,6 @@ class ObcParameters : public ImplicitSolventParameters {
};
/**---------------------------------------------------------------------------------------
Qsort/heapsort integer comparison (Simbios)
@parma a first value to compare
@param b second value to compare
@return -1, 0, 1
--------------------------------------------------------------------------------------- */
int integerComparison( const void *a, const void *b);
// ---------------------------------------------------------------------------------------
#endif // __ObcParameters_H__
platforms/reference/src/gbsa/cpuObcInterface.cpp deleted 100644 → 0
View file @ a587c652
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "cpuObcInterface.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ObcParameters.h"
#include "CpuObc.h"
using namespace OpenMM;
using namespace std;
/**---------------------------------------------------------------------------------------
Setup for Obc calculations from Gromacs
@param numberOfAtoms number of atoms
@param obcScaleFactors array of OBC scale factors (one entry each atom)
@param atomicRadii atomic radii in Angstrom (one entry each atom)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculations
@param soluteDielectric solute dielectric
@param solventDielectric solvent dielectric
@param log log reference -- if NULL, then errors/warnings
output to stderr
The method creates a CpuObc instance -- currently the OBC type II model is the
default (see paper). If the OBC type I model is desired change
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
to
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeI );
The created object is a static member of the class CpuObc;
when the force routine, cpuCalculateObcForces(), is called,
the static object is used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
extern "C" int
cpuSetObcParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* obcScaleFactors,
int includeAceApproximation,
RealOpenMM soluteDielectric, RealOpenMM solventDielectric, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\ncpuSetObcParameters: ";
// ---------------------------------------------------------------------------------------
// set log file if not NULL
if( log ){
SimTKOpenMMLog::setSimTKOpenMMLog( log );
}
// set OBC parameters (Type II)
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
obcParameters->setScaledRadiusFactors( obcScaleFactors );
obcParameters->setAtomicRadii( atomicRadii );
// dielectric constants
obcParameters->setSolventDielectric( solventDielectric );
obcParameters->setSoluteDielectric( soluteDielectric );
// ---------------------------------------------------------------------------------------
// create CpuObc instance that will calculate forces
CpuObc* cpuObc = new CpuObc( obcParameters );
// set static member for subsequent calls to calculate forces/energy
CpuImplicitSolvent::setCpuImplicitSolvent( cpuObc );
// include/do not include ACE approximation (nonpolar solvation)
cpuObc->setIncludeAceApproximation( includeAceApproximation );
// ---------------------------------------------------------------------------------------
// diagnostics
if( log ){
std::string state = cpuObc->getStateString( methodName );
(void) fprintf( log, "\n%s\nDone w/ setup\n", state.c_str() );
(void) fflush( log );
}
// ---------------------------------------------------------------------------------------
return 0;
}
/**---------------------------------------------------------------------------------------
Calculate implicit solvent forces and energy
@param atomCoordinates atom coordinates in Angstrom; format of array is
atomCoordinates[atom][3] in Angstrom
@param partialCharges partial charges
@param forces output forces in kcal/mol.A; format of array is
forces[atom][3]
@param energy output energy in kcal/mol
@param updateBornRadii if set, then Born radii are updated for current configuration;
otherwise radii correspond to configuration from previous iteration
Function calls a static method in CpuImplicitSolvent class to calculate forces/energy
--------------------------------------------------------------------------------------- */
extern "C" void
cpuCalculateImplicitSolventForces( vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges,
vector<RealVec>& forces, RealOpenMM* energy,
int updateBornRadii ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\ncpuCalculateImplicitSolventForces: ";
// ---------------------------------------------------------------------------------------
CpuImplicitSolvent::getCpuImplicitSolvent()->computeImplicitSolventForces( atomCoordinates, partialCharges,
forces, updateBornRadii );
*energy = CpuImplicitSolvent::getCpuImplicitSolvent()->getEnergy();
// printf( "\ncpuCalculateImplicitSolventForcesE=%.5e", *energy );
}
/**---------------------------------------------------------------------------------------
Retrieve the calculated energy from the static class member
The energy is calculated in cpuCalculateImplicitSolventForces()
along w/ the forces
@return the calculated energy from the static class member
--------------------------------------------------------------------------------------- */
extern "C" RealOpenMM cpuGetImplicitSolventEnergy( void ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\ncpuGetImplicitSolventEnergy: ";
// ---------------------------------------------------------------------------------------
RealOpenMM energy = CpuImplicitSolvent::getCpuImplicitSolvent()->getEnergy();
// printf( "\ncpuGetImplicitSolventEnergy E=%.5e", energy );
return energy;
}
/**---------------------------------------------------------------------------------------
Delete the Obc associated object(s)
@return 0 if static CpuObc object was set; else return -1
--------------------------------------------------------------------------------------- */
extern "C" int cpuDeleteObcParameters( void ){
return CpuImplicitSolvent::deleteCpuImplicitSolvent();
}
/**---------------------------------------------------------------------------------------
Get OBC scale factors given masses
@param numberOfAtoms number of atoms
@param masses input masses
@param scaleFactors output atomic numbers
--------------------------------------------------------------------------------------- */
extern "C" void getObcScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses,
RealOpenMM* scaleFactors ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetObcScaleFactorsGivenAtomMasses";
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
double scaleFactor = 0.8;
RealOpenMM mass = masses[atomI];
if ( mass < 1.2 && mass >= 1.0 ){ // hydrogen
scaleFactor = 0.85;
} else if( mass > 11.8 && mass < 12.2 ){ // carbon
scaleFactor = 0.72;
} else if( mass > 14.0 && mass < 15.0 ){ // nitrogen
scaleFactor = 0.79;
} else if( mass > 15.5 && mass < 16.5 ){ // oxygen
scaleFactor = 0.85;
} else if( mass > 31.5 && mass < 32.5 ){ // sulphur
scaleFactor = 0.96;
} else if( mass > 29.5 && mass < 30.5 ){ // phosphorus
scaleFactor = 0.86;
} else {
std::stringstream message;
message << methodName;
message << " Warning: mass for atom " << atomI << " mass=" << mass << "> not recognized.";
SimTKOpenMMLog::printMessage( message );
}
scaleFactors[atomI] = (RealOpenMM) scaleFactor;
}
}
/**---------------------------------------------------------------------------------------
Get OBC scale factors given atomic numbers
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param scaleFactors output atomic numbers
--------------------------------------------------------------------------------------- */
extern "C" void getObcScaleFactors( int numberOfAtoms, const int* atomicNumber, RealOpenMM* scaleFactors ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetObcScaleFactors";
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
double scaleFactor;
switch( atomicNumber[atomI] ){
case 1: // hydrogen
scaleFactor = 0.85;
break;
case 6: // carbon
scaleFactor = 0.72;
break;
case 7: // nitrogen
scaleFactor = 0.79;
break;
case 8: // oxygen
scaleFactor = 0.85;
break;
case 15: // phosphorus
scaleFactor = 0.86;
break;
case 16: // sulphur
scaleFactor = 0.85;
break;
default:
scaleFactor = 0.8;
std::stringstream message;
message << methodName;
message << " Warning: atom number=" << atomicNumber[atomI] << " for atom " << atomI << "> not handled -- useing default value.";
SimTKOpenMMLog::printMessage( message );
break;
}
scaleFactors[atomI] = (RealOpenMM) scaleFactor;
}
}
/**---------------------------------------------------------------------------------------
Get GBSA radii
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param numberOfCovalentPartners input number of covalent partners for each atom
1 for H, 2,3,4 for C, ...;
the values are only used for C, N & O
@param indexOfCovalentPartner index of covalent partner -- used only for H
e.g. if atom 22 is a H and it is bonded to atom 24,
then indexOfCovalentPartner[22] = 24
@param gbsaRadii output GBSA radii
--------------------------------------------------------------------------------------- */
extern "C" void getGbsaRadii( int numberOfAtoms, const int* atomicNumber,
const int* numberOfCovalentPartners,
const int* indexOfCovalentPartner,
RealOpenMM* gbsaRadii ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetGbsaRadii";
// ---------------------------------------------------------------------------------------
// loop over atoms
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// branch based on atomic number
double radius;
int atomIndex;
switch ( atomicNumber[atomI] ){
case 1: // H
// get index of covalent partner and validate
// radius is modified if heavy atom is N or O
atomIndex = indexOfCovalentPartner[atomI];
if( atomIndex < 0 || atomIndex >= numberOfAtoms ){
std::stringstream message;
message << methodName;
message << " Warning: atomIndex for covalent partner=" << atomIndex << " for atom " << atomI << " is invalid.";
SimTKOpenMMLog::printMessage( message );
} else if( atomicNumber[atomIndex] == 7 ){
radius = 1.15;
} else if( atomicNumber[atomIndex] == 8 ){
radius = 1.05;
} else {
radius = 1.25;
}
break;
case 3: // Li
radius = 2.432;
break;
case 6: // C
if( numberOfCovalentPartners[atomI] == 2 ){
radius = 1.825;
} else if( numberOfCovalentPartners[atomI] == 3 ){
radius = 1.875;
} else {
radius = 1.90;
}
break;
case 7: // N
if( numberOfCovalentPartners[atomI] == 4 ){
radius = 1.625;
} else if( numberOfCovalentPartners[atomI] == 1 ){
radius = 1.60;
} else {
radius = 1.7063;
}
break;
case 8: // O
if( numberOfCovalentPartners[atomI] == 1 ){
radius = 1.48;
} else {
radius = 1.535;
}
break;
case 9:
radius = 1.47;
break;
case 10:
radius = 1.39;
break;
case 11:
radius = 1.992;
break;
case 12:
radius = 1.70;
break;
case 14:
radius = 1.80;
break;
case 15:
radius = 1.87;
break;
case 16:
radius = 1.775;
break;
case 17:
radius = 1.735;
break;
case 18:
radius = 1.70;
break;
case 19:
radius = 2.123;
break;
case 20:
radius = 1.817;
break;
case 35:
radius = 1.90;
break;
case 36:
radius = 1.812;
break;
case 37:
radius = 2.26;
break;
case 53:
radius = 2.10;
break;
case 54:
radius = 1.967;
break;
case 55:
radius = 2.507;
break;
case 56:
radius = 2.188;
break;
default:
radius = 2.0;
break;
}
gbsaRadii[atomI] = (RealOpenMM) radius;
}
}
platforms/reference/src/gbsa/cpuObcInterface.h deleted 100644 → 0
View file @ a587c652
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __CpuObcInterface_H__
#define __CpuObcInterface_H__
#ifdef __cplusplus
#define externC extern "C"
#else
#define externC extern
#endif
#include "../SimTKUtilities/RealVec.h"
#include <stdio.h>
#include <vector>
/**---------------------------------------------------------------------------------------
Retrieve the calculated implicit solvation energy from the static class member
@return the calculated energy from the static class member
--------------------------------------------------------------------------------------- */
externC RealOpenMM cpuGetImplicitSolventEnergy( void );
/**---------------------------------------------------------------------------------------
Delete the Obc associated object(s)
@return 0 if static CpuObc object was set; else return -1
--------------------------------------------------------------------------------------- */
externC int cpuDeleteObcParameters( void );
/**---------------------------------------------------------------------------------------
Setup for Obc calculations from Gromacs
@param numberOfAtoms number of atoms
@param obcScaleFactors array of OBC scale factors (one entry each atom)
@param atomicRadii atomic radii in Angstrom (one entry each atom)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculations
@param soluteDielectric solute dielectric
@param solventDielectric solvent dielectric
@param log log reference -- if NULL, then errors/warnings
output to stderr
The method creates a CpuObc instance -- currently the OBC type II model is the
default (see paper). If the OBC type I model is desired change
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
to
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeI );
The created object is a static member of the class CpuObc;
when the force routine, cpuCalculateObcForces(), is called,
the static object is used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
externC
int cpuSetObcParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* obcScaleFactors,
int includeAceApproximation, RealOpenMM soluteDielectric,
RealOpenMM solventDielectric, FILE* log );
/**---------------------------------------------------------------------------------------
Calculate implicit solvent forces and energy
@param atomCoordinates atom coordinates in Angstrom; format of array is
atomCoordinates[atom][3]
@param partialCharges partial atom charges
@param forces output forces in kcal/mol.A; format of array is
forces[atom][3]
@param energy energy
@param updateBornRadii if set, then Born radii are updated for current configuration;
otherwise radii correspond to configuration from previous iteration
Function calls a static method in CpuImplicitSolvent class to calculate forces/energy
--------------------------------------------------------------------------------------- */
externC void cpuCalculateImplicitSolventForces( std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialChargesIn,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* energy,
int updateBornRadii );
/**---------------------------------------------------------------------------------------
Get OBC scale factors given masses
@param numberOfAtoms number of atoms
@param masses input masses
@param scaleFactors output atomic numbers
--------------------------------------------------------------------------------------- */
externC void getObcScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses,
RealOpenMM* scaleFactors );
/**---------------------------------------------------------------------------------------
Get OBC scale factors given atomic numbers
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param scaleFactors output atomic numbers
--------------------------------------------------------------------------------------- */
externC void getObcScaleFactors( int numberOfAtoms, const int* atomicNumber,
RealOpenMM* scaleFactors );
/**---------------------------------------------------------------------------------------
Get GBSA radii
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param numberOfCovalentPartners input number of covalent partners for each atom
1 for H, ...; only used for C, N & O
@param indexOfCovalentPartner index of covalent partner -- only used for H
e.g. if atom 22 is a H and it is bonded to atom 24
then indexOfCovalentPartner[22] = 24
@param gbsaRadii output GBSA radii
--------------------------------------------------------------------------------------- */
externC void getGbsaRadii( int numberOfAtoms, const int* atomicNumber,
const int* numberOfCovalentPartners,
const int* indexOfCovalentPartner, RealOpenMM* gbsaRadii );
#undef externC
#endif
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