Skip to content

GitLab

Commit 4d499572 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Cleaned up code

parent e9505f5b
Hide whitespace changes
Inline Side-by-side
Showing with 645 additions and 1470 deletions
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.cpp
View file @ 4d499572
......@@ -44,17 +44,12 @@ using OpenMM::RealVec;
--------------------------------------------------------------------------------------- */
CpuGBVISoftcore::CpuGBVISoftcore( ImplicitSolventParameters* gbviParameters ) : CpuImplicitSolvent( gbviParameters ){
CpuGBVISoftcore::CpuGBVISoftcore( GBVISoftcoreParameters* gbviParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::CpuGBVISoftcore";
// ---------------------------------------------------------------------------------------
_initializeGBVISoftcoreDataMembers( );
// ---------------------------------------------------------------------------------------
_gbviParameters = static_cast<GBVISoftcoreParameters*> (gbviParameters);
_gbviParameters = gbviParameters;
_switchDeriviative.resize(_gbviParameters->getNumberOfAtoms());
}
......@@ -65,36 +60,6 @@ CpuGBVISoftcore::CpuGBVISoftcore( ImplicitSolventParameters* gbviParameters ) :
--------------------------------------------------------------------------------------- */
CpuGBVISoftcore::~CpuGBVISoftcore( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::~CpuGBVISoftcore";
// ---------------------------------------------------------------------------------------
delete[] _switchDeriviative;
//if( _gbviParameters != NULL ){
// delete _gbviParameters;
//}
}
/**---------------------------------------------------------------------------------------
Initialize data members
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::_initializeGBVISoftcoreDataMembers( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::initializeDataMembers";
// ---------------------------------------------------------------------------------------
_gbviParameters = NULL;
_switchDeriviative = NULL;
}
/**---------------------------------------------------------------------------------------
......@@ -109,11 +74,7 @@ GBVISoftcoreParameters* CpuGBVISoftcore::getGBVISoftcoreParameters( void ) const
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::getGBVISoftcoreParameters";
// ---------------------------------------------------------------------------------------
return _gbviParameters;
return _gbviParameters;
}
/**---------------------------------------------------------------------------------------
......@@ -125,14 +86,7 @@ GBVISoftcoreParameters* CpuGBVISoftcore::getGBVISoftcoreParameters( void ) const
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::setGBVISoftcoreParameters( GBVISoftcoreParameters* gbviParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::setGBVISoftcoreParameters";
// ---------------------------------------------------------------------------------------
_gbviParameters = gbviParameters;
_gbviParameters = gbviParameters;
}
/**---------------------------------------------------------------------------------------
......@@ -144,37 +98,8 @@ void CpuGBVISoftcore::setGBVISoftcoreParameters( GBVISoftcoreParameters* gbviPar
--------------------------------------------------------------------------------------- */
RealOpenMM* CpuGBVISoftcore::getSwitchDeriviative( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::getSwitchDeriviative";
// ---------------------------------------------------------------------------------------
if( _switchDeriviative == NULL && _gbviParameters != NULL ){
_switchDeriviative = new RealOpenMM[_gbviParameters->getNumberOfAtoms()];
}
return _switchDeriviative;
}
/**---------------------------------------------------------------------------------------
Return switching function derivative
@return array
--------------------------------------------------------------------------------------- */
RealOpenMM* CpuGBVISoftcore::getSwitchDeriviativeConst( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::getSwitchDeriviative";
// ---------------------------------------------------------------------------------------
return _switchDeriviative;
std::vector<RealOpenMM>& CpuGBVISoftcore::getSwitchDeriviative( void ){
return _switchDeriviative;
}
/**---------------------------------------------------------------------------------------
......@@ -189,32 +114,28 @@ RealOpenMM* CpuGBVISoftcore::getSwitchDeriviativeConst( void ) const {
--------------------------------------------------------------------------------------- */
#define GBVISoftcoreDebug 0
void CpuGBVISoftcore::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 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 );
// static const char* methodName = "CpuGBVISoftcore::quinticSpline";
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
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;
}
/**---------------------------------------------------------------------------------------
......@@ -231,73 +152,59 @@ void CpuGBVISoftcore::quinticSpline( RealOpenMM x, RealOpenMM rl, RealOpenMM ru,
--------------------------------------------------------------------------------------- */
#define GBVISoftcoreDebug 0
void CpuGBVISoftcore::computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadius3, RealOpenMM bornSum,
GBVISoftcoreParameters* 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;
GBVISoftcoreParameters* gbviParameters,
RealOpenMM& bornRadius, RealOpenMM* switchDeriviative ){
static const char* methodName = "CpuGBVISoftcore::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( GBVISoftcoreDebug == 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->getQuinticUpperSplineLimit();
*switchDeriviative = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
#if( GBVISoftcoreDebug == 1 )
(void) fprintf( logFile, " Qv=%14.6e splnDrvtv=%14.6e spline[%10.3e %10.3e] ", splineValue, splineDerivative,
splineL, gbviParameters->getQuinticUpperSplineLimit() );
#endif
} else {
sum = gbviParameters->getQuinticUpperSplineLimit();
*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->getQuinticUpperSplineLimit();
*switchDeriviative = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
} else {
sum = gbviParameters->getQuinticUpperSplineLimit();
*switchDeriviative = zero;
}
} else {
sum = atomicRadius3 - bornSum;
*switchDeriviative = one;
}
bornRadius = POW( sum, minusOneThird );
}
#undef GBVISoftcoreDebug
/**---------------------------------------------------------------------------------------
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
......@@ -308,105 +215,71 @@ void CpuGBVISoftcore::computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadiu
--------------------------------------------------------------------------------------- */
#define GBVISoftcoreDebug 0
void CpuGBVISoftcore::computeBornRadii( vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii ){
return computeBornRadii( atomCoordinates, bornRadii, NULL );
}
void CpuGBVISoftcore::computeBornRadii( vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii, RealOpenMM* switchDeriviative ){
// ---------------------------------------------------------------------------------------
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;
void CpuGBVISoftcore::computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii ){
static const char* methodName = "CpuGBVISoftcore::computeBornRadii";
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
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 );
GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
const RealOpenMM* bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
// ---------------------------------------------------------------------------------------
if( switchDeriviative == NULL ){
switchDeriviative = getSwitchDeriviative();
}
GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const std::vector<RealOpenMM>& atomicRadii = gbviParameters->getAtomicRadii();
const std::vector<RealOpenMM>& scaledRadii = gbviParameters->getScaledRadii();
const std::vector<RealOpenMM>& bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
// ---------------------------------------------------------------------------------------
std::vector<RealOpenMM>& switchDeriviative = getSwitchDeriviative();
#if( GBVISoftcoreDebug == 1 )
FILE* logFile = stderr;
(void) fprintf( logFile, "\n%s\n", methodName );
#endif
// ---------------------------------------------------------------------------------------
// calculate Born radii
// calculate Born radii
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM sum = zero;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM sum = zero;
// sum over volumes
// sum over volumes
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
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 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];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
continue;
sum += bornRadiusScaleFactors[atomJ]*CpuGBVISoftcore::getVolume( r, radiusI, scaledRadii[atomJ] );
#if( GBVISoftcoreDebug == -1 )
if( atomI == 0 || atomI == 1 ){
(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( GBVISoftcoreDebug == 1 )
(void) fprintf( logFile, "%6d BornSum=%14.6e r=%14.6e r3=%14.6e (r3-sum)=%14.6e method=%d ",
atomI, sum, radiusI, POW( radiusI, minusThree ), (POW( radiusI, minusThree ) - sum), _gbviParameters->getBornRadiusScalingSoftcoreMethod() );
#endif
RealOpenMM atomicRadius3 = POW( radiusI, minusThree );
if( _gbviParameters->getBornRadiusScalingSoftcoreMethod() == GBVISoftcoreParameters::NoScaling ){
sum = atomicRadius3 - sum;
bornRadii[atomI] = POW( sum, minusOneThird );
switchDeriviative[atomI] = one;
} else if( _gbviParameters->getBornRadiusScalingSoftcoreMethod() == GBVISoftcoreParameters::QuinticSpline ){
computeBornRadiiUsingQuinticSpline( atomicRadius3, sum, gbviParameters,
bornRadii[atomI], switchDeriviative + atomI );
}
#if( GBVISoftcoreDebug == 1 )
(void) fprintf( logFile, "br=%14.6e swDrvtv=%14.6e %s\n", bornRadii[atomI], switchDeriviative[atomI], (fabs( switchDeriviative[atomI] - 1.0 ) > 1.0e-05 ? "SWWWWW" : "") );
#endif
}
if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
continue;
sum += bornRadiusScaleFactors[atomJ]*CpuGBVISoftcore::getVolume( r, radiusI, scaledRadii[atomJ] );
}
}
RealOpenMM atomicRadius3 = POW( radiusI, minusThree );
if( _gbviParameters->getBornRadiusScalingSoftcoreMethod() == GBVISoftcoreParameters::NoScaling ){
sum = atomicRadius3 - sum;
bornRadii[atomI] = POW( sum, minusOneThird );
switchDeriviative[atomI] = one;
} else if( _gbviParameters->getBornRadiusScalingSoftcoreMethod() == GBVISoftcoreParameters::QuinticSpline ){
RealOpenMM switchDeriviativeValue;
computeBornRadiiUsingQuinticSpline( atomicRadius3, sum, gbviParameters,
bornRadii[atomI], &switchDeriviativeValue );
switchDeriviative[atomI] = switchDeriviativeValue;
}
}
}
#undef GBVISoftcoreDebug
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
......@@ -421,30 +294,28 @@ if( atomI == 0 || atomI == 1 ){
RealOpenMM CpuGBVISoftcore::getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::getVolume";
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM minusThree = (RealOpenMM) -3.0;
RealOpenMM diff = (S - R);
if( FABS( diff ) < r ){
RealOpenMM diff = (S - R);
if( FABS( diff ) < r ){
RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
return (CpuGBVISoftcore::getL( r, (r + S), S ) -
CpuGBVISoftcore::getL( r, lowerBound, S ));
return (CpuGBVISoftcore::getL( r, (r + S), S ) -
CpuGBVISoftcore::getL( r, lowerBound, S ));
} else if( r <= diff ){
} else if( r <= diff ){
return CpuGBVISoftcore::getL( r, (r + S), S ) -
CpuGBVISoftcore::getL( r, (r - S), S ) +
POW( R, minusThree );
return CpuGBVISoftcore::getL( r, (r + S), S ) -
CpuGBVISoftcore::getL( r, (r - S), S ) +
POW( R, minusThree );
} else {
return zero;
}
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
......@@ -461,26 +332,24 @@ RealOpenMM CpuGBVISoftcore::getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S
RealOpenMM CpuGBVISoftcore::getL( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::getL";
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 );
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;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
RealOpenMM diff2 = (r + S)*(r - S);
return (threeHalves*xInv2)*( (fourth*rInv) - (xInv*third) + (diff2*xInv2*eighth*rInv) );
return (threeHalves*xInv2)*( (fourth*rInv) - (xInv*third) + (diff2*xInv2*eighth*rInv) );
}
/**---------------------------------------------------------------------------------------
......@@ -497,29 +366,27 @@ RealOpenMM CpuGBVISoftcore::getL( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
RealOpenMM CpuGBVISoftcore::dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGBVISoftcore::dL_dr";
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 );
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;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
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 );
}
/**---------------------------------------------------------------------------------------
......@@ -536,26 +403,24 @@ RealOpenMM CpuGBVISoftcore::dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
RealOpenMM CpuGBVISoftcore::dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::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) );
}
/**---------------------------------------------------------------------------------------
......@@ -570,21 +435,17 @@ RealOpenMM CpuGBVISoftcore::dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
RealOpenMM CpuGBVISoftcore::Sgb( RealOpenMM t ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::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 GBVISoftcoreDebug 0
/**---------------------------------------------------------------------------------------
Get GB/VI energy
......@@ -600,100 +461,71 @@ RealOpenMM CpuGBVISoftcore::Sgb( RealOpenMM t ){
RealOpenMM CpuGBVISoftcore::computeBornEnergy( const vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
static const char* methodName = "CpuGBVISoftcore::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 GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const RealOpenMM preFactor = gbviParameters->getElectricConstant();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const std::vector<RealOpenMM>& atomicRadii = gbviParameters->getAtomicRadii();
const std::vector<RealOpenMM>& gammaParameters = gbviParameters->getGammaParameters();
const GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const RealOpenMM preFactor = gbviParameters->getElectricConstant();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* gammaParameters = gbviParameters->getGammaParameters();
// ---------------------------------------------------------------------------------------
#if( GBVISoftcoreDebug == 1 )
FILE* logFile = stderr;
(void) fprintf( logFile, "\n%s\n", methodName );
(void) fflush( logFile );
#endif
// 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;
// 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
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = partialCharges[atomI];
RealOpenMM energy = zero;
RealOpenMM cavityEnergy = zero;
// self-energy term
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 = -partialChargeI2*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( GBVISoftcoreDebug == 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 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];
}
#undef GBVISoftcoreDebug
energy += two*partialChargeI*atomIEnergy;
}
energy *= preFactor;
energy -= cavityEnergy;
#define GBVISoftcoreDebug 0
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
return (conversion*energy);
}
/**---------------------------------------------------------------------------------------
......@@ -708,412 +540,190 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference
void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "CpuGBVISoftcore::computeBornForces";
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;
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
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 );
#if( GBVISoftcoreDebug == 1 || GBVISoftcoreDebug == 2 )
FILE* logFile = stderr;
(void) fprintf( logFile, "\n%s\n", methodName );
(void) fflush( logFile );
#endif
// ---------------------------------------------------------------------------------------
const GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMM* gammaParameters = gbviParameters->getGammaParameters();
const GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const std::vector<RealOpenMM>& atomicRadii = gbviParameters->getAtomicRadii();
const std::vector<RealOpenMM>& gammaParameters = gbviParameters->getGammaParameters();
const RealOpenMM preFactor = two*gbviParameters->getElectricConstant();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// constants
// set energy/forces to zero
const RealOpenMM preFactor = two*gbviParameters->getElectricConstant();
std::vector<RealOpenMM> bornForces( numberOfAtoms, 0.0 );
std::vector<RealVec> forces( numberOfAtoms );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = zero;
forces[atomI][1] = zero;
forces[atomI][2] = zero;
}
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// set energy/forces to zero
// first main loop
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;
}
vector<RealOpenMM>& bornForces = getBornForce();
bornForces.assign(numberOfAtoms, 0.0);
// ---------------------------------------------------------------------------------------
// 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)
// 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 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 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 r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
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);
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 );
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;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
if( atomI != atomJ ){
if( atomI != atomJ ){
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
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;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
}
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
// 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( GBVISoftcoreDebug == 1 )
{
double stupidFactor = three;
RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
int maxPrint = 20;
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMM* switchDeriviative = getSwitchDeriviative();
(void) fprintf( logFile, "F1: Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, 1, gbviParameters->getTau() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI];
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
RealOpenMM bF = bornForces[atomI] + (stupidFactor*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
double xx = switchDeriviative[atomI]*bF*oneThird*b2*b2;
// xx*conversion should agree w/ values pulled out of kReduceGBVISoftcoreBornForces_kernel in kForces.cu
(void) fprintf( logFile, "F1 %6d r/sclR[%14.6e %14.6e] bR=%14.6e bF=%14.6e sw=%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], xx*conversion, switchDeriviative[atomI],
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 );
int clearForces = 0;
if( clearForces ){
(void) fprintf( logFile, "Forces cleared after loop 1\n" );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = 0.0f;
forces[atomI][1] = 0.0f;
forces[atomI][2] = 0.0f;
}
}
}
#endif
#if( GBVISoftcoreDebug == 2 )
{
double stupidFactor = three;
RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
int maxPrint = 1000000;
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMM* switchDeriviative = getSwitchDeriviative();
(void) fprintf( logFile, "F1: Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, 1, gbviParameters->getTau() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI];
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
RealOpenMM bF = bornForces[atomI] + (stupidFactor*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
double xx = switchDeriviative[atomI]*bF*oneThird*b2*b2;
// xx*conversion should agree w/ values pulled out of kReduceGBVISoftcoreBornForces_kernel in kForces.cu
/*
(void) fprintf( logFile, "F1 %6d r/sclR[%14.6e %14.6e] bR=%14.6e sw=%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], bF, switchDeriviative[atomI], xx*conversion,
conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*forces[atomI][2],
atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
*/
(void) fprintf( logFile, "%6d %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
atomI, atomicRadii[atomI], scaledRadii[atomI], bornRadii[atomI], bF, xx*conversion,
conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*forces[atomI][2],
atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2], switchDeriviative[atomI] );
}
(void) fflush( logFile );
}
#endif
// second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
// ---------------------------------------------------------------------------------------
// dGpol/dBornRadius) = bornForces[]
// dBornRadius/dr = (1/3)*(bR**4)*(dV/dr)
// 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();
RealOpenMM* switchDeriviative = getSwitchDeriviative();
RealOpenMM stupidFactor = three;
const RealOpenMM* bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
const std::vector<RealOpenMM>& scaledRadii = gbviParameters->getScaledRadii();
std::vector<RealOpenMM>& switchDeriviative = getSwitchDeriviative();
const std::vector<RealOpenMM>& bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[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];
// 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;
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
bornForces[atomI] *= switchDeriviative[atomI]*oneThird*b2*b2;
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
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 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 r = SQRT( r2 );
RealOpenMM S = scaledRadii[atomJ];
RealOpenMM diff = (S - R);
// find dRb/dr, where Rb is the Born radius
RealOpenMM de = CpuGBVISoftcore::dL_dr( r, r+S, S ) + CpuGBVISoftcore::dL_dx( r, r+S, S );
if( FABS( diff ) < r ){
if( R > (r - S) ){
de -= CpuGBVISoftcore::dL_dr( r, R, S );
} else {
de -= ( CpuGBVISoftcore::dL_dr( r, (r-S), S ) + CpuGBVISoftcore::dL_dx( r, (r-S), S ) );
}
} else if( r < (S - R) ){
de -= ( CpuGBVISoftcore::dL_dr( r, r-S, S ) + CpuGBVISoftcore::dL_dx( r, r-S, S ) );
}
RealOpenMM S = scaledRadii[atomJ];
RealOpenMM diff = (S - R);
// find dRb/dr, where Rb is the Born radius
#if 0
for( int kk = 0; kk < 5; kk++ ){
RealOpenMM V1 = CpuGBVISoftcore::getVolume( r, R, S );
RealOpenMM V2 = CpuGBVISoftcore::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 = CpuGBVISoftcore::dL_drD( (double) r, r+S, S ) + CpuGBVISoftcore::dL_dxD( r, r+S, S ) - ( CpuGBVISoftcore::dL_drD( r, (r-S), S ) + CpuGBVISoftcore::dL_dxD( r, (r-S), S ) );
for( int kk = 0; kk < 5; kk++ ){
double V1 = CpuGBVISoftcore::getVolumeD( r, R, S );
double V2 = CpuGBVISoftcore::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 *= bornRadiusScaleFactors[atomJ]*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;
#if 0
if( atomI == 2613 ){
(void) fprintf( stderr, "AtomJ %5d r=%14.7e de=%14.7e bfI=%14.7e finalDe=%14.7e [%14.7e %14.7e %14.7e]\n",
atomJ, r, de, bornForces[atomI], (de*bornForces[atomI]/r),
forces[atomI][0], forces[atomI][1], forces[atomI][2] );
} else if( atomJ == 2613 ){
(void) fprintf( stderr, "AtomI %5d r=%14.7e de=%14.7e bfI=%14.7e finalDe=%14.7e [%14.7e %14.7e %14.7e]\n",
atomI, r, de, bornForces[atomI], (de*bornForces[atomI]/r),
forces[atomJ][0], forces[atomJ][1], forces[atomJ][2] );
}
#endif
}
}
}
#if( GBVISoftcoreDebug == 9 )
{
(void) fprintf( logFile, "\nPre conversion\n" );
(void) fprintf( logFile, "Atom ScaledRadii BornRadii BornForce SwitchDrv Forces\n" );
double forceSum[3] = { 0.0, 0.0, 0.0 };
RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
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, "%6d %14.6e %14.6e %14.6e %14.6e [%14.6e %14.6e %14.6e]\n",
atomI, scaledRadii[atomI], bornRadii[atomI], conversion*bornForces[atomI], switchDeriviative[atomI],
conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*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( GBVISoftcoreDebug == 1 )
{
(void) fprintf( logFile, "\nPost conversion\n" );
(void) fprintf( logFile, "Atom BornRadii BornForce SwitchDrv Forces\n" );
int maxPrint = 20;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%6d %14.6e %14.6e %14.6e 2[%14.6e %14.6e %14.6e] ttlF[%14.6e %14.6e %14.6e] %s\n",
atomI, bornRadii[atomI], conversion*bornForces[atomI], switchDeriviative[atomI],
conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*forces[atomI][2],
inputForces[atomI][0], inputForces[atomI][1], inputForces[atomI][2],
(fabs( switchDeriviative[atomI] - 1.0 ) > 1.0e-05 ? "SWWWWW" : "") );
if( atomI == maxPrint ){
atomI = numberOfAtoms - maxPrint;
if( atomI < maxPrint )atomI = numberOfAtoms;
}
}
(void) fflush( logFile );
}
#endif
#if( GBVISoftcoreDebug == 2 )
{
(void) fprintf( logFile, "\nAtom BornRadii BornForce SwitchDrv Forces Post conversion\n" );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%6d %14.6e %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], switchDeriviative[atomI] );
}
(void) fflush( logFile );
}
#endif
#undef GBVISoftcoreDebug
RealOpenMM de = CpuGBVISoftcore::dL_dr( r, r+S, S ) + CpuGBVISoftcore::dL_dx( r, r+S, S );
if( FABS( diff ) < r ){
if( R > (r - S) ){
de -= CpuGBVISoftcore::dL_dr( r, R, S );
} else {
de -= ( CpuGBVISoftcore::dL_dr( r, (r-S), S ) + CpuGBVISoftcore::dL_dx( r, (r-S), S ) );
}
} else if( r < (S - R) ){
de -= ( CpuGBVISoftcore::dL_dr( r, r-S, S ) + CpuGBVISoftcore::dL_dx( r, r-S, S ) );
}
delete[] forces;
delete[] block;
}
// de = (dG/dRb)(dRb/dr)
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string containing state
--------------------------------------------------------------------------------------- */
de *= bornRadiusScaleFactors[atomJ]*bornForces[atomI]/r;
std::string CpuGBVISoftcore::getStateString( const char* title ) const {
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;
// ---------------------------------------------------------------------------------------
}
}
// static const char* methodName = "\nCpuImplicitSolvent::getStateString";
}
// ---------------------------------------------------------------------------------------
// apply prefactor tau = (1/diel_solute - 1/diel_solvent)
std::stringstream message;
message << CpuImplicitSolvent::getStateString( title );
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];
}
return message.str();
}
/**---------------------------------------------------------------------------------------
......@@ -1132,30 +742,28 @@ std::string CpuGBVISoftcore::getStateString( const char* title ) const {
double CpuGBVISoftcore::getVolumeD( double r, double R, double S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::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 (CpuGBVISoftcore::getLD( r, (r + S), S ) -
CpuGBVISoftcore::getLD( r, lowerBound, S ));
return (CpuGBVISoftcore::getLD( r, (r + S), S ) -
CpuGBVISoftcore::getLD( r, lowerBound, S ));
} else if( r < diff ){
} else if( r < diff ){
return CpuGBVISoftcore::getLD( r, (r + S), S ) -
CpuGBVISoftcore::getLD( r, (r - S), S ) +
pow( R, minusThree );
return CpuGBVISoftcore::getLD( r, (r + S), S ) -
CpuGBVISoftcore::getLD( r, (r - S), S ) +
pow( R, minusThree );
} else {
return zero;
}
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
......@@ -1174,27 +782,25 @@ double CpuGBVISoftcore::getVolumeD( double r, double R, double S ){
double CpuGBVISoftcore::getLD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::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) );
}
/**---------------------------------------------------------------------------------------
......@@ -1215,27 +821,25 @@ double CpuGBVISoftcore::dL_drD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::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 );
}
/**---------------------------------------------------------------------------------------
......@@ -1254,24 +858,22 @@ double CpuGBVISoftcore::dL_drD( double r, double x, double S ){
double CpuGBVISoftcore::dL_dxD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "CpuGBVISoftcore::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) );
}
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.h
View file @ 4d499572
......@@ -26,11 +26,10 @@
#define __CpuGBVISoftcore_H__
#include "GBVISoftcoreParameters.h"
#include "gbsa/CpuImplicitSolvent.h"
// ---------------------------------------------------------------------------------------
class CpuGBVISoftcore : public CpuImplicitSolvent {
class CpuGBVISoftcore {
private:
......@@ -38,14 +37,9 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
GBVISoftcoreParameters* _gbviParameters;
// arrays containing switching function derivative
// vector containing switching function derivative
RealOpenMM* _switchDeriviative;
// initialize data members (more than
// one constructor, so centralize intialization here)
void _initializeGBVISoftcoreDataMembers( void );
std::vector<RealOpenMM> _switchDeriviative;
public:
......@@ -59,7 +53,7 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
CpuGBVISoftcore( ImplicitSolventParameters* gbviParameters );
CpuGBVISoftcore( GBVISoftcoreParameters* gbviParameters );
/**---------------------------------------------------------------------------------------
......@@ -98,8 +92,7 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
RealOpenMM* getSwitchDeriviative( void );
RealOpenMM* getSwitchDeriviativeConst( void ) const;
std::vector<RealOpenMM>& getSwitchDeriviative( void );
/**---------------------------------------------------------------------------------------
......@@ -111,8 +104,6 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii,
RealOpenMM* switchDeriviative = NULL );
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii );
/**---------------------------------------------------------------------------------------
......
plugins/freeEnergy/platforms/reference/src/gbsa/GBVISoftcoreParameters.cpp
View file @ 4d499572
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
* Contributors: Pande Group
*
......@@ -23,74 +22,16 @@
*/
#include <math.h>
#include <iostream>
#include <sstream>
#include <string.h>
#include "openmm/OpenMMException.h"
#include "GBVISoftcoreParameters.h"
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
// #define UseGromacsMalloc 1
#ifdef UseGromacsMalloc
extern "C" {
#include "smalloc.h"
}
#endif
/**---------------------------------------------------------------------------------------
GBVISoftcoreParameters:
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 GBVISoftcoreParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
GBVISoftcoreParameters* gb_VIParameters = new GBVISoftcoreParameters( numberOfAtoms, log );
gb_VIParameters->initializeParameters( top );
Gpu:
gb_VIParameters = new GBVISoftcoreParameters( 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 GBVISoftcoreParameters 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()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
GBVISoftcoreParameters constructor (Simbios)
......@@ -99,450 +40,249 @@ extern "C" {
--------------------------------------------------------------------------------------- */
GBVISoftcoreParameters::GBVISoftcoreParameters( int numberOfAtoms ) : ImplicitSolventParameters( numberOfAtoms ), cutoff(false), periodic(false) {
GBVISoftcoreParameters::GBVISoftcoreParameters( int numberOfAtoms ) : _numberOfAtoms(numberOfAtoms), _soluteDielectric(1.0), _solventDielectric(78.3),
_electricConstant(-0.5*ONE_4PI_EPS0), _quinticLowerLimitFactor(0.8), _bornRadiusScalingSoftcoreMethod(NoScaling),
_cutoff(false), _periodic(false) {
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::GBVISoftcoreParameters";
// ---------------------------------------------------------------------------------------
_ownScaledRadii = 0;
_scaledRadii = NULL;
_ownGammaParameters = 0;
_gammaParameters = NULL;
_ownBornRadiusScaleFactors = 0;
_bornRadiusScaleFactors = NULL;
_bornRadiusScalingSoftcoreMethod = NoScaling;
_quinticLowerLimitFactor = static_cast<RealOpenMM>(0.8);
setQuinticUpperBornRadiusLimit( static_cast<RealOpenMM>(5.0) );
_atomicRadii.resize( numberOfAtoms );
_scaledRadii.resize( numberOfAtoms );
_gammaParameters.resize( numberOfAtoms );
_bornRadiusScaleFactors.resize( numberOfAtoms );
setQuinticUpperBornRadiusLimit( static_cast<RealOpenMM>(5.0) );
}
/**---------------------------------------------------------------------------------------
GBVISoftcoreParameters destructor (Simbios)
GBVISoftcoreParameters destructor
--------------------------------------------------------------------------------------- */
GBVISoftcoreParameters::~GBVISoftcoreParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::~GBVISoftcoreParameters";
// ---------------------------------------------------------------------------------------
// 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;
delete[] _bornRadiusScaleFactors;
/*
if( getFreeArrays() ){
} */
#endif
}
/**---------------------------------------------------------------------------------------
Get the quintic spline lower limit factor
Get number of atoms
@return quintic spline lower limit factor
@return number of atoms
--------------------------------------------------------------------------------------- */
RealOpenMM GBVISoftcoreParameters::getQuinticLowerLimitFactor( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GBVISoftcoreParameters::getQuinticLowerLimitFactor:";
// ---------------------------------------------------------------------------------------
return _quinticLowerLimitFactor;
int GBVISoftcoreParameters::getNumberOfAtoms( void ) const {
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Set the quintic spline lower limit factor
Get electric constant
@param quintic spline lower limit factor
@return electric constant
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setQuinticLowerLimitFactor( RealOpenMM quinticLowerLimitFactor ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GBVISoftcoreParameters::setQuinticLowerLimitFactor:";
// ---------------------------------------------------------------------------------------
_quinticLowerLimitFactor = quinticLowerLimitFactor;
RealOpenMM GBVISoftcoreParameters::getElectricConstant( void ) const {
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Get the quintic spline upper limit
Get solvent dielectric
@return quintic spline upper limit
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM GBVISoftcoreParameters::getQuinticUpperBornRadiusLimit( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GBVISoftcoreParameters::getQuinticUpperBornRadiusLimit:";
// ---------------------------------------------------------------------------------------
return _quinticUpperBornRadiusLimit;
RealOpenMM GBVISoftcoreParameters::getSolventDielectric( void ) const {
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set the quintic spline upper limit
Set solvent dielectric
@param quintic spline upper limit
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setQuinticUpperBornRadiusLimit( RealOpenMM quinticUpperBornRadiusLimit ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GBVISoftcoreParameters::setQuinticUpperBornRadiusLimit:";
// ---------------------------------------------------------------------------------------
_quinticUpperBornRadiusLimit = quinticUpperBornRadiusLimit;
_quinticUpperSplineLimit = POW( _quinticUpperBornRadiusLimit, static_cast<RealOpenMM>(-3.0) );
void GBVISoftcoreParameters::setSolventDielectric( RealOpenMM solventDielectric ){
_solventDielectric = solventDielectric;
}
/**---------------------------------------------------------------------------------------
Get the quintic upper spline limit
Get solute dielectric
@return the quintic upper spline limit
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM GBVISoftcoreParameters::getQuinticUpperSplineLimit( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GBVISoftcoreParameters::getQuinticUpperSplineLimit:";
// ---------------------------------------------------------------------------------------
return _quinticUpperSplineLimit;
RealOpenMM GBVISoftcoreParameters::getSoluteDielectric( void ) const {
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
Set solute dielectric
@return array of atomic radii
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM* GBVISoftcoreParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
return atomicRadii;
void GBVISoftcoreParameters::setSoluteDielectric( RealOpenMM soluteDielectric ){
_soluteDielectric = soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Get the quintic spline lower limit factor
@param atomicRadii array of atomic radii
@return quintic spline lower limit factor
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
RealOpenMM GBVISoftcoreParameters::getQuinticLowerLimitFactor( void ) const {
return _quinticLowerLimitFactor;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Set the quintic spline lower limit factor
@param atomicRadii vector of atomic radii
@param quintic spline lower limit factor
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGBVISoftcoreParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
void GBVISoftcoreParameters::setQuinticLowerLimitFactor( RealOpenMM quinticLowerLimitFactor ){
_quinticLowerLimitFactor = quinticLowerLimitFactor;
}
/**---------------------------------------------------------------------------------------
Return scaled radii
If not previously set, allocate space
Get the quintic spline upper limit
@return array
@return quintic spline upper limit
--------------------------------------------------------------------------------------- */
const RealOpenMM* GBVISoftcoreParameters::getScaledRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::getScaledRadii";
// ---------------------------------------------------------------------------------------
if( _scaledRadii == NULL ){
GBVISoftcoreParameters* localThis = const_cast<GBVISoftcoreParameters* const>(this);
localThis->_scaledRadii = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownScaledRadii = true;
memset( _scaledRadii, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _scaledRadii;
RealOpenMM GBVISoftcoreParameters::getQuinticUpperBornRadiusLimit( void ) const {
return _quinticUpperBornRadiusLimit;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
Set the quintic spline upper limit
@param ownScaledRadii flag indicating whether scale factors
array should be deleted
@param quintic spline upper limit
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setOwnScaledRadii( int ownScaledRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownScaledRadii = ownScaledRadii;
void GBVISoftcoreParameters::setQuinticUpperBornRadiusLimit( RealOpenMM quinticUpperBornRadiusLimit ){
_quinticUpperBornRadiusLimit = quinticUpperBornRadiusLimit;
_quinticUpperSplineLimit = POW( _quinticUpperBornRadiusLimit, static_cast<RealOpenMM>(-3.0) );
}
/**---------------------------------------------------------------------------------------
Set scaled radii
Get the quintic upper spline limit
@param scaledRadii scaledRadii
@return the quintic upper spline limit
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setScaledRadii( RealOpenMM* scaledRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadii";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadii && _scaledRadii != scaledRadii ){
delete[] _scaledRadii;
_ownScaledRadii = false;
}
_scaledRadii = scaledRadii;
RealOpenMM GBVISoftcoreParameters::getQuinticUpperSplineLimit( void ) const {
return _quinticUpperSplineLimit;
}
#if RealOpenMMType == 0
/**---------------------------------------------------------------------------------------
Set scaled radii
Get AtomicRadii array
@param scaledRadii scaledRadii
@return array of atomic radii
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setScaledRadii( float* scaledRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setScaledRadii";
// ---------------------------------------------------------------------------------------
if( _scaledRadii == NULL ){
_scaledRadii = new RealOpenMM[getNumberOfAtoms()];
_ownScaledRadii = true;
}
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_scaledRadii[ii] = (RealOpenMM) scaledRadii[ii];
}
const std::vector<RealOpenMM>& GBVISoftcoreParameters::getAtomicRadii( void ) const {
return _atomicRadii;
}
#endif
/**---------------------------------------------------------------------------------------
Set scaled radii
Set AtomicRadii vector
@param scaledRadii scaledRadii
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setScaledRadii( const RealOpenMMVector& scaledRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadii";
void GBVISoftcoreParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
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];
}
if( atomicRadii.size() == _atomicRadii.size() ){
for( unsigned int ii = 0; ii < atomicRadii.size(); ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
} else {
std::stringstream msg;
msg << "GBVISoftcoreParameters: 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());
}
}
/**---------------------------------------------------------------------------------------
Return gamma parameters
If not previously set, allocate space
Return scaled radii
@return array
--------------------------------------------------------------------------------------- */
RealOpenMM* GBVISoftcoreParameters::getGammaParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::getGammaParameters";
// ---------------------------------------------------------------------------------------
if( _gammaParameters == NULL ){
GBVISoftcoreParameters* localThis = const_cast<GBVISoftcoreParameters* 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 GBVISoftcoreParameters::setOwnGammaParameters( int ownGammaParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownGammaParameters = ownGammaParameters;
const std::vector<RealOpenMM>& GBVISoftcoreParameters::getScaledRadii( void ) const {
return _scaledRadii;
}
/**---------------------------------------------------------------------------------------
Set gamma parameters
Set scaled radii
@param gammas gamma parameters
@param scaledRadii scaledRadii
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setGammaParameters( RealOpenMM* gammas ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setGammas";
void GBVISoftcoreParameters::setScaledRadii( const RealOpenMMVector& scaledRadii ){
// ---------------------------------------------------------------------------------------
if( _ownGammaParameters && _gammaParameters != gammas ){
delete[] _gammaParameters;
_ownGammaParameters = false;
}
_gammaParameters = gammas;
if( scaledRadii.size() == _scaledRadii.size() ){
for( unsigned int ii = 0; ii < scaledRadii.size(); ii++ ){
_scaledRadii[ii] = scaledRadii[ii];
}
} else {
std::stringstream msg;
msg << "GBVISoftcoreParameters: 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());
}
}
#if RealOpenMMType == 0
/**---------------------------------------------------------------------------------------
Set gamma parameters
Return gamma parameters
@param gammas gammas
@return array
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setGammaParameters( float* gammas ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setGammas";
// ---------------------------------------------------------------------------------------
if( _gammaParameters == NULL ){
_gammaParameters = new RealOpenMM[getNumberOfAtoms()];
_ownGammaParameters = true;
}
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_gammaParameters[ii] = (RealOpenMM) gammas[ii];
}
const std::vector<RealOpenMM>& GBVISoftcoreParameters::getGammaParameters( void ) const {
return _gammaParameters;
}
#endif
/**---------------------------------------------------------------------------------------
Set gamma parameters
......@@ -555,123 +295,31 @@ void GBVISoftcoreParameters::setGammaParameters( const RealOpenMMVector& gammas
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setGammas";
// ---------------------------------------------------------------------------------------
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];
}
if( gammas.size() == _gammaParameters.size() ){
for( unsigned int ii = 0; ii < gammas.size(); ii++ ){
_gammaParameters[ii] = gammas[ii];
}
} else {
std::stringstream msg;
msg << "GBVISoftcoreParameters: input size for gammas does not agree w/ current size: input=";
msg << gammas.size();
msg << " current size=" << _gammaParameters.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Return BornRadiusScaleFactors
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
RealOpenMM* GBVISoftcoreParameters::getBornRadiusScaleFactors( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::getBornRadiusScaleFactors";
// ---------------------------------------------------------------------------------------
if( _bornRadiusScaleFactors == NULL ){
GBVISoftcoreParameters* localThis = const_cast<GBVISoftcoreParameters* const>(this);
localThis->_bornRadiusScaleFactors = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownBornRadiusScaleFactors = true;
memset( _bornRadiusScaleFactors, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _bornRadiusScaleFactors;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
@param ownBornRadiusScaleFactors flag indicating whether Born radius scale factors
array should be deleted
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setOwnBornRadiusScaleFactors( int ownBornRadiusScaleFactorsParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownBornRadiusScaleFactors = ownBornRadiusScaleFactorsParameters;
const std::vector<RealOpenMM>& GBVISoftcoreParameters::getBornRadiusScaleFactors( void ) const {
return _bornRadiusScaleFactors;
}
/**---------------------------------------------------------------------------------------
Set BornRadiusScaleFactors
@param bornRadiusScaleFactors Born radius scale factors
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setBornRadiusScaleFactors( RealOpenMM* bornRadiusScaleFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setBornRadiusScaleFactors";
// ---------------------------------------------------------------------------------------
if( _ownBornRadiusScaleFactors && _bornRadiusScaleFactors != bornRadiusScaleFactors ){
delete[] _bornRadiusScaleFactors;
_ownBornRadiusScaleFactors = false;
}
_bornRadiusScaleFactors = bornRadiusScaleFactors;
}
#if RealOpenMMType == 0
/**---------------------------------------------------------------------------------------
Set bornRadiusScaleFactors parameters
@param bornRadiusScaleFactorss bornRadiusScaleFactorss
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setBornRadiusScaleFactors( float* bornRadiusScaleFactorss ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::setBornRadiusScaleFactors";
// ---------------------------------------------------------------------------------------
if( _bornRadiusScaleFactors == NULL ){
_bornRadiusScaleFactors = new RealOpenMM[getNumberOfAtoms()];
_ownBornRadiusScaleFactors = true;
}
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_bornRadiusScaleFactors[ii] = (RealOpenMM) bornRadiusScaleFactorss[ii];
}
return 0;
}
#endif
/**---------------------------------------------------------------------------------------
Set bornRadiusScaleFactors parameters
......@@ -684,46 +332,17 @@ void GBVISoftcoreParameters::setBornRadiusScaleFactors( const RealOpenMMVector&
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setBornRadiusScaleFactors";
// ---------------------------------------------------------------------------------------
if( _ownBornRadiusScaleFactors && _bornRadiusScaleFactors != NULL ){
delete[] _bornRadiusScaleFactors;
}
_ownBornRadiusScaleFactors = true;
_bornRadiusScaleFactors = new RealOpenMM[getNumberOfAtoms()];
for( int ii = 0; ii < (int) bornRadiusScaleFactors.size(); ii++ ){
_bornRadiusScaleFactors[ii] = bornRadiusScaleFactors[ii];
}
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string GBVISoftcoreParameters::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << ImplicitSolventParameters::getStateString( title );
std::string tab = getStringTab();
return message.str();
if( bornRadiusScaleFactors.size() == _bornRadiusScaleFactors.size() ){
for( int ii = 0; ii < (int) bornRadiusScaleFactors.size(); ii++ ){
_bornRadiusScaleFactors[ii] = bornRadiusScaleFactors[ii];
}
} else {
std::stringstream msg;
msg << "GBVISoftcoreParameters: input size for bornRadiusScaleFactors does not agree w/ current size: input=";
msg << bornRadiusScaleFactors.size();
msg << " current size=" << _bornRadiusScaleFactors.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
......@@ -735,9 +354,8 @@ std::string GBVISoftcoreParameters::getStateString( const char* title ) const {
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setUseCutoff( RealOpenMM distance ) {
cutoff = true;
cutoffDistance = distance;
_cutoff = true;
_cutoffDistance = distance;
}
/**---------------------------------------------------------------------------------------
......@@ -747,7 +365,7 @@ void GBVISoftcoreParameters::setUseCutoff( RealOpenMM distance ) {
--------------------------------------------------------------------------------------- */
bool GBVISoftcoreParameters::getUseCutoff() {
return cutoff;
return _cutoff;
}
/**---------------------------------------------------------------------------------------
......@@ -757,7 +375,7 @@ bool GBVISoftcoreParameters::getUseCutoff() {
--------------------------------------------------------------------------------------- */
RealOpenMM GBVISoftcoreParameters::getCutoffDistance() {
return cutoffDistance;
return _cutoffDistance;
}
/**---------------------------------------------------------------------------------------
......@@ -772,14 +390,16 @@ RealOpenMM GBVISoftcoreParameters::getCutoffDistance() {
void GBVISoftcoreParameters::setPeriodic( RealOpenMM* 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];
}
/**---------------------------------------------------------------------------------------
......@@ -789,7 +409,7 @@ void GBVISoftcoreParameters::setPeriodic( RealOpenMM* boxSize ) {
--------------------------------------------------------------------------------------- */
bool GBVISoftcoreParameters::getPeriodic() {
return periodic;
return _periodic;
}
/**---------------------------------------------------------------------------------------
......@@ -799,7 +419,7 @@ bool GBVISoftcoreParameters::getPeriodic() {
--------------------------------------------------------------------------------------- */
const RealOpenMM* GBVISoftcoreParameters::getPeriodicBox() {
return periodicBoxSize;
return _periodicBoxSize;
}
/**---------------------------------------------------------------------------------------
......@@ -812,23 +432,21 @@ const RealOpenMM* GBVISoftcoreParameters::getPeriodicBox() {
RealOpenMM GBVISoftcoreParameters::getTau( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGBVISoftcoreParameters::getTau:";
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;
}
return tau;
// ---------------------------------------------------------------------------------------
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;
}
return tau;
}
/**---------------------------------------------------------------------------------------
......@@ -840,15 +458,7 @@ RealOpenMM GBVISoftcoreParameters::getTau( void ) const {
--------------------------------------------------------------------------------------- */
GBVISoftcoreParameters::BornRadiusScalingSoftcoreMethod GBVISoftcoreParameters::getBornRadiusScalingSoftcoreMethod( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = " GBVISoftcoreParameters::getBornRadiusScalingSoftcoreMethod:";
// ---------------------------------------------------------------------------------------
return _bornRadiusScalingSoftcoreMethod;
return _bornRadiusScalingSoftcoreMethod;
}
/**---------------------------------------------------------------------------------------
......@@ -860,13 +470,5 @@ GBVISoftcoreParameters::BornRadiusScalingSoftcoreMethod GBVISoftcoreParameters::
--------------------------------------------------------------------------------------- */
void GBVISoftcoreParameters::setBornRadiusScalingSoftcoreMethod( BornRadiusScalingSoftcoreMethod bornRadiusScalingSoftcoreMethod ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = " GBVISoftcoreParameters::setBornRadiusScalingSoftcoreMethod:";
// ---------------------------------------------------------------------------------------
_bornRadiusScalingSoftcoreMethod = bornRadiusScalingSoftcoreMethod;
_bornRadiusScalingSoftcoreMethod = bornRadiusScalingSoftcoreMethod;
}
plugins/freeEnergy/platforms/reference/src/gbsa/GBVISoftcoreParameters.h
View file @ 4d499572
......@@ -26,11 +26,10 @@
#define __GBVISoftcoreParameters_H__
#include "SimTKUtilities/SimTKOpenMMCommon.h"
#include "gbsa/ImplicitSolventParameters.h"
// ---------------------------------------------------------------------------------------
class GBVISoftcoreParameters : public ImplicitSolventParameters {
class GBVISoftcoreParameters {
public:
......@@ -42,41 +41,41 @@ class GBVISoftcoreParameters : public ImplicitSolventParameters {
* No scaling method is applied.
*/
NoScaling = 0,
/**
* Use the method outlined in Proteins 55, 383-394 (2004), Eq. 6
*/
Tanh = 1,
/**
* Use quintic spline scaling function
*/
QuinticSpline = 2
QuinticSpline = 1
};
private:
// scaled radii
int _numberOfAtoms;
int _ownScaledRadii;
RealOpenMM* _scaledRadii;
// parameters:
// scaled radii
// gamma parameters
// BornRadiusScaleFactors parameters
// gamma parameters
int _ownGammaParameters;
RealOpenMM* _gammaParameters;
std::vector<RealOpenMM> _scaledRadii;
std::vector<RealOpenMM> _atomicRadii;
std::vector<RealOpenMM> _gammaParameters;
std::vector<RealOpenMM> _bornRadiusScaleFactors;
// BornRadiusScaleFactors parameters
int _ownBornRadiusScaleFactors;
RealOpenMM* _bornRadiusScaleFactors;
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _electricConstant;
// cutoff and periodic boundary conditions
bool cutoff;
bool periodic;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
bool _cutoff;
bool _periodic;
RealOpenMM _periodicBoxSize[3];
RealOpenMM _cutoffDistance;
// Born radii switching function params
BornRadiusScalingSoftcoreMethod _bornRadiusScalingSoftcoreMethod;
RealOpenMM _quinticLowerLimitFactor;
RealOpenMM _quinticUpperBornRadiusLimit;
RealOpenMM _quinticUpperSplineLimit;
......@@ -103,140 +102,133 @@ class GBVISoftcoreParameters : 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 );
#if RealOpenMMType == 0
void setScaledRadii( float* scaledRadii );
#endif
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 std::vector<RealOpenMM>& getScaledRadii( void ) const;
/**---------------------------------------------------------------------------------------
Get GammaParameters array
Set scaled radii
@return array of gamma values
@param vector of scaled radii
--------------------------------------------------------------------------------------- */
RealOpenMM* getGammaParameters( void ) const;
void setScaledRadii( const std::vector<RealOpenMM>& radii );
/**---------------------------------------------------------------------------------------
Set GammaParameters array
Get AtomicRadii array w/ dielectric offset applied
@param gammaParameters array of gamma parameters
@return array of atom volumes
--------------------------------------------------------------------------------------- */
void setGammaParameters( RealOpenMM* gammaParameters );
const std::vector<RealOpenMM>& getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set GammaParameters array
Set AtomicRadii array
@param gammaParameters array of gamma parameters
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
void setGammaParameters( const RealOpenMMVector& gammaParameters );
void setAtomicRadii( const RealOpenMMVector& atomicRadii );
/**---------------------------------------------------------------------------------------
Set flag indicating whether bornRadiusScaleFactor parameter array should be deleted
Get GammaParameters array
@param ownBornRadiusScaleFactors flag indicating whether bornRadiusScaleFactor parameter
array should be deleted
@return array of gamma values
--------------------------------------------------------------------------------------- */
void setOwnBornRadiusScaleFactors( int ownBornRadiusScaleFactors );
const std::vector<RealOpenMM>& getGammaParameters( void ) const;
/**---------------------------------------------------------------------------------------
Get BornRadiusScaleFactors array
Set GammaParameters array
@return array of bornRadiusScaleFactor values
@param gammaParameters array of gamma parameters
--------------------------------------------------------------------------------------- */
RealOpenMM* getBornRadiusScaleFactors( void ) const;
void setGammaParameters( const RealOpenMMVector& gammaParameters );
/**---------------------------------------------------------------------------------------
Set BornRadiusScaleFactors array
Get BornRadiusScaleFactors array
@param bornRadiusScaleFactors array of bornRadiusScaleFactor parameters
@return array of bornRadiusScaleFactor values
--------------------------------------------------------------------------------------- */
void setBornRadiusScaleFactors( RealOpenMM* bornRadiusScaleFactors );
const std::vector<RealOpenMM>& getBornRadiusScaleFactors( void ) const;
/**---------------------------------------------------------------------------------------
......@@ -248,18 +240,6 @@ class GBVISoftcoreParameters : public ImplicitSolventParameters {
void setBornRadiusScaleFactors( const RealOpenMMVector& bornRadiusScaleFactors );
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
......@@ -276,7 +256,7 @@ class GBVISoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
bool getUseCutoff();
bool getUseCutoff( void );
/**---------------------------------------------------------------------------------------
......@@ -284,7 +264,7 @@ class GBVISoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
RealOpenMM getCutoffDistance();
RealOpenMM getCutoffDistance( void );
/**---------------------------------------------------------------------------------------
......@@ -304,7 +284,7 @@ class GBVISoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
bool getPeriodic();
bool getPeriodic( void );
/**---------------------------------------------------------------------------------------
......@@ -312,7 +292,7 @@ class GBVISoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
const RealOpenMM* getPeriodicBox();
const RealOpenMM* getPeriodicBox( void );
/**---------------------------------------------------------------------------------------
......
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