Skip to content

GitLab

Commit 6718e0cc authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Code cleanup

parent 1079a319
Hide whitespace changes
Inline Side-by-side
Showing with 654 additions and 1154 deletions
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.cpp
View file @ 6718e0cc
......@@ -23,15 +23,13 @@
*/
#include <string.h>
#include <math.h>
#include <sstream>
#include <vector>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "CpuGBVISoftcore.h"
#include "../SimTKReference/ReferenceForce.h"
#include <math.h>
#include "CpuGBVISoftcore.h"
using std::vector;
using OpenMM::RealVec;
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.cpp
View file @ 6718e0cc
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
* Contributors: Pande Group
*
......@@ -26,15 +25,14 @@
#include <sstream>
#include <stdlib.h>
//#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "CpuObcSoftcore.h"
#include "../SimTKReference/ReferenceForce.h"
#include <cmath>
#include <cstdio>
#include <vector>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKReference/ReferenceForce.h"
#include "CpuObcSoftcore.h"
using std::vector;
using OpenMM::RealVec;
......@@ -46,661 +44,466 @@ using OpenMM::RealVec;
--------------------------------------------------------------------------------------- */
CpuObcSoftcore::CpuObcSoftcore( ImplicitSolventParameters* obcSoftcoreParameters ) : CpuImplicitSolvent( obcSoftcoreParameters ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::CpuObcSoftcore";
// ---------------------------------------------------------------------------------------
_initializeObcDataMembers( );
_obcSoftcoreParameters = static_cast<ObcSoftcoreParameters*> (obcSoftcoreParameters);
CpuObcSoftcore::CpuObcSoftcore( ObcSoftcoreParameters* obcSoftcoreParameters ){
_obcSoftcoreParameters = obcSoftcoreParameters;
_includeAceApproximation = 1;
_obcChain.resize(_obcSoftcoreParameters->getNumberOfAtoms());
}
/**---------------------------------------------------------------------------------------
CpuObcSoftcore destructor
CpuObcSoftcore destructor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
CpuObcSoftcore::~CpuObcSoftcore( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::~CpuObcSoftcore";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Initialize data members
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::_initializeObcDataMembers( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::initializeDataMembers";
// ---------------------------------------------------------------------------------------
_obcSoftcoreParameters = NULL;
}
/**---------------------------------------------------------------------------------------
Get ObcSoftcoreParameters reference
Get ObcSoftcoreParameters reference
@return ObcSoftcoreParameters reference
@return ObcSoftcoreParameters reference
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters* CpuObcSoftcore::getObcSoftcoreParameters( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::getObcSoftcoreParameters";
// ---------------------------------------------------------------------------------------
return _obcSoftcoreParameters;
return _obcSoftcoreParameters;
}
/**---------------------------------------------------------------------------------------
Set ObcSoftcoreParameters reference
@param ObcSoftcoreParameters reference
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::setObcSoftcoreParameters( ObcSoftcoreParameters* obcSoftcoreParameters ){
// ---------------------------------------------------------------------------------------
Set ObcSoftcoreParameters reference
// static const char* methodName = "\nCpuObcSoftcore::setObcSoftcoreParameters";
@param ObcSoftcoreParameters reference
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
_obcSoftcoreParameters = obcSoftcoreParameters;
void CpuObcSoftcore::setObcSoftcoreParameters( ObcSoftcoreParameters* obcSoftcoreParameters ){
_obcSoftcoreParameters = obcSoftcoreParameters;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcSoftcoreParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& CpuObcSoftcore::getObcChain( void ){
// ---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcSoftcoreParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
// static const char* methodName = "\nCpuObcSoftcore::getObcChain";
@return array
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
if( _obcChain.size() == 0 ){
_obcChain.resize(_obcSoftcoreParameters->getNumberOfAtoms());
}
return _obcChain;
RealOpenMMVector& CpuObcSoftcore::getObcChain( void ){
return _obcChain;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcSoftcoreParameters->getNumberOfAtoms()
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return array
@return flag
--------------------------------------------------------------------------------------- */
const vector<RealOpenMM>& CpuObcSoftcore::getObcChainConst( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::getObcChain";
// ---------------------------------------------------------------------------------------
return _obcChain;
int CpuObcSoftcore::includeAceApproximation( void ) const {
return _includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Return OBC chain temp work array of size=_obcSoftcoreParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
Set flag indicating whether AceApproximation is to be included
@return array
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& CpuObcSoftcore::getObcChainTemp( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::getImplicitSolventObcChainTemp";
// ---------------------------------------------------------------------------------------
if( _obcChainTemp.size() == 0 ){
_obcChainTemp.resize(_obcSoftcoreParameters->getNumberOfAtoms());
}
return _obcChainTemp;
void CpuObcSoftcore::setIncludeAceApproximation( int includeAceApproximation ){
_includeAceApproximation = includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Get Born radii based on papers:
J. Phys. Chem. 1996 100, 19824-19839 (HCT paper)
Proteins: Structure, Function, and Bioinformatcis 55:383-394 (2004) (OBC paper)
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::computeBornRadii( vector<RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii ){
Calculation of Born radii based on papers:
// ---------------------------------------------------------------------------------------
J. Phys. Chem. 1996 100, 19824-19839 (HCT paper)
Proteins: Structure, Function, and Bioinformatcis 55:383-394 (2004) (OBC paper)
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM two = (RealOpenMM) 2.0;
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM fourth = (RealOpenMM) 0.25;
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
static const char* methodName = "\nCpuObcSoftcore::computeBornRadii";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
void CpuObcSoftcore::computeBornRadii( vector<RealVec>& atomCoordinates, RealOpenMMVector& bornRadii ){
ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
// ---------------------------------------------------------------------------------------
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM* scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
vector<RealOpenMM>& obcChain = getObcChain();
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
const RealOpenMM* nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
// calculate Born radii
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
//FILE* logFile = SimTKOpenMMLog::getSimTKOpenMMLogFile( );
//FILE* logFile = NULL;
//FILE* logFile = fopen( "bRSoft", "w" );
const RealOpenMMVector& atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMMVector& scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
RealOpenMMVector& obcChain = getObcChain();
const RealOpenMMVector& nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
RealOpenMM radiusIInverse = one/offsetRadiusI;
RealOpenMM sum = zero;
RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
// HCT code
// ---------------------------------------------------------------------------------------
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
// calculate Born radii
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_obcSoftcoreParameters->getUseCutoff() && r > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM rInverse = one/r;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
RealOpenMM ratio = LN( (u_ij/l_ij) );
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + ( half*rInverse*ratio) + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
// this case (atom i completely inside atom j) is not considered in the original paper
// Jay Ponder and the authors of Tinker recognized this and
// worked out the details
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( radiusIInverse - l_ij);
}
sum += nonPolarScaleFactors[atomJ]*term;
}
}
}
RealOpenMM radiusIInverse = one/offsetRadiusI;
RealOpenMM sum = zero;
// OBC-specific code (Eqs. 6-8 in paper)
sum *= nonPolarScaleFactors[atomI]*half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
// HCT code
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_obcSoftcoreParameters->getUseCutoff() && r > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM rInverse = one/r;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
bornRadii[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
RealOpenMM ratio = LN( (u_ij/l_ij) );
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + ( half*rInverse*ratio) + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
// this case (atom i completely inside atom j) is not considered in the original paper
// Jay Ponder and the authors of Tinker recognized this and
// worked out the details
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( radiusIInverse - l_ij);
}
sum += nonPolarScaleFactors[atomJ]*term;
}
}
}
// OBC-specific code (Eqs. 6-8 in paper)
obcChain[atomI] = offsetRadiusI*( alphaObc - two*betaObc*sum + three*gammaObc*sum2 );
obcChain[atomI] = (one - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
sum *= nonPolarScaleFactors[atomI]*half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
bornRadii[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
obcChain[atomI] = offsetRadiusI*( alphaObc - two*betaObc*sum + three*gammaObc*sum2 );
obcChain[atomI] = (one - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
#if 0
if( logFile && atomI >= 0 ){
(void) fprintf( logFile, "RRQ %d sum %12.6e tanhS %12.6e radI %.5f %.5f born %18.10e obc %12.6e\n",
atomI, sum, tanhSum, radiusI, offsetRadiusI, bornRadii[atomI], obcChain[atomI] );
}
#endif
}
#if 0
if( logFile ){
(void) fclose( logFile );
logFile = fopen( "bRSoftJ", "w" );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%6d %18.10e %18.10e\n",
atomI, bornRadii[atomI], obcChain[atomI] );
}
(void) fclose( logFile );
}
#endif
}
return;
}
/**---------------------------------------------------------------------------------------
Get nonpolar solvation force constribution via ACE approximation
Get nonpolar solvation force constribution via ACE approximation
@param obcSoftcoreParameters parameters
@param vdwRadii Vdw radii
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
@param obcSoftcoreParameters parameters
@param vdwRadii Vdw radii
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::computeAceNonPolarForce( const ObcSoftcoreParameters* obcSoftcoreParameters,
const vector<RealOpenMM>& bornRadii, RealOpenMM* energy,
vector<RealOpenMM>& forces ) const {
const vector<RealOpenMM>& bornRadii, RealOpenMM* energy,
vector<RealOpenMM>& forces ) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM minusSix = -6.0;
// static const char* methodName = "\nCpuImplicitSolvent::computeAceNonPolarForce";
// ---------------------------------------------------------------------------------------
static const RealOpenMM minusSix = -6.0;
// compute the nonpolar solvation via ACE approximation
// ---------------------------------------------------------------------------------------
const RealOpenMM probeRadius = obcSoftcoreParameters->getProbeRadius();
const RealOpenMM surfaceAreaFactor = obcSoftcoreParameters->getPi4Asolv();
// compute the nonpolar solvation via ACE approximation
const RealOpenMMVector& atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMMVector& nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
const RealOpenMM probeRadius = obcSoftcoreParameters->getProbeRadius();
const RealOpenMM surfaceAreaFactor = obcSoftcoreParameters->getPi4Asolv();
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM* nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
// the original ACE equation is based on Eq.2 of
// 1 + 1 + pow + 3 + 1 + 2 FLOP
// M. Schaefer, C. Bartels and M. Karplus, "Solution Conformations
// and Thermodynamics of Structured Peptides: Molecular Dynamics
// Simulation with an Implicit Solvation Model", J. Mol. Biol.,
// 284, 835-848 (1998) (ACE Method)
// the original ACE equation is based on Eq.2 of
// The original equation includes the factor (atomicRadii[atomI]/bornRadii[atomI]) to the first power,
// whereas here the ratio is raised to the sixth power: (atomicRadii[atomI]/bornRadii[atomI])**6
// M. Schaefer, C. Bartels and M. Karplus, "Solution Conformations
// and Thermodynamics of Structured Peptides: Molecular Dynamics
// Simulation with an Implicit Solvation Model", J. Mol. Biol.,
// 284, 835-848 (1998) (ACE Method)
// This modification was made by Jay Ponder and is based on observations that the change yields better correlations w/
// expected values. Jay did not think it was important enough to write up, so there is
// no paper to cite.
// The original equation includes the factor (atomicRadii[atomI]/bornRadii[atomI]) to the first power,
// whereas here the ratio is raised to the sixth power: (atomicRadii[atomI]/bornRadii[atomI])**6
// This modification was made by Jay Ponder who observed it gave better correlations w/
// observed values. He did not think it was important enough to write up, so there is
// no paper to cite.
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
if( bornRadii[atomI] > 0.0 ){
RealOpenMM r = atomicRadii[atomI] + probeRadius;
RealOpenMM ratio6 = POW( atomicRadii[atomI]/bornRadii[atomI], (RealOpenMM) 6.0 );
RealOpenMM saTerm = nonPolarScaleFactors[atomI]*surfaceAreaFactor*r*r*ratio6;
*energy += saTerm;
forces[atomI] += minusSix*saTerm/bornRadii[atomI];
}
}
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
if( bornRadii[atomI] > 0.0 ){
RealOpenMM r = atomicRadii[atomI] + probeRadius;
RealOpenMM ratio6 = POW( atomicRadii[atomI]/bornRadii[atomI], static_cast<RealOpenMM>( 6.0 ) );
RealOpenMM saTerm = nonPolarScaleFactors[atomI]*surfaceAreaFactor*r*r*ratio6;
*energy += saTerm;
forces[atomI] += minusSix*saTerm/bornRadii[atomI];
}
}
}
/**---------------------------------------------------------------------------------------
Get Obc Born energy and forces
Get Obc Born energy and forces
@param bornRadii Born radii -- optional; if NULL, then ObcSoftcoreParameters
entry is used
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
The array bornRadii is also updated and the obcEnergy
The array bornRadii is also updated and the obcEnergy
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::computeBornEnergyForces( vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
RealOpenMM CpuObcSoftcore::computeBornEnergyForces( vector<RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges,
vector<RealVec>& inputForces ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::computeBornEnergyForces";
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM 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 ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
const int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
const ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
const int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// constants
// constants
const RealOpenMM preFactor = two*obcSoftcoreParameters->getElectricConstant()*(
(one/obcSoftcoreParameters->getSoluteDielectric()) -
(one/obcSoftcoreParameters->getSolventDielectric()) );
const RealOpenMM preFactor = obcSoftcoreParameters->getPreFactor();
const RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
const RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
#if 0
{
RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM* scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
RealOpenMM* obcChain = getObcChain();
FILE* logFile = fopen( "bornParameters", "w" );
(void) fprintf( logFile, "%5d dielOff=%.4e rad::hct::q::bR::Chain::coords\n", numberOfAtoms, dielectricOffset );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%5d %10.5f %10.5f %10.5f %14.7e %14.7e %14.7e %14.7e %14.7e\n", atomI,
atomicRadii[atomI], scaledRadiusFactor[atomI], partialCharges[atomI], bornRadii[atomI], obcChain[atomI],
atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
}
(void) fclose( logFile );
}
#endif
// set energy/forces to zero
// compute Born radii
RealOpenMM obcEnergy = zero;
RealOpenMMVector bornRadii( numberOfAtoms );
computeBornRadii( atomCoordinates, bornRadii );
RealOpenMM** forces = (RealOpenMM**) malloc( sizeof( RealOpenMM* )*numberOfAtoms );
RealOpenMM* block = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms*3 );
memset( block, 0, sizeof( RealOpenMM )*numberOfAtoms*3 );
RealOpenMM* blockPtr = block;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forces[ii] = blockPtr;
blockPtr += 3;
}
RealOpenMM obcEnergy = zero;
RealOpenMMVector bornForces( numberOfAtoms );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
bornForces[ii] = zero;
}
vector<RealOpenMM>& bornForces = getBornForce();
bornForces.assign(numberOfAtoms, 0.0);
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// N*( 8 + pow) ACE
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcSoftcoreParameters, bornRadii, &obcEnergy, bornForces );
}
// N*( 8 + pow) ACE
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcSoftcoreParameters, bornRadii, &obcEnergy, bornForces );
}
const RealOpenMM* nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
const RealOpenMMVector& nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// first main loop
// first main loop
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
// charges are assumed to be scaled on input so nonPolarScaleFactor is not needed
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
//RealOpenMM nonPolarScaleFactor = atomI != atomJ ? nonPolarScaleFactors[atomI]*nonPolarScaleFactors[atomJ] : nonPolarScaleFactors[atomI];
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
// Gpol *= nonPolarScaleFactor;
// charges are assumed to be scaled on input so nonPolarScaleFactor is not needed
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
if( atomI != atomJ ){
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
if( atomI != atomJ ){
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
inputForces[atomI][0] += deltaX;
inputForces[atomI][1] += deltaY;
inputForces[atomI][2] += deltaZ;
} else {
Gpol *= half;
}
inputForces[atomJ][0] -= deltaX;
inputForces[atomJ][1] -= deltaY;
inputForces[atomJ][2] -= deltaZ;
// 3 FLOP
} else {
Gpol *= half;
}
obcEnergy += Gpol;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
obcEnergy += Gpol;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
}
}
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// second main loop
// second main loop
// initialize Born radii & ObcChain temp arrays -- contain values
// used in next iteration
RealOpenMMVector& obcChain = getObcChain();
const RealOpenMMVector& atomicRadii = obcSoftcoreParameters->getAtomicRadii();
vector<RealOpenMM>& bornRadiiTemp = getBornRadiiTemp();
bornRadiiTemp.assign(numberOfAtoms, 0.0);
const RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
const RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
const RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
const RealOpenMMVector& scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
vector<RealOpenMM>& obcChainTemp = getObcChainTemp();
obcChainTemp.assign(numberOfAtoms, 0.0);
// compute factor that depends only on the outer loop index
vector<RealOpenMM>& obcChain = getObcChain();
const RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bornForces[atomI] *= bornRadii[atomI]*bornRadii[atomI]*obcChain[atomI];
}
const RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
const RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
const RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
const RealOpenMM* scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
// compute factor that depends only on the outer loop index
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bornForces[atomI] *= bornRadii[atomI]*bornRadii[atomI]*obcChain[atomI];
}
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// radius w/ dielectric offset applied
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
// radius w/ dielectric offset applied
// used to compute Born radius for next iteration
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
RealOpenMM bornSum = zero;
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 (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
// radius w/ dielectric offset applied
// radius w/ dielectric offset applied
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
if( offsetRadiusI < rScaledRadiusJ ){
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
t3 *= nonPolarScaleFactors[atomI]*nonPolarScaleFactors[atomJ];
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
t3 *= nonPolarScaleFactors[atomI]*nonPolarScaleFactors[atomJ];
RealOpenMM de = bornForces[atomI]*t3*rInverse;
RealOpenMM de = bornForces[atomI]*t3*rInverse;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
forces[atomI][0] -= deltaX;
forces[atomI][1] -= deltaY;
forces[atomI][2] -= deltaZ;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
inputForces[atomI][0] -= deltaX;
inputForces[atomI][1] -= deltaY;
inputForces[atomI][2] -= deltaZ;
forces[atomJ][0] += deltaX;
forces[atomJ][1] += deltaY;
forces[atomJ][2] += deltaZ;
inputForces[atomJ][0] += deltaX;
inputForces[atomJ][1] += deltaY;
inputForces[atomJ][2] += deltaZ;
// Born radius term
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + (half*rInverse)*LN(u_ij/l_ij) +
(fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2-u_ij2);
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( (one/offsetRadiusI) - l_ij);
}
bornSum += term;
}
}
}
// OBC-specific code (Eqs. 6-8 in paper)
bornSum *= half*offsetRadiusI;
RealOpenMM sum2 = bornSum*bornSum;
RealOpenMM sum3 = bornSum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*bornSum - betaObc*sum2 + gammaObc*sum3 );
bornRadiiTemp[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
obcChainTemp[atomI] = offsetRadiusI*( alphaObc - two*betaObc*bornSum + three*gammaObc*sum2 );
obcChainTemp[atomI] = (one - tanhSum*tanhSum)*obcChainTemp[atomI]/radiusI;
}
// cal to Joule conversion
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += forces[atomI][0];
inputForces[atomI][1] += forces[atomI][1];
inputForces[atomI][2] += forces[atomI][2];
}
setEnergy( obcEnergy );
#if 0
{
RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM* scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
RealOpenMM* obcChain = getObcChain();
//FILE* logFile = fopen( "bornParameters", "w" );
FILE* logFile = stderr;
(void) fprintf( logFile, "%5d dielOff=%.4e rad::hct::q::bR::Chain::bF::f::coords\n", numberOfAtoms, dielectricOffset );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "%5d %10.5f %10.5f q=%10.5f b[%14.7e %14.7e %14.7e] f[%14.7e %14.7e %14.7e] x[%14.7e %14.7e %14.7e]\n", atomI,
atomicRadii[atomI], scaledRadiusFactor[atomI], partialCharges[atomI], bornRadii[atomI], obcChain[atomI],
conversion*bornForces[atomI],
conversion*forces[atomI][0], conversion*forces[atomI][1], conversion*forces[atomI][2],
atomCoordinates[atomI][0], atomCoordinates[atomI][1], atomCoordinates[atomI][2] );
}
if( logFile != stderr || logFile != stdout ){
(void) fclose( logFile );
}
}
#endif
// copy new Born radii and obcChain values into permanent array
}
}
}
}
bornRadii = bornRadiiTemp;
obcChain = obcChainTemp;
return obcEnergy;;
free( (char*) block );
free( (char*) forces );
}
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.h
View file @ 6718e0cc
......@@ -30,7 +30,7 @@
// ---------------------------------------------------------------------------------------
class CpuObcSoftcore : public CpuImplicitSolvent {
class CpuObcSoftcore {
private:
......@@ -40,13 +40,12 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
// arrays containing OBC chain derivative
std::vector<RealOpenMM> _obcChain;
std::vector<RealOpenMM> _obcChainTemp;
RealOpenMMVector _obcChain;
// initialize data members (more than
// one constructor, so centralize intialization here)
// flag to signal whether ACE approximation
// is to be included
void _initializeObcDataMembers( void );
int _includeAceApproximation;
public:
......@@ -60,7 +59,7 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
CpuObcSoftcore( ImplicitSolventParameters* obcSoftcoreParameters );
CpuObcSoftcore( ObcSoftcoreParameters* obcSoftcoreParameters );
/**---------------------------------------------------------------------------------------
......@@ -92,26 +91,45 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return array
@return flag
--------------------------------------------------------------------------------------- */
int includeAceApproximation( void ) const;
/**---------------------------------------------------------------------------------------
Set flag indicating whether AceApproximation is to be included
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void setIncludeAceApproximation( int includeAceApproximation );
/**---------------------------------------------------------------------------------------
Get energy
std::vector<RealOpenMM>& getObcChain( void );
const std::vector<RealOpenMM>& getObcChainConst( void ) const;
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM getEnergy( void ) const;
/**---------------------------------------------------------------------------------------
Return OBC chain temp work array of size=_implicitSolventParameters->getNumberOfAtoms()
Return OBC chain derivative: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getObcChainTemp( void );
RealOpenMMVector& getObcChain( void );
/**---------------------------------------------------------------------------------------
......@@ -124,7 +142,7 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii );
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMMVector& bornRadii );
/**---------------------------------------------------------------------------------------
......@@ -139,22 +157,23 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
void computeAceNonPolarForce( const ObcSoftcoreParameters* obcSoftcoreParameters,
const std::vector<RealOpenMM>& bornRadii, RealOpenMM* energy,
std::vector<RealOpenMM>& forces ) const;
const RealOpenMMVector& bornRadii, RealOpenMM* energy,
RealOpenMMVector& forces ) const;
/**---------------------------------------------------------------------------------------
Get Born energy and forces based on OBC
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@return energy
--------------------------------------------------------------------------------------- */
void computeBornEnergyForces( std::vector<RealOpenMM>& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& forces );
RealOpenMM computeBornEnergyForces( std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges, std::vector<OpenMM::RealVec>& forces );
};
// ---------------------------------------------------------------------------------------
......
plugins/freeEnergy/platforms/reference/src/gbsa/ObcSoftcoreParameters.cpp
View file @ 6718e0cc
......@@ -26,764 +26,425 @@
#include <string.h>
#include <sstream>
#include "openmm/OpenMMException.h"
#include "ObcSoftcoreParameters.h"
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
// #define UseGromacsMalloc 1
#ifdef UseGromacsMalloc
extern "C" {
#include "smalloc.h"
}
#endif
const std::string ObcSoftcoreParameters::ParameterFileName = std::string( "params.agb" );
/**---------------------------------------------------------------------------------------
ObcSoftcoreParameters:
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 ObcSoftcoreParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
ObcSoftcoreParameters* obcParameters = new ObcSoftcoreParameters( numberOfAtoms, log );
obcParameters->initializeParameters( top );
Gpu:
obcParameters = new ObcSoftcoreParameters( 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 ObcSoftcoreParameters destructor does not free arrays
obcParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuObc::obcVdwRadii )->getData() );
obcParameters->setVolume( getBrookStreamWrapperAtIndex( GpuObc::obcVolume )->getData() );
obcParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuObc::obcGpolFixed )->getData() );
obcParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuObc::obcBornRadii )->getData() );
obcParameters->setFreeArrays( false );
obcParameters->initializeParameters( top );
Issues:
Tinker's atom radii are used.
The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
only tested for generic proteins
see mapGmxAtomNameToTinkerAtomNumber()
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters constructor
@param numberOfAtoms number of atoms
@param obcType OBC type (Eq. 7 or 8 in paper)
/**---------------------------------------------------------------------------------------
ObcSoftcoreParameters constructor (Simbios)
@param numberOfAtoms number of atoms
@param obcType OBC type (Eq. 7 or 8 in paper)
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters::ObcSoftcoreParameters( int numberOfAtoms, ObcSoftcoreParameters::ObcType obcType ) : ImplicitSolventParameters( numberOfAtoms ), cutoff(false), periodic(false) {
ObcSoftcoreParameters::ObcSoftcoreParameters( int numberOfAtoms, ObcSoftcoreParameters::ObcType obcType ) :
_numberOfAtoms(numberOfAtoms),
_obcType(obcType), _dielectricOffset(0.009), _nonPolarPreFactor(2.25936),
_soluteDielectric(1.0), _solventDielectric(78.3), _probeRadius(0.14),
_electricConstant(-0.5*ONE_4PI_EPS0), _pi4Asolv( 28.3919551),
_cutoff(false), _periodic(false) {
// ---------------------------------------------------------------------------------------
_atomicRadii.resize( numberOfAtoms );
_scaledRadiusFactors.resize( numberOfAtoms );
_nonPolarScaleFactors.resize( numberOfAtoms );
// static const char* methodName = "\nObcSoftcoreParameters::ObcSoftcoreParameters";
// ---------------------------------------------------------------------------------------
_obcType = obcType;
_dielectricOffset = 0.009f;
_ownScaledRadiusFactors = 0;
_scaledRadiusFactors = NULL;
_ownNonPolarScaleFactors = 0;
_nonPolarScaleFactors = NULL;
_nonPolarPreFactor = (RealOpenMM) 2.25936;
setObcTypeParameters( obcType );
setObcTypeParameters( obcType );
}
/**---------------------------------------------------------------------------------------
ObcSoftcoreParameters destructor (Simbios)
ObcSoftcoreParameters destructor
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters::~ObcSoftcoreParameters( ){
}
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::~ObcSoftcoreParameters";
// ---------------------------------------------------------------------------------------
// in GPU runs, arrays may be 'owned' by BrookStreamWrapper -- hence they should not
// be freed here, i.e., _freeArrays should be 'false'
#ifdef UseGromacsMalloc
/*
if( _freeArrays ){
if( _vdwRadii != NULL ){
save_free( "_vdwRadii", __FILE__, __LINE__, _vdwRadii );
}
} */
#else
if( _ownScaledRadiusFactors ){
delete[] _scaledRadiusFactors;
}
if( _ownNonPolarScaleFactors ){
delete[] _nonPolarScaleFactors;
}
/**---------------------------------------------------------------------------------------
/*
if( getFreeArrays() ){
Get number of atoms
} */
@return number of atoms
#endif
--------------------------------------------------------------------------------------- */
int ObcSoftcoreParameters::getNumberOfAtoms( void ) const {
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Get OBC type
Get OBC type
@return OBC type
@return OBC type
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters::ObcType ObcSoftcoreParameters::getObcType( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::getObcType:";
// ---------------------------------------------------------------------------------------
return _obcType;
return _obcType;
}
/**---------------------------------------------------------------------------------------
Set OBC type specific parameters
Set OBC type specific parameters
@param obcType OBC type (ObcTypeI or ObcTypeII -- Eq. 7 or 8)
@param obcType OBC type (ObcTypeI or ObcTypeII -- Eq. 7 or 8)
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setObcTypeParameters( ObcSoftcoreParameters::ObcType obcType ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::setObcTypeParameters:";
// ---------------------------------------------------------------------------------------
if( obcType == ObcTypeI ){
_alphaObc = 0.8f;
_betaObc = 0.0f;
_gammaObc = 2.91f;
} else {
_alphaObc = 1.0f;
_betaObc = 0.8f;
_gammaObc = 4.85f;
}
_obcType = obcType;
if( obcType == ObcTypeI ){
_alphaObc = 0.8f;
_betaObc = 0.0f;
_gammaObc = 2.91f;
} else {
_alphaObc = 1.0f;
_betaObc = 0.8f;
_gammaObc = 4.85f;
}
_obcType = obcType;
}
/**---------------------------------------------------------------------------------------
Get dielectricOffset
Get dielectricOffset
@return _dielectricOffset
@return _dielectricOffset
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcSoftcoreParameters::getDielectricOffset( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::getDielectricOffset:";
// ---------------------------------------------------------------------------------------
return _dielectricOffset;
return _dielectricOffset;
}
/**---------------------------------------------------------------------------------------
Get alpha OBC (Eqs. 6 & 7) in Proteins paper
Get alpha OBC (Eqs. 6 & 7) in Proteins paper
@return alphaObc
@return alphaObc
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcSoftcoreParameters::getAlphaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::getAlphaObc:";
// ---------------------------------------------------------------------------------------
return _alphaObc;
return _alphaObc;
}
/**---------------------------------------------------------------------------------------
Get beta OBC (Eqs. 6 & 7) in Proteins paper
Get beta OBC (Eqs. 6 & 7) in Proteins paper
@return betaObc
@return betaObc
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcSoftcoreParameters::getBetaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::getBetaObc:";
// ---------------------------------------------------------------------------------------
return _betaObc;
return _betaObc;
}
/**---------------------------------------------------------------------------------------
Get gamma OBC (Eqs. 6 & 7) in Proteins paper
Get gamma OBC (Eqs. 6 & 7) in Proteins paper
@return gammaObc
@return gammaObc
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcSoftcoreParameters::getGammaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::getGammaObc:";
// ---------------------------------------------------------------------------------------
return _gammaObc;
return _gammaObc;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
Get solvent dielectric
@return array of atomic radii
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM* ObcSoftcoreParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
// if dielectric offset applied, then unapply
return atomicRadii;
RealOpenMM ObcSoftcoreParameters::getSolventDielectric( void ) const {
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Set solvent dielectric
@param atomicRadii array of atomic radii
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcSoftcoreParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
void ObcSoftcoreParameters::setSolventDielectric( RealOpenMM solventDielectric ){
_solventDielectric = solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
Get solute dielectric
@param atomicRadii vector of atomic radii
@return soluteDielectric
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nObcSoftcoreParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters::setAtomicRadii( atomicRadii );
RealOpenMM ObcSoftcoreParameters::getSoluteDielectric( void ) const {
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Return OBC scale factors
If not previously set, allocate space
Set solute dielectric
@return array
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcSoftcoreParameters::getScaledRadiusFactors( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::getScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _scaledRadiusFactors == NULL ){
ObcSoftcoreParameters* localThis = const_cast<ObcSoftcoreParameters* const>(this);
localThis->_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownScaledRadiusFactors = true;
memset( _scaledRadiusFactors, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _scaledRadiusFactors;
void ObcSoftcoreParameters::setSoluteDielectric( RealOpenMM soluteDielectric ){
_soluteDielectric = soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
Get electric constant
@param ownScaledRadiusFactors flag indicating whether scale factors
array should be deleted
@return electricConstant
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setOwnScaleFactors( int ownScaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownScaledRadiusFactors = ownScaledRadiusFactors;
RealOpenMM ObcSoftcoreParameters::getElectricConstant( void ) const {
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Set OBC scale factors
Get probe radius
@param scaledRadiusFactors scaledRadiusFactors
@return probeRadius
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != scaledRadiusFactors ){
delete[] _scaledRadiusFactors;
_ownScaledRadiusFactors = false;
}
_scaledRadiusFactors = scaledRadiusFactors;
RealOpenMM ObcSoftcoreParameters::getProbeRadius( void ) const {
return _probeRadius;
}
#if RealOpenMMType == 0
/**---------------------------------------------------------------------------------------
Set OBC scale factors
Set probe radius
@param scaledRadiusFactors scaledRadiusFactors
@param probeRadius probe radius
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setScaledRadiusFactors( float* scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _scaledRadiusFactors == NULL ){
_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
_ownScaledRadiusFactors = true;
}
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_scaledRadiusFactors[ii] = (RealOpenMM) scaledRadiusFactors[ii];
}
void ObcSoftcoreParameters::setProbeRadius( RealOpenMM probeRadius ){
_probeRadius = probeRadius;
}
#endif
/**---------------------------------------------------------------------------------------
Set OBC scale factors
Get pi*4*Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049*1000.0; (Still)
((RealOpenMM) M_PI)*4.0f*0.0054*1000.0; (OBC)
@param scaledRadiusFactors scaledRadiusFactors
@return pi4Asolv
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != NULL ){
delete[] _scaledRadiusFactors;
}
_ownScaledRadiusFactors = true;
_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
for( int ii = 0; ii < (int) scaledRadiusFactors.size(); ii++ ){
_scaledRadiusFactors[ii] = scaledRadiusFactors[ii];
}
RealOpenMM ObcSoftcoreParameters::getPi4Asolv( void ) const {
return _pi4Asolv;
}
/**---------------------------------------------------------------------------------------
Map Gmx atom name to Tinker atom number (Simbios)
@param atomName atom name (CA, HA, ...); upper and lower case should both work
@param log if set, then print error messages to log file
@return Tinker atom number if atom name is valid; else return -1
--------------------------------------------------------------------------------------- */
int ObcSoftcoreParameters::mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const {
Get AtomicRadii array
// ---------------------------------------------------------------------------------------
@return array of atomic radii
static int mapCreated = 0;
static int atomNameMap[26];
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
// set up atomNameMap array on first call to this method
// atomNameMap[ii] = Tinker atom number
// where ii = (the ASCII index - 65) of the first character in the
// input atom name; name may be lower case
if( !mapCreated ){
mapCreated = 1;
for( int ii = 0; ii < 26; ii++ ){
atomNameMap[ii] = -1;
}
// H
atomNameMap[7] = 1;
// C
atomNameMap[2] = 6;
// N
atomNameMap[13] = 7;
// O
atomNameMap[14] = 8;
// S
atomNameMap[18] = 16;
}
const RealOpenMMVector& ObcSoftcoreParameters::getAtomicRadii( void ) const {
return _atomicRadii;
}
// map first letter in atom name to Tinker atom number
/**---------------------------------------------------------------------------------------
int firstAsciiValue = ((int) atomName[0]) - 65;
Set AtomicRadii array
// check for lower case
@param atomicRadii vector of atomic radii
if( firstAsciiValue > 25 ){
firstAsciiValue -= 32;
}
--------------------------------------------------------------------------------------- */
// validate
void ObcSoftcoreParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii ){
if( firstAsciiValue < 0 || firstAsciiValue > 25 ){
if( log != NULL ){
(void) fprintf( log, "Atom name=<%s> unrecognized.", atomName );
}
(void) fprintf( stderr, "Atom name=<%s> unrecognized.", atomName );
return -1;
}
return atomNameMap[firstAsciiValue];
if( atomicRadii.size() == _atomicRadii.size() ){
for( unsigned int ii = 0; ii < atomicRadii.size(); ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
} else {
std::stringstream msg;
msg << "ObcSoftcoreParameters: 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());
}
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string ObcSoftcoreParameters::getStateString( const char* title ) const {
Return OBC scale factors
// ---------------------------------------------------------------------------------------
@return array
// static const char* methodName = "\nObcSoftcoreParameters::getStateString";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
const RealOpenMMVector& ObcSoftcoreParameters::getScaledRadiusFactors( void ) const {
return _scaledRadiusFactors;
}
std::stringstream message;
message << ImplicitSolventParameters::getStateString( title );
/**---------------------------------------------------------------------------------------
std::string tab = getStringTab();
Set OBC scale factors
if( getObcType() == ObcTypeI ){
message << tab << "OBC type: Type I";
} else {
message << tab << "OBC type: Type II";
}
message << tab << "Alpha: " << getAlphaObc();
message << tab << "Beta: " << getBetaObc();
message << tab << "Gamma: " << getGammaObc();
@param scaledRadiusFactors scaledRadiusFactors
return message.str();
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors ){
if( scaledRadiusFactors.size() == _scaledRadiusFactors.size() ){
for( unsigned int ii = 0; ii < scaledRadiusFactors.size(); ii++ ){
_scaledRadiusFactors[ii] = scaledRadiusFactors[ii];
}
} else {
std::stringstream msg;
msg << "ObcSoftcoreParameters: input size for scaled radius factors does not agree w/ current size: input=";
msg << scaledRadiusFactors.size();
msg << " current size=" << _scaledRadiusFactors.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
Set the force to use a cutoff.
@param distance the cutoff distance
@param distance the cutoff distance
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setUseCutoff( RealOpenMM distance ) {
cutoff = true;
cutoffDistance = distance;
_cutoff = true;
_cutoffDistance = distance;
}
/**---------------------------------------------------------------------------------------
Get whether to use a cutoff.
Get whether to use a cutoff.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
bool ObcSoftcoreParameters::getUseCutoff() {
return cutoff;
return _cutoff;
}
/**---------------------------------------------------------------------------------------
Get the cutoff distance.
Get the cutoff distance.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
RealOpenMM ObcSoftcoreParameters::getCutoffDistance() {
return cutoffDistance;
return _cutoffDistance;
}
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param boxSize the X, Y, and Z widths of the periodic box
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::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);
Get whether to use periodic boundary conditions.
assert(boxSize[0] >= 2.0*_cutoffDistance);
assert(boxSize[1] >= 2.0*_cutoffDistance);
assert(boxSize[2] >= 2.0*_cutoffDistance);
--------------------------------------------------------------------------------------- */
bool ObcSoftcoreParameters::getPeriodic() {
return periodic;
_periodic = true;
_periodicBoxSize[0] = boxSize[0];
_periodicBoxSize[1] = boxSize[1];
_periodicBoxSize[2] = boxSize[2];
}
/**---------------------------------------------------------------------------------------
Get the periodic box dimension
Get whether to use periodic boundary conditions.
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcSoftcoreParameters::getPeriodicBox() {
return periodicBoxSize;
bool ObcSoftcoreParameters::getPeriodic() {
return _periodic;
}
/**---------------------------------------------------------------------------------------
Return non-polar scale factors
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcSoftcoreParameters::getNonPolarScaleFactors( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "ObcSoftcoreParameters::getNonPolarScaleFactors";
Get the periodic box dimension
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
if( _nonPolarScaleFactors == NULL ){
ObcSoftcoreParameters* localThis = const_cast<ObcSoftcoreParameters* const>(this);
localThis->_nonPolarScaleFactors = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownNonPolarScaleFactors = true;
memset( _nonPolarScaleFactors, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _nonPolarScaleFactors;
const RealOpenMM* ObcSoftcoreParameters::getPeriodicBox() {
return _periodicBoxSize;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
@param ownNonPolarScaleFactors flag indicating whether scale factors
array should be deleted
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setOwnNonPolarScaleFactors( int ownNonPolarScaleFactors ){
Return non-polar scale factors
// ---------------------------------------------------------------------------------------
@return array
// static const char* methodName = "\nCpuObc::setOwnScaleFactors";
--------------------------------------------------------------------------------------- */
// ---------------------------------------------------------------------------------------
_ownNonPolarScaleFactors = ownNonPolarScaleFactors;
const RealOpenMMVector& ObcSoftcoreParameters::getNonPolarScaleFactors( void ) const {
return _nonPolarScaleFactors;
}
/**---------------------------------------------------------------------------------------
Set non-polar scale factors
Set non-polar scale factors
@param nonPolarScaleFactors nonPolarScaleFactors
@param nonPolarScaleFactors nonPolarScaleFactors
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setNonPolarScaleFactors( const RealOpenMMVector& nonPolarScaleFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setNonPolarScaleFactors";
// ---------------------------------------------------------------------------------------
if( _ownNonPolarScaleFactors ){
delete[] _nonPolarScaleFactors;
}
_ownNonPolarScaleFactors = true;
_nonPolarScaleFactors = new RealOpenMM[nonPolarScaleFactors.size()];
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_nonPolarScaleFactors[ii] = nonPolarScaleFactors[ii];
}
}
#if RealOpenMMType == 0
/**---------------------------------------------------------------------------------------
Set non-polar scale factors
@param nonPolarScaleFactors nonPolarScaleFactors
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setNonPolarScaleFactors( float* nonPolarScaleFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setNonPolarScaleFactors";
// ---------------------------------------------------------------------------------------
if( _nonPolarScaleFactors == NULL ){
_nonPolarScaleFactors = new RealOpenMM[getNumberOfAtoms()];
_ownNonPolarScaleFactors = true;
}
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_nonPolarScaleFactors[ii] = (RealOpenMM) nonPolarScaleFactors[ii];
}
if( nonPolarScaleFactors.size() == _nonPolarScaleFactors.size() ){
for( unsigned int ii = 0; ii < nonPolarScaleFactors.size(); ii++ ){
_nonPolarScaleFactors[ii] = nonPolarScaleFactors[ii];
}
} else {
std::stringstream msg;
msg << "ObcSoftcoreParameters: input size for non-polar scale factors does not agree w/ current size: input=";
msg << nonPolarScaleFactors.size();
msg << " current size=" << _nonPolarScaleFactors.size();
throw OpenMM::OpenMMException(msg.str());
}
}
#endif
/**---------------------------------------------------------------------------------------
Set OBC scale factors
Set OBC scale factors
@param nonPolarScaleFactors nonPolarScaleFactors
@param nonPolarScaleFactors nonPolarScaleFactors
--------------------------------------------------------------------------------------- */
--------------------------------------------------------------------------------------- */
void ObcSoftcoreParameters::setNonPolarPrefactor( RealOpenMM nonPolarPreFactor ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setNonPolarScaleFactors";
// ---------------------------------------------------------------------------------------
_nonPolarPreFactor = nonPolarPreFactor;
_nonPolarPreFactor = nonPolarPreFactor;
}
plugins/freeEnergy/platforms/reference/src/gbsa/ObcSoftcoreParameters.h
View file @ 6718e0cc
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
* Contributors: Pande Group
*
......@@ -26,11 +25,10 @@
#define __ObcSoftcoreParameters_H__
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "gbsa/ImplicitSolventParameters.h"
// ---------------------------------------------------------------------------------------
class ObcSoftcoreParameters : public ImplicitSolventParameters {
class ObcSoftcoreParameters {
public:
......@@ -38,37 +36,42 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
enum ObcType { ObcTypeI, ObcTypeII };
static const std::string ParameterFileName;
private:
int _numberOfAtoms;
// OBC constants & parameters
RealOpenMM _dielectricOffset;
RealOpenMM _alphaObc;
RealOpenMM _betaObc;
RealOpenMM _gammaObc;
RealOpenMM _probeRadius;
RealOpenMM _pi4Asolv;
ObcType _obcType;
RealOpenMM _nonPolarPreFactor;
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _electricConstant;
// scaling factors for nonpolar term
int _ownNonPolarScaleFactors;
RealOpenMM* _nonPolarScaleFactors;
RealOpenMMVector _atomicRadii;
RealOpenMMVector _nonPolarScaleFactors;
// scaled radius factors (S_kk in HCT paper)
int _ownScaledRadiusFactors;
RealOpenMM* _scaledRadiusFactors;
RealOpenMMVector _scaledRadiusFactors;
// cutoff and periodic boundary conditions
bool cutoff;
bool periodic;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
bool _cutoff;
bool _periodic;
RealOpenMM _periodicBoxSize[3];
RealOpenMM _cutoffDistance;
public:
......@@ -90,6 +93,97 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
~ObcSoftcoreParameters( );
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int getNumberOfAtoms( void ) const;
/**---------------------------------------------------------------------------------------
Get electric constant
@return electric constant
--------------------------------------------------------------------------------------- */
RealOpenMM getElectricConstant( void ) const;
/**---------------------------------------------------------------------------------------
Get probe radius (Simbios)
@return probeRadius
--------------------------------------------------------------------------------------- */
RealOpenMM getProbeRadius( void ) const;
/**---------------------------------------------------------------------------------------
Set probe radius (Simbios)
@param probeRadius probe radius
--------------------------------------------------------------------------------------- */
void setProbeRadius( RealOpenMM probeRadius );
/**---------------------------------------------------------------------------------------
Get pi4Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049f*1000.0f; (Simbios)
@return pi4Asolv
--------------------------------------------------------------------------------------- */
RealOpenMM getPi4Asolv( void ) const;
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSolventDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void setSolventDielectric( RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSoluteDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void setSoluteDielectric( RealOpenMM soluteDielectric );
/**---------------------------------------------------------------------------------------
Get OBC type
......@@ -158,33 +252,18 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
const RealOpenMM* getScaledRadiusFactors( void ) const;
const RealOpenMMVector& getScaledRadiusFactors( void ) const;
/**---------------------------------------------------------------------------------------
Return OBC scale factors
Set OBC scale factors
@return array
@param input vector of radius factors
--------------------------------------------------------------------------------------- */
void setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors );
#if RealOpenMMType == 0
void setScaledRadiusFactors( float* scaledRadiusFactors );
#endif
void setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors );
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors arra should be deleted
@param ownScaledRadiusFactors flag indicating whether scale factors
array should be deleted
--------------------------------------------------------------------------------------- */
void setOwnScaleFactors( int ownScaledRadiusFactors );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
......@@ -193,17 +272,7 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
RealOpenMM* getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
--------------------------------------------------------------------------------------- */
void setAtomicRadii( RealOpenMM* atomicRadii );
const RealOpenMMVector& getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
......@@ -215,31 +284,6 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
void setAtomicRadii( const RealOpenMMVector& atomicRadii );
/**---------------------------------------------------------------------------------------
Map Gmx atom name to Tinker atom number (Simbios)
@param atomName atom name (CA, HA, ...); upper and lower case should both work
@param log if set, then print error messages to log file
return Tinker atom number if atom name is valid; else return -1
--------------------------------------------------------------------------------------- */
int mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
......@@ -256,7 +300,7 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
bool getUseCutoff();
bool getUseCutoff( void );
/**---------------------------------------------------------------------------------------
......@@ -264,7 +308,7 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
RealOpenMM getCutoffDistance();
RealOpenMM getCutoffDistance( void );
/**---------------------------------------------------------------------------------------
......@@ -284,7 +328,7 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
bool getPeriodic();
bool getPeriodic( void );
/**---------------------------------------------------------------------------------------
......@@ -292,29 +336,17 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
--------------------------------------------------------------------------------------- */
const RealOpenMM* getPeriodicBox();
const RealOpenMM* getPeriodicBox( void );
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
@param ownNonPolarScaleFactors flag indicating whether scale factors
array should be deleted
--------------------------------------------------------------------------------------- */
void setOwnNonPolarScaleFactors( int ownNonPolarScaleFactors );
/**---------------------------------------------------------------------------------------
Return non-polar scale factors
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* getNonPolarScaleFactors( void ) const;
const RealOpenMMVector& getNonPolarScaleFactors( void ) const;
/**---------------------------------------------------------------------------------------
......@@ -338,19 +370,6 @@ class ObcSoftcoreParameters : public ImplicitSolventParameters {
};
/**---------------------------------------------------------------------------------------
Qsort/heapsort integer comparison (Simbios)
@parma a first value to compare
@param b second value to compare
@return -1, 0, 1
--------------------------------------------------------------------------------------- */
int integerComparison( const void *a, const void *b);
// ---------------------------------------------------------------------------------------
#endif // __ObcSoftcoreParameters_H__
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