/* Portions copyright (c) 2006 Stanford University and Simbios. * Contributors: Pande Group * * Permission is hereby granted, free of charge, to any person obtaining * a copy of this software and associated documentation files (the * "Software"), to deal in the Software without restriction, including * without limitation the rights to use, copy, modify, merge, publish, * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject * to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include "../SimTKUtilities/SimTKOpenMMCommon.h" #include "../SimTKUtilities/SimTKOpenMMLog.h" #include "../SimTKUtilities/SimTKOpenMMUtilities.h" #include "CpuImplicitSolvent.h" #include using namespace OpenMM; using namespace std; //#define UseGromacsMalloc 1 // Replacement new/delete w/ Gromac's smalloc() and sfree() // static data member created by call to cpuSetImplicitSolventParameters() // stores parameter settings, ... // used to calculate GBSA forces/energy CpuImplicitSolvent* CpuImplicitSolvent::_cpuImplicitSolvent = NULL; /**--------------------------------------------------------------------------------------- CpuImplicitSolvent constructor @param implicitSolventParameters implicitSolventParameters object --------------------------------------------------------------------------------------- */ CpuImplicitSolvent::CpuImplicitSolvent( ImplicitSolventParameters* implicitSolventParameters ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::CpuImplicitSolvent(2)"; // --------------------------------------------------------------------------------------- _initializeDataMembers( ); _implicitSolventParameters = implicitSolventParameters; } /**--------------------------------------------------------------------------------------- CpuImplicitSolvent destructor --------------------------------------------------------------------------------------- */ CpuImplicitSolvent::~CpuImplicitSolvent( ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::~CpuImplicitSolvent"; // --------------------------------------------------------------------------------------- delete _implicitSolventParameters; } /**--------------------------------------------------------------------------------------- Initialize data members --------------------------------------------------------------------------------------- */ void CpuImplicitSolvent::_initializeDataMembers( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::_initializeDataMembers"; // --------------------------------------------------------------------------------------- _includeAceApproximation = 0; _forceConversionFactor = (RealOpenMM) 1.0; _energyConversionFactor = (RealOpenMM) 1.0; _forceCallIndex = 0; _implicitSolventEnergy = (RealOpenMM) 0.0; } /**--------------------------------------------------------------------------------------- Return number of atoms @return number of atoms --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::getNumberOfAtoms( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getNumberOfAtoms"; // --------------------------------------------------------------------------------------- return _implicitSolventParameters->getNumberOfAtoms(); } /**--------------------------------------------------------------------------------------- Return energy @return energy --------------------------------------------------------------------------------------- */ RealOpenMM CpuImplicitSolvent::getEnergy( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getEnergy"; // --------------------------------------------------------------------------------------- return _implicitSolventEnergy; } /**--------------------------------------------------------------------------------------- Delete static _cpuImplicitSolvent object if set @return SimTKOpenMMCommon::DefaultReturn if _cpuImplicitSolvent was set; otherwise return SimTKOpenMMCommon::ErrorReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::deleteCpuImplicitSolvent( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::deleteCpuImplicitSolvent"; // --------------------------------------------------------------------------------------- if( _cpuImplicitSolvent != NULL ){ delete _cpuImplicitSolvent; _cpuImplicitSolvent = NULL; return SimTKOpenMMCommon::DefaultReturn; } else { return SimTKOpenMMCommon::ErrorReturn; } } /**--------------------------------------------------------------------------------------- Set static member _cpuImplicitSolvent @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::setCpuImplicitSolvent( CpuImplicitSolvent* cpuImplicitSolvent ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::setCpuImplicitSolvent"; // --------------------------------------------------------------------------------------- _cpuImplicitSolvent = cpuImplicitSolvent; return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Get static member cpuImplicitSolvent @return static member cpuImplicitSolvent --------------------------------------------------------------------------------------- */ CpuImplicitSolvent* CpuImplicitSolvent::getCpuImplicitSolvent( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getCpuImplicitSolvent"; // --------------------------------------------------------------------------------------- return _cpuImplicitSolvent; } /**--------------------------------------------------------------------------------------- Set energy @param energy @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::setEnergy( RealOpenMM energy ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::setEnergy"; // --------------------------------------------------------------------------------------- _implicitSolventEnergy = energy; return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Return ImplicitSolventParameters @return ImplicitSolventParameters --------------------------------------------------------------------------------------- */ ImplicitSolventParameters* CpuImplicitSolvent::getImplicitSolventParameters( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getImplicitSolventParameters"; // --------------------------------------------------------------------------------------- return _implicitSolventParameters; } /**--------------------------------------------------------------------------------------- Set ImplicitSolventParameters @param ImplicitSolventParameters @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::setImplicitSolventParameters( ImplicitSolventParameters* implicitSolventParameters ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::setImplicitSolventParameters"; // --------------------------------------------------------------------------------------- _implicitSolventParameters = implicitSolventParameters; return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Return flag signalling whether AceApproximation for nonpolar term is to be included @return flag --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::includeAceApproximation( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::includeAceApproximation"; // --------------------------------------------------------------------------------------- return _includeAceApproximation; } /**--------------------------------------------------------------------------------------- Set flag indicating whether AceApproximation is to be included @param includeAceApproximation new includeAceApproximation value @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::setIncludeAceApproximation( int includeAceApproximation ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::setImplicitSolventParameters"; // --------------------------------------------------------------------------------------- _includeAceApproximation = includeAceApproximation; return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Return ForceConversionFactor for units @return ForceConversionFactor --------------------------------------------------------------------------------------- */ RealOpenMM CpuImplicitSolvent::getForceConversionFactor( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getForceConversionFactor"; // --------------------------------------------------------------------------------------- return _forceConversionFactor; } /**--------------------------------------------------------------------------------------- Set ForceConversionFactor @param ForceConversionFactor (units conversion) @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::setForceConversionFactor( RealOpenMM forceConversionFactor ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::setForceConversionFactor"; // --------------------------------------------------------------------------------------- _forceConversionFactor = forceConversionFactor; return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Return EnergyConversionFactor for units @return EnergyConversionFactor --------------------------------------------------------------------------------------- */ RealOpenMM CpuImplicitSolvent::getEnergyConversionFactor( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getEnergyConversionFactor"; // --------------------------------------------------------------------------------------- return _energyConversionFactor; } /**--------------------------------------------------------------------------------------- Set EnergyConversionFactor @param EnergyConversionFactor (units conversion) @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::setEnergyConversionFactor( RealOpenMM energyConversionFactor ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::setEnergyConversionFactor"; // --------------------------------------------------------------------------------------- _energyConversionFactor = energyConversionFactor; return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Return ForceCallIndex -- number of times forces have been calculated @return ForceCallIndex --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::getForceCallIndex( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getForceCallIndex"; // --------------------------------------------------------------------------------------- return _forceCallIndex; } /**--------------------------------------------------------------------------------------- Increment ForceCallIndex @return incremented forceCallIndex --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::incrementForceCallIndex( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::incrementForceCallIndex"; // --------------------------------------------------------------------------------------- _forceCallIndex++; return _forceCallIndex; } /**--------------------------------------------------------------------------------------- Return bornForce, a work array of size _implicitSolventParameters->getNumberOfAtoms()*sizeof( RealOpenMM ) On first call, memory for array is allocated if not set @return array --------------------------------------------------------------------------------------- */ vector& CpuImplicitSolvent::getBornForce( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getImplicitSolventBornForce"; // --------------------------------------------------------------------------------------- if( _bornForce.size() == 0 ){ _bornForce.resize(_implicitSolventParameters->getNumberOfAtoms()); } return _bornForce; } /**--------------------------------------------------------------------------------------- Return Born radii: size = _implicitSolventParameters->getNumberOfAtoms() On first call, memory for array is allocated if it is not set @return array --------------------------------------------------------------------------------------- */ vector& CpuImplicitSolvent::getBornRadii( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getBornRadii"; // --------------------------------------------------------------------------------------- if( _bornRadii.size() == 0 ){ _bornRadii.resize(_implicitSolventParameters->getNumberOfAtoms(), 0.0); } return _bornRadii; } /**--------------------------------------------------------------------------------------- Return Born radii: size = _implicitSolventParameters->getNumberOfAtoms() @return array --------------------------------------------------------------------------------------- */ const vector& CpuImplicitSolvent::getBornRadiiConst( void ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getBornRadii"; // --------------------------------------------------------------------------------------- return _bornRadii; } /**--------------------------------------------------------------------------------------- Return Born radii temp work array of size=_implicitSolventParameters->getNumberOfAtoms() On first call, memory for array is allocated if not set @return array --------------------------------------------------------------------------------------- */ vector& CpuImplicitSolvent::getBornRadiiTemp( void ){ // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getImplicitSolventBornRadiiTemp"; // --------------------------------------------------------------------------------------- if( _tempBornRadii.size() == 0 ){ _tempBornRadii.resize(_implicitSolventParameters->getNumberOfAtoms(), 0.0); } return _tempBornRadii; } /**--------------------------------------------------------------------------------------- Compute Born radii @param atomCoordinates atomic coordinates @param bornRadii output array of Born radii @param obcChain output array of Obc chain derivatives @return SimTKOpenMMCommon::DefaultReturn or SimTKOpenMMCommon::ErrorReturn if problems encountered --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::computeBornRadii( vector& atomCoordinates, vector& bornRadii ){ // --------------------------------------------------------------------------------------- static const char* methodName = "\nCpuImplicitSolvent::computeBornRadii"; // --------------------------------------------------------------------------------------- std::stringstream message; message << methodName; message << " Error: calling from base class."; SimTKOpenMMLog::printError( message ); return SimTKOpenMMCommon::ErrorReturn; } /**--------------------------------------------------------------------------------------- Get Born energy and forces @param atomCoordinates atomic coordinates @param partialCharges partial charges @param forces forces (output) @param updateBornRadii if set, then Born radii are updated for current configuration; otherwise radii correspond to configuration from previous iteration @return SimTKOpenMMCommon::DefaultReturn; abort if cpuImplicitSolvent is not set --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::computeImplicitSolventForces( vector& atomCoordinates, const RealOpenMM* partialCharges, vector& forces, int updateBornRadii ){ // --------------------------------------------------------------------------------------- static const char* methodName = "\nCpuImplicitSolvent::computeImplicitSolventForces"; // --------------------------------------------------------------------------------------- int callId = incrementForceCallIndex(); // check parameters have been initialized ImplicitSolventParameters* implicitSolventParameters = getImplicitSolventParameters(); if( !implicitSolventParameters ){ std::stringstream message; message << methodName; message << " implicitSolventParameters has not been initialized!"; SimTKOpenMMLog::printError( message ); return SimTKOpenMMCommon::ErrorReturn; } // check to see if Born radii have been previously calculated // if not, then calculate; // logic here assumes that the radii are intitialized to zero // and then once computed, always greater than zero. // after first iteration Born radii are updated in force calculation (computeBornEnergyForces()) // unless updateBornRadii is set vector& bornRadii = getBornRadii(); if( updateBornRadii || bornRadii[0] < (RealOpenMM) 0.0001 || callId == 1 ){ computeBornRadii( atomCoordinates, bornRadii ); } // compute forces computeBornEnergyForces( getBornRadii(), atomCoordinates, partialCharges, forces ); return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Get Born energy and forces based on J. Phys. Chem. A V101 No 16, p. 3005 (Simbios) @param bornRadii Born radii -- optional; if NULL, then ImplicitSolventParameters entry is used @param atomCoordinates atomic coordinates @param partialCharges partial charges @param forces forces @return SimTKOpenMMCommon::ErrorReturn since the call should be implemented in a derived class --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::computeBornEnergyForces( vector& bornRadii, vector& atomCoordinates, const RealOpenMM* partialCharges, vector& forces ){ // --------------------------------------------------------------------------------------- static const char* methodName = "\nCpuImplicitSolvent::computeBornEnergyForces"; // --------------------------------------------------------------------------------------- std::stringstream message; message << methodName; message << " Error: calling from base class."; SimTKOpenMMLog::printError( message ); return SimTKOpenMMCommon::ErrorReturn; } /**--------------------------------------------------------------------------------------- Get nonpolar solvation force constribution via ACE approximation @param implicitSolventParameters parameters @param vdwRadii Vdw radii @param bornRadii Born radii @param energy energy (output): value is incremented from input value @param forces forces: values are incremented from input values @return SimTKOpenMMCommon::DefaultReturn --------------------------------------------------------------------------------------- */ int CpuImplicitSolvent::computeAceNonPolarForce( const ImplicitSolventParameters* implicitSolventParameters, const vector& bornRadii, RealOpenMM* energy, vector& forces ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::computeAceNonPolarForce"; static const RealOpenMM minusSix = -6.0; // --------------------------------------------------------------------------------------- // compute the nonpolar solvation via ACE approximation const RealOpenMM probeRadius = implicitSolventParameters->getProbeRadius(); const RealOpenMM surfaceAreaFactor = implicitSolventParameters->getPi4Asolv(); const RealOpenMM* atomicRadii = implicitSolventParameters->getAtomicRadii(); int numberOfAtoms = implicitSolventParameters->getNumberOfAtoms(); // 1 + 1 + pow + 3 + 1 + 2 FLOP // the original ACE equation is based on Eq.2 of // M. Schaefer, C. Bartels and M. Karplus, "Solution Conformations // and Thermodynamics of Structured Peptides: Molecular Dynamics // Simulation with an Implicit Solvation Model", J. Mol. Biol., // 284, 835-848 (1998) (ACE Method) // The original equation includes the factor (atomicRadii[atomI]/bornRadii[atomI]) to the first power, // whereas here the ratio is raised to the sixth power: (atomicRadii[atomI]/bornRadii[atomI])**6 // This modification was made by Jay Ponder who observed it gave better correlations w/ // observed values. He did not think it was important enough to write up, so there is // no paper to cite. for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ if( bornRadii[atomI] > 0.0 ){ RealOpenMM r = atomicRadii[atomI] + probeRadius; RealOpenMM ratio6 = POW( atomicRadii[atomI]/bornRadii[atomI], (RealOpenMM) 6.0 ); RealOpenMM saTerm = surfaceAreaFactor*r*r*ratio6; *energy += saTerm; forces[atomI] += minusSix*saTerm/bornRadii[atomI]; } } return SimTKOpenMMCommon::DefaultReturn; } /**--------------------------------------------------------------------------------------- Get string w/ state @param title title (optional) @return string containing state --------------------------------------------------------------------------------------- */ std::string CpuImplicitSolvent::getStateString( const char* title ) const { // --------------------------------------------------------------------------------------- // static const char* methodName = "\nCpuImplicitSolvent::getStateString"; // --------------------------------------------------------------------------------------- std::stringstream message; if( title ){ message << title; } message << "\nImplicitSolvent info:"; message << "\nForce conversion=" << getForceConversionFactor() << " Energy conversion=" << getEnergyConversionFactor(); message << "\nInclude ACE approximation="; if( includeAceApproximation() ){ message << "Y"; } else { message << "N"; } // get parameters message << getImplicitSolventParameters()->getStateString( NULL ); return message.str(); }