Skip to content

GitLab

Commit f3e77e33 authored by Peter Eastman's avatar Peter Eastman
Browse files

Moved GBSA reference code into the reference platform directory

parent 0aad1354
Hide whitespace changes
Inline Side-by-side
Showing with 8753 additions and 0 deletions
platforms/reference/src/Gbsa/Gbsa.vcproj 0 → 100644
View file @ f3e77e33
<?xml version="1.0" encoding="UTF-8"?>
<VisualStudioProject
ProjectType="Visual C++"
Version="9.00"
Name="Gbsa"
ProjectGUID="{519D78EF-DC26-47E3-B240-994A447774E5}"
RootNamespace="Gbsa"
Keyword="Win32Proj"
TargetFrameworkVersion="131072"
>
<Platforms>
<Platform
Name="Win32"
/>
</Platforms>
<ToolFiles>
</ToolFiles>
<Configurations>
<Configuration
Name="Debug|Win32"
OutputDirectory="Debug"
IntermediateDirectory="Debug"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
Optimization="0"
AdditionalIncludeDirectories="..\gromacs\include;..\simtk;..\SimTKUtilities"
PreprocessorDefinitions="WIN32;_DEBUG;_LIB"
MinimalRebuild="true"
BasicRuntimeChecks="3"
RuntimeLibrary="3"
UsePrecompiledHeader="0"
WarningLevel="3"
DebugInformationFormat="4"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
OutputFile="$(OutDir)\$(ProjectName).lib"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
<Configuration
Name="Release|Win32"
OutputDirectory="Release"
IntermediateDirectory="Release"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
AdditionalIncludeDirectories="..\gromacs\include;..\simtk;..\SimTKUtilities"
PreprocessorDefinitions="WIN32;NDEBUG;_LIB;"
RuntimeLibrary="2"
UsePrecompiledHeader="0"
WarningLevel="3"
DebugInformationFormat="3"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
OutputFile="$(OutDir)\Gbsa.lib"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
<Configuration
Name="Debug_FAHCORE|Win32"
OutputDirectory="$(ConfigurationName)"
IntermediateDirectory="$(ConfigurationName)"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
Optimization="0"
AdditionalIncludeDirectories="..\gromacs\include;..\simtk;..\SimTKUtilities"
PreprocessorDefinitions="WIN32;_DEBUG;_LIB;"
MinimalRebuild="true"
BasicRuntimeChecks="3"
RuntimeLibrary="3"
UsePrecompiledHeader="0"
WarningLevel="3"
DebugInformationFormat="4"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
OutputFile="$(OutDir)/Gbsa.lib"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
<Configuration
Name="Release_FAHCORE|Win32"
OutputDirectory="$(ConfigurationName)"
IntermediateDirectory="$(ConfigurationName)"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
Optimization="0"
AdditionalIncludeDirectories="..\gromacs\include;..\simtk;../../corewrap;..\SimTKUtilities"
PreprocessorDefinitions="WIN32;_LIB;FULLINDIRECT"
MinimalRebuild="true"
BasicRuntimeChecks="3"
RuntimeLibrary="2"
UsePrecompiledHeader="0"
WarningLevel="3"
DebugInformationFormat="4"
ForcedIncludeFiles="swindirect.h"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
OutputFile="$(OutDir)/Gbsa.lib"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
<Configuration
Name="DebugNoGPU_FAH|Win32"
OutputDirectory="$(ConfigurationName)"
IntermediateDirectory="$(ConfigurationName)"
ConfigurationType="4"
>
<Tool
Name="VCPreBuildEventTool"
/>
<Tool
Name="VCCustomBuildTool"
/>
<Tool
Name="VCXMLDataGeneratorTool"
/>
<Tool
Name="VCWebServiceProxyGeneratorTool"
/>
<Tool
Name="VCMIDLTool"
/>
<Tool
Name="VCCLCompilerTool"
Optimization="0"
AdditionalIncludeDirectories="..\gromacs\include;..\simtk;..\SimTKUtilities"
PreprocessorDefinitions="WIN32;_DEBUG;_LIB;NoGPU"
MinimalRebuild="true"
BasicRuntimeChecks="3"
RuntimeLibrary="3"
UsePrecompiledHeader="0"
WarningLevel="3"
DebugInformationFormat="4"
/>
<Tool
Name="VCManagedResourceCompilerTool"
/>
<Tool
Name="VCResourceCompilerTool"
/>
<Tool
Name="VCPreLinkEventTool"
/>
<Tool
Name="VCLibrarianTool"
OutputFile="$(OutDir)\$(ProjectName).lib"
/>
<Tool
Name="VCALinkTool"
/>
<Tool
Name="VCXDCMakeTool"
/>
<Tool
Name="VCBscMakeTool"
/>
<Tool
Name="VCFxCopTool"
/>
<Tool
Name="VCPostBuildEventTool"
/>
</Configuration>
</Configurations>
<References>
</References>
<Files>
<Filter
Name="Header Files"
Filter="h;hpp;hxx;hm;inl;inc;xsd"
UniqueIdentifier="{93995380-89BD-4b04-88EB-625FBE52EBFB}"
>
<File
RelativePath=".\CpuGbsa.h"
>
</File>
<File
RelativePath=".\CpuGrycuk.h"
>
</File>
<File
RelativePath=".\cpuGrycukInterface.h"
>
</File>
<File
RelativePath=".\CpuImplicitSolvent.h"
>
</File>
<File
RelativePath=".\CpuObc.h"
>
</File>
<File
RelativePath=".\cpuObcInterface.h"
>
</File>
<File
RelativePath=".\gmx_modifications.h"
>
</File>
<File
RelativePath=".\gmxGbsa.h"
>
</File>
<File
RelativePath=".\gmxObc.h"
>
</File>
<File
RelativePath=".\gromacsCpuObcInterface.h"
>
</File>
<File
RelativePath=".\GrycukParameters.h"
>
</File>
<File
RelativePath=".\ImplicitSolventAtom.h"
>
</File>
<File
RelativePath=".\ImplicitSolventParameters.h"
>
</File>
<File
RelativePath=".\ObcParameters.h"
>
</File>
</Filter>
<Filter
Name="Resource Files"
Filter="rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx"
UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}"
>
</Filter>
<Filter
Name="Source Files"
Filter="cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx"
UniqueIdentifier="{4FC737F1-C7A5-4376-A066-2A32D752A2FF}"
>
<File
RelativePath=".\CpuGrycuk.cpp"
>
</File>
<File
RelativePath=".\cpuGrycukInterface.cpp"
>
</File>
<File
RelativePath=".\CpuImplicitSolvent.cpp"
>
</File>
<File
RelativePath=".\CpuObc.cpp"
>
</File>
<File
RelativePath=".\cpuObcInterface.cpp"
>
</File>
<File
RelativePath=".\gromacsCpuGrycukInterface.cpp"
>
</File>
<File
RelativePath=".\gromacsCpuObcInterface.cpp"
>
</File>
<File
RelativePath=".\GrycukParameters.cpp"
>
</File>
<File
RelativePath=".\ImplicitSolventParameters.cpp"
>
</File>
<File
RelativePath=".\ObcParameters.cpp"
>
</File>
</Filter>
</Files>
<Globals>
</Globals>
</VisualStudioProject>
platforms/reference/src/Gbsa/GbsaAtomParameter.cpp 0 → 100644
View file @ f3e77e33
/* 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 <math.h>
// #include "UtilitiesSimTk.h"
#include "GbsaAtomParameter.h"
/**---------------------------------------------------------------------------------------
GbsaAtomParameter:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in Gbsa equation gPol
(4) list of atoms that should be excluded in calculating
force -- nonbonded atoms (1-2, and 1-3 atoms)
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
GbsaParameters constructor (Simbios)
@param parameterLineTokens tokens from parameter file
--------------------------------------------------------------------------------------- */
GbsaAtomParameter::GbsaAtomParameter( const StringVector& parameterLineTokens ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::GbsaAtomParameter";
// ---------------------------------------------------------------------------------------
if( parameterLineTokens.size() < 2 ){
// XXX
return;
}
setTypeId( parameterLineTokens[0] );
setVdwRadius( parameterLineTokens[1] );
}
/**---------------------------------------------------------------------------------------
GbsaAtomParameter destructor (Simbios)
--------------------------------------------------------------------------------------- */
GbsaAtomParameter::~GbsaAtomParameter( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::~GbsaAtomParameter";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Get type id (Simbios)
@return type id
--------------------------------------------------------------------------------------- */
std::string GbsaAtomParameter::getTypeId( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::getTypeId:";
// ---------------------------------------------------------------------------------------
return _typeId;
}
/**---------------------------------------------------------------------------------------
Set type id (Simbios)
@param typeId type id
@return 0
--------------------------------------------------------------------------------------- */
int GbsaAtomParameter::setTypeId( std::string typeId ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::setTypeId:";
// ---------------------------------------------------------------------------------------
_typeId = typeId;
return 0;
}
/**---------------------------------------------------------------------------------------
Get vdw radius (Simbios)
@return vdw radius
--------------------------------------------------------------------------------------- */
Real GbsaAtomParameter::getVdwRadius( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::getVdwRadius:";
// ---------------------------------------------------------------------------------------
return _vdwRadius;
}
/**---------------------------------------------------------------------------------------
Set vdw radius (Simbios)
@param vdwRadius new vdw radius
@return 0
--------------------------------------------------------------------------------------- */
int GbsaAtomParameter::setVdwRadius( Real vdwRadius ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::setVdwRadius:";
// ---------------------------------------------------------------------------------------
_vdwRadius = vdwRadius;
return 0;
}
/**---------------------------------------------------------------------------------------
Set vdw radius (Simbios)
@param vdwRadius new vdw radius (string)
@return 0
--------------------------------------------------------------------------------------- */
int GbsaAtomParameter::setVdwRadius( const std::string& vdwRadius ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaAtomParameter::setVdwRadius(s):";
// ---------------------------------------------------------------------------------------
// bool SimTKOpenMMUtilities::isValidReal( std::string stringToCheck );
return setVdwRadius( (Real) atof( vdwRadius.c_str() ) );
}
platforms/reference/src/Gbsa/GbsaAtomParameter.h 0 → 100755
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __GbsaAtomParameter_H__
#define __GbsaAtomParameter_H__
// STL includes
#include <vector>
#include <map>
#include <set>
#include <string>
#include "CommonSimTk.h"
#include "RealTypeSimTk.h"
//#ifndef Real
//#define Real float
//#endif
/*
class StringComparisonForMap : public std::binary_function< const std::string&, const std::string&, bool > {
public:
bool operator()( const std::string& str1, const std::string& str2 ) const {
return strcmp( str1.c_str(), str2.c_str() ) < 0;
}
};
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
*/
// ---------------------------------------------------------------------------------------
class GbsaAtomParameter {
private:
// type id -- int Gromacs 'amber99_10', for example
std::string _typeId;
// vdw radius
Real _vdwRadius;
public:
/**---------------------------------------------------------------------------------------
GbsaParameters constructor (Simbios)
@param parameterLineTokens tokens from parameter file
--------------------------------------------------------------------------------------- */
GbsaAtomParameter( const StringVector& parameterLineTokens );
/**---------------------------------------------------------------------------------------
GbsaAtomParameter destructor (Simbios)
--------------------------------------------------------------------------------------- */
~GbsaAtomParameter( );
/**---------------------------------------------------------------------------------------
Get type id (Simbios)
@return type id
--------------------------------------------------------------------------------------- */
std::string getTypeId( void ) const;
/**---------------------------------------------------------------------------------------
Set type id (Simbios)
@param typeId type id
@return 0
--------------------------------------------------------------------------------------- */
int setTypeId( std::string typeId );
/**---------------------------------------------------------------------------------------
Get vdw radius (Simbios)
@return vdw radius
--------------------------------------------------------------------------------------- */
Real getVdwRadius( void ) const;
/**---------------------------------------------------------------------------------------
Set vdw radius (Simbios)
@param vdwRadius new vdw radius
@return 0
--------------------------------------------------------------------------------------- */
int setVdwRadius( Real vdwRadius );
/**---------------------------------------------------------------------------------------
Set vdw radius (Simbios)
@param vdwRadius new vdw radius
@return 0
--------------------------------------------------------------------------------------- */
int setVdwRadius( const std::string& vdwRadius );
};
typedef std::map<std::string, GbsaAtomParameter*, StringComparisonForMap> GbsaAtomParameterMap;
typedef GbsaAtomParameterMap::iterator GbsaAtomParameterMapI;
typedef GbsaAtomParameterMap::const_iterator GbsaAtomParameterMapCI;
// ---------------------------------------------------------------------------------------
#endif // __GbsaAtomParameter_H__
platforms/reference/src/Gbsa/GbsaParameters.cpp 0 → 100644
View file @ f3e77e33
/* 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.
*/
// generic c includes
#include <math.h>
#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMGromacsUtilities.h"
#include "GbsaAtomParameter.h"
#include "GbsaParameters.h"
#define UseGromacsMalloc 1
#ifdef UseGromacsMalloc
extern "C" {
#include "smalloc.h"
}
#endif
/**---------------------------------------------------------------------------------------
GbsaParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in Gbsa 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 GbsaParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
GbsaParameters* gbsaParameters = new GbsaParameters( numberOfAtoms, log );
gbsaParameters->initializeParameters( top );
Gpu:
gbsaParameters = new GbsaParameters( 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 GbsaParameters destructor does not free arrays
gbsaParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuGbsa::gbsaVdwRadii )->getData() );
gbsaParameters->setVolume( getBrookStreamWrapperAtIndex( GpuGbsa::gbsaVolume )->getData() );
gbsaParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuGbsa::gbsaGpolFixed )->getData() );
gbsaParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuGbsa::gbsaBornRadii )->getData() );
gbsaParameters->setFreeArrays( false );
gbsaParameters->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()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
GbsaParameters constructor (Simbios)
@param numberOfAtoms number of atoms
@param log file descriptor (can be null)
--------------------------------------------------------------------------------------- */
GbsaParameters::GbsaParameters( int numberOfAtoms, FILE* log ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::GbsaParameters";
// ---------------------------------------------------------------------------------------
_numberOfAtoms = numberOfAtoms;
_log = log;
_gbsaBondsArray = NULL;
_exclusionWorkArray = NULL;
_vdwRadii = NULL;
_volume = NULL;
_gPolFixed = NULL;
_bornRadii = NULL;
// see comments in ~GbsaParameters for explanation
_freeArrays = true;
}
/**---------------------------------------------------------------------------------------
GbsaParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
GbsaParameters::~GbsaParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::~GbsaParameters";
// ---------------------------------------------------------------------------------------
// 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 );
}
if( _volume != NULL ){
save_free( "_volume", __FILE__, __LINE__, _volume );
}
if( _gPolFixed != NULL ){
save_free( "_gPolFixed", __FILE__, __LINE__, _gPolFixed );
}
if( _bornRadii != NULL ){
save_free( "_bornRadii", __FILE__, __LINE__, _bornRadii );
}
}
if( _exclusionWorkArray != NULL ){
save_free( "_exclusionWorkArray", __FILE__, __LINE__, _exclusionWorkArray );
}
#else
if( _freeArrays ){
if( _vdwRadii != NULL ){
delete[] _vdwRadii;
}
if( _volume != NULL ){
delete[] _volume;
}
if( _gPolFixed != NULL ){
delete[] _gPolFixed;
}
if( _bornRadii != NULL ){
delete[] _bornRadii;
}
}
if( _exclusionWorkArray != NULL ){
delete[] _exclusionWorkArray;
}
#endif
if( _gbsaBondsArray != NULL ){
freeBondArray( _numberOfAtoms, _gbsaBondsArray );
}
}
/**---------------------------------------------------------------------------------------
Initialize Gbsa Parameters (Simbios)
--------------------------------------------------------------------------------------- */
void GbsaParameters::initializeConstants( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::initializeConstants:";
// ---------------------------------------------------------------------------------------
// GBSA parameters from '97 paper
_innerDielectric = 1.0f;
_solventDielectric = 78.3f;
_kcalA_To_kJNm = 4.184f*0.1f;
_phi = -0.09f;
_P1 = 0.073f;
_P2 = 0.921f;
_P3 = 6.211f;
_P4 = 15.236f;
_P4_2 = 0.5f*sqrt( _P4 );
_P5 = 1.254f;
_probeRadius = 0.14f;
_electricConstant = -166.02691f;
// ---------------------------------------------------------------------------------------
_resetPreFactor();
_P5Inverse = 1.0f/_P5;
_piP5 = ((Real) M_PI)*_P5;
_pi4Asolv = ((Real) M_PI)*4.0f*0.0049f*1000.0f;
}
/**---------------------------------------------------------------------------------------
Reset prefactor (Simbios)
called when _electricConstant, _innerDielectric, or _solventDielectric are modified
--------------------------------------------------------------------------------------- */
void GbsaParameters::_resetPreFactor( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::_resetPreFactor:";
// ---------------------------------------------------------------------------------------
if( _innerDielectric != 0.0f && _solventDielectric != 0.0f ){
_preFactor = 2.0f*_electricConstant*( (1.0f/_innerDielectric) - (1.0f/_solventDielectric) );
} else {
_preFactor = 0.0;
}
}
/**---------------------------------------------------------------------------------------
Set solvent dielectric (Simbios)
@param solventDielectric solvent dielectric
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::setSolventDielectric( Real solventDielectric ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::setSolventDielectric:";
// ---------------------------------------------------------------------------------------
_solventDielectric = solventDielectric;
_resetPreFactor();
return 0;
}
/**---------------------------------------------------------------------------------------
Set inner dielectric (Simbios)
@param innerDielectric: inner dielectric
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::setInnerDielectric( Real innerDielectric ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::setInnerDielectric:";
// ---------------------------------------------------------------------------------------
_innerDielectric = innerDielectric;
_resetPreFactor();
return 0;
}
/**---------------------------------------------------------------------------------------
Set electric constant (Simbios)
@param electricConstant electric constant
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::setElectricConstant( Real electricConstant ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::setElectricConstant:";
// ---------------------------------------------------------------------------------------
_electricConstant = electricConstant;
_resetPreFactor();
return 0;
}
/**---------------------------------------------------------------------------------------
Initialize Gbsa Parameters (Simbios)
@param top GMX topology data struct -- used to get harmonic angle
@return 0 if no errors
--------------------------------------------------------------------------------------- */
int GbsaParameters::initializeParameters( const t_topology* top ){
// ---------------------------------------------------------------------------------------
bool printParameters = false;
bool debugOn = false;
// bool debugOn = false;
static const char* methodName = "\nGbsaParameters::initializeParameters:";
// ---------------------------------------------------------------------------------------
FILE* log = getLog();
initializeConstants();
// find stretch bond entries via F_BONDS & F_SHAKE (for H's) entries
// Settle 'bonds' handled separtely
gbsaBonds** gbsaBondsArray = (gbsaBonds**) allocateBondsArray( getNumberOfAtoms() );
int errors = addStretchBonds( getNumberOfAtoms(), gbsaBondsArray, top, F_BONDS, 3, 1, log );
errors += addStretchBonds( getNumberOfAtoms(), gbsaBondsArray, top, F_SHAKE, 3, 1, log );
errors += addSettleStretchBonds( getNumberOfAtoms(), gbsaBondsArray, top, F_SETTLE, 2, 1, log );
if( errors > 0 ){
exit(-1);
} else if( log != NULL ){
(void) fprintf( log, "\nNo 1-2 bond errors." );
(void) fflush( log );
}
// ---------------------------------------------------------------------------------------
// find angle bond entries based on stretch bond entries
errors = addAngleBonds( getNumberOfAtoms(), gbsaBondsArray, top, F_ANGLES, log );
setGbsaBondsArray( gbsaBondsArray );
if( errors > 0 ){
printBondsArray( getNumberOfAtoms(), top, gbsaBondsArray, "Bonds", log );
exit(-1);
} else if( log != NULL ){
(void) fprintf( log, "\nNo 1-3 bond errors." );
(void) fflush( log );
}
// set exclusion array
setGbsaExclusions( top, debugOn, log );
if( log != NULL ){
(void) fprintf( log, "%s exclusions set", methodName );
(void) fflush( log );
}
// print bonds and report if errors
if( debugOn && log ){
printBondsArray( getNumberOfAtoms(), top, gbsaBondsArray, "Bonds", log );
}
// ---------------------------------------------------------------------------------------
// allocate memory for work arrays, if not already allocated
// when Gbsa is run on GPU, the memory for these arrays may
// have been allocated when the BrookStreamWrappers
// were created (GpuGbsa constructor) and hence do not
// need to be allocated here
unsigned int arraySz = getNumberOfAtoms()*sizeof( Real );
#ifdef UseGromacsMalloc
if( _vdwRadii == NULL ){
_vdwRadii = (Real*) save_malloc( "_vdwRadii", __FILE__, __LINE__, arraySz );
}
if( _volume == NULL ){
_volume = (Real*) save_malloc( "_volume", __FILE__, __LINE__, arraySz );
}
if( _gPolFixed == NULL ){
_gPolFixed = (Real*) save_malloc( "_gPolFixed", __FILE__, __LINE__, arraySz );
}
if( _bornRadii == NULL ){
_bornRadii = (Real*) save_malloc( "_bornRadii", __FILE__, __LINE__, arraySz );
}
#else
if( _vdwRadii == NULL ){
_vdwRadii = new Real[getNumberOfAtoms()];
}
if( _volume == NULL ){
_volume = new Real[getNumberOfAtoms()];
}
if( _gPolFixed == NULL ){
_gPolFixed = new Real[getNumberOfAtoms()];
}
if( _bornRadii == NULL ){
_bornRadii = new Real[getNumberOfAtoms()];
}
#endif
memset( _vdwRadii, 0, arraySz );
memset( _volume, 0, arraySz );
memset( _gPolFixed, 0, arraySz );
memset( _bornRadii, 0, arraySz );
// compute atomic volumes (eqs 6 & 7) in J. Phys. Chem. A V101 No 16, p. 3005
// In Tinker: ksolv.f ~ line 340 (solvtyp.eq.'STILL')
// get Macromodel vdW radii and volumes (excluding overlap)
// and calculate fixed GBSA terms
std::string parameterFileName = "Params_macromodel_agb.txt";
fprintf( log, "\nCalling getMacroModelAtomicRadii %s", parameterFileName.c_str() );
getMacroModelAtomicRadii( getNumberOfAtoms(), parameterFileName, top, top->atoms.atomname, _vdwRadii, log );
/*
getMacroModelAtomicRadii( getNumberOfAtoms(), gbsaBondsArray, top->atoms.atomname, _bornRadii, log );
char*** names = top->atoms.atomname;
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
Real diff = abs( _bornRadii[ii] - _vdwRadii[ii] );
fprintf( log, "\n%d %s %.3f %3f %.3f", ii, (*(names[ii])), diff, _vdwRadii[ii], _bornRadii[ii] );
}
fprintf( log, "\nDone getMacroModelAtomicRadii %s", parameterFileName.c_str() );
fflush( log );
exit(0);
*/
getMacroModelAtomicVolumes( getNumberOfAtoms(), top, gbsaBondsArray, _vdwRadii, _volume, log );
getFixedGBSA_GPol( getNumberOfAtoms(), top, gbsaBondsArray, _vdwRadii, _volume, _gPolFixed, log );
// ---------------------------------------------------------------------------------------
// print parameters
if( printParameters && log != NULL ){
printParameterInfo( getNumberOfAtoms(), top,
"GmxParameters.txt", log );
}
// ---------------------------------------------------------------------------------------
if( log != NULL ){
(void) fprintf( log, "%s done", methodName );
(void) fflush( log );
}
return 0;
}
/**---------------------------------------------------------------------------------------
Set Gbsa Exclusions (Simbios)
@param top GMX topology data struct -- used to get harmonic angle
@param printOn print diagnostics
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::setGbsaExclusions( const t_topology* top, bool printOn, FILE* log ){
// ---------------------------------------------------------------------------------------
static const int collectionBinMax = 500;
int collectionBin[collectionBinMax];
bool saveExclusionIndices = true;
// static const char* methodName = "\nGbsaParameters::setGbsaExclusions";
// ---------------------------------------------------------------------------------------
int numberOfAtoms = getNumberOfAtoms();
gbsaBonds** gbsaBondsArray = getGbsaBondsArray();
// excludedAtoms: used to mark which atoms are to excluded
#ifdef UseGromacsMalloc
int* excludedAtoms = (int*) save_malloc( "excludedAtoms", __FILE__, __LINE__, (numberOfAtoms + 2)*sizeof( int ) );
#else
int* excludedAtoms = new int[numberOfAtoms+2];
#endif
for( int ii = 0; ii < numberOfAtoms + 2; ii++ ){
excludedAtoms[ii] = -1;
}
// gather indices to be excluded for atom ii and then
// set values in excludeIndices array
// excludedAtoms is used to insure each atom is only
// excluded once
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( gbsaBondsArray[ii] != NULL ){
int numberOfIndicesToExclude =
collectGbsaBondIndices( gbsaBondsArray[ii], collectionBinMax, collectionBin, ii, log );
gbsaBondsArray[ii]->minExclusionIndex = numberOfAtoms + 1;
gbsaBondsArray[ii]->maxExclusionIndex = -1;
// SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( ii, top, MAX_BUFFER_SZ, buffer, numberOfAtoms, ATOM_ID_STRING_TAB );
// (void) fprintf( log, "\n[%d %s] CalcExcnt=%d ", ii, buffer, numberOfIndicesToExclude );
for( int jj = 0; jj < numberOfIndicesToExclude; jj++ ){
int excludeIndex = collectionBin[jj];
if( excludedAtoms[excludeIndex] != ii ){
excludedAtoms[excludeIndex] = ii;
if( saveExclusionIndices && gbsaBondsArray[ii]->numberOfExclusions < MAX_BOND_EXCLUSIONS ){
gbsaBondsArray[ii]->exclusions[gbsaBondsArray[ii]->numberOfExclusions++] = excludeIndex;
// (void) fprintf( log, "[%d %d] ", excludeIndex, gbsaBondsArray[ii]->numberOfExclusions );
}
if( gbsaBondsArray[ii]->minExclusionIndex > excludeIndex ){
gbsaBondsArray[ii]->minExclusionIndex = excludeIndex;
}
if( gbsaBondsArray[ii]->maxExclusionIndex < excludeIndex ){
gbsaBondsArray[ii]->maxExclusionIndex = excludeIndex;
}
}
}
}
}
// free excludedAtoms array
#ifdef UseGromacsMalloc
save_free( "excludedAtoms", __FILE__, __LINE__, excludedAtoms );
#else
delete[] excludedAtoms;
#endif
// ---------------------------------------------------------------------------------------
return 0;
}
/**---------------------------------------------------------------------------------------
Collect Gbsa bond indices for a particular atom (Simbios)
@param gbsaBond gbsaBond
@param collectionBin array of indices to be excluded
@param excludeIndex exclude index
@return number of collected indices (includes excludeIndex at end of array)
--------------------------------------------------------------------------------------- */
int GbsaParameters::collectGbsaBondIndices( gbsaBonds* gbsaBond, int collectionBinMax,
int* collectionBin, int excludeIndex, FILE* log ){
// ---------------------------------------------------------------------------------------
static const int gbsaStretchBondArrayMax = 500;
gbsaStretchBond* gbsaStretchBondArray[gbsaStretchBondArrayMax];
static const char* methodName = "\nGbsaParameters:;collectGbsaBondIndices";
// ---------------------------------------------------------------------------------------
// first collect all stretch bonds from the stretchBonds and angleBonds data structs
int numberOfStretchBonds = 0;
gbsaStretchBond* stretchBonds = gbsaBond->stretchBonds;
while( stretchBonds != NULL && numberOfStretchBonds < gbsaStretchBondArrayMax ){
gbsaStretchBondArray[numberOfStretchBonds++] = stretchBonds;
stretchBonds = stretchBonds->nextBond;
}
/*
(void) fprintf( log, "\n excludeIndex=%d str=%d", excludeIndex, numberOfStretchBonds );
(void) fflush( log );
*/
gbsaAngleBond* angleBonds = gbsaBond->angleBonds;
while( angleBonds != NULL && numberOfStretchBonds < (gbsaStretchBondArrayMax-1) ){
if( angleBonds->stretchBondI != NULL ){
gbsaStretchBondArray[numberOfStretchBonds++] = angleBonds->stretchBondI;
}
if( angleBonds->stretchBondJ != NULL ){
gbsaStretchBondArray[numberOfStretchBonds++] = angleBonds->stretchBondJ;
}
angleBonds = angleBonds->nextBond;
}
if( numberOfStretchBonds >= gbsaStretchBondArrayMax ){
(void) fprintf( log ? log : stderr, "%s gbsaStretchBondArrayMax=%d too small", methodName, gbsaStretchBondArrayMax );
return -1;
}
/*
(void) fprintf( log, " str2=%d", numberOfStretchBonds );
(void) fflush( log );
*/
// ---------------------------------------------------------------------------------------
// now collect all indices including the exclude index
int collectionBinSize = 0;
for( int ii = 0; ii < numberOfStretchBonds && collectionBinSize < (collectionBinMax+2); ii++ ){
/*
(void) fprintf( log, " %s", (gbsaStretchBondArray[ii] == NULL) ? "NULL" : "Ok" );
(void) fflush( log );
*/
if( gbsaStretchBondArray[ii]->atomI != excludeIndex ){
collectionBin[collectionBinSize++] = gbsaStretchBondArray[ii]->atomI;
}
if( gbsaStretchBondArray[ii]->atomJ != excludeIndex ){
collectionBin[collectionBinSize++] = gbsaStretchBondArray[ii]->atomJ;
}
/*
(void) fprintf( log, " collect=[%d %d %d]", collectionBinSize, gbsaStretchBondArray[ii]->atomI, gbsaStretchBondArray[ii]->atomJ );
(void) fflush( log );
*/
}
if( collectionBinSize < collectionBinMax ){
collectionBin[collectionBinSize++] = excludeIndex;
} else {
(void) fprintf( log ? log : stderr, "%s collectionBinMax=%d is too small.", methodName, collectionBinMax );
}
return collectionBinSize;
}
/**---------------------------------------------------------------------------------------
Allocate memory for bond array (Simbios)
@param maxAtoms max number of atoms
array entries are initialized to zero
free array[0] and then array when done
@returns ptr to allocated array or NULL if out of memory
--------------------------------------------------------------------------------------- */
gbsaBonds** GbsaParameters::allocateBondsArray( int maxAtoms ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGbsaParameters::allocateBondsArray";
// ---------------------------------------------------------------------------------------
#ifdef UseGromacsMalloc
gbsaBonds** gbsaBondsArray = (gbsaBonds**) save_malloc( "gbsaBondsArray", __FILE__, __LINE__, maxAtoms*sizeof( gbsaBonds* ) );
gbsaBonds* gbsaBondsBlock = (gbsaBonds* ) save_malloc( "gbsaBondsBlock", __FILE__, __LINE__, maxAtoms*sizeof( gbsaBonds ) );
#else
gbsaBonds** gbsaBondsArray = new gbsaBonds*[maxAtoms];
gbsaBonds* gbsaBondsBlock = new gbsaBonds[maxAtoms];
#endif
if( gbsaBondsArray == NULL || gbsaBondsBlock == NULL ){
(void) fprintf( stderr, "%s apparently out of memory: maxAtoms=%d", methodName, maxAtoms );
(void) fflush( stderr );
return NULL;
}
memset( gbsaBondsBlock, 0, sizeof( gbsaBonds )*maxAtoms );
for( int ii = 0; ii < maxAtoms; ii++ ){
gbsaBondsArray[ii] = gbsaBondsBlock++;
}
return gbsaBondsArray;
}
/**---------------------------------------------------------------------------------------
Deallocate memory for bond array (Simbios)
@param maxAtoms max number of atoms
@param gbsaBondsArray array to be freed
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::freeBondArray( int maxAtoms, gbsaBonds** gbsaBondsArray ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::freeBondArray";
// ---------------------------------------------------------------------------------------
for( int ii = 0; ii < maxAtoms; ii++ ){
if( gbsaBondsArray[ii] != NULL ){
gbsaStretchBond* nextStretch;
gbsaStretchBond* stretchBonds = gbsaBondsArray[ii]->stretchBonds;
while( stretchBonds != NULL ){
nextStretch = stretchBonds->nextBond;
#ifdef UseGromacsMalloc
save_free( "stretchBonds", __FILE__, __LINE__, stretchBonds );
#else
delete stretchBonds;
#endif
stretchBonds = nextStretch;
}
gbsaAngleBond* angleBonds = gbsaBondsArray[ii]->angleBonds;
gbsaAngleBond* nextAngle;
while( angleBonds != NULL ){
nextAngle = angleBonds->nextBond;
#ifdef UseGromacsMalloc
save_free( "angleBonds", __FILE__, __LINE__, angleBonds );
#else
delete angleBonds;
#endif
angleBonds = nextAngle;
}
}
}
#ifdef UseGromacsMalloc
save_free( "gbsaBondsArray[0]", __FILE__, __LINE__, gbsaBondsArray[0] );
save_free( "gbsaBondsArray", __FILE__, __LINE__, gbsaBondsArray );
#else
delete[] gbsaBondsArray[0];
delete[] gbsaBondsArray;
#endif
return 0;
}
/**---------------------------------------------------------------------------------------
Print bond array (Simbios)
@param numberOfAtoms number of atoms
@param top GMX topology data struct
@param gbsaBonds bond array
@param title title string (optional)
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::printBondsArray( int numberOfAtoms, const t_topology* top,
gbsaBonds** gbsaBondsArray, const char* title, FILE* log ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::printBondsArray";
// ---------------------------------------------------------------------------------------
if( !log ){
return 0;
}
// print title (if set) and bonds
if( title != NULL ){
(void) fprintf( log, "\n%s", title );
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( gbsaBondsArray[ii] != NULL ){
printBond( numberOfAtoms, top, ii, gbsaBondsArray[ii], log );
}
}
(void) fflush( log );
return 0;
}
/**---------------------------------------------------------------------------------------
Print bond (Simbios)
@param top GMX topology data struct -- used to get harmonic angle
@param atomIndex atom index
@param gbsaBond bond data struct
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::printBond( int numberOfAtoms, const t_topology* top, int atomIndex,
const gbsaBonds* gbsaBond, FILE* log ) const {
// ---------------------------------------------------------------------------------------
int stretchCount;
int angleCount;
char buffer[4][MAX_BUFFER_SZ];
// static const char* methodName = "\nGbsaParameters::printBond";
// ---------------------------------------------------------------------------------------
// print bond
if( log == NULL ){
return 0;
}
if( gbsaBond != NULL ){
getBondCounts( gbsaBond, &stretchCount, &angleCount );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( atomIndex, top, MAX_BUFFER_SZ, buffer[0], numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "\n%d %s StretchCount=%d AngleCount=%d ExclRange=[%d %d] ExclWidth=%d Save=[%d %d] width=%d",
atomIndex, buffer[0], stretchCount, angleCount,
gbsaBond->minExclusionIndex, gbsaBond->maxExclusionIndex,
gbsaBond->maxExclusionIndex - gbsaBond->minExclusionIndex,
gbsaBond->startExclusionIndex,
gbsaBond->stopExclusionIndex,
gbsaBond->stopExclusionIndex - gbsaBond->startExclusionIndex
);
if( gbsaBond->numberOfExclusions ){
(void) fprintf( log, "\nExclusions: %d [ ", gbsaBond->numberOfExclusions );
for( int ii = 0; ii < gbsaBond->numberOfExclusions; ii++ ){
(void) fprintf( log, "%d ", gbsaBond->exclusions[ii] );
}
(void) fprintf( log, "]", gbsaBond->numberOfExclusions );
}
gbsaStretchBond* stretchBonds = gbsaBond->stretchBonds;
int count = 0;
while( stretchBonds != NULL ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBonds->atomJ, top, MAX_BUFFER_SZ, buffer[0],
numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "\n %d %s %.6f", count++, buffer[0], stretchBonds->bondLength );
stretchBonds = stretchBonds->nextBond;
}
gbsaAngleBond* angleBonds = gbsaBond->angleBonds;
count = 0;
while( angleBonds != NULL ){
gbsaStretchBond* stretchBondI = angleBonds->stretchBondI;
gbsaStretchBond* stretchBondJ = angleBonds->stretchBondJ;
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBondI->atomI, top, MAX_BUFFER_SZ, buffer[0],
numberOfAtoms, ATOM_ID_STRING_TAB );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBondI->atomJ, top, MAX_BUFFER_SZ, buffer[1],
numberOfAtoms, ATOM_ID_STRING_TAB );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBondJ->atomI, top, MAX_BUFFER_SZ, buffer[2],
numberOfAtoms, ATOM_ID_STRING_TAB );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBondJ->atomJ, top, MAX_BUFFER_SZ, buffer[3],
numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "\n %d [ %s %s %s %s] [%.4f %.4f] pvt=%d vol=%d Ang=%.2f Ord=[%d %d %d]",
count++, buffer[0], buffer[1], buffer[2], buffer[3],
stretchBondI->bondLength,
stretchBondJ->bondLength,
angleBonds->pivotAtomIndex,
angleBonds->volumeIndex,
angleBonds->harmonicAngleWidth,
angleBonds->orderedIndices[0],
angleBonds->orderedIndices[1],
angleBonds->orderedIndices[2]
);
angleBonds = angleBonds->nextBond;
}
}
(void) fflush( log );
return 0;
}
/**---------------------------------------------------------------------------------------
Print parameter info (Simbios)
@param numberOfAtoms number of atoms
@param top GMX topology data struct -- used to get harmonic angle
@param parameterInfoFileName parameterInfo FileName
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::printParameterInfo( int numberOfAtoms, const t_topology* top,
const char* parameterInfoFileName, FILE* log ){
// ---------------------------------------------------------------------------------------
char buffer[4][MAX_BUFFER_SZ];
static const char* methodName = "\nGbsaParameters::printParameterInfo";
// ---------------------------------------------------------------------------------------
gbsaBonds** gbsaBondsArray = getGbsaBondsArray();
FILE* parameterInfoFile = NULL;
#ifdef WIN32
fopen_s( &parameterInfoFile, parameterInfoFileName, "w" );
#else
parameterInfoFile = fopen( parameterInfoFileName, "w" );
#endif
if( parameterInfoFile != NULL ){
if( log != NULL ){
(void) fprintf( log, "%s opened file=<%s>.", methodName, parameterInfoFileName );
(void) fflush( log );
}
} else {
if( log != NULL ){
(void) fprintf( log, "%s could not open file=<%s> -- abort output.",
methodName, parameterInfoFile );
(void) fflush( log );
}
return -1;
}
(void) fprintf( parameterInfoFile, "# %d atoms\n", numberOfAtoms );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( ii, top, MAX_BUFFER_SZ, buffer[0], numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( parameterInfoFile, "%d %s\n", ii, buffer[0] );
}
(void) fflush( parameterInfoFile );
(void) fprintf( parameterInfoFile, "StretchBond\n" );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( gbsaBondsArray[ii] != NULL ){
gbsaStretchBond* stretchBonds = gbsaBondsArray[ii]->stretchBonds;
while( stretchBonds != NULL ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBonds->atomI, top, MAX_BUFFER_SZ, buffer[0],
numberOfAtoms, ATOM_ID_STRING_TAB );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( stretchBonds->atomJ, top, MAX_BUFFER_SZ, buffer[1],
numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( parameterInfoFile, "%d %d %.6e %s %s\n",
stretchBonds->atomI, stretchBonds->atomJ, stretchBonds->bondLength,
buffer[0], buffer[1] );
stretchBonds = stretchBonds->nextBond;
}
}
}
(void) fflush( parameterInfoFile );
(void) fprintf( parameterInfoFile, "AngleBond\n" );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( gbsaBondsArray[ii] != NULL ){
gbsaAngleBond* angleBonds = gbsaBondsArray[ii]->angleBonds;
while( angleBonds != NULL ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( angleBonds->orderedIndices[0], top, MAX_BUFFER_SZ, buffer[0],
numberOfAtoms, ATOM_ID_STRING_TAB );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( angleBonds->orderedIndices[1], top, MAX_BUFFER_SZ, buffer[1],
numberOfAtoms, ATOM_ID_STRING_TAB );
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( angleBonds->orderedIndices[2], top, MAX_BUFFER_SZ, buffer[2],
numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( parameterInfoFile, "%d %d %d %.6e %s %s %s\n",
angleBonds->orderedIndices[0],
angleBonds->orderedIndices[1],
angleBonds->orderedIndices[2],
angleBonds->harmonicAngleWidth, buffer[0], buffer[1], buffer[2] );
angleBonds = angleBonds->nextBond;
}
}
}
(void) fclose( parameterInfoFile );
return 0;
}
/**---------------------------------------------------------------------------------------
Add angle bonds using stretch bond lists (Simbios)
Also load in angle width
@param maxAtoms max number of atoms
@param gbsaBondsArray array of gbsaBonds
@param top GMX topology data struct -- used to get harmonic angle
@param idefArrayIndex parameter index for GMX iparams data struct
@param log if set, then print error messages to log file
@return 0 if no errors
--------------------------------------------------------------------------------------- */
int GbsaParameters::addAngleBonds( int maxAtoms, gbsaBonds** gbsaBondsArray,
const t_topology* top, int idefArrayIndex, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGbsaParameters::addAngleBonds";
// ---------------------------------------------------------------------------------------
// loop over atoms
// for each atomI loop over stretch bonds
// for each atomJ in stretch bond, check if atomI/atomJ != ii
// create new angle bond w/ atomJ as the pivot atom
for( int ii = 0; ii < maxAtoms; ii++ ){
gbsaStretchBond* stretchBonds = gbsaBondsArray[ii]->stretchBonds;
while( stretchBonds != NULL ){
int atomJ = stretchBonds->atomJ;
gbsaStretchBond* stretchBondsAtomJ = gbsaBondsArray[atomJ]->stretchBonds;
while( stretchBondsAtomJ != NULL ){
if( stretchBondsAtomJ->atomI != ii && stretchBondsAtomJ->atomJ != ii ){
addStretchBondsToAngleBondList( stretchBonds, stretchBondsAtomJ,
atomJ, gbsaBondsArray, ii, log );
if( stretchBondsAtomJ->atomI != atomJ ){
addStretchBondsToAngleBondList( stretchBonds, stretchBondsAtomJ,
atomJ, gbsaBondsArray, stretchBondsAtomJ->atomI, log );
}
if( stretchBondsAtomJ->atomJ != atomJ ){
addStretchBondsToAngleBondList( stretchBonds, stretchBondsAtomJ,
atomJ, gbsaBondsArray, stretchBondsAtomJ->atomJ, log );
}
}
stretchBondsAtomJ = stretchBondsAtomJ->nextBond;
}
stretchBonds = stretchBonds->nextBond;
}
}
// printBondsArray( maxAtoms, top, gbsaBondsArray, "PreAngleWidth", log );
// add angle 'width'
// match atom indices w/ GMX indices and then use GMX harmonic.rA value
// to set angle width
t_iatom* atoms = top->idef.il[idefArrayIndex].iatoms;
t_iparams* params = top->idef.iparams;
int offset = 4;
// loop over GMX angles
// set angle width in each of the 3 atoms
// print warning message if a bond was not found
for( int ii = 0; ii < top->idef.il[idefArrayIndex].nr; ii += offset ){
int type = atoms[ii];
Real angleWidth = params[ type ].harmonic.rA*SimTKOpenMMCommon::DegreeToRadians ;
for( int jj = 0; jj < 3; jj++ ){
if( jj != 1 ){
gbsaAngleBond* angleBond = findAngleBond( gbsaBondsArray[atoms[ii + jj + 1]],
atoms[ii + 1],
atoms[ii + 2],
atoms[ii + 3] );
if( angleBond != NULL && angleBond->harmonicAngleWidth <= 0.0f ){
angleBond->harmonicAngleWidth = angleWidth;
} else if( log != NULL ){
char buffer[3][MAX_BUFFER_SZ];
for( int kk = 0; kk < 3; kk++ ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( atoms[ii + kk + 1], top, MAX_BUFFER_SZ, buffer[kk],
maxAtoms, ATOM_ID_STRING_TAB );
}
(void) fprintf( log, "%s Warning No angle bond for [%s %s %s] [%d %d %d] in atm=%d",
methodName, buffer[0], buffer[1], buffer[2], atoms[ii + 1],
atoms[ii + 2], atoms[ii + 3], atoms[ii + jj + 1] );
}
}
}
}
return 0;
}
/**---------------------------------------------------------------------------------------
Add stretch bonds to angle data struct (Simbios)
@param stretchBondI stretch bond I
@param stretchBondJ stretch bond J
@param pivotAtomIndex index of shared atom
@param gbsaBonds gbsaBond
@param saveIndex index to save????
@param log if set, then print error messages to log file
@return allocated gbsaAngleBond
--------------------------------------------------------------------------------------- */
gbsaAngleBond* GbsaParameters::addStretchBondsToAngleBondList( gbsaStretchBond* stretchBondI,
gbsaStretchBond* stretchBondJ,
int pivotAtomIndex, gbsaBonds** gbsaBonds,
int saveIndex, FILE* log ){
// ---------------------------------------------------------------------------------------
int orderedIndices[3];
static const char* methodName = "\nGbsaParameters::addStretchBondsToAngleBondList";
// ---------------------------------------------------------------------------------------
int duplicateCount =
sortAtomAngleIndices( stretchBondI->atomI, stretchBondI->atomJ,
stretchBondJ->atomI, stretchBondJ->atomJ,
orderedIndices );
if( duplicateCount != 1 && log ){
(void) fprintf( log, "%s duplicate count from sortAtomAngleIndices is not 1 (=%d) ind=[%d %d %d %d]",
methodName, duplicateCount, stretchBondI->atomI, stretchBondI->atomJ,
stretchBondJ->atomI, stretchBondJ->atomJ );
}
gbsaAngleBond* newAngleBond;
if( findAngleBondGivenOrderedArray( gbsaBonds[saveIndex], orderedIndices ) == NULL ){
#ifdef UseGromacsMalloc
newAngleBond = (gbsaAngleBond*) save_malloc( "gbsaAngleBond", __FILE__, __LINE__, sizeof( gbsaAngleBond ) );
#else
newAngleBond = new gbsaAngleBond;
#endif
memset( newAngleBond, 0, sizeof( gbsaAngleBond ) );
// gbsaAngleBond* newAngleBond;
// snew( newAngleBond, 1 );
newAngleBond->stretchBondI = stretchBondI;
newAngleBond->stretchBondJ = stretchBondJ;
newAngleBond->pivotAtomIndex = pivotAtomIndex;
memcpy( newAngleBond->orderedIndices, orderedIndices, sizeof(int)*3 );
int volumeIndex = -1;
for( int ii = 0; ii < 3 && volumeIndex < 0; ii++ ){
if( orderedIndices[ii] != saveIndex && orderedIndices[ii] != pivotAtomIndex ){
volumeIndex = orderedIndices[ii];
}
}
newAngleBond->volumeIndex = volumeIndex;
// push current entry onto list
newAngleBond->nextBond = gbsaBonds[saveIndex]->angleBonds;
gbsaBonds[saveIndex]->angleBonds = newAngleBond;
gbsaBonds[saveIndex]->angleBondCount++;
} else {
newAngleBond = NULL;
}
return newAngleBond;
}
/**---------------------------------------------------------------------------------------
Find angle bond with atom indices that match input angle indices (Simbios)
@param gbsaBond gbsaBond
@param atomI index of atomI
@param atomJ index of atomJ
@param atomK index of atomK
@return allocated gbsaAngleBond
--------------------------------------------------------------------------------------- */
gbsaAngleBond* GbsaParameters::findAngleBond( const gbsaBonds* gbsaBond, int atomI,
int atomJ, int atomK ) const {
// ---------------------------------------------------------------------------------------
int orderedAtomIndices[3];
// ---------------------------------------------------------------------------------------
// sort atom indices and them compare to all available bonds
// return bond when all indices match or NULL if no match
sortAtomAngleIndices( atomI, atomI, atomJ, atomK, orderedAtomIndices );
return findAngleBondGivenOrderedArray( gbsaBond, orderedAtomIndices );
}
/**---------------------------------------------------------------------------------------
Find angle bond with atom indices that match input ordered angle indices (Simbios)
@param gbsaBond gbsaBond
@param orderedAtomIndices array of ordered indices
@return gbsaAngleBond if found
--------------------------------------------------------------------------------------- */
gbsaAngleBond* GbsaParameters::findAngleBondGivenOrderedArray( const gbsaBonds* gbsaBond,
int* orderedAtomIndices ) const {
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// sort atom indices and them compare to all available bonds
// return bond when all indices match or NULL if no match
gbsaAngleBond* angleBond = gbsaBond->angleBonds;
while( angleBond != NULL ){
if( angleBond->orderedIndices[0] == orderedAtomIndices[0] &&
angleBond->orderedIndices[1] == orderedAtomIndices[1] &&
angleBond->orderedIndices[2] == orderedAtomIndices[2] ){
return angleBond;
}
angleBond = angleBond->nextBond;
}
return NULL;
}
/**---------------------------------------------------------------------------------------
Sort atom indices
@param atomI index of atomI
@param atomJ index of atomJ
@param atomK index of atomK
@param atomL index of atomL
@param orderedIndices output array of ordered indices assumed to be of size 3
@return count (should be 3, if 4, then all the atom indices were distinct)
--------------------------------------------------------------------------------------- */
int GbsaParameters::sortAtomAngleIndices( int atomI, int atomJ,
int atomK, int atomL,
int* orderedIndices ) const {
// ---------------------------------------------------------------------------------------
int orderedAtomIndices[4];
// ---------------------------------------------------------------------------------------
// sort atom indices and then remove index with more than 1 entry
// input should be of the form: (i, j, k, k ), i.e., one duplicate
orderedAtomIndices[0] = atomI;
orderedAtomIndices[1] = atomJ;
orderedAtomIndices[2] = atomK;
orderedAtomIndices[3] = atomL;
size_t numberToSort = 4;
qsort( orderedAtomIndices, numberToSort, sizeof( int ), &integerComparison );
orderedIndices[0] = orderedAtomIndices[0];
int count = 1;
int duplicateCount = 0;
for( int ii = 1; ii < 4; ii++ ){
if( orderedAtomIndices[ii] != orderedAtomIndices[ii-1] && count < 3 ){
orderedIndices[count++] = orderedAtomIndices[ii];
} else {
duplicateCount++;
}
}
return duplicateCount;
}
/**---------------------------------------------------------------------------------------
Get bond counts (Simbios)
@param gbsaBond gbsaBonds to check
@param stretchCount stretch count on return
@param angleCount angle count on return
@return 0
--------------------------------------------------------------------------------------- */
int GbsaParameters::getBondCounts( const gbsaBonds* gbsaBond, int* stretchCount, int* angleCount ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\ngetBondCounts";
// ---------------------------------------------------------------------------------------
*stretchCount = getStretchBondCount( gbsaBond );
*angleCount = getAngleBondCount( gbsaBond );
return 0;
}
/**---------------------------------------------------------------------------------------
Get stretch bond count (Simbios)
@param gbsaBond gbsaBonds to check
@return stretch bond count
--------------------------------------------------------------------------------------- */
int GbsaParameters::getStretchBondCount( const gbsaBonds* gbsaBond ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GbsaParameters::getStretchBondCount";
// ---------------------------------------------------------------------------------------
gbsaStretchBond* stretchBonds = gbsaBond->stretchBonds;
int stretchCount = 0;
while( stretchBonds != NULL ){
stretchCount++;
stretchBonds = stretchBonds->nextBond;
}
return stretchCount;
}
/**---------------------------------------------------------------------------------------
Get angle bond count (Simbios)
@param gbsaBond gbsaBonds to check
@return angle bond count
--------------------------------------------------------------------------------------- */
int GbsaParameters::getAngleBondCount( const gbsaBonds* gbsaBond ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GbsaParameters::getAngleBondCount";
// ---------------------------------------------------------------------------------------
gbsaAngleBond* angleBonds = gbsaBond->angleBonds;
int angleCount = 0;
while( angleBonds != NULL ){
angleCount++;
angleBonds = angleBonds->nextBond;
}
return angleCount;
}
/**---------------------------------------------------------------------------------------
Add stretch (1-2) bonds (Simbios)
@param maxAtoms max number of atoms
@param gbsaBondsArray array of gbsaBonds data structs
@param top Gromacs t_topolgy struct
@param idefArrayIndex index to bond parameters (F_BONDS, F_SHAKE, ... )
@param offset number of entries for each bond block
@param atomIndexOffset offset into block for atom indices
@param log if set, then print error messages to log file
@return 0 if no errors or
return x, where x is the number of errors encountered
--------------------------------------------------------------------------------------- */
int GbsaParameters::addStretchBonds( int maxAtoms, gbsaBonds** gbsaBondsArray,
const t_topology* top, int idefArrayIndex,
int offset, int atomIndexOffset, FILE* log ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GbsaParameters::addStretchBonds";
// ---------------------------------------------------------------------------------------
// load 1-2 bonds
t_iatom* atoms = top->idef.il[idefArrayIndex].iatoms;
t_iparams* params = top->idef.iparams;
int errors = 0;
for( int ii = 0; ii < top->idef.il[idefArrayIndex].nr; ii += offset ){
int type = (int) atoms[ ii ];
int atomI = (int) atoms[ ii + atomIndexOffset ];
int atomJ = (int) atoms[ ii + atomIndexOffset + 1 ];
// validate indices
if( atomI >= maxAtoms || atomI < 0 ){
if( log ){
(void) fprintf( log, "\nAtom indexI=%d (ii=%d) too large: max=%d", atomI, ii, maxAtoms );
}
errors++;
atomI = -1;
}
if( atomJ >= maxAtoms || atomJ < 0 ){
if( log ){
(void) fprintf( log, "\nAtom indexJ=%d (ii=%d) too large: max=%d", atomJ, ii, maxAtoms );
}
errors++;
atomJ = -1;
}
// add new stretch bond, if atomJ and/or atomI are not in stretch bond lists
Real bondLength = params[type].harmonic.rA;
if( atomI >= 0 && atomJ >= 0 ){
if( !isAtomInStretchBondList( atomJ, gbsaBondsArray[atomI], log ) ){
addStretchBond( atomJ, atomI, gbsaBondsArray[atomI], bondLength, log );
}
if( !isAtomInStretchBondList( atomI, gbsaBondsArray[atomJ], log ) ){
addStretchBond( atomI, atomJ, gbsaBondsArray[atomJ], bondLength, log );
}
}
}
return errors;
}
/**---------------------------------------------------------------------------------------
Add SETTLE stretch (1-2) bonds (Simbios)
@param maxAtoms max number of atoms
@param gbsaBondsArray array of gbsaBonds data structs
@param top Gromacs t_topolgy struct
@param idefArrayIndex index to bond parameters (F_BONDS, F_SHAKE, ... )
@param offset number of entries for each bond block
@param atomIndexOffset offset into block for atom indices
@param log if set, then print error messages to log file
@return 0 if no errors or
return x, where x is the number of errors encountered
--------------------------------------------------------------------------------------- */
int GbsaParameters::addSettleStretchBonds( int maxAtoms, gbsaBonds** gbsaBondsArray,
const t_topology* top, int idefArrayIndex,
int offset, int atomIndexOffset, FILE* log ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "GbsaParameters::addStretchBonds";
// ---------------------------------------------------------------------------------------
// load bonds via SETTLE parameter
t_iatom* atoms = top->idef.il[idefArrayIndex].iatoms;
t_iparams* params = top->idef.iparams;
int errors = 0;
// fprintf( log, "\naddBond idx=%d nr=%d off=%d %d", idefArrayIndex, top->idef.il[idefArrayIndex].nr, offset,
// atomIndexOffset );
// fflush( log );
for( int ii = 0; ii < top->idef.il[idefArrayIndex].nr; ii += offset ){
int type = (int) atoms[ ii ];
int atomI = (int) atoms[ ii + atomIndexOffset ];
Real OHBondLength = params[type].harmonic.rA;
Real HHBondLength = params[type].harmonic.krA;
Real bondAngle = 1.0 - ( (HHBondLength*HHBondLength)/ (2.0*OHBondLength*OHBondLength));
bondAngle = ACOS( bondAngle );
//fprintf( log, "\n%d type=%d [%d %d %d] OH=%.3f HH=%.4f angle=%.3f (deg)",
// ii, type, atomI, atomI+1, atomI+2, OHBondLength, HHBondLength, bondAngle*180.0f/M_PI );
//fflush( log );
// validate indices
if( (atomI + 2) >= maxAtoms || atomI < 0 ){
if( log ){
(void) fprintf( log, "\nAtom indexI=%d (ii=%d) too large: max=%d", atomI, ii, maxAtoms );
(void) fflush( log );
}
errors++;
atomI = -1;
}
// add new stretch bond, if atomJ and/or atomI are not in stretch bond lists
if( atomI >= 0 ){
gbsaStretchBond* oxygenStretchBonds[2] = { NULL , NULL };
int hits = 0;
for( int jj = 1; jj <= 2; jj++ ){
if( !isAtomInStretchBondList( atomI+jj, gbsaBondsArray[atomI], log ) ){
hits++;
oxygenStretchBonds[jj-1] = addStretchBond( atomI+jj, atomI, gbsaBondsArray[atomI], OHBondLength, log );
}
if( !isAtomInStretchBondList( atomI, gbsaBondsArray[atomI+jj], log ) ){
addStretchBond( atomI, atomI+jj, gbsaBondsArray[atomI+jj], OHBondLength, log );
}
}
// SETTLE values
// HOH angle ~ 104.5 deg
// OH distance=0.09572
// HH distance=0.15139
// HH distance= sqrt( (1-cos(1.82))*2*OH distance*OH distance)
if( hits == 2 && oxygenStretchBonds[0] && oxygenStretchBonds[1] ){
for( int jj = 0; jj < 3; jj++ ){
gbsaAngleBond* gbsaAngleBond = addStretchBondsToAngleBondList( oxygenStretchBonds[0], oxygenStretchBonds[1],
atomI, gbsaBondsArray, atomI + jj, log );
gbsaAngleBond->harmonicAngleWidth = bondAngle;
}
}
}
}
return errors;
}
/**---------------------------------------------------------------------------------------
Check if atom is in bond list (as atom j) (Simbios)
@param atomIndex index of atom to be searched for
@param bond gbsaBond data struct
@param log if set, then print error messages to log file
@return 0 if atom is in StretchBondList; 1 otherwise
--------------------------------------------------------------------------------------- */
int GbsaParameters::isAtomInStretchBondList( int atomIndex, const gbsaBonds* bond, FILE* log ) const {
// check if atomIndex is in stretch bond list
gbsaStretchBond* stretchBond = bond->stretchBonds;
while( stretchBond != NULL ){
if( stretchBond->atomJ == atomIndex ){
return 1;
}
stretchBond = stretchBond->nextBond;
}
return 0;
}
/**---------------------------------------------------------------------------------------
Add a stretch bonds (Simbios)
@param atomIndexI index of atom I
@param atomIndexJ index of atom J
@param bonds bond data struct
@param bondLength bond length
@param log if set, then print error messages to log file
@return gbsaStretchBond
--------------------------------------------------------------------------------------- */
gbsaStretchBond* GbsaParameters::addStretchBond( int atomIndexJ, int atomIndexI,
gbsaBonds* bonds, Real bondLength, FILE* log ){
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
// allocate memory, set fields and add to stretchBonds list
#ifdef UseGromacsMalloc
gbsaStretchBond* newStretchBond = (gbsaStretchBond*) save_malloc( "gbsaStretchBond", __FILE__, __LINE__, sizeof( gbsaStretchBond ) );
#else
gbsaStretchBond* newStretchBond = new gbsaStretchBond;
#endif
newStretchBond->atomI = atomIndexI;
newStretchBond->atomJ = atomIndexJ;
newStretchBond->bondLength = bondLength;
// push current entry onto list
newStretchBond->nextBond = bonds->stretchBonds;
bonds->stretchBonds = newStretchBond;
bonds->stretchBondCount++;
return newStretchBond;
}
/**---------------------------------------------------------------------------------------
Assign standard radii for GB/SA methods other than ACE;
taken from Macromodel and OPLS-AA, except for hydrogens (Simbios)
Logic follows that in Tinker's ksolv.f
Currently only works for standard amino acid atoms
If invalid atom name is encountered, a message is printed to log file and the
radius for that atom is set to 1.0f
@param numberOfAtoms number of atoms
@param gbsaBondsArray array of gbsaBonds
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
@return 0 always
--------------------------------------------------------------------------------------- */
int GbsaParameters::getMacroModelAtomicRadii( int numberOfAtoms, gbsaBonds** gbsaBondsArray,
char*** atomNames, Real* radii, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGbsaParameters::getMacroModelAtomicRadii";
// ---------------------------------------------------------------------------------------
// loop over atoms
// get Tinker atom number from atom name
// using atom number and bonds12 array, set atom radius
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
int tinkerAtomNumber = mapGmxAtomNameToTinkerAtomNumber( (*(atomNames[atomI])), log );
int stretchBondCount = gbsaBondsArray[atomI]->stretchBondCount;
Real radius;
int bondedAtom;
int bondedTinkerAtomNumber;
switch( tinkerAtomNumber ){
// H
case 1:
// NH?
bondedAtom = gbsaBondsArray[atomI]->stretchBonds->atomJ;
bondedTinkerAtomNumber = mapGmxAtomNameToTinkerAtomNumber( (*(atomNames[bondedAtom])), log );
/*
(void) fprintf( log, "\n H j=%d Name=<%s>", bondedAtom, (*(atomNames[bondedAtom])) );
(void) fprintf( log, " tinkerNo=%d", bondedTinkerAtomNumber );
(void) fflush( log );
*/
if( bondedTinkerAtomNumber == 7 ){
radius = 1.15f;
} else if( bondedTinkerAtomNumber == 8 ){
radius = 1.05f;
} else {
radius = 1.25f;
}
break;
// ?
case 3:
radius = 1.432f;
break;
// C
case 6:
// C-terminal
if( stretchBondCount == 3 ){
radius = 1.875f;
} else if( stretchBondCount == 2 ){
radius = 1.825f;
} else {
radius = 1.90f;
}
break;
// N
case 7:
if( stretchBondCount == 4 ){
radius = 1.625f;
} else if( stretchBondCount == 1 ){
radius = 1.60f;
} else {
radius = 1.7063f;
}
break;
// O
case 8:
if( stretchBondCount == 1 ){
radius = 1.48f;
} else {
radius = 1.535f;
}
break;
case 9:
radius = 1.47f;
break;
case 10:
radius = 1.39f;
break;
case 11:
radius = 1.992f;
break;
case 12:
radius = 1.70f;
break;
case 14:
radius = 1.80f;
break;
case 15:
radius = 1.87f;
break;
case 16:
radius = 1.775f;
break;
case 17:
radius = 1.735f;
break;
case 18:
radius = 1.70f;
break;
case 19:
radius = 2.123f;
break;
case 20:
radius = 1.817f;
break;
case 35:
radius = 1.90f;
break;
case 36:
radius = 1.812f;
break;
case 37:
radius = 2.26f;
break;
case 53:
radius = 2.10f;
break;
case 54:
radius = 1.967f;
break;
case 55:
radius = 2.507f;
break;
case 56:
radius = 2.188f;
break;
default:
radius = 1.0f;
(void) fprintf( log ? log : stderr, "\nWarning: %s Tinker atom number=%d unrecognized -- name=<%s>.",
methodName, tinkerAtomNumber, (*(atomNames[bondedAtom])) );
break;
}
// convert from A to nm
radii[atomI] = 0.1f*radius;
}
return 0;
}
/**---------------------------------------------------------------------------------------
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 GbsaParameters::mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const {
// ---------------------------------------------------------------------------------------
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;
}
// map first letter in atom name to Tinker atom number
int firstAsciiValue = ((int) atomName[0]) - 65;
// check for lower case
if( firstAsciiValue > 25 ){
firstAsciiValue -= 32;
}
// validate
if( firstAsciiValue < 0 || firstAsciiValue > 25 ){
if( log != NULL ){
(void) fprintf( log, "Atom name=<%s> unrecognized.", atomName );
}
(void) fprintf( stderr, "Atom name=<%s> unrecognized.", atomName );
return -1;
}
return atomNameMap[firstAsciiValue];
}
/**---------------------------------------------------------------------------------------
Get atomic volumes minus subvolumes that lie inside directly bonded atoms
Eqs 6 & 7 in J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
Use Macromodel radii
@param numberOfAtoms number of atoms
@param top GMX t_topology struct
@param gbsaBondsArray array of gbsaBonds
@param vdwRadius array of Macromodel radii for each atom
@param volume output array of volumes
@param log if set, then print error messages to log file
@return 0 always; atomic volumes in array 'volume'
--------------------------------------------------------------------------------------- */
int GbsaParameters::getMacroModelAtomicVolumes( int numberOfAtoms, const t_topology* top,
gbsaBonds** gbsaBondsArray, Real* vdwRadius,
Real* volume, FILE* log ){
// ---------------------------------------------------------------------------------------
char buffer[MAX_BUFFER_SZ];
bool debugOn = false;
// static const char* methodName = "\nGbsaParameters::getMacroModelAtomicVolumes";
// ---------------------------------------------------------------------------------------
debugOn = debugOn && log != NULL;
Real pi_3 = ((Real) M_PI)/3.0f;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// include pi/3 after finished summing over h_ij
Real volume_I = 4.0*POW( vdwRadius[atomI], 3.0 );
Real atomI_R = vdwRadius[atomI];
Real atomI_R2 = atomI_R*atomI_R;
// loop over directly bonded neighbors
gbsaStretchBond* stretchBond = gbsaBondsArray[atomI]->stretchBonds;
while( stretchBond != NULL ){
int atomJ = stretchBond->atomJ;
Real bondLength = stretchBond->bondLength*1.01f;
Real atomJ_R = vdwRadius[atomJ];
Real ratio = (atomJ_R*atomJ_R - atomI_R2 - bondLength*bondLength)/(2.0f*atomI_R*bondLength);
Real h_IJ = atomI_R*( 1.0f + ratio );
volume_I -= h_IJ*h_IJ*( 3.0f*atomI_R - h_IJ );
if( debugOn ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( atomJ, top, MAX_BUFFER_SZ, buffer, numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "\n atmJ=%d %s bndL=%.3f jR=%.3f h=%.3f offset=%.3f",
atomJ, buffer, bondLength, atomJ_R, h_IJ, (h_IJ*h_IJ*( 3.0f*atomI_R - h_IJ )) );
}
stretchBond = stretchBond->nextBond;
}
volume[atomI] = volume_I*pi_3;
if( debugOn ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( atomI, top, MAX_BUFFER_SZ, buffer, numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "\n%s %.3f %.3f %d Bonded: ", buffer,
10.0f*vdwRadius[atomI], 1000.0f*volume[atomI], gbsaBondsArray[atomI]->stretchBondCount );
stretchBond = gbsaBondsArray[atomI]->stretchBonds;
while( stretchBond != NULL ){
(void) fprintf( log, " %d", stretchBond->atomJ );
stretchBond = stretchBond->nextBond;
}
}
}
return 0;
}
/**---------------------------------------------------------------------------------------
Calculate first two terms in Gpol equation (nearest and next-nearest neighbors)
Eqs 6 & 7 in J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
Use Macromodel radii
@param numberOfAtoms number of atoms
@param top GMX t_topology struct
@param gbsaBondsArray array of gbsa bonds
@param vdwRadius array of Macromodel radii for each atom
@param volume array of atomic volumes
@param gPolFixed output array of fixed GBSA values for GPol
@param log if set, then print error messages to log file
@return 0 always; fixed value for G_pol in array 'gPolFixed'
--------------------------------------------------------------------------------------- */
int GbsaParameters::getFixedGBSA_GPol( int numberOfAtoms, const t_topology* top,
gbsaBonds** gbsaBondsArray, Real* vdwRadius,
Real* volume, Real* gPolFixed, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGbsaParameters::getFixedGBSA_GPol";
// ---------------------------------------------------------------------------------------
// GBSA parameters from '97 paper
Real phi = getPhi();
Real P1 = getP1();
Real P2 = getP2();
Real P3 = getP3();
Real electricConstant = getElectricConstant();
// ---------------------------------------------------------------------------------------
// 1st term
// convert to nm
Real P1_Phi = 0.1f*(P1 + phi);
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
gPolFixed[atomI] = electricConstant/(vdwRadius[atomI]+P1_Phi);
}
// stretch term (1-2 bonds)
char idBuffer[MAX_BUFFER_SZ];
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// loop over directly bonded neighbors
Real volumeR4Sum = 0.0f;
gbsaStretchBond* stretchBond = gbsaBondsArray[atomI]->stretchBonds;
while( stretchBond != NULL ){
int atomJ = stretchBond->atomJ;
Real volumeJ_R = volume[atomJ];
Real bondLength = stretchBond->bondLength;
Real r4 = POW( bondLength, 4.0 );
if( r4 > 0.0 ){
volumeR4Sum += volumeJ_R/r4;
} else if( log != NULL ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( atomI, top, MAX_BUFFER_SZ, idBuffer, numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "%s %s 1-2 Bonded atom atomIndex=%d has bond length=0.",
methodName, idBuffer, atomJ );
}
stretchBond = stretchBond->nextBond;
}
gPolFixed[atomI] += P2*volumeR4Sum;
}
// bend term (1-3 bonds)
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// loop over angles
Real volumeR4Sum = 0.0f;
gbsaAngleBond* angleBond = gbsaBondsArray[atomI]->angleBonds;
while( angleBond != NULL ){
gbsaStretchBond* gbsaStretchBond1 = angleBond->stretchBondI;
gbsaStretchBond* gbsaStretchBond2 = angleBond->stretchBondJ;
Real bondLength1 = gbsaStretchBond1->bondLength;
Real bondLength2 = gbsaStretchBond2->bondLength;
Real cosOfAngle = cos( angleBond->harmonicAngleWidth );
Real distance2 = bondLength1*bondLength1 +
bondLength2*bondLength2 -
2.0f*bondLength1*bondLength2*cosOfAngle;
Real r4 = distance2*distance2;
if( r4 > 0.0 ){
volumeR4Sum += volume[angleBond->volumeIndex]/r4;
} else if( log != NULL ){
SimTKOpenMMGromacsUtilities::getAtomIdStringGivenAtomIndex( atomI, top, MAX_BUFFER_SZ, idBuffer, numberOfAtoms, ATOM_ID_STRING_TAB );
(void) fprintf( log, "%s %s angle-bonded atom has distance 0.",
methodName, idBuffer );
}
angleBond = angleBond->nextBond;
}
gPolFixed[atomI] += P3*volumeR4Sum;
}
return 0;
}
/**---------------------------------------------------------------------------------------
Get exclusions for specific atom (Simbios)
@param numberOfAtoms number of atoms
@param gbsaBondsArray array of gbsa bonds
@param atomI atom index of atom for which exclusions are to be set
@param exclusionWorkArray exclusionWorkArray[j] = 1, if atom j is to be excluded
@param value may be null on input in which space is allocated
@param previousIndex previousIndex -- if < 0, then iniitialize all entries to
@param 0
@param log if set, then print error messages to log file
@return 0
Abort if exclusionWorkArray == NULL
--------------------------------------------------------------------------------------- */
int GbsaParameters::getExclusionsForAtom( int numberOfAtoms, const gbsaBonds** gbsaBondsArray,
int atomI, int* exclusionWorkArray, int previousIndex,
FILE* log ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGbsaParameters::getExclusionsForAtom";
// ---------------------------------------------------------------------------------------
// allocate memory if exclusionWorkArray is null
if( exclusionWorkArray == NULL ){
(void) fprintf( log ? log : stderr, "%s error exclusionWorkArray is NULL.", methodName );
(void) fflush( log ? log : stderr );
exit(-1);
return -1;
}
// set exclusion factors
if( previousIndex < 0 ){
if( previousIndex == -1 ){
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionWorkArray[ii] = -1;
}
} else {
memset( exclusionWorkArray, 0, numberOfAtoms*sizeof( int ) );
}
} else {
setExclusionValue( gbsaBondsArray[previousIndex], -1, exclusionWorkArray, log );
}
setExclusionValue( gbsaBondsArray[atomI], atomI, exclusionWorkArray, log );
return 0;
}
/**---------------------------------------------------------------------------------------
Set exclusions for specific atom (Simbios)
@param gbsaBonds gbsa bond
@param setValue set value
@param exclusionWorkArray exclusionWorkArray[j] = 1, if atom j is to be excluded
@param value may be null on input in which space is allocated
@param log if set, then print error messages to log file
@return array of Born radii
--------------------------------------------------------------------------------------- */
int GbsaParameters::setExclusionValue( const gbsaBonds* gbsaBonds, int setValue,
int* exclusionWorkArray, FILE* log ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::setExclusionsValue";
// ---------------------------------------------------------------------------------------
gbsaStretchBond* stretchBond = gbsaBonds->stretchBonds;
while( stretchBond != NULL ){
exclusionWorkArray[stretchBond->atomI] = setValue;
exclusionWorkArray[stretchBond->atomJ] = setValue;
stretchBond = stretchBond->nextBond;
}
gbsaAngleBond* angleBond = gbsaBonds->angleBonds;
while( angleBond != NULL ){
exclusionWorkArray[angleBond->stretchBondI->atomI] = setValue;
exclusionWorkArray[angleBond->stretchBondI->atomJ] = setValue;
exclusionWorkArray[angleBond->stretchBondJ->atomI] = setValue;
exclusionWorkArray[angleBond->stretchBondJ->atomJ] = setValue;
angleBond = angleBond->nextBond;
}
return 0;
}
/**---------------------------------------------------------------------------------------
Get array for accessing exclusion entries (Simbios)
(work array)
@return array of size _numberOfAtoms (int)
--------------------------------------------------------------------------------------- */
int* GbsaParameters::getExclusionWorkArray( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::getExclusionWorkArray";
// ---------------------------------------------------------------------------------------
if( _exclusionWorkArray == NULL && _numberOfAtoms > 0 ){
#ifdef UseGromacsMalloc
_exclusionWorkArray = (int*) save_malloc( "_exclusionWorkArray", __FILE__, __LINE__, sizeof( int )*_numberOfAtoms );
#else
_exclusionWorkArray = new int[_numberOfAtoms];
#endif
memset( _exclusionWorkArray, 0, _numberOfAtoms*sizeof( int ) );
}
return _exclusionWorkArray;
}
/**---------------------------------------------------------------------------------------
Print state to log file (Simbios)
@param title title (optional)
@param log print state to log file
@return 0 always
--------------------------------------------------------------------------------------- */
int GbsaParameters::logState( const char* title, FILE* log ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGbsaParameters::logState";
// ---------------------------------------------------------------------------------------
if( log == NULL ){
return 0;
}
if( title ){
(void) fprintf( log, "%s", title );
}
(void) fprintf( log, "\nGbsaParameter info:\n" );
(void) fprintf( log, "\n Solvent dielectric: %12.4f", _solventDielectric );
(void) fprintf( log, "\n Inner dielectric: %12.4f", _innerDielectric );
(void) fprintf( log, "\n Electric constant: %12.4f", _electricConstant );
(void) fprintf( log, "\n Gpol prefactor: %12.4f", _preFactor );
(void) fprintf( log, "\n Probe radius: %12.4f", _probeRadius );
(void) fprintf( log, "\n Phi: %12.4f", _phi );
(void) fprintf( log, "\n P1: %12.4f", _P1 );
(void) fprintf( log, "\n P2: %12.4f", _P2 );
(void) fprintf( log, "\n P3: %12.4f", _P3 );
(void) fprintf( log, "\n P4: %12.4f", _P4 );
(void) fprintf( log, "\n P5: %12.4f", _P5 );
(void) fprintf( log, "\n Number of atoms: %12d", _numberOfAtoms );
(void) fprintf( log, "\n Free arrays: %12d", _freeArrays );
return 0;
}
/* ---------------------------------------------------------------------------------------
Override C++ new w/ Gromac's smalloc/sfree (Simbios)
@param size bytes to allocate
@return ptr to allocated memory
--------------------------------------------------------------------------------------- */
void* GbsaParameters::operator new( size_t size ){
void *ptr;
smalloc(ptr, (int) size);
/*
(void) fprintf( stdout, "\nGbsaParameters new called -- size=%u", size );
(void) fflush( stdout );
*/
return ptr;
}
/* ---------------------------------------------------------------------------------------
Override C++ delete w/ Gromac's sfree (Simbios)
@param ptr ptr to block to free
--------------------------------------------------------------------------------------- */
void GbsaParameters::operator delete( void *ptr ){
/*
(void) fprintf( stdout, "\nGbsaParameters delete called." );
(void) fflush( stdout );
*/
sfree( ptr );
}
/* ---------------------------------------------------------------------------------------
Override C++ new w/ Gromac's smalloc/sfree (Simbios)
@param size bytes to allocate
@return ptr to allocated memory
--------------------------------------------------------------------------------------- */
void* GbsaParameters::operator new[]( size_t size ){
void *ptr;
smalloc(ptr, (int) size);
/*
(void) fprintf( stdout, "\nGbsaParameters new[] called -- size=%u", size );
(void) fflush( stdout );
*/
return ptr;
}
/* ---------------------------------------------------------------------------------------
Override C++ delete w/ Gromac's sfree (Simbios)
@param ptr ptr to block to free
--------------------------------------------------------------------------------------- */
void GbsaParameters::operator delete[]( void *ptr ){
/*
(void) fprintf( stdout, "\nGbsaParameters delete[] called." );
(void) fflush( stdout );
*/
sfree( ptr );
}
/**---------------------------------------------------------------------------------------
Get Vdw radii from parameter file (Simbios)
@param numberOfAtoms number of atoms
@param parameterFileName parameterFileName
@param top Gromacs topology data struct
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
@return 0 always
--------------------------------------------------------------------------------------- */
int GbsaParameters::getMacroModelAtomicRadii( int numberOfAtoms, const std::string parameterFileName,
const t_topology* top, char*** atomNames,
Real* radii, FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nGbsaParameters::getMacroModelAtomicRadii";
// ---------------------------------------------------------------------------------------
StringVector* fileContents = readFile( parameterFileName );
GbsaAtomParameterMap* parameterMap = parseParameterFile( *fileContents );
delete fileContents;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
char* atomTypeC = *(top->atoms.atomtype[ii]);
std::string atomType = atomTypeC;
GbsaAtomParameterMapCI entry = parameterMap->find( atomType );
if( entry != parameterMap->end() ){
GbsaAtomParameter* gbsaAtomParameter = (*entry).second;
radii[ii] = (Real) gbsaAtomParameter->getVdwRadius();
radii[ii] *= 0.1;
std::stringstream message;
message << " type found for atom=<" << (*(atomNames[ii])) << "> type=<" << atomType << "> r=" << radii[ii];
(void) fprintf( log, "%s\n" , message.str().c_str() );
} else {
radii[ii] = (Real) 0.0001f;
std::stringstream message;
message << methodName;
message << " no type found for atom=<" << (*(atomNames[ii])) << "> type=<" << atomType << ">";
(void) fprintf( log, "%s\n" , message.str().c_str() );
}
}
delete parameterMap;
return 0;
}
/**---------------------------------------------------------------------------------------
Read file (Simbios)
@param fileName file name
@return vector of strings, one string per line; return NULL if file is not found
--------------------------------------------------------------------------------------- */
StringVector* GbsaParameters::readFile( const std::string& fileName ){
// ---------------------------------------------------------------------------------------
static int bufferSize = 2048;
static char buffer[2048];
static const std::string methodName = "\nGbsaParameters::readParameterFile";
// ---------------------------------------------------------------------------------------
// open file
FILE* file = NULL;
#ifdef WIN32
fopen_s( &file, fileName.c_str(), "r" );
#else
file = fopen( fileName.c_str(), "r" );
#endif
if( file != NULL ){
std::stringstream message;
message << methodName.c_str() << " opened file=<" << fileName.c_str() << ">.";
// AmoebaLog::printMessage( message );
} else {
std::stringstream message;
message << methodName.c_str() << " could not open file=<" << fileName.c_str() << ">.";
(void) printf( "\n%s\n", message.str().c_str() );
(void) fflush( stdout );
// AmoebaLog::printMessage( message );
exit(-1);
return NULL;
}
// loop over all lines in file, checking for end-of-file
int lineNumber = 0;
StringVector* fileContents = new StringVector();
while( !feof( file ) ){
// read next line
int bufferLen;
lineNumber++;
if( fgets( buffer, bufferSize, file ) != NULL ){
bufferLen = (int) strlen( buffer );
if( bufferLen > 0 ){
buffer[bufferLen-1] = '\0';
}
// (void) fprintf( log, "%s", buffer );
// (void) fflush( log );
fileContents->push_back( buffer );
}
}
// done
(void) fclose( file );
if( file ){
//std::stringstream message;
//message << methodName.c_str() << " read " << lineNumber << " lines from file=<" << fileName->c_str() << ">.";
// AmoebaLog::printMessage( message );
}
return fileContents;
}
/**---------------------------------------------------------------------------------------
Parse parameter file (Simbios)
@param fileContents file contents
@return vector of strings, one string per line
--------------------------------------------------------------------------------------- */
GbsaAtomParameterMap* GbsaParameters::parseParameterFile( const StringVector& fileContents ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nGbsaParameters::parseParameterFile";
static const std::string delimiter = " ";
// ---------------------------------------------------------------------------------------
GbsaAtomParameterMap* map = new GbsaAtomParameterMap();
for( StringVectorCI ii = fileContents.begin(); ii != fileContents.end(); ii++ ){
std::string line = *ii;
StringVector tokenVector;
tokenizeString( line, tokenVector, delimiter, 0 );
if( tokenVector[0] != "#" && tokenVector[0] != "@" ){
GbsaAtomParameter* gbsaAtomParameter = new GbsaAtomParameter( tokenVector );
(*map)[tokenVector[0]] = gbsaAtomParameter;
}
}
return map;
}
/**---------------------------------------------------------------------------------------
Tokenize a string (static method) (Simbios)
@param line string to tokenize
@param tokenVector upon return vector of tokens
@param delimiter token delimter
@param clearTokenVector if true, clear tokenVector
@return 1
--------------------------------------------------------------------------------------- */
int GbsaParameters::tokenizeString( const std::string& line, StringVector& tokenVector,
const std::string& delimiter, int clearTokenVector ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nGbsaParameters::tokenizeString";
// ---------------------------------------------------------------------------------------
char *ptr_c;
if( clearTokenVector ){
tokenVector.clear();
}
char* lineBuffer = (char*) malloc( 2048*sizeof( char) );
#ifdef WIN32
(void) sprintf_s( lineBuffer, 2048, "%s", line.c_str() );
#else
(void) sprintf( lineBuffer, "%s", line.c_str() );
#endif
while( (ptr_c = strsep( &lineBuffer, delimiter.c_str() )) != NULL ){
if( *ptr_c ){
tokenVector.push_back( std::string( ptr_c ) );
}
}
free( lineBuffer );
return 0;
}
/**---------------------------------------------------------------------------------------
Replacement of sorts for strtok() (static method) (Simbios)
Used to parse parameter file lines
Should be moved to Utilities file
@param lineBuffer string to tokenize
@param delimiter token delimter
@return number of args; if return value equals maxTokens, then more tokens than allocated
--------------------------------------------------------------------------------------- */
char* GbsaParameters::strsep( char** lineBuffer, const char* delimiter ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nTinkerParameterSet::strsep"
char *s;
const char *spanp;
int c, sc;
char *tok;
// ---------------------------------------------------------------------------------------
s = *lineBuffer;
if( s == NULL ){
return (NULL);
}
for( tok = s;; ){
c = *s++;
spanp = delimiter;
do {
if( (sc = *spanp++) == c ){
if( c == 0 ){
s = NULL;
} else {
s[-1] = 0;
}
*lineBuffer = s;
return( tok );
}
} while( sc != 0 );
}
}
/**---------------------------------------------------------------------------------------
Qsort/heapsort integer comparison (Simbios)
@param a first value to compare
@param b second value to compare
@return 1, 0, -1 based on comparison
--------------------------------------------------------------------------------------- */
int integerComparison( const void *a, const void *b){
int diff = *(int *) a - *(int *) b;
if( diff < 0 ){
return -1;
} else if( diff > 0 ){
return 1;
} else {
return 0;
}
}
platforms/reference/src/Gbsa/GbsaParameters.h 0 → 100755
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __GbsaParameters_H__
#define __GbsaParameters_H__
// Gromacs data structs
#ifndef M_PI
#define M_PI 3.1415926535
#endif
#include "typedefs.h"
// ---------------------------------------------------------------------------------------
// buffer size of atom id strings
#define MAX_BUFFER_SZ 128
#include "GbsaAtomParameter.h"
struct _gbsaStretchBondStruct {
int atomI;
int atomJ;
Real bondLength;
struct _gbsaStretchBondStruct* nextBond;
};
typedef struct _gbsaStretchBondStruct gbsaStretchBond;
struct _gbsaAngleBondStruct {
gbsaStretchBond *stretchBondI;
gbsaStretchBond *stretchBondJ;
int pivotAtomIndex;
int volumeIndex;
int orderedIndices[3];
Real harmonicAngleWidth;
struct _gbsaAngleBondStruct* nextBond;
};
typedef struct _gbsaAngleBondStruct gbsaAngleBond;
#define MAX_BOND_EXCLUSIONS 20
struct _gbsaBondStruct {
int stretchBondCount;
gbsaStretchBond* stretchBonds;
int angleBondCount;
gbsaAngleBond* angleBonds;
int numberOfExclusions;
int exclusions[MAX_BOND_EXCLUSIONS];
int minExclusionIndex;
int maxExclusionIndex;
int startExclusionIndex;
int stopExclusionIndex;
};
typedef struct _gbsaBondStruct gbsaBonds;
/*
class GbsaAtomInfo {
private:
static const int MaxBondExclusions = 20;
int _stretchBondCount;
gbsaStretchBond* _stretchBonds;
int _angleBondCount;
gbsaAngleBond* _angleBonds;
int _numberOfExclusions;
int _exclusions[MaxBondExclusions];
int _minExclusionIndex;
int _maxExclusionIndex;
int _startExclusionIndex;
int _stopExclusionIndex;
public:
// constructor/destructor
GbsaAtomInfo( );
~GbsaAtomInfo( );
// accessors
int getStretchBondCount( void ){ return _stretchBondCount; };
gbsaStretchBond* getGbsaStretchBond( void ){ return _stretchBonds; };
int getAngleBondCount( void ){ return _angleBondCount; };
gbsaAngleBond* getGbsaAngleBond( void ){ return _angleBonds; };
int getNumberOfExclusions( void ){ return _numberOfExclusions; };
int* get_exclusions( void ){ return _exclusions; };
int getMinExclusionIndex( void ){ return _minExclusionIndex; };
int getMaxExclusionIndex( void ){ return _maxExclusionIndex; };
int getStartExclusionIndex( void ){ return _startExclusionIndex; };
int getStopExclusionIndex( void ){ return _stopExclusionIndex; };
};
*/
class GbsaParameters {
private:
// GBSA constants & parameters; parameters from '97 paper
Real _solventDielectric;
Real _innerDielectric;
Real _kcalA_To_kJNm;
Real _phi;
Real _P1;
Real _P2;
Real _P3;
Real _P4;
Real _P4_2;
Real _P5;
Real _electricConstant;
Real _probeRadius;
Real _preFactor;
Real _P5Inverse;
Real _piP5;
Real _pi4Asolv;
// ---------------------------------------------------------------------------------------
int _numberOfAtoms;
FILE* _log;
bool _freeArrays;
Real* _vdwRadii;
Real* _volume;
Real* _gPolFixed;
Real* _bornRadii;
// bond and exclusion info
gbsaBonds** _gbsaBondsArray;
// work array
int* _exclusionWorkArray;
// currently not exposed
int setGbsaBondsArray( gbsaBonds** gbsaBondsArray ){ _gbsaBondsArray = gbsaBondsArray; return 0; };
/**---------------------------------------------------------------------------------------
Reset prefactor (Simbios)
called when _electricConstant, _innerDielectric, or _solventDielectric are modified
--------------------------------------------------------------------------------------- */
void _resetPreFactor( void );
public:
// constructor/destructor
GbsaParameters( int numberOfAtoms, FILE* log );
~GbsaParameters( );
// override of new/delete
static void* operator new( size_t size );
static void operator delete( void *p );
static void* operator new[]( size_t size );
static void operator delete[]( void *p );
// accessors
int getNumberOfAtoms( void ){ return _numberOfAtoms; };
FILE* getLog( void ){ return _log; };
int setLog( FILE* log ){ _log = log; return 0; };
// array accessors
Real* getVdwRadii( void ){ return _vdwRadii; };
Real* getVolume( void ){ return _volume; };
Real* getGPolFixed( void ){ return _gPolFixed; };
Real* getBornRadii( void ){ return _bornRadii; };
int setVdwRadii( Real* vdwRadii ){ _vdwRadii = vdwRadii; return 0; };
int setVolume( Real* volume ){ _volume = volume; return 0; };
int setGPolFixed( Real* gPolFixed ){ _gPolFixed = gPolFixed; return 0; };
int setBornRadii( Real* bornRadii ){ _bornRadii = bornRadii; return 0; };
// flag signalling whether arrays should be deleted
// if destructor is called
int setFreeArrays( bool freeArrays ){ _freeArrays = freeArrays; return 0; };
// work array
int* getExclusionWorkArray( void );
// bond info
gbsaBonds** getGbsaBondsArray( void ){ return _gbsaBondsArray; };
// get accessor for constants
Real getSolventDielectric( void ){ return _solventDielectric; };
Real getInnerDielectric( void ){ return _innerDielectric; };
Real getKcalA_To_kJNm( void ){ return _kcalA_To_kJNm; };
Real getPhi( void ){ return _phi; };
Real getP1( void ){ return _P1; };
Real getP2( void ){ return _P2; };
Real getP3( void ){ return _P3; };
Real getP4( void ){ return _P4; };
Real getP4_2( void ){ return _P4_2; };
Real getP5( void ){ return _P5; };
Real getElectricConstant( void ){ return _electricConstant; };
Real getProbeRadius( void ){ return _probeRadius; };
Real getPreFactor( void ){ return _preFactor; };
Real getP5Inverse( void ){ return _P5Inverse; };
Real getPiP5( void ){ return _piP5; };
Real getPi4Asolv( void ){ return _pi4Asolv; };
/**---------------------------------------------------------------------------------------
Set solvent dielectric (Simbios)
@param solventDielectric solvent dielectric
@return 0
--------------------------------------------------------------------------------------- */
int setSolventDielectric( Real solventDielectric );
/**---------------------------------------------------------------------------------------
Set inner dielectric (Simbios)
@param innerDielectric inner dielectric
@return 0
--------------------------------------------------------------------------------------- */
int setInnerDielectric( Real innerDielectric );
/**---------------------------------------------------------------------------------------
Set electric constant (Simbios)
@param electricConstant electric constant
@return 0
--------------------------------------------------------------------------------------- */
int setElectricConstant( Real electricConstant );
/**---------------------------------------------------------------------------------------
Initialize Gbsa Constants (Simbios)
--------------------------------------------------------------------------------------- */
void initializeConstants( void );
/**---------------------------------------------------------------------------------------
Initialize Gbsa Parameters (Simbios)
@param top GMX topology data struct
@return 0
--------------------------------------------------------------------------------------- */
int initializeParameters( const t_topology* top );
/**---------------------------------------------------------------------------------------
Set Gbsa Exclusions (Simbios)
@param top GMX topology data struct -- used to get harmonic angle
@param printOn print diagnostics
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int setGbsaExclusions( const t_topology* top, bool printOn, FILE* log );
/**---------------------------------------------------------------------------------------
Allocate memory for bond array (Simbios)
@param maxAtoms max number of atoms
array entries are intialized to zero
free array[0] and then array when done
@return ptr to allocated array or NULL if out of memory
--------------------------------------------------------------------------------------- */
gbsaBonds** allocateBondsArray( int maxAtoms );
/**---------------------------------------------------------------------------------------
Deallocate memory for bond array (Simbios)
@param maxAtoms max number of atoms
@param gbsaBondsArray array to be freed
@return 0
--------------------------------------------------------------------------------------- */
int freeBondArray( int maxAtoms, gbsaBonds** gbsaBondsArray );
/**---------------------------------------------------------------------------------------
Print bond array (Simbios)
@param numberOfAtoms number of atoms
@param top GMX topology data struct -- used to get harmonic angle
@param gbsaBondsArray bond array
@param title title string (optional)
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int printBondsArray( int numberOfAtoms, const t_topology* top,
gbsaBonds** gbsaBondsArray,
const char* title, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Print bond (Simbios)
@param numberOfAtoms number of atoms
@param top GMX topology data struct -- used to get harmonic angle
@param atomIndex atom index
@param gbsaBond bond data struct
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int printBond( int numberOfAtoms, const t_topology* top, int atomIndex,
const gbsaBonds* gbsaBond, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Print parameter info (Simbios)
@param numberOfAtoms number of atoms
@param top GMX topology data struct -- used to get harmonic angle
@param parameterInfoFileName parameterInfo FileName
@param log file descriptor
@return 0
--------------------------------------------------------------------------------------- */
int printParameterInfo( int numberOfAtoms, const t_topology* top,
const char* parameterInfoFileName, FILE* log );
/**---------------------------------------------------------------------------------------
Add angle bonds using stretch bond lists (Simbios)
Also load in angle width
@param maxAtoms max number of atoms
@param gbsaBondsArray array of gbsaBonds
@param top GMX topology data struct -- used to get harmonic angle
@param idefArrayIndex parameter index for GMX iparams data struct
@param log if set, then print error messages to log file
@return 0 if no errors or
--------------------------------------------------------------------------------------- */
int addAngleBonds( int maxAtoms, gbsaBonds** gbsaBondsArray,
const t_topology* top, int idefArrayIndex, FILE* log );
/**---------------------------------------------------------------------------------------
Add stretch bonds to angle data struct (Simbios)
@param stretchBondI stretch bond I
@param stretchBondJ stretch bond J
@param pivotAtomIndex index of shared atom
@param gbsaBonds gbsaBond
@param saveIndex index to save????
@param log if set, then print error messages to log file
@return allocated gbsaAngleBond
--------------------------------------------------------------------------------------- */
gbsaAngleBond* addStretchBondsToAngleBondList( gbsaStretchBond* stretchBondI,
gbsaStretchBond* stretchBondJ,
int pivotAtomIndex,
gbsaBonds** gbsaBonds, int saveIndex, FILE* log );
/**---------------------------------------------------------------------------------------
Find angle bond with atom indices that match input angle indices (Simbios)
@param gbsaBond gbsaBond
@param atomI index of atomI
@param atomJ index of atomJ
@param atomK index of atomK
@return allocated gbsaAngleBond
--------------------------------------------------------------------------------------- */
gbsaAngleBond* findAngleBond( const gbsaBonds* gbsaBond, int atomI, int atomJ, int atomK ) const;
/**---------------------------------------------------------------------------------------
Find angle bond with atom indices that match input ordered angle indices (Simbios)
gbsaBond gbsaBond
orderedAtomIndices array of ordered indices
return gbsaAngleBond if found
--------------------------------------------------------------------------------------- */
gbsaAngleBond* findAngleBondGivenOrderedArray( const gbsaBonds* gbsaBond, int* orderedAtomIndices ) const;
/**---------------------------------------------------------------------------------------
Sort atom indices
@param atomI index of atomI
@param atomJ index of atomJ
@param atomK index of atomK
@param atomL index of atomL
@param orderedIndices output array of ordered indices assumed to be of size 3
@return count (should be 3, if 4, then all the atom indices were distinct)
--------------------------------------------------------------------------------------- */
int sortAtomAngleIndices( int atomI, int atomJ, int atomK, int atomL,
int* orderedIndices ) const;
/**---------------------------------------------------------------------------------------
Get bond counts (Simbios)
@param gbsaBond gbsaBonds to check
@param stretchCount stretch count on return
@param angleCount angle count on return
@return 0
--------------------------------------------------------------------------------------- */
int getBondCounts( const gbsaBonds* gbsaBond, int* stretchCount, int* angleCount ) const;
/**---------------------------------------------------------------------------------------
Get stretch bond count (Simbios)
@param gbsaBond gbsaBonds to check
@return stretch bond count
--------------------------------------------------------------------------------------- */
int getStretchBondCount( const gbsaBonds* gbsaBond ) const;
/**---------------------------------------------------------------------------------------
Get angle bond count (Simbios)
@param gbsaBond gbsaBonds to check
@return angle bond count
--------------------------------------------------------------------------------------- */
int getAngleBondCount( const gbsaBonds* gbsaBond ) const;
/**---------------------------------------------------------------------------------------
Add stretch bonds (Simbios)
@param maxAtoms max number of atoms
@param gbsaBondsArray array of gbsaBonds data structs
@param top Gromacs t_topolgy struct
@param idefArrayIndex index to bond parameters (F_BONDS, F_SHAKE, ... )
@param offset number of entries for each bond block
@param atomIndexOffset offset into block for atom indices
@param log if set, then print error messages to log file
@return 0 if no errors or
return x, where x is the number of errors encountered
--------------------------------------------------------------------------------------- */
int addStretchBonds( int maxAtoms, gbsaBonds** gbsaBondsArray,
const t_topology* top, int idefArrayIndex, int offset,
int atomIndexOffset, FILE* log );
/**---------------------------------------------------------------------------------------
Add SETTLE stretch bonds (Simbios)
@param maxAtoms max number of atoms
@param gbsaBondsArray array of gbsaBonds data structs
@param top Gromacs t_topolgy struct
@param idefArrayIndex index to bond parameters (F_BONDS, F_SHAKE, ... )
@param offset number of entries for each bond block
@param atomIndexOffset offset into block for atom indices
@param log if set, then print error messages to log file
@return 0 if no errors or
return x, where x is the number of errors encountered
--------------------------------------------------------------------------------------- */
int addSettleStretchBonds( int maxAtoms, gbsaBonds** gbsaBondsArray,
const t_topology* top, int idefArrayIndex, int offset,
int atomIndexOffset, FILE* log );
/**---------------------------------------------------------------------------------------
Check if atom is in bond list (as atom j) (Simbios)
@param atomIndex index of atom to be searched for
@param bond gbsaBond data struct
@param log if set, then print error messages to log file
@return 0 if atom is in StretchBondList; 1 otherwise
--------------------------------------------------------------------------------------- */
int isAtomInStretchBondList( int atomIndex, const gbsaBonds* bond, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Add a stretch bonds (Simbios)
@param atomIndexI index of atom I
@param atomIndexJ index of atom J
@param bonds gpu bond data struct
@param bondLength bond length
@param log if set, then print error messages to log file
@return gbsaStretchBond
--------------------------------------------------------------------------------------- */
gbsaStretchBond* addStretchBond( int atomIndexJ, int atomIndexI,
gbsaBonds* bonds, Real bondLength, FILE* log );
/**---------------------------------------------------------------------------------------
Assign standard radii for GB/SA methods other than ACE;
taken from Macromodel and OPLS-AA, except for hydrogens (Simbios)
Logic based on logic in Tinker's ksolv.f
Currently only works for standard amino acid atoms
If invalid atom name is encountered, a message is printed to log file and the
radius for that atom is set to 1.0f
@param numberOfAtoms number of atoms
@param gbsaBondsArray array of gbsaBonds
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
@return 0 always
--------------------------------------------------------------------------------------- */
int getMacroModelAtomicRadii( int numberOfAtoms, gbsaBonds** gbsaBondsArray,
char*** atomNames, Real* radii, FILE* log );
/**---------------------------------------------------------------------------------------
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 atomic volumes minus subvolumes that lie inside directly bonded atoms
Eqs 6 & 7 in J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
Use Macromodel radii
@param numberOfAtoms number of atoms
@param top GMX t_topology struct
@param gbsaBondsArray array of gbsaBonds
@param vdwRadius array of Macromodel radii for each atom
@param volume output array of volumes
@param log if set, then print error messages to log file
@return 0 always; atomic volumes in array 'volume'
--------------------------------------------------------------------------------------- */
int getMacroModelAtomicVolumes( int numberOfAtoms, const t_topology* top,
gbsaBonds** gbsaBondsArray, Real* vdwRadius,
Real* volume, FILE* log );
/**---------------------------------------------------------------------------------------
Calculate first two terms in Gpol equation (nearest and next-nearest neighbors)
Eqs 6 & 7 in J. Phys. Chem. A V101 No 16, p. 3005 (Simbios)
Use Macromodel radii
@param numberOfAtoms number of atoms
@param top GMX t_topology struct
@param gbsaBondsArray array of gpu bonds
@param vdwRadius array of Macromodel radii for each atom
@param volume array of atomic volumes
@param gPolFixed output array of fixed GBSA values for GPol
@param log if set, then print error messages to log file
return 0 always; fixed value for G_pol in array 'gPolFixed'
--------------------------------------------------------------------------------------- */
int getFixedGBSA_GPol( int numberOfAtoms, const t_topology* top,
gbsaBonds** gbsaBondsArray,
Real* vdwRadius, Real* volume, Real* gPolFixed, FILE* log );
/**---------------------------------------------------------------------------------------
Get exclusions for specific atom (Simbios)
@param numberOfAtoms number of atoms
@param gbsaBondsArray array of gpu bonds
@param atomI atom index of atom for which exclusions are to be set
@param exclusionWorkArray exclusionWorkArray[j] = 1, if atom j is to be excluded
value may be null on input in which space is allocated
@param previousIndex previousIndex -- if < 0, then iniitialize all entries to
0
@param log if set, then print error messages to log file
@return 0 if no problems
--------------------------------------------------------------------------------------- */
int getExclusionsForAtom( int numberOfAtoms, const gbsaBonds** gbsaBondsArray,
int atomI, int* exclusionWorkArray, int previousIndex, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Set exclusions for specific atom (Simbios)
@param gbsaBonds gpu bond
@param setValue set value
@param exclusionWorkArray exclusionWorkArray[j] = 1, if atom j is to be excluded
value may be null on input in which space is allocated
@param log if set, then print error messages to log file
@return array of Born radii
--------------------------------------------------------------------------------------- */
int setExclusionValue( const gbsaBonds* gbsaBonds, int setValue, int* exclusionWorkArray, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Collect Gbsa Bond Indices for a particular atom (Simbios)
@param gbsaBond gbsaBond
@param collectionBinMax size of collectionBin array
@param collectionBin array of indices to be excluded
@param excludeIndex exclude index
@param log if set, then print error messages to log file
@return number of collected indices (includes excludeIndex at end of array)
--------------------------------------------------------------------------------------- */
int collectGbsaBondIndices( gbsaBonds* gbsaBond, int collectionBinMax,
int* collectionBin, int excludeIndex, FILE* log );
/**---------------------------------------------------------------------------------------
Print state to log file (Simbios)
@param title title (optional)
@param log print state to log file
@return 0 always;
--------------------------------------------------------------------------------------- */
int logState( const char* title, FILE* log ) const;
/**---------------------------------------------------------------------------------------
Get Vdw radii from parameter file (Simbios)
@param numberOfAtoms number of atoms
@param parameterFileName parameterFileName
@param top Gromacs topology data struct
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
@return 0 always
--------------------------------------------------------------------------------------- */
int getMacroModelAtomicRadii( int numberOfAtoms, const std::string parameterFileName,
const t_topology* top, char*** atomNames,
Real* radii, FILE* log );
/**---------------------------------------------------------------------------------------
Read file (Simbios)
@param fileName file name
@return vector of strings, one string per line
--------------------------------------------------------------------------------------- */
static StringVector* readFile( const std::string& fileName );
/**---------------------------------------------------------------------------------------
Parse parameter file (Simbios)
@param fileContents file contents
@return vector of strings, one string per line
--------------------------------------------------------------------------------------- */
static GbsaAtomParameterMap* parseParameterFile( const StringVector& fileContents );
/**---------------------------------------------------------------------------------------
Tokenize a string (static method) (Simbios)
@param line string to tokenize
@param tokenVector upon return vector of tokens
@param delimiter token delimter
@param clearTokenVector if true, clear tokenVector
@return 1
--------------------------------------------------------------------------------------- */
static int tokenizeString( const std::string& line, StringVector& tokenVector,
const std::string& delimiter, int clearTokenVector );
/**---------------------------------------------------------------------------------------
Replacement of sorts for strtok() (static method) (Simbios)
Used to parse parameter file lines
Should be moved to Utilities file
@param lineBuffer string to tokenize
@param delimiter token delimter
@return number of args; if return value equals maxTokens, then more tokens than allocated
--------------------------------------------------------------------------------------- */
static char* strsep( char** lineBuffer, const char* delimiter );
};
/**---------------------------------------------------------------------------------------
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 // __GbsaParameters_H__
platforms/reference/src/Gbsa/GrycukParameters.cpp 0 → 100644
View file @ f3e77e33
/* 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 <math.h>
#include <sstream>
#include "GrycukParameters.h"
#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
// #define UseGromacsMalloc 1
#ifdef UseGromacsMalloc
extern "C" {
#include "smalloc.h"
}
#endif
const std::string GrycukParameters::ParameterFileName = std::string( "params.agb" );
/**---------------------------------------------------------------------------------------
GrycukParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in Grycuk 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 GrycukParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
GrycukParameters* grycukParameters = new GrycukParameters( numberOfAtoms, log );
grycukParameters->initializeParameters( top );
Gpu:
grycukParameters = new GrycukParameters( 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 GrycukParameters destructor does not free arrays
grycukParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuGrycuk::grycukVdwRadii )->getData() );
grycukParameters->setVolume( getBrookStreamWrapperAtIndex( GpuGrycuk::grycukVolume )->getData() );
grycukParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuGrycuk::grycukGpolFixed )->getData() );
grycukParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuGrycuk::grycukBornRadii )->getData() );
grycukParameters->setFreeArrays( false );
grycukParameters->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()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
GrycukParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
GrycukParameters::GrycukParameters( int numberOfAtoms ) : ImplicitSolventParameters( numberOfAtoms ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGrycukParameters::GrycukParameters";
// ---------------------------------------------------------------------------------------
_ownScaledRadiusFactors = 0;
_scaledRadiusFactors = NULL;
}
/**---------------------------------------------------------------------------------------
GrycukParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
GrycukParameters::~GrycukParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGrycukParameters::~GrycukParameters";
// ---------------------------------------------------------------------------------------
// in GPU runs, arrays may be 'owned' by BrookStreamWrapper -- hence they should not
// be freed here, i.e., _freeArrays should be 'false'
#ifdef UseGromacsMalloc
/*
if( _freeArrays ){
if( _vdwRadii != NULL ){
save_free( "_vdwRadii", __FILE__, __LINE__, _vdwRadii );
}
} */
#else
if( _ownScaledRadiusFactors ){
delete[] _scaledRadiusFactors;
}
/*
if( getFreeArrays() ){
} */
#endif
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atomic radii
--------------------------------------------------------------------------------------- */
RealOpenMM* GrycukParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
// if dielectric offset applied, then unapply
return atomicRadii;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int GrycukParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGrycukParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
return ImplicitSolventParameters::setAtomicRadii( atomicRadii, units );
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int GrycukParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii, int units ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGrycukParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
return ImplicitSolventParameters::setAtomicRadii( atomicRadii, units );
}
/**---------------------------------------------------------------------------------------
Return Grycuk scale factors
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* GrycukParameters::getScaledRadiusFactors( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGrycuk::getScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _scaledRadiusFactors == NULL ){
GrycukParameters* localThis = const_cast<GrycukParameters* const>(this);
localThis->_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownScaledRadiusFactors = true;
memset( _scaledRadiusFactors, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _scaledRadiusFactors;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
@param ownScaledRadiusFactors flag indicating whether scale factors
array should be deleted
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int GrycukParameters::setOwnScaleFactors( int ownScaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGrycuk::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownScaledRadiusFactors = ownScaledRadiusFactors;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set Grycuk scale factors
@param scaledRadiusFactors scaledRadiusFactors
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int GrycukParameters::setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGrycuk::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != scaledRadiusFactors ){
delete[] _scaledRadiusFactors;
_ownScaledRadiusFactors = false;
}
_scaledRadiusFactors = scaledRadiusFactors;
return SimTKOpenMMCommon::DefaultReturn;
}
#if RealOpenMMType == 2
/**---------------------------------------------------------------------------------------
Set Grycuk scale factors
@param scaledRadiusFactors scaledRadiusFactors
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int GrycukParameters::setScaledRadiusFactors( float* scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGrycuk::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _scaledRadiusFactors == NULL ){
_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
_ownScaledRadiusFactors = true;
}
for( int ii = 0; ii < getNumberOfAtoms(); ii++ ){
_scaledRadiusFactors[ii] = (RealOpenMM) scaledRadiusFactors[ii];
}
return SimTKOpenMMCommon::DefaultReturn;
}
#endif
/**---------------------------------------------------------------------------------------
Set Grycuk scale factors
@param scaledRadiusFactors scaledRadiusFactors
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int GrycukParameters::setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuGrycuk::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];
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
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 GrycukParameters::mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const {
// ---------------------------------------------------------------------------------------
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;
}
// map first letter in atom name to Tinker atom number
int firstAsciiValue = ((int) atomName[0]) - 65;
// check for lower case
if( firstAsciiValue > 25 ){
firstAsciiValue -= 32;
}
// validate
if( firstAsciiValue < 0 || firstAsciiValue > 25 ){
if( log != NULL ){
(void) fprintf( log, "Atom name=<%s> unrecognized.", atomName );
}
(void) fprintf( stderr, "Atom name=<%s> unrecognized.", atomName );
return -1;
}
return atomNameMap[firstAsciiValue];
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string GrycukParameters::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nGrycukParameters::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << ImplicitSolventParameters::getStateString( title );
// std::string tab = getStringTab();
return message.str();
}
/**---------------------------------------------------------------------------------------
Return zero value if all parameters set; else return nonzero
@return ready status
--------------------------------------------------------------------------------------- */
int GrycukParameters::isNotReady( void ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nGrycukParameters::isNotReady";
// ---------------------------------------------------------------------------------------
int isReady = ImplicitSolventParameters::isNotReady();
int errors = 0;
std::stringstream message;
message << methodName;
const RealOpenMM* scaledRadiusFactors = getScaledRadiusFactors();
if( scaledRadiusFactors == NULL || scaledRadiusFactors[0] <= 0.0 ){
errors++;
message << "\n scaledRadiusFactors is not set";
}
// check scale factors are in correct units
RealOpenMM average, stdDev, maxValue, minValue;
int minIndex, maxIndex;
SimTKOpenMMUtilities::getArrayStatistics( getNumberOfAtoms(), scaledRadiusFactors, &average,
&stdDev, &minValue, &minIndex,
&maxValue, &maxIndex );
if( average < 0.3 || average > 2.0 || minValue < 0.1 ){
errors++;
message << "\n scale factors for atomic radii appear not to be set correctly -- radii should be in Angstroms";
message << "\n average radius=" << average << " min radius=" << minValue << " at atom index=" << minIndex;
}
if( errors ){
message << std::endl;
SimTKOpenMMLog::printMessage( message );
}
errors += isReady;
return errors;
}
platforms/reference/src/Gbsa/GrycukParameters.h 0 → 100644
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __GrycukParameters_H__
#define __GrycukParameters_H__
#include "SimTKOpenMMCommon.h"
#include "ImplicitSolventParameters.h"
// ---------------------------------------------------------------------------------------
class GrycukParameters : public ImplicitSolventParameters {
public:
static const std::string ParameterFileName;
private:
// scaled radius factors (S_kk in HCT paper)
int _ownScaledRadiusFactors;
RealOpenMM* _scaledRadiusFactors;
public:
/**---------------------------------------------------------------------------------------
GrycukParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
GrycukParameters( int numberOfAtoms );
/**---------------------------------------------------------------------------------------
GrycukParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
~GrycukParameters( );
// override of new/delete
//static void* operator new( size_t size );
//static void operator delete( void *p );
//static void* operator new[]( size_t size );
//static void operator delete[]( void *p );
/**---------------------------------------------------------------------------------------
Return Grycuk scale factors
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* getScaledRadiusFactors( void ) const;
/**---------------------------------------------------------------------------------------
Return Grycuk scale factors
@return array
--------------------------------------------------------------------------------------- */
int setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors );
#if RealOpenMMType == 2
int setScaledRadiusFactors( float* scaledRadiusFactors );
#endif
int 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
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setOwnScaleFactors( int ownScaledRadiusFactors );
/**---------------------------------------------------------------------------------------
Assign standard radii for GB/SA methods other than ACE;
taken from Macromodel and OPLS-AA, except for hydrogens (Simbios)
Logic based on logic in Tinker's ksolv.f
Currently only works for standard amino acid atoms
If invalid atom name is encountered, a message is printed to log file and the
radius for that atom is set to 1.0f
@param numberOfAtoms number of atoms
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int getMacroModelAtomicRadii( int numberOfAtoms,
char*** atomNames, RealOpenMM* radii, FILE* log );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
@return array of atom volumes
--------------------------------------------------------------------------------------- */
RealOpenMM* getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setAtomicRadii( RealOpenMM* atomicRadii, int units = SimTKOpenMMCommon::MdUnits );
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setAtomicRadii( const RealOpenMMVector& atomicRadii, int units = SimTKOpenMMCommon::MdUnits );
/**---------------------------------------------------------------------------------------
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;
/**---------------------------------------------------------------------------------------
Return zero value if all parameters set; else return nonzero
@return ready status
--------------------------------------------------------------------------------------- */
int isNotReady( void ) const;
};
/**---------------------------------------------------------------------------------------
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 // __GrycukParameters_H__
platforms/reference/src/Gbsa/ImplicitSolventParameters.cpp 0 → 100644
View file @ f3e77e33
/* 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 "SimTKOpenMMCommon.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
#include "ImplicitSolventParameters.h"
#include <sstream>
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in ImplicitSolvent equation gPol
(4) list of atoms that should be excluded in calculating
force -- nonbonded atoms (1-2, and 1-3 atoms)
Implementation:
Slightly different sequence of calls when running on CPU vs GPU.
Difference arise because the CPU-side data arrays for the Brook
streams are allocated by the BrookStreamWrapper objects. These
arrays are then used by ImplicitSolventParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
ImplicitSolventParameters* implicitSolventParameters = new ImplicitSolventParameters( numberOfAtoms, log );
implicitSolventParameters->initializeParameters( top );
Gpu:
implicitSolventParameters = new ImplicitSolventParameters( gpu->natoms, log );
// set arrays for cpu using stream data field;
// initializeParameters() only allocates space for arrays if they are not set (==NULL)
// also set flag so that ImplicitSolventParameters destructor does not free arrays
implicitSolventParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventVdwRadii )->getData() );
implicitSolventParameters->setVolume( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventVolume )->getData() );
implicitSolventParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventGpolFixed )->getData() );
implicitSolventParameters->setAtomicRadii( getBrookStreamWrapperAtIndex( GpuImplicitSolvent::implicitSolventAtomicRadii )->getData() );
implicitSolventParameters->setFreeArrays( false );
implicitSolventParameters->initializeParameters( top );
Issues:
Tinker's atom radii are used.
The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
only tested for generic proteins
see mapGmxAtomNameToTinkerAtomNumber()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters::ImplicitSolventParameters( int numberOfAtoms ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::ImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
_numberOfAtoms = numberOfAtoms;
_ownAtomicRadii = false;
_atomicRadii = NULL;
// see comments in ~ImplicitSolventParameters for explanation
_freeArrays = false;
_initializeImplicitSolventConstants();
}
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters::~ImplicitSolventParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::~ImplicitSolventParameters";
// ---------------------------------------------------------------------------------------
if( _ownAtomicRadii ){
delete[] _atomicRadii;
}
}
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::getNumberOfAtoms( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getNumberOfAtoms:";
// ---------------------------------------------------------------------------------------
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getSolventDielectric( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getSolventDielectric:";
// ---------------------------------------------------------------------------------------
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setSolventDielectric( RealOpenMM solventDielectric ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setSolventDielectric:";
// ---------------------------------------------------------------------------------------
_solventDielectric = solventDielectric;
_resetPreFactor();
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getSoluteDielectric( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getSoluteDielectric:";
// ---------------------------------------------------------------------------------------
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setSoluteDielectric( RealOpenMM soluteDielectric ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setSoluteDielectric:";
// ---------------------------------------------------------------------------------------
_soluteDielectric = soluteDielectric;
_resetPreFactor();
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get electric constant (Simbios)
@return electricConstant
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getElectricConstant( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getElectricConstant:";
// ---------------------------------------------------------------------------------------
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Set electric constant (Simbios)
@param electricConstant electric constant
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setElectricConstant( RealOpenMM electricConstant ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setElectricConstant:";
// ---------------------------------------------------------------------------------------
_electricConstant = electricConstant;
_resetPreFactor();
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get probe radius (Simbios)
@return probeRadius
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getProbeRadius( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getProbeRadius:";
// ---------------------------------------------------------------------------------------
return _probeRadius;
}
/**---------------------------------------------------------------------------------------
Set probe radius (Simbios)
@param probeRadius probe radius
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setProbeRadius( RealOpenMM probeRadius ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setProbeRadius:";
// ---------------------------------------------------------------------------------------
_probeRadius = probeRadius;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get pi*4*Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049*1000.0; (Still)
((RealOpenMM) M_PI)*4.0f*0.0054*1000.0; (OBC)
@return pi4Asolv
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getPi4Asolv( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getPi4Asolv:";
// ---------------------------------------------------------------------------------------
return _pi4Asolv;
}
/**---------------------------------------------------------------------------------------
Get prefactor
@return prefactor
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getPreFactor( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getPreFactor:";
// ---------------------------------------------------------------------------------------
return _preFactor;
}
/**---------------------------------------------------------------------------------------
Set pi4Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049*1000.0; (Still)
((RealOpenMM) M_PI)*4.0f*0.0054*1000.0; (OBC)
@param pi4Asolv probe radius
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setPi4Asolv( RealOpenMM pi4Asolv ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setPi4Asolv:";
// ---------------------------------------------------------------------------------------
_pi4Asolv = pi4Asolv;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get conversion factor for kcal/A -> kJ/nm (Simbios)
@return kcalA_To_kJNm factor
--------------------------------------------------------------------------------------- */
RealOpenMM ImplicitSolventParameters::getKcalA_To_kJNm( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getKcalA_To_kJNm:";
// ---------------------------------------------------------------------------------------
return _kcalA_To_kJNm;
}
/**---------------------------------------------------------------------------------------
Set KcalA_To_kJNm
@param kcalA_To_kJNm probe radius
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setKcalA_To_kJNm( RealOpenMM kcalA_To_kJNm ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setKcalA_To_kJNm:";
// ---------------------------------------------------------------------------------------
_kcalA_To_kJNm = kcalA_To_kJNm;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atomic radii
--------------------------------------------------------------------------------------- */
RealOpenMM* ImplicitSolventParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
if( _atomicRadii == NULL ){
ImplicitSolventParameters* localThis = const_cast<ImplicitSolventParameters* const>(this);
localThis->_atomicRadii = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownAtomicRadii = true;
memset( localThis->_atomicRadii, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _atomicRadii;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
if( _ownAtomicRadii && atomicRadii != _atomicRadii ){
delete[] _atomicRadii;
_ownAtomicRadii = false;
}
_atomicRadii = atomicRadii;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii, int units ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nImplicitSolventParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
int numberOfAtoms = getNumberOfAtoms();
if( _ownAtomicRadii ){
delete[] _atomicRadii;
_atomicRadii = new RealOpenMM[numberOfAtoms];
} else if( _atomicRadii == NULL ){
_atomicRadii = new RealOpenMM[numberOfAtoms];
_ownAtomicRadii = true;
}
if( numberOfAtoms != (int) atomicRadii.size() ){
std::stringstream message;
message << methodName;
message << " number of object atoms=" << numberOfAtoms << " does not match number in input vector=" << atomicRadii.size();
SimTKOpenMMLog::printWarning( message );
numberOfAtoms = numberOfAtoms < (int) atomicRadii.size() ? numberOfAtoms : (int) atomicRadii.size();
}
// force kcal/A units
if( units == SimTKOpenMMCommon::MdUnits ){
RealOpenMM ten = (RealOpenMM) 10.0;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
_atomicRadii[ii] = ten*atomicRadii[ii];
}
} else {
for( int ii = 0; ii < numberOfAtoms; ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nImplicitSolventParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
int numberOfAtoms = getNumberOfAtoms();
if( _ownAtomicRadii ){
delete[] _atomicRadii;
_atomicRadii = new RealOpenMM[numberOfAtoms];
} else if( _atomicRadii == NULL ){
_atomicRadii = new RealOpenMM[numberOfAtoms];
_ownAtomicRadii = true;
}
// force kcal/A units
if( units == SimTKOpenMMCommon::MdUnits ){
RealOpenMM ten = (RealOpenMM) 10.0;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
_atomicRadii[ii] = ten*atomicRadii[ii];
}
} else {
for( int ii = 0; ii < numberOfAtoms; ii++ ){
_atomicRadii[ii] = atomicRadii[ii];
}
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setOwnAtomicRadii( int ownAtomicRadii ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nImplicitSolventParameters::setOwnAtomicRadii:";
// ---------------------------------------------------------------------------------------
_ownAtomicRadii = ownAtomicRadii;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get free array flag -- if set then work arrays are freed when destructor is called
@return freeArrays flag
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::getFreeArrays( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setFreeArrays:";
// ---------------------------------------------------------------------------------------
return _freeArrays;
}
/**---------------------------------------------------------------------------------------
Set free array flag -- if set then work arrays are freed when destructor is called
@param freeArrays flag
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::setFreeArrays( int freeArrays ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::setFreeArrays:";
// ---------------------------------------------------------------------------------------
_freeArrays = freeArrays;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Initialize ImplicitSolvent Parameters (Simbios)
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::_initializeImplicitSolventConstants( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::_initializeImplicitSolventConstants:";
// ---------------------------------------------------------------------------------------
_soluteDielectric = (RealOpenMM) 1.0;
_solventDielectric = (RealOpenMM) 78.3;
_kcalA_To_kJNm = (RealOpenMM) 0.4184;
_probeRadius = (RealOpenMM) 1.4;
_electricConstant = (RealOpenMM) -166.02691;
// _pi4Asolv = (RealOpenMM) PI_M*4.0*0.0049*1000.0;
//_pi4Asolv = (RealOpenMM) PI_M*19.6;
// _pi4Asolv = (RealOpenMM) PI_M*4.0*0.0054;
_pi4Asolv = (RealOpenMM) (PI_M*0.0216);
_resetPreFactor();
}
/**---------------------------------------------------------------------------------------
Reset prefactor (Simbios)
called when _electricConstant, _soluteDielectric, or _solventDielectric are modified
--------------------------------------------------------------------------------------- */
void ImplicitSolventParameters::_resetPreFactor( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::_resetPreFactor:";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
if( getSoluteDielectric() != zero && getSolventDielectric() != zero ){
_preFactor = two*getElectricConstant()*( (one/getSoluteDielectric()) - (one/getSolventDielectric()) );
} else {
_preFactor = zero;
}
}
/**---------------------------------------------------------------------------------------
Get state (Simbios)
@param title title (optional)
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
std::string ImplicitSolventParameters::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
if( title ){
message << title;
}
std::string tab = getStringTab();
message << "\nImplicitSolventParameters :";
message << tab << "Solute dielectric: " << getSoluteDielectric();
message << tab << "Solvent dielectric: " << getSolventDielectric();
message << tab << "Electric constant: " << getElectricConstant();
message << tab << "Probe radius: " << getProbeRadius();
message << tab << "Number of atoms: " << getNumberOfAtoms();
message << tab << "Free arrays: " << getFreeArrays();
return message.str();
}
/**---------------------------------------------------------------------------------------
Return zero value if all parameters set; else return nonzero
@return ready status
--------------------------------------------------------------------------------------- */
int ImplicitSolventParameters::isNotReady( void ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nImplicitSolventParameters::isNotReady";
// ---------------------------------------------------------------------------------------
int errors = 0;
int warning = 0;
std::stringstream message;
message << methodName;
if( getNumberOfAtoms() <= 0 ){
errors++;
message << "\n number of atoms=" << getNumberOfAtoms() << " is invalid.";
}
RealOpenMM* atomicRadii = getAtomicRadii();
if( atomicRadii == NULL ){
errors++;
message << "\n scaledRadiusFactors is not set";
}
// check radii are in correct units
RealOpenMM average, stdDev, maxValue, minValue;
int minIndex, maxIndex;
SimTKOpenMMUtilities::getArrayStatistics( getNumberOfAtoms(), atomicRadii, &average,
&stdDev, &minValue, &minIndex,
&maxValue, &maxIndex );
if( average < 1.0 || average > 10.0 || minValue < 0.5 ){
errors++;
message << "\n atomic radii appear not to be set correctly -- radii should be in Angstroms";
message << "\n average radius=" << average << " min radius=" << minValue << " at atom index=" << minIndex;
}
if( getPreFactor() == 0.0 ){
errors++;
message << "\n prefactor is not set.";
}
if( getSolventDielectric() <= 0.0 ){
errors++;
message << "\n solvent dielectric=" << getSolventDielectric() << " is invalid.";
}
if( getSolventDielectric() >= 1000.0 ){
warning++;
message << "\n Warning: solvent dielectric=" << getSolventDielectric() << " is large.";
}
if( getSoluteDielectric() <= 0.0 ){
errors++;
message << "\n solute dielectric=" << getSoluteDielectric() << " is invalid.";
}
if( getSoluteDielectric() >= 1000.0 ){
warning++;
message << "\n Warning: solute dielectric=" << getSoluteDielectric() << " is large.";
}
if( getElectricConstant() >= 0.0 ){
errors++;
message << "\n electric constant=" << getElectricConstant() << " is invalid.";
}
if( getElectricConstant() >= 1000.0 ){
warning++;
message << "\n Warning: electric constant=" << getElectricConstant() << " is large.";
}
if( getProbeRadius() <= 0.0 ){
errors++;
message << "\n probe radius=" << getProbeRadius() << " is invalid.";
}
if( getProbeRadius() >= 10.0 ){
warning++;
message << "\n Warning: probe radius=" << getProbeRadius() << " is large.";
}
if( getPi4Asolv() <= 0.0 ){
errors++;
message << "\n Pi4Asolv=" << getPi4Asolv() << " is invalid.";
}
if( getPi4Asolv() >= 1000.0 ){
warning++;
message << "\n Warning: Pi4Asolv=" << getPi4Asolv() << " is large.";
}
if( errors || warning ){
message << std::endl;
SimTKOpenMMLog::printMessage( message );
}
return errors;
}
/**---------------------------------------------------------------------------------------
Get string tab -- centralized
@return tab string
--------------------------------------------------------------------------------------- */
std::string ImplicitSolventParameters::getStringTab( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getStringTab";
// ---------------------------------------------------------------------------------------
return std::string( "\n " );
}
platforms/reference/src/Gbsa/ImplicitSolventParameters.h 0 → 100755
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ImplicitSolventParameters_H__
#define __ImplicitSolventParameters_H__
#include "SimTKOpenMMRealType.h"
#include <string>
// ---------------------------------------------------------------------------------------
class ImplicitSolventParameters {
protected:
// parameters common to implicit solvent models
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _kcalA_To_kJNm;
RealOpenMM _electricConstant;
RealOpenMM _probeRadius;
RealOpenMM _pi4Asolv;
RealOpenMM _preFactor;
// ---------------------------------------------------------------------------------------
// atom count
int _numberOfAtoms;
// flag signalling whether arrays are to be freed
int _freeArrays;
// atomic radii
int _ownAtomicRadii;
RealOpenMM* _atomicRadii;
/**---------------------------------------------------------------------------------------
Initialize ImplicitSolvent Parameters (Simbios)
--------------------------------------------------------------------------------------- */
void _initializeImplicitSolventConstants( void );
/**---------------------------------------------------------------------------------------
Set KcalA_To_kJNm
@param kcalA_To_kJNm probe radius
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setKcalA_To_kJNm( RealOpenMM kcalA_To_kJNm );
/**---------------------------------------------------------------------------------------
Set free array flag -- if set then work arrays are freed when destructor is called
@param freeArrays flag
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setFreeArrays( int freeArrays );
/**---------------------------------------------------------------------------------------
Reset prefactor (Simbios)
called when _electricConstant, _soluteDielectric, or _solventDielectric are modified
--------------------------------------------------------------------------------------- */
void _resetPreFactor( void );
public:
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
ImplicitSolventParameters( int numberOfAtoms );
/**---------------------------------------------------------------------------------------
ImplicitSolventParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
~ImplicitSolventParameters( );
// override of new/delete
//static void* operator new( size_t size );
//static void operator delete( void *p );
//static void* operator new[]( size_t size );
//static void operator delete[]( void *p );
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int getNumberOfAtoms( void ) const;
/**---------------------------------------------------------------------------------------
Get free array flag -- if set then work arrays are freed when destructor is called
@return freeArrays flag
--------------------------------------------------------------------------------------- */
int getFreeArrays( void ) const;
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSolventDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setSolventDielectric( RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSoluteDielectric( void ) const;
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setSoluteDielectric( RealOpenMM soluteDielectric );
/**---------------------------------------------------------------------------------------
Get conversion factor for kcal/A -> kJ/nm (Simbios)
@return kcalA_To_kJNm factor
--------------------------------------------------------------------------------------- */
RealOpenMM getKcalA_To_kJNm( void ) const;
/**---------------------------------------------------------------------------------------
Get electric constant (Simbios)
@return electricConstant
--------------------------------------------------------------------------------------- */
RealOpenMM getElectricConstant( void ) const;
/**---------------------------------------------------------------------------------------
Set electric constant (Simbios)
@param electricConstant electric constant
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setElectricConstant( RealOpenMM electricConstant );
/**---------------------------------------------------------------------------------------
Get probe radius (Simbios)
@return probeRadius
--------------------------------------------------------------------------------------- */
RealOpenMM getProbeRadius( void ) const;
/**---------------------------------------------------------------------------------------
Set probe radius (Simbios)
@param probeRadius probe radius
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setProbeRadius( RealOpenMM probeRadius );
/**---------------------------------------------------------------------------------------
Get pi4Asolv: used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049f*1000.0f; (Simbios)
@return pi4Asolv
--------------------------------------------------------------------------------------- */
RealOpenMM getPi4Asolv( void ) const;
/**---------------------------------------------------------------------------------------
Get prefactor
@returnpreFactor
--------------------------------------------------------------------------------------- */
RealOpenMM getPreFactor( void ) const;
/**---------------------------------------------------------------------------------------
Set values used in ACE approximation for nonpolar term
((RealOpenMM) M_PI)*4.0f*0.0049f*1000.0f; (Simbios)
@param pi4Asolv see above
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setPi4Asolv( RealOpenMM pi4Asolv );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atom volumes
--------------------------------------------------------------------------------------- */
virtual RealOpenMM* getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
virtual int setAtomicRadii( RealOpenMM* atomicRadii );
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
@param units units flag SimTKOpenMMCommon::MdUnits or
SimTKOpenMMCommon::KcalAngUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
virtual int setAtomicRadii( const RealOpenMMVector& atomicRadii,
int units = SimTKOpenMMCommon::MdUnits );
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
virtual int setAtomicRadii( RealOpenMM* atomicRadii, int units );
/**---------------------------------------------------------------------------------------
Set flag indicating whether AtomicRadii array is owned by
object; if flag is set, then when the object is deleted,
the array will be freed
@param ownAtomicRadii flag indicating whether array of atomic radii
should be freed
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setOwnAtomicRadii( int ownAtomicRadii );
/**---------------------------------------------------------------------------------------
Print state to log file (Simbios)
@param title title (optional)
@param log print state to log file
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
virtual std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Get string tab -- centralized
@return tab string
--------------------------------------------------------------------------------------- */
std::string getStringTab( void ) const;
/**---------------------------------------------------------------------------------------
Return nonzero value errors
@return ready status
--------------------------------------------------------------------------------------- */
virtual int isNotReady( void ) const;
};
// ---------------------------------------------------------------------------------------
#endif // __ImplicitSolventParameters_H__
platforms/reference/src/Gbsa/Makefile 0 → 100644
View file @ f3e77e33
# ---------------------------------------------------------------------------------------
include ../makeConfig/Top.mk
# ---------------------------------------------------------------------------------------
# target is gmxGbsa library
# ld reference is in ../gromacs/vc7.mk
# include file reference is in ../gromacs/gmxgpu.h
STATIC_LIBRARY := gmxGbsa
# ---------------------------------------------------------------------------------------
BaseFiles := $(GbsaBaseFiles)
IncludeDirectories := $(SimTkDirectory) $(GromacsIncludeDirectory)
# ---------------------------------------------------------------------------------------
include $(MakeDirectory)/BaseMake.mk
# ---------------------------------------------------------------------------------------
platforms/reference/src/Gbsa/ObcParameters.cpp 0 → 100644
View file @ f3e77e33
/* 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 <math.h>
#include <sstream>
#include "ObcParameters.h"
#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
// #define UseGromacsMalloc 1
#ifdef UseGromacsMalloc
extern "C" {
#include "smalloc.h"
}
#endif
const std::string ObcParameters::ParameterFileName = std::string( "params.agb" );
/**---------------------------------------------------------------------------------------
ObcParameters:
Calculates for each atom
(1) the van der Waal radii
(2) volume
(3) fixed terms in Obc equation gPol
(4) list of atoms that should be excluded in calculating
force -- nonbonded atoms (1-2, and 1-3 atoms)
Implementation:
Slightly different sequence of calls when running on CPU vs GPU.
Difference arise because the CPU-side data arrays for the Brook
streams are allocated by the BrookStreamWrapper objects. These
arrays are then used by ObcParameters when initializing the
the values (vdwRadii, volume, ...) to be used in the calculation.
Cpu:
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, log );
obcParameters->initializeParameters( top );
Gpu:
obcParameters = new ObcParameters( gpu->natoms, log );
// set arrays for cpu using stream data field;
// initializeParameters() only allocates space for arrays if they are not set (==NULL)
// also set flag so that ObcParameters destructor does not free arrays
obcParameters->setVdwRadii( getBrookStreamWrapperAtIndex( GpuObc::obcVdwRadii )->getData() );
obcParameters->setVolume( getBrookStreamWrapperAtIndex( GpuObc::obcVolume )->getData() );
obcParameters->setGPolFixed( getBrookStreamWrapperAtIndex( GpuObc::obcGpolFixed )->getData() );
obcParameters->setBornRadii( getBrookStreamWrapperAtIndex( GpuObc::obcBornRadii )->getData() );
obcParameters->setFreeArrays( false );
obcParameters->initializeParameters( top );
Issues:
Tinker's atom radii are used.
The logic for mapping the Gromacs atom names to Tinker type may be incomplete;
only tested for generic proteins
see mapGmxAtomNameToTinkerAtomNumber()
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
ObcParameters constructor (Simbios)
@param numberOfAtoms number of atoms
@param obcType OBC type (Eq. 7 or 8 in paper)
--------------------------------------------------------------------------------------- */
ObcParameters::ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType ) : ImplicitSolventParameters( numberOfAtoms ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::ObcParameters";
// ---------------------------------------------------------------------------------------
_obcType = obcType;
_dielectricOffset = 0.09f;
_ownScaledRadiusFactors = 0;
_scaledRadiusFactors = NULL;
setObcTypeParameters( obcType );
}
/**---------------------------------------------------------------------------------------
ObcParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
ObcParameters::~ObcParameters( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::~ObcParameters";
// ---------------------------------------------------------------------------------------
// in GPU runs, arrays may be 'owned' by BrookStreamWrapper -- hence they should not
// be freed here, i.e., _freeArrays should be 'false'
#ifdef UseGromacsMalloc
/*
if( _freeArrays ){
if( _vdwRadii != NULL ){
save_free( "_vdwRadii", __FILE__, __LINE__, _vdwRadii );
}
} */
#else
if( _ownScaledRadiusFactors ){
delete[] _scaledRadiusFactors;
}
/*
if( getFreeArrays() ){
} */
#endif
}
/**---------------------------------------------------------------------------------------
Get OBC type
@return OBC type
--------------------------------------------------------------------------------------- */
ObcParameters::ObcType ObcParameters::getObcType( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getObcType:";
// ---------------------------------------------------------------------------------------
return _obcType;
}
/**---------------------------------------------------------------------------------------
Set OBC type specific parameters
@param obcType OBC type (ObcTypeI or ObcTypeII -- Eq. 7 or 8)
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ObcParameters::setObcTypeParameters( ObcParameters::ObcType obcType ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::setObcTypeParameters:";
// ---------------------------------------------------------------------------------------
if( obcType == ObcTypeI ){
_alphaObc = 0.8f;
_betaObc = 0.0f;
_gammaObc = 2.91f;
} else {
_alphaObc = 1.0f;
_betaObc = 0.8f;
_gammaObc = 4.85f;
}
_obcType = obcType;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get dielectricOffset
@return _dielectricOffset
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getDielectricOffset( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getDielectricOffset:";
// ---------------------------------------------------------------------------------------
return _dielectricOffset;
}
/**---------------------------------------------------------------------------------------
Get alpha OBC (Eqs. 6 & 7) in Proteins paper
@return alphaObc
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getAlphaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getAlphaObc:";
// ---------------------------------------------------------------------------------------
return _alphaObc;
}
/**---------------------------------------------------------------------------------------
Get beta OBC (Eqs. 6 & 7) in Proteins paper
@return betaObc
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getBetaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getBetaObc:";
// ---------------------------------------------------------------------------------------
return _betaObc;
}
/**---------------------------------------------------------------------------------------
Get gamma OBC (Eqs. 6 & 7) in Proteins paper
@return gammaObc
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getGammaObc( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getGammaObc:";
// ---------------------------------------------------------------------------------------
return _gammaObc;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atomic radii
--------------------------------------------------------------------------------------- */
RealOpenMM* ObcParameters::getAtomicRadii( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nImplicitSolventParameters::getAtomicRadii:";
// ---------------------------------------------------------------------------------------
RealOpenMM* atomicRadii = ImplicitSolventParameters::getAtomicRadii();
// if dielectric offset applied, then unapply
return atomicRadii;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ObcParameters::setAtomicRadii( RealOpenMM* atomicRadii, int units ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
return ImplicitSolventParameters::setAtomicRadii( atomicRadii, units );
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ObcParameters::setAtomicRadii( const RealOpenMMVector& atomicRadii, int units ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nObcParameters::setAtomicRadii:";
// ---------------------------------------------------------------------------------------
return ImplicitSolventParameters::setAtomicRadii( atomicRadii, units );
}
/**---------------------------------------------------------------------------------------
Return OBC scale factors
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcParameters::getScaledRadiusFactors( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::getScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _scaledRadiusFactors == NULL ){
ObcParameters* localThis = const_cast<ObcParameters* const>(this);
localThis->_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
localThis->_ownScaledRadiusFactors = true;
memset( _scaledRadiusFactors, 0, sizeof( RealOpenMM )*getNumberOfAtoms() );
}
return _scaledRadiusFactors;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether scale factors array should be deleted
@param ownScaledRadiusFactors flag indicating whether scale factors
array should be deleted
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int ObcParameters::setOwnScaleFactors( int ownScaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setOwnScaleFactors";
// ---------------------------------------------------------------------------------------
_ownScaledRadiusFactors = ownScaledRadiusFactors;
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Set OBC scale factors
@param scaledRadiusFactors scaledRadiusFactors
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int ObcParameters::setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != scaledRadiusFactors ){
delete[] _scaledRadiusFactors;
_ownScaledRadiusFactors = false;
}
_scaledRadiusFactors = scaledRadiusFactors;
return SimTKOpenMMCommon::DefaultReturn;
}
#if RealOpenMMType == 2
/**---------------------------------------------------------------------------------------
Set OBC scale factors
@param scaledRadiusFactors scaledRadiusFactors
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int ObcParameters::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];
}
return SimTKOpenMMCommon::DefaultReturn;
}
#endif
/**---------------------------------------------------------------------------------------
Set OBC scale factors
@param scaledRadiusFactors scaledRadiusFactors
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int ObcParameters::setScaledRadiusFactors( const RealOpenMMVector& scaledRadiusFactors ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::setScaledRadiusFactors";
// ---------------------------------------------------------------------------------------
if( _ownScaledRadiusFactors && _scaledRadiusFactors != NULL ){
delete[] _scaledRadiusFactors;
}
_ownScaledRadiusFactors = true;
_scaledRadiusFactors = new RealOpenMM[getNumberOfAtoms()];
for( int ii = 0; ii < (int) scaledRadiusFactors.size(); ii++ ){
_scaledRadiusFactors[ii] = scaledRadiusFactors[ii];
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
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 ObcParameters::mapGmxAtomNameToTinkerAtomNumber( const char* atomName, FILE* log ) const {
// ---------------------------------------------------------------------------------------
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;
}
// map first letter in atom name to Tinker atom number
int firstAsciiValue = ((int) atomName[0]) - 65;
// check for lower case
if( firstAsciiValue > 25 ){
firstAsciiValue -= 32;
}
// validate
if( firstAsciiValue < 0 || firstAsciiValue > 25 ){
if( log != NULL ){
(void) fprintf( log, "Atom name=<%s> unrecognized.", atomName );
}
(void) fprintf( stderr, "Atom name=<%s> unrecognized.", atomName );
return -1;
}
return atomNameMap[firstAsciiValue];
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string
--------------------------------------------------------------------------------------- */
std::string ObcParameters::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nObcParameters::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << ImplicitSolventParameters::getStateString( title );
std::string tab = getStringTab();
if( getObcType() == ObcTypeI ){
message << tab << "OBC type: Type I";
} else {
message << tab << "OBC type: Type II";
}
message << tab << "Alpha: " << getAlphaObc();
message << tab << "Beta: " << getBetaObc();
message << tab << "Gamma: " << getGammaObc();
return message.str();
}
/**---------------------------------------------------------------------------------------
Return zero value if all parameters set; else return nonzero
@return ready status
--------------------------------------------------------------------------------------- */
int ObcParameters::isNotReady( void ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nObcParameters::isNotReady";
// ---------------------------------------------------------------------------------------
int isReady = ImplicitSolventParameters::isNotReady();
int errors = 0;
std::stringstream message;
message << methodName;
const RealOpenMM* scaledRadiusFactors = getScaledRadiusFactors();
if( scaledRadiusFactors == NULL || scaledRadiusFactors[0] <= 0.0 ){
errors++;
message << "\n scaledRadiusFactors is not set";
}
// check scale factors are in correct units
RealOpenMM average, stdDev, maxValue, minValue;
int minIndex, maxIndex;
SimTKOpenMMUtilities::getArrayStatistics( getNumberOfAtoms(), scaledRadiusFactors, &average,
&stdDev, &minValue, &minIndex,
&maxValue, &maxIndex );
if( average < 0.3 || average > 2.0 || minValue < 0.1 ){
errors++;
message << "\n scale factors for atomic radii appear not to be set correctly -- radii should be in Angstroms";
message << "\n average radius=" << average << " min radius=" << minValue << " at atom index=" << minIndex;
}
if( errors ){
message << std::endl;
SimTKOpenMMLog::printMessage( message );
}
errors += isReady;
return errors;
}
platforms/reference/src/Gbsa/ObcParameters.h 0 → 100755
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ObcParameters_H__
#define __ObcParameters_H__
#include "SimTKOpenMMCommon.h"
#include "ImplicitSolventParameters.h"
// ---------------------------------------------------------------------------------------
class ObcParameters : public ImplicitSolventParameters {
public:
// OBC types
enum ObcType { ObcTypeI, ObcTypeII };
static const std::string ParameterFileName;
private:
// OBC constants & parameters
RealOpenMM _dielectricOffset;
RealOpenMM _alphaObc;
RealOpenMM _betaObc;
RealOpenMM _gammaObc;
ObcType _obcType;
// scaled radius factors (S_kk in HCT paper)
int _ownScaledRadiusFactors;
RealOpenMM* _scaledRadiusFactors;
/**---------------------------------------------------------------------------------------
Set solvent dielectric (Simbios)
@param dielectricOffset solvent dielectric
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setDielectricOffset( RealOpenMM dielectricOffset );
public:
/**---------------------------------------------------------------------------------------
ObcParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType = ObcTypeII );
/**---------------------------------------------------------------------------------------
ObcParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
~ObcParameters( );
// override of new/delete
//static void* operator new( size_t size );
//static void operator delete( void *p );
//static void* operator new[]( size_t size );
//static void operator delete[]( void *p );
/**---------------------------------------------------------------------------------------
Get OBC type
@return OBC type
--------------------------------------------------------------------------------------- */
ObcParameters::ObcType getObcType( void ) const;
/**---------------------------------------------------------------------------------------
Set OBC type specific parameters
@param obcType OBC type (ObcTypeI or ObcTypeII -- Eq. 7 or 8)
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setObcTypeParameters( ObcParameters::ObcType obcType );
/**---------------------------------------------------------------------------------------
Get alpha OBC (Eqs. 6 & 7) in Proteins paper
@return alphaObc
--------------------------------------------------------------------------------------- */
RealOpenMM getAlphaObc( void ) const;
/**---------------------------------------------------------------------------------------
Get beta OBC (Eqs. 6 & 7) in Proteins paper
@return betaObc
--------------------------------------------------------------------------------------- */
RealOpenMM getBetaObc( void ) const;
/**---------------------------------------------------------------------------------------
Get gamma OBC (Eqs. 6 & 7) in Proteins paper
@return gammaObc
--------------------------------------------------------------------------------------- */
RealOpenMM getGammaObc( void ) const;
/**---------------------------------------------------------------------------------------
Get solvent dielectric (Simbios)
@return dielectricOffset dielectric offset
--------------------------------------------------------------------------------------- */
RealOpenMM getDielectricOffset( void ) const;
/**---------------------------------------------------------------------------------------
Return OBC scale factors
@return array
--------------------------------------------------------------------------------------- */
const RealOpenMM* getScaledRadiusFactors( void ) const;
/**---------------------------------------------------------------------------------------
Set OBC scale factors
@return array
--------------------------------------------------------------------------------------- */
int setScaledRadiusFactors( RealOpenMM* scaledRadiusFactors );
#if RealOpenMMType == 2
int setScaledRadiusFactors( float* scaledRadiusFactors );
#endif
int 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
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setOwnScaleFactors( int ownScaledRadiusFactors );
/**---------------------------------------------------------------------------------------
Assign standard radii for GB/SA methods other than ACE;
taken from Macromodel and OPLS-AA, except for hydrogens (Simbios)
Logic based on logic in Tinker's ksolv.f
Currently only works for standard amino acid atoms
If invalid atom name is encountered, a message is printed to log file and the
radius for that atom is set to 1.0f
@param numberOfAtoms number of atoms
@param atomNames array of atom names from GMX top data struct
@param radii array to store Macromodel radii for each atom
@param log if set, then print error messages to log file
@return SimTKOpenMMCommon::DefaultReturn always
--------------------------------------------------------------------------------------- */
int getMacroModelAtomicRadii( int numberOfAtoms,
char*** atomNames, RealOpenMM* radii, FILE* log );
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
@return array of atom volumes
--------------------------------------------------------------------------------------- */
RealOpenMM* getAtomicRadii( void ) const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii array of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setAtomicRadii( RealOpenMM* atomicRadii, int units = SimTKOpenMMCommon::MdUnits );
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
@param units units flag: SimTKOpenMMCommon::KcalAngUnits or
SimTKOpenMMCommon::MdUnits
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
int setAtomicRadii( const RealOpenMMVector& atomicRadii, int units = SimTKOpenMMCommon::MdUnits );
/**---------------------------------------------------------------------------------------
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;
/**---------------------------------------------------------------------------------------
Return zero value if all parameters set; else return nonzero
@return ready status
--------------------------------------------------------------------------------------- */
int isNotReady( void ) const;
};
/**---------------------------------------------------------------------------------------
Qsort/heapsort integer comparison (Simbios)
@parma a first value to compare
@param b second value to compare
@return -1, 0, 1
--------------------------------------------------------------------------------------- */
int integerComparison( const void *a, const void *b);
// ---------------------------------------------------------------------------------------
#endif // __ObcParameters_H__
platforms/reference/src/Gbsa/ObcReadme.txt 0 → 100644
View file @ f3e77e33
The OBC Generalized Born model is based on the following papers:
J. Phys. Chem. 1996 100, 19824-19839 (HCT paper)
Proteins: Structure, Function, and Bioinformatics 55:383-394 (2004) (OBC paper)
---------------------------------------------------------------------------------------
The two main interface methods perform the following tasks:
(1) set the GBSA parameters; this method is called once prior to execution of
the main simulation loop
(2) calculate the GBSA forces/energy given the atom coordinates and charges
for each simulation step
The two methods and helper methods are in the file cpuObcInterface.cpp
---------------------------------------------------------------------------------------
The setup call is
cpuSetObcParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* obcScaleFactors,
int includeAceApproximation,
RealOpenMM soluteDielectric, RealOpenMM solventDielectric, FILE* log );
The parameters include the following:
(1) number of atoms
(2) the OBC atomic radii (akin to VDW radii, but optimized for GBSA calculations)
(3) the OBC scale factors for each atom
(4) a flag indicating whether the nonpolar ACE approximation is to be included
in calculating the forces and energy: I do not have a reference for this term
(5-6) the solute and solvent dielectric (typically 1.0 & 78.3)
(7) a log file reference; the reference may be set to NULL,
if logging is unwanted. Note: if the reference is NULL, warnings, ... will be output to stderr
cpuSetObcParameters() creates a static CpuObc object which is then referenced when the
GBSA forces and energy are to be calculated
The default OBC model is the Type II model: Eq. 8 in the 2004 OBC paper;
for the Type I model, replace the call
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
with
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeI );
---------------------------------------------------------------------------------------
The OBC scale factors may be obtained via calls to either of the following 'helper' methods
getObcScaleFactors( int numberOfAtoms, const int* atomicNumber, RealOpenMM* scaleFactors )
or
getObcScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses, RealOpenMM* scaleFactors )
Documentation for the inputs to these routines is contained in the file, cpuObcInterface.cpp.
The scale factors are those used in Tinker5, and are based on the element type
of each atom (H, C, N, O, ...).
---------------------------------------------------------------------------------------
The OBC radii may be obtained via the helper method:
getGbsaRadii( int numberOfAtoms, const int* atomicNumber,
const int* numberOfCovalentPartners, const int* indexOfCovalentPartner,
RealOpenMM* gbsaRadii );
Documentation for the inputs to this routine is contained in the file, cpuObcInterface.cpp.
The radii are based on the atom type (methyl C, carbonyl C, ...). The value for a hydrogen
depends on its covalent heavy atom partner. As for the OBC scale factors, the values
hardwired into the method are the same as those used in Tinker5.
---------------------------------------------------------------------------------------
The energy/forces are calculated via the call
cpuCalculateImplicitSolventForces( RealOpenMM** atomCoordinates, const RealOpenMM* partialCharges,
RealOpenMM** forces, RealOpenMM* energy, int updateBornRadii );
The coordinate and force arrays should have the layout: atomCoordinates[atomIndex][3] &
forces[atomIndex][3]
If updateBornRadii = 0, then the Born radii from the previous iteration are used; this reduces the
number of O(N**2) loops from 3 to 2.
If updateBornRadii != 0, then the Born radii are calculated for the input configuration.
For MD runs, setting updateBornRadii = 0 is relatively safe since the radii change little
from one iteration to the next. However, for minimization techniques updateBornRadii
should probably be set to a nonzero value since the configurations can change significantly
between adjacent calls.
---------------------------------------------------------------------------------------
Notes:
(1) The type (float or double) for floating-point numbers is set in the file SimTKOpenMMRealType.h.
If 'RealOpenMMType' is set to 1 in the file, the calculations are performed in
single-precision; otherwise they are performed in double-precision
(2) The units are Angstroms for the spatial dimensions, kcal/mol for energy, and kcal/mol.A for the forces
(3) The code in gromacsCpuObcInterface.cpp provides an example of the parameter setup, ... used for
interfacing w/ Gromacs; the details are somewhat different from the code in cpuObcInterface.cpp
(4) Output from Tinker for a singles ubiquitin configuration are in the directory 'Examples'.
The xyz file gives the coordinates and atom types (Amber99 -- note the charge for atom type 268
(line 2008 in the Tinker amber99.prm file) was changed from -0.2737 to -0.2774 to agree w/
the value used in the version of Gromacs we have access to.)
The key file specifies that OBC should be used for the implicit solvent.
The results are in the files 111UBQ.gbsaAce and 111UBQ.gbsaNoAce. The first file includes the
ACE approximation, while the second does not. On the first line of each file, the number of atoms
and GBSA energy are reported. The remaining lines (one for each atom) give the coordinates (3 entries),
Born radius, partial charge, GBSA radius (rsolv in Tinker), the OBC scaling factor, and the forces (3). The
last three fields give the atomic number, Tinker Amber99 atom name and type.
---------------------------------------------------------------------------------------
platforms/reference/src/Gbsa/asciiTable 0 → 100644
View file @ f3e77e33
/* Ascii table for reference
Binary Dec Hex Glyph
0010 0000 32 20 SP
0010 0001 33 21 !
0010 0010 34 22 "
0010 0011 35 23 #
0010 0100 36 24 $
0010 0101 37 25 %
0010 0110 38 26 &
0010 0111 39 27 '
0010 1000 40 28 (
0010 1001 41 29 )
0010 1010 42 2A *
0010 1011 43 2B +
0010 1100 44 2C ,
0010 1101 45 2D -
0010 1110 46 2E .
0010 1111 47 2F /
0011 0000 48 30 0
0011 0001 49 31 1
0011 0010 50 32 2
0011 0011 51 33 3
0011 0100 52 34 4
0011 0101 53 35 5
0011 0110 54 36 6
0011 0111 55 37 7
0011 1000 56 38 8
0011 1001 57 39 9
0011 1010 58 3A :
0011 1011 59 3B ;
0011 1100 60 3C <
0011 1101 61 3D =
0011 1110 62 3E >
0011 1111 63 3F ?
0100 0000 64 40 @
0100 0001 65 41 A
0100 0010 66 42 B
0100 0011 67 43 C
0100 0100 68 44 D
0100 0101 69 45 E
0100 0110 70 46 F
0100 0111 71 47 G
0100 1000 72 48 H
0100 1001 73 49 I
0100 1010 74 4A J
0100 1011 75 4B K
0100 1100 76 4C L
0100 1101 77 4D M
0100 1110 78 4E N
0100 1111 79 4F O
0101 0000 80 50 P
0101 0001 81 51 Q
0101 0010 82 52 R
0101 0011 83 53 S
0101 0100 84 54 T
0101 0101 85 55 U
0101 0110 86 56 V
0101 0111 87 57 W
0101 1000 88 58 X
0101 1001 89 59 Y
0101 1010 90 5A Z
0101 1011 91 5B [
0101 1100 92 5C \
0101 1101 93 5D ]
0101 1110 94 5E ^
0101 1111 95 5F _
0110 0000 96 60 x
0110 0001 97 61 a
0110 0010 98 62 b
0110 0011 99 63 c
0110 0100 100 64 d
0110 0101 101 65 e
0110 0110 102 66 f
0110 0111 103 67 g
0110 1000 104 68 h
0110 1001 105 69 i
0110 1010 106 6A j
0110 1011 107 6B k
0110 1100 108 6C l
0110 1101 109 6D m
0110 1110 110 6E n
0110 1111 111 6F o
0111 0000 112 70 p
0111 0001 113 71 q
0111 0010 114 72 r
0111 0011 115 73 s
0111 0100 116 74 t
0111 0101 117 75 u
0111 0110 118 76 v
0111 0111 119 77 w
0111 1000 120 78 x
0111 1001 121 79 y
0111 1010 122 7A z
0111 1011 123 7B {
0111 1100 124 7C |
0111 1101 125 7D }
0111 1110 126 7E ~
*/
platforms/reference/src/Gbsa/checklist 0 → 100644
View file @ f3e77e33
In gbsaCpu:
cpuSetGbsaParameters:
cpuGbsaLog setting
ACE approximation
force/energy conversion factors
cpuCalculateGbsaForces
cpuGbsaLog setting
In CpuGbsa::computeGbsaForces
printSampleOutput = 10
platforms/reference/src/Gbsa/cpuGrycukInterface.cpp 0 → 100644
View file @ f3e77e33
/* 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 "cpuGrycukInterface.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
// #include "SimTKOpenMMGromacsUtilities.h"
#include "GrycukParameters.h"
#include "CpuGrycuk.h"
/**---------------------------------------------------------------------------------------
Setup for Grycuk calculations from Gromacs
@param numberOfAtoms number of atoms
@param grycukScaleFactors array of Grycuk scale factors (one entry each atom)
@param atomicRadii atomic radii in Angstrom (one entry each atom)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculations
@param soluteDielectric solute dielectric
@param solventDielectric solvent dielectric
@param log log reference -- if NULL, then errors/warnings
output to stderr
The method creates a CpuGrycuk instance -- currently the Grycuk type II model is the
default (see paper). If the Grycuk type I model is desired change
GrycukParameters* grycukParameters = new GrycukParameters( numberOfAtoms, GrycukParameters::GrycukTypeII );
to
GrycukParameters* grycukParameters = new GrycukParameters( numberOfAtoms, GrycukParameters::GrycukTypeI );
The created object is a static member of the class CpuGrycuk;
when the force routine, cpuCalculateGrycukForces(), is called,
the static object is used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
extern "C" int
cpuSetGrycukParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* grycukScaleFactors,
int includeAceApproximation,
RealOpenMM soluteDielectric, RealOpenMM solventDielectric, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\ncpuSetGrycukParameters: ";
// ---------------------------------------------------------------------------------------
// set log file if not NULL
if( log ){
SimTKOpenMMLog::setSimTKOpenMMLog( log );
}
// set Grycuk parameters
GrycukParameters* grycukParameters = new GrycukParameters( numberOfAtoms );
grycukParameters->setScaledRadiusFactors( grycukScaleFactors );
grycukParameters->setAtomicRadii( atomicRadii, SimTKOpenMMCommon::KcalAngUnits );
// dielectric constants
grycukParameters->setSolventDielectric( solventDielectric );
grycukParameters->setSoluteDielectric( soluteDielectric );
// ---------------------------------------------------------------------------------------
// create CpuGrycuk instance that will calculate forces
CpuGrycuk* cpuGrycuk = new CpuGrycuk( grycukParameters );
// set static member for subsequent calls to calculate forces/energy
CpuImplicitSolvent::setCpuImplicitSolvent( cpuGrycuk );
// set base file name, ...
//cpuGrycuk->readInfoFile( "CpuImplicitSolventInfo" );
// include/do not include ACE approximation (nonpolar solvation)
cpuGrycuk->setIncludeAceApproximation( includeAceApproximation );
// ---------------------------------------------------------------------------------------
// diagnostics
if( log ){
std::string state = cpuGrycuk->getStateString( methodName );
(void) fprintf( log, "\n%s\nDone w/ setup\n", state.c_str() );
(void) fflush( log );
}
// ---------------------------------------------------------------------------------------
return 0;
}
/**---------------------------------------------------------------------------------------
Get Grycuk scale factors given masses
@param numberOfAtoms number of atoms
@param masses input masses
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
extern "C" int getGrycukScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses,
RealOpenMM* scaleFactors ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetGrycukScaleFactorsGivenAtomMasses";
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
double scaleFactor = 0.8;
RealOpenMM mass = masses[atomI];
if ( mass < 1.2 && mass >= 1.0 ){ // hydrogen
scaleFactor = 0.85;
} else if( mass > 11.8 && mass < 12.2 ){ // carbon
scaleFactor = 0.72;
} else if( mass > 14.0 && mass < 15.0 ){ // nitrogen
scaleFactor = 0.79;
} else if( mass > 15.5 && mass < 16.5 ){ // oxygen
scaleFactor = 0.85;
} else if( mass > 31.5 && mass < 32.5 ){ // sulphur
scaleFactor = 0.96;
} else if( mass > 29.5 && mass < 30.5 ){ // phosphorus
scaleFactor = 0.86;
} else {
std::stringstream message;
message << methodName;
message << " Warning: mass for atom " << atomI << " mass=" << mass << "> not recognized.";
SimTKOpenMMLog::printMessage( message );
}
scaleFactors[atomI] = (RealOpenMM) scaleFactor;
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get Grycuk scale factors given atomic numbers
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
extern "C" int getGrycukScaleFactors( int numberOfAtoms, const int* atomicNumber, RealOpenMM* scaleFactors ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetGrycukScaleFactors";
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
double scaleFactor;
switch( atomicNumber[atomI] ){
case 1: // hydrogen
scaleFactor = 0.85;
break;
case 6: // carbon
scaleFactor = 0.72;
break;
case 7: // nitrogen
scaleFactor = 0.79;
break;
case 8: // oxygen
scaleFactor = 0.85;
break;
case 15: // phosphorus
scaleFactor = 0.86;
break;
case 16: // sulphur
scaleFactor = 0.85;
break;
default:
scaleFactor = 0.8;
std::stringstream message;
message << methodName;
message << " Warning: atom number=" << atomicNumber[atomI] << " for atom " << atomI << "> not handled -- useing default value.";
SimTKOpenMMLog::printMessage( message );
break;
}
scaleFactors[atomI] = (RealOpenMM) scaleFactor;
}
return SimTKOpenMMCommon::DefaultReturn;
}
platforms/reference/src/Gbsa/cpuGrycukInterface.h 0 → 100644
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __CpuGrycukInterface_H__
#define __CpuGrycukInterface_H__
#ifdef __cplusplus
#define externC extern "C"
#else
#define externC extern
#endif
#include "SimTKOpenMMRealType.h"
#include <stdio.h>
/**---------------------------------------------------------------------------------------
Delete the Grycuk associated object(s)
@return 0 if static CpuGrycuk object was set; else return -1
--------------------------------------------------------------------------------------- */
externC int cpuDeleteGrycukParameters( void );
/**---------------------------------------------------------------------------------------
Setup for Grycuk calculations from Gromacs
@param numberOfAtoms number of atoms
@param grycukScaleFactors array of Grycuk scale factors (one entry each atom)
@param atomicRadii atomic radii in Angstrom (one entry each atom)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculations
@param soluteDielectric solute dielectric
@param solventDielectric solvent dielectric
@param log log reference -- if NULL, then errors/warnings
output to stderr
The method creates a CpuGrycuk instance -- currently the Grycuk type II model is the
default (see paper). If the Grycuk type I model is desired change
GrycukParameters* grycukParameters = new GrycukParameters( numberOfAtoms, GrycukParameters::GrycukTypeII );
to
GrycukParameters* grycukParameters = new GrycukParameters( numberOfAtoms, GrycukParameters::GrycukTypeI );
The created object is a static member of the class CpuGrycuk;
when the force routine, cpuCalculateGrycukForces(), is called,
the static object is used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
externC
int cpuSetGrycukParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* grycukScaleFactors,
int includeAceApproximation, RealOpenMM soluteDielectric,
RealOpenMM solventDielectric, FILE* log );
/**---------------------------------------------------------------------------------------
Get Grycuk scale factors given masses
@param numberOfAtoms number of atoms
@param masses input masses
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
externC int getGrycukScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses,
RealOpenMM* scaleFactors );
/**---------------------------------------------------------------------------------------
Get Grycuk scale factors given atomic numbers
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
externC int getGrycukScaleFactors( int numberOfAtoms, const int* atomicNumber,
RealOpenMM* scaleFactors );
#undef externC
#endif
platforms/reference/src/Gbsa/cpuObcInterface.cpp 0 → 100644
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "cpuObcInterface.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
// #include "SimTKOpenMMGromacsUtilities.h"
#include "ObcParameters.h"
#include "CpuObc.h"
/**---------------------------------------------------------------------------------------
Setup for Obc calculations from Gromacs
@param numberOfAtoms number of atoms
@param obcScaleFactors array of OBC scale factors (one entry each atom)
@param atomicRadii atomic radii in Angstrom (one entry each atom)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculations
@param soluteDielectric solute dielectric
@param solventDielectric solvent dielectric
@param log log reference -- if NULL, then errors/warnings
output to stderr
The method creates a CpuObc instance -- currently the OBC type II model is the
default (see paper). If the OBC type I model is desired change
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
to
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeI );
The created object is a static member of the class CpuObc;
when the force routine, cpuCalculateObcForces(), is called,
the static object is used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
extern "C" int
cpuSetObcParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* obcScaleFactors,
int includeAceApproximation,
RealOpenMM soluteDielectric, RealOpenMM solventDielectric, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\ncpuSetObcParameters: ";
// ---------------------------------------------------------------------------------------
// set log file if not NULL
if( log ){
SimTKOpenMMLog::setSimTKOpenMMLog( log );
}
// set OBC parameters (Type II)
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
obcParameters->setScaledRadiusFactors( obcScaleFactors );
obcParameters->setAtomicRadii( atomicRadii, SimTKOpenMMCommon::KcalAngUnits );
// dielectric constants
obcParameters->setSolventDielectric( solventDielectric );
obcParameters->setSoluteDielectric( soluteDielectric );
// ---------------------------------------------------------------------------------------
// create CpuObc instance that will calculate forces
CpuObc* cpuObc = new CpuObc( obcParameters );
// set static member for subsequent calls to calculate forces/energy
CpuImplicitSolvent::setCpuImplicitSolvent( cpuObc );
// set base file name, ...
//cpuObc->readInfoFile( "CpuImplicitSolventInfo" );
// include/do not include ACE approximation (nonpolar solvation)
cpuObc->setIncludeAceApproximation( includeAceApproximation );
// ---------------------------------------------------------------------------------------
// diagnostics
if( log ){
std::string state = cpuObc->getStateString( methodName );
(void) fprintf( log, "\n%s\nDone w/ setup\n", state.c_str() );
(void) fflush( log );
}
// ---------------------------------------------------------------------------------------
return 0;
}
/**---------------------------------------------------------------------------------------
Calculate implicit solvent forces and energy
@param atomCoordinates atom coordinates in Angstrom; format of array is
atomCoordinates[atom][3] in Angstrom
@param partialCharges partial charges
@param forces output forces in kcal/mol.A; format of array is
forces[atom][3]
@param energy output energy in kcal/mol
@param updateBornRadii if set, then Born radii are updated for current configuration;
otherwise radii correspond to configuration from previous iteration
Function calls a static method in CpuImplicitSolvent class to calculate forces/energy
@return result from CpuImplicitSolvent::computeImplicitSolventForces
--------------------------------------------------------------------------------------- */
extern "C" int
cpuCalculateImplicitSolventForces( RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges,
RealOpenMM** forces, RealOpenMM* energy,
int updateBornRadii ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\ncpuCalculateImplicitSolventForces: ";
// ---------------------------------------------------------------------------------------
int status = CpuImplicitSolvent::computeImplicitSolventForces( atomCoordinates, partialCharges,
forces, updateBornRadii );
*energy = CpuImplicitSolvent::getCpuImplicitSolvent()->getEnergy();
// printf( "\ncpuCalculateImplicitSolventForcesE=%.5e", *energy );
return status;
}
/**---------------------------------------------------------------------------------------
Retrieve the calculated energy from the static class member
The energy is calculated in cpuCalculateImplicitSolventForces()
along w/ the forces
@return the calculated energy from the static class member
--------------------------------------------------------------------------------------- */
extern "C" float cpuGetObcEnergy( void ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\ncpuGetObcEnergy: ";
// ---------------------------------------------------------------------------------------
RealOpenMM energy = CpuImplicitSolvent::getCpuImplicitSolvent()->getEnergy();
// printf( "\ncpuGetObcEnergyE=%.5e", energy );
return (float) energy;
}
/**---------------------------------------------------------------------------------------
Delete the Obc associated object(s)
@return 0 if static CpuObc object was set; else return -1
--------------------------------------------------------------------------------------- */
extern "C" int cpuDeleteObcParameters( void ){
return CpuImplicitSolvent::deleteCpuImplicitSolvent();
}
/**---------------------------------------------------------------------------------------
Get OBC scale factors given masses
@param numberOfAtoms number of atoms
@param masses input masses
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
extern "C" int getObcScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses,
RealOpenMM* scaleFactors ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetObcScaleFactorsGivenAtomMasses";
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
double scaleFactor = 0.8;
RealOpenMM mass = masses[atomI];
if ( mass < 1.2 && mass >= 1.0 ){ // hydrogen
scaleFactor = 0.85;
} else if( mass > 11.8 && mass < 12.2 ){ // carbon
scaleFactor = 0.72;
} else if( mass > 14.0 && mass < 15.0 ){ // nitrogen
scaleFactor = 0.79;
} else if( mass > 15.5 && mass < 16.5 ){ // oxygen
scaleFactor = 0.85;
} else if( mass > 31.5 && mass < 32.5 ){ // sulphur
scaleFactor = 0.96;
} else if( mass > 29.5 && mass < 30.5 ){ // phosphorus
scaleFactor = 0.86;
} else {
std::stringstream message;
message << methodName;
message << " Warning: mass for atom " << atomI << " mass=" << mass << "> not recognized.";
SimTKOpenMMLog::printMessage( message );
}
scaleFactors[atomI] = (RealOpenMM) scaleFactor;
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get OBC scale factors given atomic numbers
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
extern "C" int getObcScaleFactors( int numberOfAtoms, const int* atomicNumber, RealOpenMM* scaleFactors ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetObcScaleFactors";
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
double scaleFactor;
switch( atomicNumber[atomI] ){
case 1: // hydrogen
scaleFactor = 0.85;
break;
case 6: // carbon
scaleFactor = 0.72;
break;
case 7: // nitrogen
scaleFactor = 0.79;
break;
case 8: // oxygen
scaleFactor = 0.85;
break;
case 15: // phosphorus
scaleFactor = 0.86;
break;
case 16: // sulphur
scaleFactor = 0.85;
break;
default:
scaleFactor = 0.8;
std::stringstream message;
message << methodName;
message << " Warning: atom number=" << atomicNumber[atomI] << " for atom " << atomI << "> not handled -- useing default value.";
SimTKOpenMMLog::printMessage( message );
break;
}
scaleFactors[atomI] = (RealOpenMM) scaleFactor;
}
return SimTKOpenMMCommon::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Get GBSA radii
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param numberOfCovalentPartners input number of covalent partners for each atom
1 for H, 2,3,4 for C, ...;
the values are only used for C, N & O
@param indexOfCovalentPartner index of covalent partner -- used only for H
e.g. if atom 22 is a H and it is bonded to atom 24,
then indexOfCovalentPartner[22] = 24
@param gbsaRadii output GBSA radii
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
extern "C" int getGbsaRadii( int numberOfAtoms, const int* atomicNumber,
const int* numberOfCovalentPartners,
const int* indexOfCovalentPartner,
RealOpenMM* gbsaRadii ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\ngetGbsaRadii";
// ---------------------------------------------------------------------------------------
// loop over atoms
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// branch based on atomic number
double radius;
int atomIndex;
switch ( atomicNumber[atomI] ){
case 1: // H
// get index of covalent partner and validate
// radius is modified if heavy atom is N or O
atomIndex = indexOfCovalentPartner[atomI];
if( atomIndex < 0 || atomIndex >= numberOfAtoms ){
std::stringstream message;
message << methodName;
message << " Warning: atomIndex for covalent partner=" << atomIndex << " for atom " << atomI << " is invalid.";
SimTKOpenMMLog::printMessage( message );
} else if( atomicNumber[atomIndex] == 7 ){
radius = 1.15;
} else if( atomicNumber[atomIndex] == 8 ){
radius = 1.05;
} else {
radius = 1.25;
}
break;
case 3: // Li
radius = 2.432;
break;
case 6: // C
if( numberOfCovalentPartners[atomI] == 2 ){
radius = 1.825;
} else if( numberOfCovalentPartners[atomI] == 3 ){
radius = 1.875;
} else {
radius = 1.90;
}
break;
case 7: // N
if( numberOfCovalentPartners[atomI] == 4 ){
radius = 1.625;
} else if( numberOfCovalentPartners[atomI] == 1 ){
radius = 1.60;
} else {
radius = 1.7063;
}
break;
case 8: // O
if( numberOfCovalentPartners[atomI] == 1 ){
radius = 1.48;
} else {
radius = 1.535;
}
break;
case 9:
radius = 1.47;
break;
case 10:
radius = 1.39;
break;
case 11:
radius = 1.992;
break;
case 12:
radius = 1.70;
break;
case 14:
radius = 1.80;
break;
case 15:
radius = 1.87;
break;
case 16:
radius = 1.775;
break;
case 17:
radius = 1.735;
break;
case 18:
radius = 1.70;
break;
case 19:
radius = 2.123;
break;
case 20:
radius = 1.817;
break;
case 35:
radius = 1.90;
break;
case 36:
radius = 1.812;
break;
case 37:
radius = 2.26;
break;
case 53:
radius = 2.10;
break;
case 54:
radius = 1.967;
break;
case 55:
radius = 2.507;
break;
case 56:
radius = 2.188;
break;
default:
radius = 2.0;
break;
}
gbsaRadii[atomI] = (RealOpenMM) radius;
}
return SimTKOpenMMCommon::DefaultReturn;
}
platforms/reference/src/Gbsa/cpuObcInterface.h 0 → 100644
View file @ f3e77e33
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __CpuObcInterface_H__
#define __CpuObcInterface_H__
#ifdef __cplusplus
#define externC extern "C"
#else
#define externC extern
#endif
#include "SimTKOpenMMRealType.h"
#include <stdio.h>
/**---------------------------------------------------------------------------------------
Delete the Obc associated object(s)
@return 0 if static CpuObc object was set; else return -1
--------------------------------------------------------------------------------------- */
externC int cpuDeleteObcParameters( void );
/**---------------------------------------------------------------------------------------
Setup for Obc calculations from Gromacs
@param numberOfAtoms number of atoms
@param obcScaleFactors array of OBC scale factors (one entry each atom)
@param atomicRadii atomic radii in Angstrom (one entry each atom)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculations
@param soluteDielectric solute dielectric
@param solventDielectric solvent dielectric
@param log log reference -- if NULL, then errors/warnings
output to stderr
The method creates a CpuObc instance -- currently the OBC type II model is the
default (see paper). If the OBC type I model is desired change
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
to
ObcParameters* obcParameters = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeI );
The created object is a static member of the class CpuObc;
when the force routine, cpuCalculateObcForces(), is called,
the static object is used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
externC
int cpuSetObcParameters( int numberOfAtoms, RealOpenMM* atomicRadii, RealOpenMM* obcScaleFactors,
int includeAceApproximation, RealOpenMM soluteDielectric,
RealOpenMM solventDielectric, FILE* log );
/**---------------------------------------------------------------------------------------
Calculate implicit solvent forces and energy
@param atomCoordinates atom coordinates in Angstrom; format of array is
atomCoordinates[atom][3]
@param partialCharges partial atom charges
@param forces output forces in kcal/mol.A; format of array is
forces[atom][3]
@param energy energy
@param updateBornRadii if set, then Born radii are updated for current configuration;
otherwise radii correspond to configuration from previous iteration
Function calls a static method in CpuImplicitSolvent class to calculate forces/energy
@return result from CpuImplicitSolvent::computeImplicitSolventForces
--------------------------------------------------------------------------------------- */
externC int cpuCalculateImplicitSolventForces( RealOpenMM** atomCoordinates,
const RealOpenMM* partialChargesIn,
RealOpenMM** forces, RealOpenMM* energy,
int updateBornRadii );
/**---------------------------------------------------------------------------------------
Get OBC scale factors given masses
@param numberOfAtoms number of atoms
@param masses input masses
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
externC int getObcScaleFactorsGivenAtomMasses( int numberOfAtoms, const RealOpenMM* masses,
RealOpenMM* scaleFactors );
/**---------------------------------------------------------------------------------------
Get OBC scale factors given atomic numbers
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param scaleFactors output atomic numbers
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
externC int getObcScaleFactors( int numberOfAtoms, const int* atomicNumber,
RealOpenMM* scaleFactors );
/**---------------------------------------------------------------------------------------
Get GBSA radii
@param numberOfAtoms number of atoms
@param atomicNumber input atomic number for each atom
@param numberOfCovalentPartners input number of covalent partners for each atom
1 for H, ...; only used for C, N & O
@param indexOfCovalentPartner index of covalent partner -- only used for H
e.g. if atom 22 is a H and it is bonded to atom 24
then indexOfCovalentPartner[22] = 24
@param gbsaRadii output GBSA radii
@return SimTKOpenMMCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
externC int getGbsaRadii( int numberOfAtoms, const int* atomicNumber,
const int* numberOfCovalentPartners,
const int* indexOfCovalentPartner, RealOpenMM* gbsaRadii );
#undef externC
#endif
platforms/reference/src/Gbsa/gbsaCpu.cpp 0 → 100644
View file @ f3e77e33
/* 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 "CpuGbsa.h"
// 'gmxGbsa.h' has declarations for cpuSetGbsaParameters() and
// and cpuCalculateGbsaForces()
#include "gmxGbsa.h"
/**---------------------------------------------------------------------------------------
Example calling sequence:
float** forces = NULL;
float energy;
cpuSetGbsaParameters( mdatoms->nr, top, stdlog );
// if 'forces' array is initially null, then memory will
// be allocated by cpuCalculateGbsaForces()
// format for 'forces' array:
// forces[3][numberAtoms]
forces = cpuCalculateGbsaForces( x, mdatoms->chargeA, &energy,
forces, stdlog );
--------------------------------------------------------------------------------------- */
/**---------------------------------------------------------------------------------------
Setup for Gbsa calculations
@param numberOfAtoms number of atoms
@param top Gromacs t_topology (as in md.c)
@param log log reference (stdlog in md.c)
@param includeAceApproximation if true, then include nonpolar
ACE term in calculataions
@param innerDielectric nonsolvent dielectric
@param solventDielectric solvent dielectric
function creates a CpuGbsa instance; the created object is a static member of the
class CpuGbsa; when the force routine (a static method) is called, the object is
used to compute the forces and energy
@return 0
--------------------------------------------------------------------------------------- */
extern "C" int
cpuSetGbsaParameters( int numberOfAtoms, const t_topology* top, FILE* log,
bool includeAceApproximation,
float innerDielectric, float solventDielectric ){
// ---------------------------------------------------------------------------------------
static bool printGbsaParameters = false;
static const char* methodName = "\ncpuSetGbsaParameters: ";
// static bool debug = true;
static bool debug = false;
// ---------------------------------------------------------------------------------------
// toggle whether logging enabled
FILE* cpuGbsaLog = debug ? log : NULL;
CpuGbsa* cpuGbsa = new CpuGbsa( numberOfAtoms, top, cpuGbsaLog );
// include/not include ACE approximation (nonpolar solvation)
cpuGbsa->setIncludeAceApproximation( includeAceApproximation );
// cpuGbsa->setIncludeAceApproximation( false );
// cpuGbsa->setIncludeAceApproximation( true );
// default values are 1.0/1.0 -> gives values in 10.0*kcal/A (force) and 100*kcal/A (energy)
// if values are 0.41868/0.41868 -> gives values in Gromacs units kJ/nm
if( !debug ){
cpuGbsa->setForceConversionFactor( 0.41868f );
cpuGbsa->setEnergyConversionFactor( 0.41868f );
}
// dielectrics
cpuGbsa->getGbsaParameters()->setSolventDielectric( solventDielectric );
cpuGbsa->getGbsaParameters()->setInnerDielectric( innerDielectric );
// set static member for subsequent calls to calculate forces/energy
CpuGbsa::setCpuGbsa( cpuGbsa );
// ---------------------------------------------------------------------------------------
if( log != NULL ){
cpuGbsa->logState( methodName, log );
(void) fprintf( log, "%s done", methodName );
(void) fflush( log );
}
// ---------------------------------------------------------------------------------------
return 0;
}
/**---------------------------------------------------------------------------------------
Calculate Gbsa forces and energy \n
@param atomCoordinates Gromacs atom coordinates ('x' in md.c)
@param partialCharges Gromacs charges ('mdatoms->chargeA' in md.c)
@param energy output energy
@param forces output forces (format currently forces[3][numberOfAtoms])
if 'forces' reference is NULL upon input, then \n
memory is allocated for array; responsibilty of callee \n
to free \n
delete forces[0]; \n
delete forces; \n
Note memory is allocated as single block of size \n
3*numberOfAtoms*sizeof( float ), so only free \n
'forces[0]' \n
@param log log reference (stdlog in md.c)
Function calls a static method in CpuGbsa class to calculate forces/energy
Logging can be toggled by setting cpuGbsaLog to input 'log' or NULL
@return force array; energy is returned in *energy
--------------------------------------------------------------------------------------- */
extern "C" float**
cpuCalculateGbsaForces( const rvec *atomCoordinates, const float* partialCharges,
float** forces, FILE* log ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\ncpuCalculateGbsaForces: ";
// ---------------------------------------------------------------------------------------
// FILE* cpuGbsaLog = NULL;
FILE* cpuGbsaLog = log;
// (void) fprintf( log, "%s start", methodName );
// (void) fflush( log );
forces = CpuGbsa::computeGbsaForces( atomCoordinates, partialCharges, forces, cpuGbsaLog );
// *energy = CpuGbsa::getCpuGbsa()->getEnergy();
// (void) fprintf( log, "%s done E=%.4f", methodName, *energy );
// (void) fflush( log );
return forces;
}
/**---------------------------------------------------------------------------------------
Retrieve the calculated Gbsa energy from the static class member
@return the calculated Gbsa energy from the static class member
--------------------------------------------------------------------------------------- */
extern "C" float cpuGetEnergy( void ){
return CpuGbsa::getCpuGbsa()->getEnergy();
}
/**---------------------------------------------------------------------------------------
Delete the Gbsa associated object(s)
@return 0 if static CpuGbsa object was set; else return -1
--------------------------------------------------------------------------------------- */
extern "C" int cpuDeleteGbsaParameters( void ){
return CpuGbsa::deleteCpuGbsa();
}
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