Skip to content

GitLab

Commit 167ae8a0 authored by peastman's avatar peastman
Browse files

Merge pull request #1271 from peastman/gbvi 

Deleted GBVIForce
parents b20c20fe 14e78600
Hide whitespace changes
Inline Side-by-side
Showing with 3 additions and 5560 deletions
CMakeLists.txt
View file @ 167ae8a0
......@@ -87,7 +87,7 @@ ENDIF(${CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT})
# The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
SET(OPENMM_SOURCE_SUBDIRS . openmmapi olla libraries/jama libraries/quern libraries/lepton libraries/sfmt libraries/lbfgs libraries/hilbert libraries/csha1 platforms/reference serialization libraries/validate libraries/irrxml libraries/vecmath)
SET(OPENMM_SOURCE_SUBDIRS . openmmapi olla libraries/jama libraries/quern libraries/lepton libraries/sfmt libraries/lbfgs libraries/hilbert libraries/csha1 platforms/reference serialization libraries/irrxml libraries/vecmath)
IF(WIN32)
SET(OPENMM_SOURCE_SUBDIRS ${OPENMM_SOURCE_SUBDIRS} libraries/pthreads)
ELSE(WIN32)
......@@ -316,14 +316,14 @@ ENDIF (MSVC)
IF(OPENMM_BUILD_SHARED_LIB)
ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_LINK_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_VALIDATE_BUILDING_SHARED_LIBRARY -DPTHREAD_BUILDING_SHARED_LIBRARY")
SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_LINK_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DPTHREAD_BUILDING_SHARED_LIBRARY")
ENDIF(OPENMM_BUILD_SHARED_LIB)
SET(OPENMM_BUILD_STATIC_LIB OFF CACHE BOOL "Whether to build static OpenMM libraries")
IF(OPENMM_BUILD_STATIC_LIB)
ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
SET(EXTRA_COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_USE_STATIC_LIBRARIES -DLEPTON_USE_STATIC_LIBRARIES -DPTW32_STATIC_LIB")
SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_LINK_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_BUILDING_STATIC_LIBRARY -DOPENMMM_VALIDATE_BUILDING_STATIC_LIBRARY -DOPENMM_VALIDATE_BUILDING_STATIC_LIBRARY -DPTHREAD_BUILDING_STATIC_LIBRARY")
SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_LINK_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_BUILDING_STATIC_LIBRARY -DPTHREAD_BUILDING_STATIC_LIBRARY")
ENDIF(OPENMM_BUILD_STATIC_LIB)
IF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
......
libraries/validate/include/ValidateOpenMM.h deleted 100644 → 0
View file @ b20c20fe
#ifndef VALIDATE_OPENMM_H_
#define VALIDATE_OPENMM_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "OpenMM.h"
#include "../../../platforms/reference/include/ReferencePlatform.h"
#include "openmm/Context.h"
#include "openmm/System.h"
#include "ValidateWindowsIncludes.h"
// free-energy plugin includes
//#define INCLUDE_FREE_ENERGY_PLUGIN
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
#include "../../../plugins/freeEnergy/openmmapi/include/OpenMMFreeEnergy.h"
#include "../../../plugins/freeEnergy/openmmapi/include/openmm/freeEnergyKernels.h"
//#include "../../../plugins/freeEnergy/platforms/reference/include/ReferenceFreeEnergyPlatform.h"
#include "../../../plugins/freeEnergy/platforms/reference/include/ReferenceFreeEnergyKernelFactory.h"
#endif
#include <sstream>
#include <typeinfo>
#include <limits>
#include <cstdlib>
#include <cstring>
#include <cstdio>
namespace OpenMM {
typedef std::map< std::string, int > StringIntMap;
typedef StringIntMap::iterator StringIntMapI;
typedef StringIntMap::const_iterator StringIntMapCI;
typedef std::map< std::string, double > StringDoubleMap;
typedef StringDoubleMap::iterator StringDoubleMapI;
typedef StringDoubleMap::const_iterator StringDoubleMapCI;
typedef std::map< std::string, unsigned int > StringUIntMap;
typedef StringUIntMap::iterator StringUIntMapI;
typedef StringUIntMap::const_iterator StringUIntMapCI;
typedef std::vector< std::string > StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
typedef std::vector< int > IntVector;
typedef IntVector::iterator IntVectorI;
typedef IntVector::const_iterator IntVectorCI;
typedef std::map< std::string, std::vector<std::string> > StringStringVectorMap;
typedef StringStringVectorMap::iterator StringStringVectorMapI;
typedef StringStringVectorMap::const_iterator StringStringVectorMapCI;
/**
* Base class w/ common functionality
*/
class ValidateOpenMM {
public:
ValidateOpenMM();
~ValidateOpenMM();
// force names
static const std::string HARMONIC_BOND_FORCE;
static const std::string HARMONIC_ANGLE_FORCE;
static const std::string PERIODIC_TORSION_FORCE;
static const std::string RB_TORSION_FORCE;
static const std::string NB_FORCE;
static const std::string NB_SOFTCORE_FORCE;
static const std::string NB_EXCEPTION_FORCE;
static const std::string NB_EXCEPTION_SOFTCORE_FORCE;
static const std::string GBSA_OBC_FORCE;
static const std::string GBSA_OBC_SOFTCORE_FORCE;
static const std::string GBVI_FORCE;
static const std::string GBVI_SOFTCORE_FORCE;
static const std::string CM_MOTION_REMOVER;
static const std::string ANDERSEN_THERMOSTAT;
static const std::string CUSTOM_BOND_FORCE;
static const std::string CUSTOM_EXTERNAL_FORCE;
static const std::string CUSTOM_NONBONDED_FORCE;
/**
* Return true if input number is nan or infinity
*
* @param number number to test
*
* @return true if number is nan or infinity
*/
static int isNanOrInfinity( double number );
/**
* Get force name
*
* @param force OpenMM system force
*
* @return force name or "NA" if force is not recognized
*/
std::string getForceName(const Force& force ) const;
/**
* Copy force
*
* @param force OpenMM system force to copy
*
* @return force or NULL if not recognized
*/
Force* copyForce(const Force& force) const;
/**
* Get copy of input system, but omit forces
*
* @param systemToCopy system to copy
*
* @return copy of system but w/o forces
*/
System* copySystemExcludingForces( const System& systemToCopy ) const;
/**
*
* Set the velocities/positions of context2 to those of context1
*
* @param context1 context1
* @param context2 context2
*
* @return 0
*/
void synchContexts( const Context& context1, Context& context2 ) const;
/**
*
* Get log FILE* reference
*
* @return log
*
*/
FILE* getLog() const;
/**
*
* Set log FILE* reference
*
* @param log log
*
*/
void OPENMM_VALIDATE_EXPORT setLog( FILE* log );
/**---------------------------------------------------------------------------------------
Copy constraints
@param systemToCopy system whose constraints are to be copied
@param system system to add constraints to
@param log log file pointer -- may be NULL
--------------------------------------------------------------------------------------- */
void copyConstraints( const System& systemToCopy, System* system, FILE* log = NULL ) const;
/**---------------------------------------------------------------------------------------
Get force dependencies
@param forceName force to check if there exist any dependencies
@param returnVector vector of forces the input force is dependent on (example: GBSAOBC force requires Nonbonded force since
on Cuda platofrm they are computed in same loop and hence are inseparable)
--------------------------------------------------------------------------------------- */
void getForceDependencies( std::string forceName, StringVector& returnVector ) const;
/**
* Write masses to parameter file
*
* @param filePtr file to write masses to
* @param system write masses in system
*/
void writeMasses( FILE* filePtr, const System& system ) const;
/**
* Write constraints to parameter file
*
* @param filePtr file to write constraints to
* @param system write constraints in system
*
*/
void writeConstraints( FILE* filePtr, const System& system ) const;
/**
* Write harmonicBondForce parameters to file
*
* @param filePtr file to write forces to
* @param harmonicBondForce write harmonicBondForce parameters
*
*/
void writeHarmonicBondForce( FILE* filePtr, const HarmonicBondForce& harmonicBondForce ) const;
/**
* Write harmonicAngleForce parameters to file
*
* @param filePtr file to write forces to
* @param harmonicAngleForce write harmonicAngleForce parameters
*
*/
void writeHarmonicAngleForce( FILE* filePtr, const HarmonicAngleForce& harmonicAngleForce ) const;
/**
* Write rbTorsionForce parameters to file
*
* @param filePtr file to write forces to
* @param rbTorsionForce write rbTorsionForce parameters
*
*/
void writeRbTorsionForce( FILE* filePtr, const RBTorsionForce& rbTorsionForce ) const;
/**
* Write periodicTorsionForce parameters to file
*
* @param filePtr file to write forces to
* @param periodicTorsionForce write periodicTorsionForce parameters
*
*/
void writePeriodicTorsionForce( FILE* filePtr, const PeriodicTorsionForce& periodicTorsionForce ) const;
/**
* Write nonbonded parameters to file
*
* @param filePtr file to write forces to
* @param nonbondedForce write nonbondedForce parameters
*
*/
void writeNonbondedForce( FILE* filePtr, const NonbondedForce & nonbondedForce ) const;
/**
* Write GBSAOBCForce parameters to file
*
* @param filePtr file to write forces to
* @param gbsaObcForce write gbsaObcForce parameters
*
*/
void writeGbsaObcForce( FILE* filePtr, const GBSAOBCForce& gbsaObcForce ) const;
/**
* Write GBSA GB/VI Force parameters to file
*
* @param filePtr file to write forces to
* @param gbviObcForce write gbviObcForce parameters
*
*/
void writeGBVIForce( FILE* filePtr, const GBVIForce& gbviForce ) const;
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
/**
* Write nonbonded softcore parameters to file
*
* @param filePtr file to write forces to
* @param nonbondedForce write nonbondedForce parameters
*/
void writeNonbondedSoftcoreForce( FILE* filePtr, const NonbondedSoftcoreForce & nonbondedSoftcoreForce ) const;
/**
* Write GBSAOBCSoftcoreForce parameters to file
*
* @param filePtr file to write forces to
* @param gbsaObcForce write gbsaObcForce parameters
*
*/
void writeGbsaObcSoftcoreForce( FILE* filePtr, const GBSAOBCSoftcoreForce& gbsaObcForce ) const;
/**
* Write GBSA GB/VI softcore force parameters to file
*
* @param filePtr file to write forces to
* @param gbviObcForce write gbviObcForce parameters
*
*/
void writeGBVISoftcoreForce( FILE* filePtr, const GBVISoftcoreForce& gbviSoftcoreForce ) const;
#endif
/**
* Write coordinates, velocities, ... to file
*
* @param filePtr file to write Vec3 entries to
* @param vect3Array write array of Vec3
*
*/
void writeVec3( FILE* filePtr, const std::vector<Vec3>& vect3Array ) const;
/**
* Write context info to file (positions, velocities, forces, energies)
*
* @param filePtr file to write entries to
* @param context write context positions, velocities, forces, energies to file
*
*/
void writeContext( FILE* filePtr, const Context& context ) const;
/**
* Write integrator
*
* @param filePtr file to write integrator info to
* @param integrator write integrator info (time step, seed, ... as applicable)
*
*/
void writeIntegrator( FILE* filePtr, const Integrator& integrator ) const;
/**
* Write parameter file
* @param context context whose entries are to be written to file
* @param parameterFileName file name
*
*/
void writeParameterFile( const Context& context, const std::string& parameterFileName ) const;
private:
FILE* _log;
// map of force dependencies (e.g., GBSAObc requires NB force on CudaPlatform)
StringStringVectorMap _forceDependencies;
};
} // namespace OpenMM
#endif /*VALIDATE_OPENMM_H_*/
libraries/validate/include/ValidateOpenMMForces.h deleted 100644 → 0
View file @ b20c20fe
#ifndef VALIDATE_OPENMM_FORCES_H_
#define VALIDATE_OPENMM_FORCES_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "ValidateOpenMM.h"
namespace OpenMM {
typedef std::map< int, int > MapIntInt;
typedef MapIntInt::iterator MapIntIntI;
typedef MapIntInt::const_iterator MapIntIntCI;
// Helper class for ValidateOpenMMForces class used to store force results and facitilitate comparisons of
// resulting forces
class ForceValidationResult {
public:
ForceValidationResult(const Context& context1, const Context& context2, StringUIntMap& forceNamesMap);
~ForceValidationResult();
/**
* Get potential energy at specified platform index (0 || 1)
*
* @return potential energy for spercifed platform
*
* @throws OpenMMException if energyIndex is not 0 or 1
*/
double getPotentialEnergy(int energyIndex) const;
/**
* Get array of forces at specified platform index (0 || 1)
*
* @return array of force norms
*
* @throws OpenMMException if forceIndex is not 0 or 1
*/
std::vector<double> getForceNorms(int forceIndex) const;
/**
* Get array of forces at platform index (0 || 1)
*
* @return force array
*
* @throws OpenMMException if forceIndex is not 0 or 1
*/
std::vector<Vec3> getForces(int forceIndex) const;
/**
* Get maximum delta in force norm
*
* @param maxIndex return atom index of entry with maximum delta norm (optional)
*
* @return max delta in norm of forces
*/
double getMaxDeltaForceNorm(int* maxIndex = NULL) const;
/**
* Get maximum relative delta in force norm
*
* @param maxIndex return atom index of entry w/ maximum relative delta norm (optional)
*
* @return max relative delta in norm of forces
*/
double getMaxRelativeDeltaForceNorm(int* maxIndex = NULL) const;
/**
* Get maximum dot product between forces
*
* @param maxIndex return atom index of entry w/ maximum dot product between forces (optional)
*
* @return max dot product between forces
*/
double getMaxDotProduct(int* maxIndex = NULL) const;
/**
* Get name of force associated w/ computed results
*
* @return force name(s); if more than one force active in computation,
* then names are concatenated and separated by '::' (e.g., 'NB_FORCE::GBSA_OBC_FORCE')
*/
std::string getForceName() const;
/**
* Get platform name
*
* @param index index of platform (0 or 1)
*
* @return platform name
*
* @throws OpenMMException if index is not 0 or 1
*/
std::string getPlatformName(int index) const;
/**
* Register index of two entries that differ by a specified tolerance
*
* @param index inconsistent index
*
*/
void registerInconsistentForceIndex(int index, int value = 1);
/**
* Clear list of entries that differ by a specified tolerance
*
*/
void clearInconsistentForceIndexList();
/**
* Get list of entries that differ by a specified tolerance
*
*/
void getInconsistentForceIndexList(std::vector<int>& inconsistentIndices) const;
/**
* Get number of entries in inconsistent index list
*
*/
int getNumberOfInconsistentForceEntries() const;
/**
* Return true if nans were detected
*
* @return true if nans were detected
*/
int nansDetected() const;
/**
* Determine if force norms are valid
*
* @param tolerance tolerance
*/
void compareForceNorms(double tolerance);
/**
* Determine if forces are valid
*
* @param tolerance tolerance
*/
void compareForces(double tolerance);
private:
// computed potential energies and forces fror two platforms
double _potentialEnergies[2];
std::vector<Vec3> _forces[2];
// platform and force names
std::string _platforms[2];
std::vector<std::string> _forceNames;
// force norms and stat entries
std::vector<double> _norms[2];
std::vector<double> _normStatVectors[2];
// map of indicies w/ inconsistent force entries
std::map<int, int> _inconsistentForceIndicies;
// if set, then nans detected
int _nansDetected;
/**
* Calculate norms of vectors
*
*/
void _calculateNorms();
/**
* Calculate norms of specified vector
*
*/
void _calculateNormOfForceVector(int forceIndex);
// stat indices
static const int STAT_AVG = 0;
static const int STAT_STD = 1;
static const int STAT_MIN = 2;
static const int STAT_ID1 = 3;
static const int STAT_MAX = 4;
static const int STAT_ID2 = 5;
static const int STAT_CNT = 6;
static const int STAT_END = 7;
/**
* Find vector stats
*
*/
void _findStatsForDouble(const std::vector<double>& array, std::vector<double>& statVector) const;
};
// Class used to compare forces/potential energies on two platforms
class ValidateOpenMMForces : public ValidateOpenMM {
public:
OPENMM_VALIDATE_EXPORT ValidateOpenMMForces();
OPENMM_VALIDATE_EXPORT ~ValidateOpenMMForces();
/**
* Validate force/energy by comparing the results between the forces/energies computed on user-provided context platform
* with Reference platform
*
* @param context context reference
* @param summaryString output summary string of results of comparison (optional)
*
* @return number of inconsistent entries
*/
int OPENMM_VALIDATE_EXPORT compareWithReferencePlatform(Context& context, std::string* summaryString = NULL);
/**
* Validate force/energy by comparing the results between the forces/energies computed on two different platforms
*
* @param context context reference
* @param compareForces indices of force to be tested
* @param platform1 first platform to compute forces
* @param platform2 second platform to compute forces
*
* @return ForceValidationResult reference containing results of force/energy computations
* on the two input platforms
*/
ForceValidationResult* compareForce(Context& context, std::vector<int>& compareForces,
Platform& platform1, Platform& platform2) const;
/**
* Compare individual forces by comparing calculations across two platforms (platform associated w/ input context and
* comparisonPlatform)
*
* @param context context reference
* @param platform comparsion platform reference
* @param forceValidationResults output vector of ForceValidationResult ptrs (user is responsible for deleting
* individual ForceValidationResult objects)
*/
void compareOpenMMForces(Context& context, Platform& comparisonPlatform, std::vector<ForceValidationResult*>& forceValidationResults) const;
/**
* Determine if results are consistent
*
* @param forceValidationResults vector of ForceValidationResult ptrs to check if forces are consistent
*/
void checkForInconsistentForceEntries(std::vector<ForceValidationResult*>& forceValidationResults) const;
/**
* Get total number of force entries that are inconsistent
*
* @param forceValidationResults vector of ForceValidationResult ptrs to check if forces are consistent
*/
int getTotalNumberOfInconsistentForceEntries(std::vector<ForceValidationResult*>& forceValidationResults) const;
/**
* Get summary string of results
*
* @param forceValidationResults vector of ForceValidationResult ptrs
*/
std::string getSummary(std::vector<ForceValidationResult*>& forceValidationResults) const;
/**
* Set force tolerance
*
* @param tolerance force tolerance
*/
void setForceTolerance(double tolerance);
/**
* Get force tolerance
*
* @return force tolerance
*/
double getForceTolerance() const;
/*
* Get force tolerance for specified force
*
* @param forceName name of force
*
* @return force tolerance
*
* */
double getForceTolerance(const std::string& forceName) const;
/*
* Get max errors to print in summary string
*
* @return max errors to print
*
* */
int getMaxErrorsToPrint() const;
/*
* Set max errors to print in summary string
*
* @param maxErrorsToPrint max errors to print
*
* */
void setMaxErrorsToPrint(int maxErrorsToPrint);
/*
* Return true if force is not to be validated (Andersen thermostat, CM motion remover, ...)
*
* @param forceName force name
*
* @return true if force is not currently validated
**/
int isExcludedForce(std::string forceName) const;
private:
// initialize class entries
void _initialize();
/*
* Format output line
*
* @param tab tab
* @param description description
* @param value value
*
* @return string containing contents
*
* */
std::string _getLine(const std::string& tab,
const std::string& description,
const std::string& value) const;
std::vector<ForceValidationResult*> _forceValidationResults;
// max errors to print
int _maxErrorsToPrint;
// tolerence
double _forceTolerance;
// map of force tolerances to type (name)
StringDoubleMap _forceTolerances;
// forces to be excluded from validation
StringIntMap _forcesToBeExcluded;
};
} // namespace OpenMM
#endif /*VALIDATE_OPENMM_FORCES_H_*/
libraries/validate/include/ValidateWindowsIncludes.h deleted 100644 → 0
View file @ b20c20fe
#ifndef OPENMM_VALIDATE_WINDOW_INCLUDE_H_
#define OPENMM_VALIDATE_WINDOW_INCLUDE_H_
/*
* Shared libraries are messy in Visual Studio. We have to distinguish three
* cases:
* (1) this header is being used to build the OpenMMValidate shared library
* (dllexport)
* (2) this header is being used by a *client* of the OpenMMValidate shared
* library (dllimport)
* (3) we are building the OpenMMValidate static library, or the client is
* being compiled with the expectation of linking with the
* OpenMMValidate static library (nothing special needed)
* In the CMake script for building this library, we define one of the symbols
* OpenMMValidate_BUILDING_{SHARED|STATIC}_LIBRARY
* Client code normally has no special symbol defined, in which case we'll
* assume it wants to use the shared library. However, if the client defines
* the symbol OPENMM_VALIDATE_USE_STATIC_LIBRARIES we'll suppress the dllimport so
* that the client code can be linked with static libraries. Note that
* the client symbol is not library dependent, while the library symbols
* affect only the OpenMMValidate library, meaning that other libraries can
* be clients of this one. However, we are assuming all-static or all-shared.
*/
#ifdef _MSC_VER
// We don't want to hear about how sprintf is "unsafe".
#pragma warning(disable:4996)
#if defined(OPENMM_VALIDATE_BUILDING_SHARED_LIBRARY)
#define OPENMM_VALIDATE_EXPORT __declspec(dllexport)
// Keep MS VC++ quiet about lack of dll export of private members.
#pragma warning(disable:4251)
#elif defined(OPENMM_VALIDATE_BUILDING_STATIC_LIBRARY) || defined(OPENMM_VALIDATE_USE_STATIC_LIBRARIES)
#define OPENMM_VALIDATE_EXPORT
#else
#define OPENMM_VALIDATE_EXPORT __declspec(dllimport) // i.e., a client of a shared library
#endif
#else
#define OPENMM_VALIDATE_EXPORT // Linux, Mac
#endif
#endif // OPENMM_VALIDATE_WINDOW_INCLUDE_H_
libraries/validate/src/ValidateOpenMM.cpp deleted 100644 → 0
View file @ b20c20fe
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "ValidateOpenMM.h"
using namespace OpenMM;
// fixed force names
const std::string ValidateOpenMM::HARMONIC_BOND_FORCE = "HarmonicBond";
const std::string ValidateOpenMM::HARMONIC_ANGLE_FORCE = "HarmonicAngle";
const std::string ValidateOpenMM::PERIODIC_TORSION_FORCE = "PeriodicTorsion";
const std::string ValidateOpenMM::RB_TORSION_FORCE = "RbTorsion";
const std::string ValidateOpenMM::NB_FORCE = "Nb";
const std::string ValidateOpenMM::NB_SOFTCORE_FORCE = "NbSoftcore";
const std::string ValidateOpenMM::NB_EXCEPTION_FORCE = "NbException";
const std::string ValidateOpenMM::NB_EXCEPTION_SOFTCORE_FORCE = "NbSoftcoreException";
const std::string ValidateOpenMM::GBSA_OBC_FORCE = "GBSAOBC";
const std::string ValidateOpenMM::GBSA_OBC_SOFTCORE_FORCE = "GBSAOBCSoftcore";
const std::string ValidateOpenMM::GBVI_FORCE = "GBVI";
const std::string ValidateOpenMM::GBVI_SOFTCORE_FORCE = "GBVISoftcore";
const std::string ValidateOpenMM::CM_MOTION_REMOVER = "CMMotionRemover";
const std::string ValidateOpenMM::ANDERSEN_THERMOSTAT = "AndersenThermostat";
const std::string ValidateOpenMM::CUSTOM_BOND_FORCE = "CustomBond";
const std::string ValidateOpenMM::CUSTOM_EXTERNAL_FORCE = "CustomExternal";
const std::string ValidateOpenMM::CUSTOM_NONBONDED_FORCE = "CustomNonBonded";
ValidateOpenMM::ValidateOpenMM() {
_log = NULL;
// force dependencies
// these may need to specialized depending on platform since
// currently apply to Cuda platform, but may not apply to OpenCL platform
std::vector< std::string > nbVector;
nbVector.push_back( NB_FORCE );
_forceDependencies[GBSA_OBC_FORCE] = nbVector;
_forceDependencies[GBVI_FORCE] = nbVector;
std::vector< std::string > nbSoftcoreVector;
nbSoftcoreVector.push_back( NB_SOFTCORE_FORCE );
_forceDependencies[GBSA_OBC_SOFTCORE_FORCE] = nbSoftcoreVector;
_forceDependencies[GBVI_SOFTCORE_FORCE] = nbSoftcoreVector;
}
ValidateOpenMM::~ValidateOpenMM() {
}
int ValidateOpenMM::isNanOrInfinity( double number ){
return (number != number || number == std::numeric_limits<double>::infinity() || number == -std::numeric_limits<double>::infinity()) ? 1 : 0;
}
FILE* ValidateOpenMM::getLog() const {
return _log;
}
void ValidateOpenMM::setLog( FILE* log ){
_log = log;
}
std::string ValidateOpenMM::getForceName(const Force& force) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ValidateForce::getForceName";
// ---------------------------------------------------------------------------------------
std::string forceName = "NA";
// bond force
try {
const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
return HARMONIC_BOND_FORCE;
} catch( std::bad_cast ){
}
// angle force
try {
const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
return HARMONIC_ANGLE_FORCE;
} catch( std::bad_cast ){
}
// periodic torsion force
try {
const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
return PERIODIC_TORSION_FORCE;
} catch( std::bad_cast ){
}
// RB torsion force
try {
const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
return RB_TORSION_FORCE;
} catch( std::bad_cast ){
}
// nonbonded force
try {
const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
return NB_FORCE;
} catch( std::bad_cast ){
}
// GBSA OBC
try {
const GBSAOBCForce& obcForce = dynamic_cast<const GBSAOBCForce&>(force);
return GBSA_OBC_FORCE;
} catch( std::bad_cast ){
}
// GB/VI
try {
const GBVIForce& obcForce = dynamic_cast<const GBVIForce&>(force);
return GBVI_FORCE;
} catch( std::bad_cast ){
}
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
// free energy plugin forces
// nonbonded softcore
try {
const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
return NB_SOFTCORE_FORCE;
} catch( std::bad_cast ){
}
// GBSA OBC softcore
try {
const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
return GBSA_OBC_SOFTCORE_FORCE;
} catch( std::bad_cast ){
}
// GB/VI softcore
try {
const GBVISoftcoreForce& gbviForce = dynamic_cast<const GBVISoftcoreForce&>(force);
return GBVI_SOFTCORE_FORCE;
} catch( std::bad_cast ){
}
#endif
// CMMotionRemover
try {
const CMMotionRemover& cMMotionRemover = dynamic_cast<const CMMotionRemover&>(force);
return CM_MOTION_REMOVER;
} catch( std::bad_cast ){
}
// AndersenThermostat
try {
const AndersenThermostat & andersenThermostat = dynamic_cast<const AndersenThermostat&>(force);
return ANDERSEN_THERMOSTAT;
} catch( std::bad_cast ){
}
// CustomBondForce
try {
const CustomBondForce & customBondForce = dynamic_cast<const CustomBondForce&>(force);
return CUSTOM_BOND_FORCE;
} catch( std::bad_cast ){
}
// CustomExternalForce
try {
const CustomExternalForce& customExternalForce = dynamic_cast<const CustomExternalForce&>(force);
return CUSTOM_EXTERNAL_FORCE;
} catch( std::bad_cast ){
}
// CustomNonbondedForce
try {
const CustomNonbondedForce& customNonbondedForce = dynamic_cast<const CustomNonbondedForce&>(force);
return CUSTOM_NONBONDED_FORCE;
} catch( std::bad_cast ){
}
return forceName;
}
Force* ValidateOpenMM::copyForce(const Force& force) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ValidateForce::copyForce";
// ---------------------------------------------------------------------------------------
// bond force
try {
const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
return new HarmonicBondForce( harmonicBondForce );
} catch( std::bad_cast ){
}
// angle force
try {
const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
return new HarmonicAngleForce( harmonicAngleForce );
} catch( std::bad_cast ){
}
// periodic torsion force
try {
const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
return new PeriodicTorsionForce( periodicTorsionForce );
} catch( std::bad_cast ){
}
// RB torsion force
try {
const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
return new RBTorsionForce( rBTorsionForce );
} catch( std::bad_cast ){
}
// nonbonded force
try {
const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
return new NonbondedForce( nbForce );
} catch( std::bad_cast ){
}
// GBSA OBC
try {
const GBSAOBCForce& obcForce = dynamic_cast<const GBSAOBCForce&>(force);
return new GBSAOBCForce( obcForce );
} catch( std::bad_cast ){
}
// GB/VI
try {
const GBVIForce& gbviForce = dynamic_cast<const GBVIForce&>(force);
return new GBVIForce( gbviForce );
} catch( std::bad_cast ){
}
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
// free energy plugin forces
// nonbonded softcore
try {
const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
return new NonbondedSoftcoreForce( nbForce );
} catch( std::bad_cast ){
}
// GBSA OBC softcore
try {
const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
return new GBSAOBCSoftcoreForce( obcForce );
} catch( std::bad_cast ){
}
// GB/VI softcore
try {
const GBVISoftcoreForce& gbviForce = dynamic_cast<const GBVISoftcoreForce&>(force);
return new GBVISoftcoreForce( gbviForce );
} catch( std::bad_cast ){
}
#endif
// CMMotionRemover
try {
const CMMotionRemover& cMMotionRemover = dynamic_cast<const CMMotionRemover&>(force);
return new CMMotionRemover( cMMotionRemover );
} catch( std::bad_cast ){
}
// AndersenThermostat
try {
const AndersenThermostat & andersenThermostat = dynamic_cast<const AndersenThermostat&>(force);
return new AndersenThermostat( andersenThermostat );
} catch( std::bad_cast ){
}
// CustomBondForce
try {
const CustomBondForce & customBondForce = dynamic_cast<const CustomBondForce&>(force);
return new CustomBondForce( customBondForce );
} catch( std::bad_cast ){
}
// CustomExternalForce
try {
const CustomExternalForce& customExternalForce = dynamic_cast<const CustomExternalForce&>(force);
return new CustomExternalForce( customExternalForce );
} catch( std::bad_cast ){
}
// CustomNonbondedForce
try {
const CustomNonbondedForce& customNonbondedForce = dynamic_cast<const CustomNonbondedForce&>(force);
return new CustomNonbondedForce( customNonbondedForce );
} catch( std::bad_cast ){
}
return NULL;
}
/**
* Get copy of input system, but leave out forces
*
* @param systemToCopy system to copy
* @return system w/o forces
*/
System* ValidateOpenMM::copySystemExcludingForces( const System& systemToCopy ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ValidateOpenMM::copySystemExcludingForces";
// ---------------------------------------------------------------------------------------
// create new system and add masses and box dimensions
System* system = new System();
for(int ii = 0; ii < systemToCopy.getNumParticles(); ii++ ){
double mass = systemToCopy.getParticleMass( ii );
system->addParticle( mass );
}
Vec3 a, b, c;
systemToCopy.getDefaultPeriodicBoxVectors( a, b, c );
system->setDefaultPeriodicBoxVectors( a, b, c );
copyConstraints( systemToCopy, system );
return system;
}
/**---------------------------------------------------------------------------------------
Set the velocities/positions of context2 to those of context1
@param context1 context1
@param context2 context2
@return 0
--------------------------------------------------------------------------------------- */
void ValidateOpenMM::synchContexts( const Context& context1, Context& context2 ) const {
// ---------------------------------------------------------------------------------------
//static const char* methodName = "ValidateOpenMM::synchContexts: ";
// ---------------------------------------------------------------------------------------
const State state = context1.getState(State::Positions | State::Velocities);
const std::vector<Vec3>& positions = state.getPositions();
const std::vector<Vec3>& velocities = state.getVelocities();
context2.setPositions( positions );
context2.setVelocities( velocities );
return;
}
/**---------------------------------------------------------------------------------------
Copy constraints
@param systemToCopy system whose constraints are to be copied
@param system system to add constraints to
@param log log file pointer -- may be NULL
--------------------------------------------------------------------------------------- */
void ValidateOpenMM::copyConstraints( const System& systemToCopy, System* system, FILE* log ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "copyConstraints";
// ---------------------------------------------------------------------------------------
for( int ii = 0; ii < systemToCopy.getNumConstraints(); ii++ ){
int particle1, particle2;
double distance;
systemToCopy.getConstraintParameters( ii, particle1, particle2, distance );
system->addConstraint( particle1, particle2, distance );
}
}
/**---------------------------------------------------------------------------------------
Get force dependencies
@param forceName force to check if there exist any dependencies
@param returnVector vector of forces the input force is dependent on (example: GBSAOBC force requires Nonbonded force)
--------------------------------------------------------------------------------------- */
void ValidateOpenMM::getForceDependencies( std::string forceName, StringVector& returnVector ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "getForceDependencies";
// ---------------------------------------------------------------------------------------
StringStringVectorMapCI dependencyVector = _forceDependencies.find( forceName );
if( dependencyVector != _forceDependencies.end() ){
StringVector dependices = (*dependencyVector).second;
for( unsigned int ii = 0; ii < dependices.size(); ii++ ){
returnVector.push_back( dependices[ii] );
}
}
}
/**
* Write masses and box dimensions to parameter file
*/
void ValidateOpenMM::writeMasses( FILE* filePtr, const System& system ) const {
(void) fprintf( filePtr, "Masses %d\n", system.getNumParticles() );
for(int ii = 0; ii < system.getNumParticles(); ii++ ){
double mass = system.getParticleMass( ii );
(void) fprintf( filePtr, "%8d %14.7e\n", ii, mass );
}
Vec3 a, b, c;
system.getDefaultPeriodicBoxVectors( a, b, c);
(void) fprintf( filePtr, "Box %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n",
a[0], a[1], a[2], b[0], b[1], b[2], c[0], c[1], c[2] );
}
/**
* Write constraints to parameter file
*/
void ValidateOpenMM::writeConstraints( FILE* filePtr, const System& system ) const {
(void) fprintf( filePtr, "Constraints %d\n", system.getNumConstraints() );
for(int ii = 0; ii < system.getNumConstraints(); ii++ ){
int particle1, particle2;
double distance;
system.getConstraintParameters( ii, particle1, particle2, distance );
(void) fprintf( filePtr, "%8d %8d %8d %14.7e\n", ii, particle1, particle2, distance );
}
}
/**
* Write harmonicBondForce parameters to file
*/
void ValidateOpenMM::writeHarmonicBondForce( FILE* filePtr, const HarmonicBondForce& harmonicBondForce ) const {
(void) fprintf( filePtr, "HarmonicBondForce %d\n", harmonicBondForce.getNumBonds() );
for(int ii = 0; ii < harmonicBondForce.getNumBonds(); ii++ ){
int particle1, particle2;
double length, k;
harmonicBondForce.getBondParameters( ii, particle1, particle2, length, k );
(void) fprintf( filePtr, "%8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, length, k);
}
}
/**
* Write harmonicAngleForce parameters to file
*/
void ValidateOpenMM::writeHarmonicAngleForce( FILE* filePtr, const HarmonicAngleForce& harmonicAngleForce ) const {
(void) fprintf( filePtr, "HarmonicAngleForce %d\n", harmonicAngleForce.getNumAngles() );
for(int ii = 0; ii < harmonicAngleForce.getNumAngles(); ii++ ){
int particle1, particle2, particle3;
double angle, k;
harmonicAngleForce.getAngleParameters( ii, particle1, particle2, particle3, angle, k );
(void) fprintf( filePtr, "%8d %8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, particle3, angle, k);
}
}
/**
* Write rbTorsionForce parameters to file
*/
void ValidateOpenMM::writeRbTorsionForce( FILE* filePtr, const RBTorsionForce& rbTorsionForce ) const {
(void) fprintf( filePtr, "RBTorsionForce %d\n", rbTorsionForce.getNumTorsions() );
for(int ii = 0; ii < rbTorsionForce.getNumTorsions(); ii++ ){
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
rbTorsionForce.getTorsionParameters( ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
(void) fprintf( filePtr, "%8d %8d %8d %8d %8d %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 );
}
}
/**
* Write periodicTorsionForce parameters to file
*/
void ValidateOpenMM::writePeriodicTorsionForce( FILE* filePtr, const PeriodicTorsionForce& periodicTorsionForce ) const {
(void) fprintf( filePtr, "PeriodicTorsionForce %d\n", periodicTorsionForce.getNumTorsions() );
for(int ii = 0; ii < periodicTorsionForce.getNumTorsions(); ii++ ){
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
periodicTorsionForce.getTorsionParameters( ii, particle1, particle2, particle3, particle4, periodicity, phase, k );
(void) fprintf( filePtr, "%8d %8d %8d %8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, particle3, particle4, periodicity, phase, k );
}
}
/**
* Write GBSAOBCForce parameters to file
*/
void ValidateOpenMM::writeGbsaObcForce( FILE* filePtr, const GBSAOBCForce& gbsaObcForce ) const {
(void) fprintf( filePtr, "GBSAOBCForce %d\n", gbsaObcForce.getNumParticles() );
for(int ii = 0; ii < gbsaObcForce.getNumParticles(); ii++ ){
double charge, radius, scalingFactor;
gbsaObcForce.getParticleParameters( ii, charge, radius, scalingFactor );
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e\n", ii, charge, radius, scalingFactor );
}
(void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbsaObcForce.getSoluteDielectric() );
(void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbsaObcForce.getSolventDielectric() );
}
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
/**
* Write GBSAOBCSoftcoreForce parameters to file
*/
void ValidateOpenMM::writeGbsaObcSoftcoreForce( FILE* filePtr, const GBSAOBCSoftcoreForce& gbsaObcForce ) const {
(void) fprintf( filePtr, "GBSAOBCSoftcoreForce %d\n", gbsaObcForce.getNumParticles() );
for(int ii = 0; ii < gbsaObcForce.getNumParticles(); ii++ ){
double charge, radius, scalingFactor;
double nonPolarScalingFactor;
gbsaObcForce.getParticleParameters( ii, charge, radius, scalingFactor, nonPolarScalingFactor );
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e %14.7e\n", ii, charge, radius, scalingFactor, nonPolarScalingFactor );
}
(void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbsaObcForce.getSoluteDielectric() );
(void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbsaObcForce.getSolventDielectric() );
}
#endif
/**
* Write GBSA GB/VI Force parameters to file
*/
void ValidateOpenMM::writeGBVIForce( FILE* filePtr, const GBVIForce& gbviForce ) const {
(void) fprintf( filePtr, "GBVIForce %d\n", gbviForce.getNumParticles() );
for(int ii = 0; ii < gbviForce.getNumParticles(); ii++ ){
double charge, radius, gamma;
gbviForce.getParticleParameters( ii, charge, radius, gamma );
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e\n", ii, charge, radius, gamma );
}
(void) fprintf( filePtr, "GBVIBonds %d\n", gbviForce.getNumBonds() );
for(int ii = 0; ii < gbviForce.getNumBonds(); ii++ ){
int atomI, atomJ;
double bondLength;
gbviForce.getBondParameters( ii, atomI, atomJ, bondLength );
(void) fprintf( filePtr, "%8d %8d %8d %14.7e\n", ii, atomI, atomJ, bondLength );
}
(void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbviForce.getSoluteDielectric() );
(void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbviForce.getSolventDielectric() );
}
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
/**
* Write GBSA GB/VI softcore force parameters to file
*/
void ValidateOpenMM::writeGBVISoftcoreForce( FILE* filePtr, const GBVISoftcoreForce& gbviSoftcoreForce ) const {
(void) fprintf( filePtr, "GBVISoftcoreForce %d\n", gbviSoftcoreForce.getNumParticles() );
for(int ii = 0; ii < gbviSoftcoreForce.getNumParticles(); ii++ ){
double charge, radius, gamma;
double nonPolarScalingFactor;
gbviSoftcoreForce.getParticleParameters( ii, charge, radius, gamma, nonPolarScalingFactor );
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e %14.7e\n", ii, charge, radius, gamma, nonPolarScalingFactor );
}
(void) fprintf( filePtr, "GBVISoftcoreBonds %d\n", gbviSoftcoreForce.getNumBonds() );
for(int ii = 0; ii < gbviSoftcoreForce.getNumBonds(); ii++ ){
int atomI, atomJ;
double bondLength;
gbviSoftcoreForce.getBondParameters( ii, atomI, atomJ, bondLength );
(void) fprintf( filePtr, "%8d %8d %8d %14.7e\n", ii, atomI, atomJ, bondLength );
}
(void) fprintf( filePtr, "SoluteDielectric %14.7e\n", gbviSoftcoreForce.getSoluteDielectric() );
(void) fprintf( filePtr, "SolventDielectric %14.7e\n", gbviSoftcoreForce.getSolventDielectric() );
(void) fprintf( filePtr, "BornRadiusScalingMethod %d\n", gbviSoftcoreForce.getBornRadiusScalingMethod() );
(void) fprintf( filePtr, "QuinticLowerLimitFactor %14.7e\n", gbviSoftcoreForce.getQuinticLowerLimitFactor() );
(void) fprintf( filePtr, "QuinticUpperBornRadiusLimit %14.7e\n", gbviSoftcoreForce.getQuinticUpperBornRadiusLimit() );
}
#endif
/**
* Write coordinates to file
*/
void ValidateOpenMM::writeVec3( FILE* filePtr, const std::vector<Vec3>& vect3Array ) const {
for( unsigned int ii = 0; ii < vect3Array.size(); ii++ ){
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e\n", ii,
vect3Array[ii][0], vect3Array[ii][1], vect3Array[ii][2] );
}
}
/**
* Write context info to file
*/
void ValidateOpenMM::writeContext( FILE* filePtr, const Context& context ) const {
State state = context.getState( State::Positions | State::Velocities | State::Forces | State::Energy );
std::vector<Vec3> positions = state.getPositions();
std::vector<Vec3> velocities = state.getVelocities();
std::vector<Vec3> forces = state.getForces();
(void) fprintf( filePtr, "Positions %zu\n", positions.size() );
writeVec3( filePtr, positions );
(void) fprintf( filePtr, "Velocities %zu\n", velocities.size() );
writeVec3( filePtr, velocities );
(void) fprintf( filePtr, "Forces %zu\n", forces.size() );
writeVec3( filePtr, forces );
(void) fprintf( filePtr, "KineticEnergy %14.7e\n", state.getKineticEnergy() );
(void) fprintf( filePtr, "PotentialEnergy %14.7e\n", state.getPotentialEnergy() );
}
/**
* Write nonbonded parameters to file
*/
void ValidateOpenMM::writeNonbondedForce( FILE* filePtr, const NonbondedForce & nonbondedForce ) const {
// charge and vdw parameters
(void) fprintf( filePtr, "NonbondedForce %d\n", nonbondedForce.getNumParticles() );
for(int ii = 0; ii < nonbondedForce.getNumParticles(); ii++ ){
double charge, sigma, epsilon;
nonbondedForce.getParticleParameters( ii, charge, sigma, epsilon );
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e\n", ii, charge, sigma, epsilon );
}
// cutoff, dielectric, Ewald tolerance
(void) fprintf( filePtr, "CutoffDistance %14.7e\n", nonbondedForce.getCutoffDistance() );
(void) fprintf( filePtr, "RFDielectric %14.7e\n", nonbondedForce.getReactionFieldDielectric() );
(void) fprintf( filePtr, "EwaldRTolerance %14.7e\n", nonbondedForce.getEwaldErrorTolerance() );
// cutoff mode
std::string nonbondedForceMethod;
switch( nonbondedForce.getNonbondedMethod() ){
case NonbondedForce::NoCutoff:
nonbondedForceMethod = "NoCutoff";
break;
case NonbondedForce::CutoffNonPeriodic:
nonbondedForceMethod = "CutoffNonPeriodic";
break;
case NonbondedForce::CutoffPeriodic:
nonbondedForceMethod = "CutoffPeriodic";
break;
case NonbondedForce::Ewald:
nonbondedForceMethod = "Ewald";
break;
case NonbondedForce::PME:
nonbondedForceMethod = "PME";
break;
default:
nonbondedForceMethod = "Unknown";
}
(void) fprintf( filePtr, "NonbondedForceMethod %s\n", nonbondedForceMethod.c_str() );
(void) fprintf( filePtr, "NonbondedForceExceptions %d\n", nonbondedForce.getNumExceptions() );
for(int ii = 0; ii < nonbondedForce.getNumExceptions(); ii++ ){
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbondedForce.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon );
(void) fprintf( filePtr, "%8d %8d %8d %14.7e %14.7e %14.7e\n", ii, particle1, particle2, chargeProd, sigma, epsilon );
}
}
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
/**
* Write nonbonded softcore parameters to file
*/
void ValidateOpenMM::writeNonbondedSoftcoreForce( FILE* filePtr, const NonbondedSoftcoreForce & nonbondedSoftcoreForce ) const {
(void) fprintf( filePtr, "NonbondedSoftcoreForce %d\n", nonbondedSoftcoreForce.getNumParticles() );
for(int ii = 0; ii < nonbondedSoftcoreForce.getNumParticles(); ii++ ){
double charge, sigma, epsilon, softCoreLJ;
nonbondedSoftcoreForce.getParticleParameters( ii, charge, sigma, epsilon, softCoreLJ );
(void) fprintf( filePtr, "%8d %14.7e %14.7e %14.7e %14.7e\n", ii, charge, sigma, epsilon, softCoreLJ );
}
(void) fprintf( filePtr, "CutoffDistance %14.7e\n", nonbondedSoftcoreForce.getCutoffDistance() );
(void) fprintf( filePtr, "RFDielectric %14.7e\n", nonbondedSoftcoreForce.getReactionFieldDielectric() );
(void) fprintf( filePtr, "EwaldRTolerance %14.7e\n", nonbondedSoftcoreForce.getEwaldErrorTolerance() );
std::string nonbondedSoftcoreForceMethod;
switch( nonbondedSoftcoreForce.getNonbondedMethod() ){
case NonbondedSoftcoreForce::NoCutoff:
nonbondedSoftcoreForceMethod = "NoCutoff";
break;
/*
case NonbondedSoftcoreForce::CutoffNonPeriodic:
nonbondedSoftcoreForceMethod = "CutoffNonPeriodic";
break;
case NonbondedSoftcoreForce::CutoffPeriodic:
nonbondedSoftcoreForceMethod = "CutoffPeriodic";
break;
case NonbondedSoftcoreForce::Ewald:
nonbondedSoftcoreForceMethod = "Ewald";
break;
case NonbondedSoftcoreForce::PME:
nonbondedSoftcoreForceMethod = "PME";
break;
*/
default:
nonbondedSoftcoreForceMethod = "Unknown";
}
(void) fprintf( filePtr, "NonbondedSoftcoreForceMethod %s\n", nonbondedSoftcoreForceMethod.c_str() );
(void) fprintf( filePtr, "NonbondedSoftcoreForceExceptions %d\n", nonbondedSoftcoreForce.getNumExceptions() );
for(int ii = 0; ii < nonbondedSoftcoreForce.getNumExceptions(); ii++ ){
int particle1, particle2;
double chargeProd, sigma, epsilon, softCoreLJ;
nonbondedSoftcoreForce.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon, softCoreLJ );
(void) fprintf( filePtr, "%8d %8d %8d %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, chargeProd, sigma, epsilon, softCoreLJ );
}
return;
}
#endif
void ValidateOpenMM::writeIntegrator( FILE* filePtr, const Integrator& integrator ) const {
// LangevinIntegrator
try {
const LangevinIntegrator langevinIntegrator = dynamic_cast<const LangevinIntegrator&>(integrator);
(void) fprintf( filePtr, "Integrator LangevinIntegrator\n" );
(void) fprintf( filePtr, "StepSize %14.7e\n", langevinIntegrator.getStepSize() );
(void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", langevinIntegrator.getConstraintTolerance() );
(void) fprintf( filePtr, "Temperature %14.7e\n", langevinIntegrator.getTemperature() );
(void) fprintf( filePtr, "Friction %14.7e\n", langevinIntegrator.getFriction() );
(void) fprintf( filePtr, "RandomNumberSeed %d\n", langevinIntegrator.getRandomNumberSeed() );
return;
} catch( std::bad_cast ){
}
// VariableLangevinIntegrator
try {
const VariableLangevinIntegrator& langevinIntegrator = dynamic_cast<const VariableLangevinIntegrator&>(integrator);
(void) fprintf( filePtr, "Integrator VariableLangevinIntegrator\n" );
(void) fprintf( filePtr, "StepSize %14.7e\n", langevinIntegrator.getStepSize() );
(void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", langevinIntegrator.getConstraintTolerance() );
(void) fprintf( filePtr, "Temperature %14.7e\n", langevinIntegrator.getTemperature() );
(void) fprintf( filePtr, "Friction %14.7e\n", langevinIntegrator.getFriction() );
(void) fprintf( filePtr, "RandomNumberSeed %d\n", langevinIntegrator.getRandomNumberSeed() );
(void) fprintf( filePtr, "ErrorTolerance %14.7e\n", langevinIntegrator.getErrorTolerance() );
return;
} catch( std::bad_cast ){
}
// VerletIntegrator
try {
const VerletIntegrator& verletIntegrator = dynamic_cast<const VerletIntegrator&>(integrator);
(void) fprintf( filePtr, "Integrator VerletIntegrator\n" );
(void) fprintf( filePtr, "StepSize %14.7e\n", verletIntegrator.getStepSize() );
(void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", verletIntegrator.getConstraintTolerance() );
return;
} catch( std::bad_cast ){
}
// VariableVerletIntegrator
try {
const VariableVerletIntegrator & variableVerletIntegrator = dynamic_cast<const VariableVerletIntegrator&>(integrator);
(void) fprintf( filePtr, "Integrator VariableVerletIntegrator\n" );
(void) fprintf( filePtr, "StepSize %14.7e\n", variableVerletIntegrator.getStepSize() );
(void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", variableVerletIntegrator.getConstraintTolerance() );
(void) fprintf( filePtr, "ErrorTolerance %14.7e\n", variableVerletIntegrator.getErrorTolerance() );
return;
} catch( std::bad_cast ){
}
// BrownianIntegrator
try {
const BrownianIntegrator& brownianIntegrator = dynamic_cast<const BrownianIntegrator&>(integrator);
(void) fprintf( filePtr, "Integrator BrownianIntegrator\n" );
(void) fprintf( filePtr, "StepSize %14.7e\n", brownianIntegrator.getStepSize() );
(void) fprintf( filePtr, "ConstraintTolerance %14.7e\n", brownianIntegrator.getConstraintTolerance() );
(void) fprintf( filePtr, "Temperature %14.7e\n", brownianIntegrator.getTemperature() );
(void) fprintf( filePtr, "Friction %14.7e\n", brownianIntegrator.getFriction() );
(void) fprintf( filePtr, "RandomNumberSeed %d\n", brownianIntegrator.getRandomNumberSeed() );
return;
} catch( std::bad_cast ){
}
if( getLog() ){
(void) fprintf( getLog(), "Integrator not recognized." );
(void) fflush( getLog() );
}
return;
}
/**
* Write parameter file: Custom forces not implemented
* Mesage is sent to stderr if a force is not recognized
*
*/
void ValidateOpenMM::writeParameterFile( const Context& context, const std::string& parameterFileName ) const {
// open file
FILE* filePtr = fopen( parameterFileName.c_str(), "w" );
const System& system = context.getSystem();
const Integrator& integrator = context.getIntegrator();
// (void) fprintf( filePtr, "Version %s\n", versionString.c_str() );
(void) fprintf( filePtr, "Particles %8d\n", system.getNumParticles() );
writeMasses( filePtr, system );
(void) fprintf( filePtr, "NumberOfForces %8d\n", system.getNumForces() );
// print active forces and relevant parameters
for(int i = 0; i < system.getNumForces(); ++i) {
int hit = 0;
const Force& force = system.getForce(i);
// bond
if( !hit ){
try {
const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
writeHarmonicBondForce( filePtr, harmonicBondForce );
hit++;
} catch( std::bad_cast ){
}
}
// angle
if( !hit ){
try {
const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
writeHarmonicAngleForce( filePtr, harmonicAngleForce );
hit++;
} catch( std::bad_cast ){
}
}
// PeriodicTorsionForce
if( !hit ){
try {
const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
writePeriodicTorsionForce( filePtr, periodicTorsionForce );
hit++;
} catch( std::bad_cast ){
}
}
// RBTorsionForce
if( !hit ){
try {
const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
writeRbTorsionForce( filePtr, rBTorsionForce );
hit++;
} catch( std::bad_cast ){
}
}
// nonbonded
if( !hit ){
try {
const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
writeNonbondedForce( filePtr, nbForce );
hit++;
} catch( std::bad_cast ){
}
}
// nonbonded softcore
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
if( !hit ){
try {
const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
writeNonbondedSoftcoreForce( filePtr, nbForce );
hit++;
} catch( std::bad_cast ){
}
}
#endif
// GBSA OBC
if( !hit ){
try {
const GBSAOBCForce& obcForce = dynamic_cast<const GBSAOBCForce&>(force);
writeGbsaObcForce( filePtr, obcForce );
hit++;
} catch( std::bad_cast ){
}
}
// GBSA OBC softcore
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
if( !hit ){
try {
const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
writeGbsaObcSoftcoreForce( filePtr, obcForce );
hit++;
} catch( std::bad_cast ){
}
}
#endif
// GB/VI
if( !hit ){
try {
const GBVIForce& obcForce = dynamic_cast<const GBVIForce&>(force);
writeGBVIForce( filePtr, obcForce );
hit++;
} catch( std::bad_cast ){
}
}
// GB/VI softcore
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
if( !hit ){
try {
const GBVISoftcoreForce& obcForce = dynamic_cast<const GBVISoftcoreForce&>(force);
writeGBVISoftcoreForce( filePtr, obcForce );
hit++;
} catch( std::bad_cast ){
}
}
#endif
// COM
if( !hit ){
try {
const CMMotionRemover& cMMotionRemover = dynamic_cast<const CMMotionRemover&>(force);
(void) fprintf( filePtr, "CMMotionRemover %d\n", cMMotionRemover.getFrequency() );
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
char buffer[1024];
(void) sprintf( buffer, " %2d force not recognized.\n", i );
(void) fprintf( stderr, "%s\n", buffer );
// throw OpenMMException( buffer );
}
}
// constraints
writeConstraints( filePtr, system );
// context
writeContext( filePtr, context );
// integrator
writeIntegrator( filePtr, integrator );
// close file
(void) fclose( filePtr );
}
libraries/validate/src/ValidateOpenMMForces.cpp deleted 100644 → 0
View file @ b20c20fe
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include <iostream>
#include <sstream>
#include "openmm/OpenMMException.h"
#include "ValidateOpenMMForces.h"
using namespace OpenMM;
ForceValidationResult::ForceValidationResult( const Context& context1, const Context& context2, StringUIntMap& forceNamesMap ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::ForceValidationResult";
// ---------------------------------------------------------------------------------------
_nansDetected = 0;
// calculate forces/energies
int types = State::Forces | State::Energy;
State state1 = context1.getState( types );
State state2 = context2.getState( types );
_potentialEnergies[0] = state1.getPotentialEnergy();
_potentialEnergies[1] = state2.getPotentialEnergy();
if( ValidateOpenMM::isNanOrInfinity( _potentialEnergies[0] ) || ValidateOpenMM::isNanOrInfinity( _potentialEnergies[1] ) ){
_nansDetected++;
}
_forces[0] = state1.getForces();
_forces[1] = state2.getForces();
_platforms[0] = context1.getPlatform().getName();
_platforms[1] = context2.getPlatform().getName();
for( StringUIntMapI ii = forceNamesMap.begin(); ii != forceNamesMap.end(); ii++ ){
_forceNames.push_back( (*ii).first );
}
_calculateNorms();
}
ForceValidationResult::~ForceValidationResult( ){
}
void ForceValidationResult::_calculateNorms(){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::_calculateNorms";
// ---------------------------------------------------------------------------------------
for( int jj = 0; jj < 2; jj++ ){
if( _norms[jj].size() < 1 ){
_calculateNormOfForceVector( jj );
_findStatsForDouble( _norms[jj], _normStatVectors[jj] );
}
}
}
/**---------------------------------------------------------------------------------------
Find stats for double array
@param array array
@param statVector vector of stats
@return 0
--------------------------------------------------------------------------------------- */
void ForceValidationResult::_findStatsForDouble( const std::vector<double>& array, std::vector<double>& statVector ) const {
// ---------------------------------------------------------------------------------------
//static const char* methodName = "ForceValidationResult::_findStatsForDouble";
// ---------------------------------------------------------------------------------------
statVector.resize( STAT_END );
double avgValue = 0.0;
double stdValue = 0.0;
double minValue = 1.0e+30;
double maxValue = -1.0e+30;
int minValueIndex = 0;
int maxValueIndex = 0;
for( unsigned int ii = 0; ii < array.size(); ii++ ){
double norm = array[ii];
avgValue += norm;
stdValue += norm*norm;
if( norm > maxValue ){
maxValue = norm;
maxValueIndex = ii;
}
if( norm < minValue ){
minValue = norm;
minValueIndex = ii;
}
}
double count = static_cast<double>(array.size());
double iCount = count > 0.0 ? 1.0/count : 0.0;
statVector[STAT_AVG] = avgValue*iCount;
statVector[STAT_STD] = stdValue - avgValue*avgValue*count;
if( count > 1.0 ){
statVector[STAT_STD] = std::sqrt( stdValue/( count - 1.0 ) );
}
statVector[STAT_MIN] = minValue;
statVector[STAT_ID1] = static_cast<double>(minValueIndex);
statVector[STAT_MAX] = maxValue;
statVector[STAT_ID2] = static_cast<double>(maxValueIndex);
statVector[STAT_CNT] = count;
return;
}
void ForceValidationResult::_calculateNormOfForceVector( int forceIndex ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::_calculateNormOfForceVector";
// ---------------------------------------------------------------------------------------
for( unsigned int ii = 0; ii < _forces[forceIndex].size(); ii++ ){
Vec3 f1 = _forces[forceIndex][ii];
_norms[forceIndex].push_back( std::sqrt( (f1[0]*f1[0]) + (f1[1]*f1[1]) + (f1[2]*f1[2]) ) );
if( ValidateOpenMM::isNanOrInfinity( f1[0] ) || ValidateOpenMM::isNanOrInfinity( f1[1] ) || ValidateOpenMM::isNanOrInfinity( f1[2] ) ){
_nansDetected++;
}
}
return;
}
double ForceValidationResult::getPotentialEnergy( int energyIndex ) const {
if( energyIndex == 0 || energyIndex == 1 ){
return _potentialEnergies[energyIndex];
} else {
char buffer[1024];
(void) sprintf( buffer, "getPotentialEnergy: energyIndex=%d is inconsistent.", energyIndex );
throw OpenMMException( buffer );
}
}
std::vector<Vec3> ForceValidationResult::getForces( int forceIndex ) const {
if( forceIndex == 0 || forceIndex == 1 ){
return _forces[forceIndex];
} else {
char buffer[1024];
(void) sprintf( buffer, "getForces: forceIndex=%d is inconsistent.", forceIndex );
throw OpenMMException( buffer );
}
}
std::vector<double> ForceValidationResult::getForceNorms( int forceIndex ) const {
if( forceIndex == 0 || forceIndex == 1 ){
return _norms[forceIndex];
} else {
char buffer[1024];
(void) sprintf( buffer, "getForceNorms: forceIndex=%d is inconsistent.", forceIndex );
throw OpenMMException( buffer );
}
}
int ForceValidationResult::nansDetected() const {
return _nansDetected;
}
void ForceValidationResult::registerInconsistentForceIndex( int index, int value ){
_inconsistentForceIndicies[index] = value;
}
void ForceValidationResult::clearInconsistentForceIndexList(){
_inconsistentForceIndicies.clear();
}
void ForceValidationResult::getInconsistentForceIndexList( std::vector<int>& inconsistentIndices ) const {
for( MapIntIntCI ii = _inconsistentForceIndicies.begin(); ii != _inconsistentForceIndicies.end(); ii++ ){
inconsistentIndices.push_back( (*ii).first );
}
}
int ForceValidationResult::getNumberOfInconsistentForceEntries() const {
return static_cast<int>(_inconsistentForceIndicies.size() );
}
double ForceValidationResult::getMaxDeltaForceNorm( int* maxIndex ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::calculateNorm";
// ---------------------------------------------------------------------------------------
double maxDelta = -1.0e+90;
int maxDeltaIndex = -1;
for( unsigned int ii = 0; ii < _forces[0].size(); ii++ ){
Vec3 f1 = _forces[0][ii];
double normF1 = _norms[0][ii];
Vec3 f2 = _forces[1][ii];
double normF2 = _norms[1][ii];
double delta = std::sqrt( (f1[0]-f2[0])*(f1[0]-f2[0]) + (f1[1]-f2[1])*(f1[1]-f2[1]) + (f1[2]-f2[2])*(f1[2]-f2[2]) );
if( maxDelta < delta ){
maxDelta = delta;
maxDeltaIndex = static_cast<int>(ii);
}
}
if( maxIndex ){
*maxIndex = maxDeltaIndex;
}
return maxDelta;
}
double ForceValidationResult::getMaxRelativeDeltaForceNorm( int* maxIndex ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::getMaxRelativeDeltaForceNorm";
// ---------------------------------------------------------------------------------------
double maxRelativeDelta = -1.0e+90;
int maxRelativeDeltaIndex = -1;
for( unsigned int ii = 0; ii < _forces[0].size(); ii++ ){
Vec3 f1 = _forces[0][ii];
double normF1 = _norms[0][ii];
Vec3 f2 = _forces[1][ii];
double normF2 = _norms[1][ii];
double delta = std::sqrt( (f1[0]-f2[0])*(f1[0]-f2[0]) + (f1[1]-f2[1])*(f1[1]-f2[1]) + (f1[2]-f2[2])*(f1[2]-f2[2]) );
double forceSum = 0.5*(normF1 + normF2);
if( forceSum > 0.0 ){
delta /= forceSum;
}
if( maxRelativeDelta < delta ){
maxRelativeDelta = delta;
maxRelativeDeltaIndex = static_cast<int>(ii);
}
}
if( maxIndex ){
*maxIndex = maxRelativeDeltaIndex;
}
return maxRelativeDelta;
}
double ForceValidationResult::getMaxDotProduct( int* maxIndex ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::getMaxDotProduct";
// ---------------------------------------------------------------------------------------
double maxDotProduct = -1.0e+90;
int maxDotProductIndex = -1;
for( unsigned int ii = 0; ii < _forces[0].size(); ii++ ){
Vec3 f1 = _forces[0][ii];
double normF1 = _norms[0][ii];
Vec3 f2 = _forces[1][ii];
double normF2 = _norms[1][ii];
double dotProduct = f1[0]*f2[0] + f1[1]*f2[1] + f1[2]*f2[2];
if( (normF1*normF2) > 0.0 ){
dotProduct /= (normF1*normF2);
dotProduct = 1.0 - dotProduct;
} else {
dotProduct = 0.0;
}
if( maxDotProduct < dotProduct ){
maxDotProduct = dotProduct;
maxDotProductIndex = static_cast<int>(ii);
}
}
if( maxIndex ){
*maxIndex = maxDotProductIndex;
}
return maxDotProduct;
}
std::string ForceValidationResult::getForceName() const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::getForceName";
// ---------------------------------------------------------------------------------------
std::stringstream forceName;
for( unsigned int ii = 0; ii < _forceNames.size(); ii++ ){
forceName << _forceNames[ii];
if( ii < (_forceNames.size()-1) )forceName << "::";
}
return forceName.str();
}
std::string ForceValidationResult::getPlatformName( int index ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ForceValidationResult::getPlatformName";
// ---------------------------------------------------------------------------------------
if( index == 0 || index == 1 ){
return _platforms[index];
} else {
char buffer[1024];
(void) sprintf( buffer, "getPlatformName: index=%d is inconsistent.", index );
throw OpenMMException( buffer );
}
}
void ForceValidationResult::compareForces( double tolerance ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = " ForceValidationResult::compareForces";
// ---------------------------------------------------------------------------------------
std::vector<Vec3> forces1 = getForces( 0 );
std::vector<double> forceNorms1 = getForceNorms( 0 );
std::vector<Vec3> forces2 = getForces( 1 );
std::vector<double> forceNorms2 = getForceNorms( 1 );
clearInconsistentForceIndexList();
for( unsigned int jj = 0; jj < forces1.size(); jj++ ){
if( ValidateOpenMM::isNanOrInfinity( forceNorms1[jj] ) || ValidateOpenMM::isNanOrInfinity( forceNorms2[jj] ) ){
registerInconsistentForceIndex( jj );
} else {
double delta = std::sqrt( (forces1[jj][0]-forces2[jj][0])*(forces1[jj][0]-forces2[jj][0]) +
(forces1[jj][1]-forces2[jj][1])*(forces1[jj][1]-forces2[jj][1]) +
(forces1[jj][2]-forces2[jj][2])*(forces1[jj][2]-forces2[jj][2]) );
double sum = 0.5*( fabs( forceNorms1[jj] ) + fabs( forceNorms2[jj] ) );
double diff = delta > 0.0 ? delta/sum : 0.0;
if( diff > tolerance && sum > tolerance ){
registerInconsistentForceIndex( jj );
}
}
}
}
void ForceValidationResult::compareForceNorms( double tolerance ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "ForceValidationResult::compareForceNorms";
// ---------------------------------------------------------------------------------------
std::vector<double> forceNorms0 = getForceNorms( 0 );
std::vector<double> forceNorms1 = getForceNorms( 1 );
clearInconsistentForceIndexList();
for( unsigned int jj = 0; jj < forceNorms0.size(); jj++ ){
if( ValidateOpenMM::isNanOrInfinity( forceNorms0[jj] ) || ValidateOpenMM::isNanOrInfinity( forceNorms1[jj] ) ){
registerInconsistentForceIndex( jj );
} else {
double sum = 0.5*( fabs( forceNorms0[jj] ) + fabs( forceNorms1[jj] ) );
double diff = sum > 0.0 ? fabs( forceNorms0[jj] - forceNorms1[jj] )/sum : 0.0;
if( diff > tolerance && sum > tolerance ){
registerInconsistentForceIndex( jj );
}
}
}
}
ValidateOpenMMForces::ValidateOpenMMForces() {
_initialize();
}
void ValidateOpenMMForces::_initialize(){
_forceTolerance = 1.0e-02;
_maxErrorsToPrint = 25;
// force tolerances by name
_forceTolerances[HARMONIC_BOND_FORCE] = _forceTolerance;
_forceTolerances[HARMONIC_ANGLE_FORCE] = _forceTolerance;
//_forceTolerances[PERIODIC_TORSION_FORCE] = _forceTolerance;
_forceTolerances[PERIODIC_TORSION_FORCE] = 0.3;
//_forceTolerances[RB_TORSION_FORCE] = _forceTolerance;
_forceTolerances[RB_TORSION_FORCE] = 0.3;
_forceTolerances[NB_FORCE] = _forceTolerance;
_forceTolerances[NB_SOFTCORE_FORCE] = _forceTolerance;
_forceTolerances[GBSA_OBC_FORCE] = _forceTolerance;
_forceTolerances[GBSA_OBC_SOFTCORE_FORCE] = _forceTolerance;
_forceTolerances[GBVI_FORCE] = _forceTolerance;
_forceTolerances[GBVI_SOFTCORE_FORCE] = _forceTolerance;
// forces to be excluded from validation
_forcesToBeExcluded[CM_MOTION_REMOVER] = 1;
_forcesToBeExcluded[ANDERSEN_THERMOSTAT] = 1;
}
ValidateOpenMMForces::~ValidateOpenMMForces(){
for( unsigned int ii = 0; ii < _forceValidationResults.size(); ii++ ){
delete _forceValidationResults[ii];
}
_forceValidationResults.resize( 0 );
}
double ValidateOpenMMForces::getForceTolerance() const {
return _forceTolerance;
}
void ValidateOpenMMForces::setForceTolerance( double forceTolerance ){
_forceTolerance = forceTolerance;
}
int ValidateOpenMMForces::compareWithReferencePlatform( Context& context, std::string* summaryString ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "ValidateOpenMMForces::compareWithReferencePlatform";
// ---------------------------------------------------------------------------------------
ReferencePlatform referencePlatform;
#ifdef INCLUDE_FREE_ENERGY_PLUGIN
ReferenceFreeEnergyKernelFactory* factory = new ReferenceFreeEnergyKernelFactory();
referencePlatform.registerKernelFactory(CalcNonbondedSoftcoreForceKernel::Name(), factory);
referencePlatform.registerKernelFactory(CalcGBSAOBCSoftcoreForceKernel::Name(), factory);
referencePlatform.registerKernelFactory(CalcGBVISoftcoreForceKernel::Name(), factory);
#endif
compareOpenMMForces( context, referencePlatform, _forceValidationResults );
checkForInconsistentForceEntries( _forceValidationResults );
if( summaryString ){
*summaryString = getSummary( _forceValidationResults );
}
return getTotalNumberOfInconsistentForceEntries( _forceValidationResults );
}
void ValidateOpenMMForces::compareOpenMMForces( Context& context, Platform& platform,
std::vector<ForceValidationResult*>& forceValidationResults ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "ValidateOpenMMForces::compareOpenMMForces";
// ---------------------------------------------------------------------------------------
// loop over forces
const System& system = context.getSystem();
std::vector<int> compareAllForces;
for( int ii = 0; ii < system.getNumForces(); ii++ ){
std::vector<int> compareForces;
std::string forceName = getForceName( system.getForce( ii ) );
if( isExcludedForce( forceName ) == 0 ){
compareForces.push_back( ii );
compareAllForces.push_back( ii );
ForceValidationResult* forceValidationResult = compareForce(context, compareForces, platform, context.getPlatform() );
forceValidationResults.push_back( forceValidationResult );
} else if( getLog() ){
(void) fprintf( getLog(), "Force %s is not being validated.\n", forceName.c_str() );
}
}
ForceValidationResult* forceValidationResult = compareForce(context, compareAllForces, platform, context.getPlatform() );
forceValidationResults.push_back( forceValidationResult );
}
ForceValidationResult* ValidateOpenMMForces::compareForce(Context& context, std::vector<int>& compareForces,
Platform& platform1, Platform& platform2 ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "ValidateOpenMMForces::compareForce";
// ---------------------------------------------------------------------------------------
// note if platforms are identical
if( getLog() && platform1.getName().compare( platform2.getName() ) == 0 ){
(void) fprintf( getLog(), "Note: Platforms to compares %s are identical.\n", platform1.getName().c_str() );
(void) fflush( getLog() );
}
const System& system = context.getSystem();
// collect systemForceNameMap[forceName] = index in system
// systemForceNameIndex[index] = force name
StringIntMap systemForceNameMap;
StringVector systemForceNameIndex;
systemForceNameIndex.resize( system.getNumForces() );
for( int ii = 0; ii < system.getNumForces(); ii++ ){
std::string forceName = getForceName( system.getForce( ii ) );
if( forceName.compare( "NA" ) == 0 ){
std::stringstream message;
message << "Force at index=" << ii << " not found -- aborting!";
std::cerr << message.str() << std::endl;
throw OpenMM::OpenMMException(message.str());
}
systemForceNameMap[forceName] = ii;
systemForceNameIndex[ii] = forceName;
}
// diagnostics
if( 0 && getLog() ){
for( StringIntMapI ii = systemForceNameMap.begin(); ii != systemForceNameMap.end(); ii++ ){
int index = (*ii).second;
(void) fprintf( getLog(), " System force map %s index=%d reverse map=%s\n", (*ii).first.c_str(), index, systemForceNameIndex[index].c_str() );
}
for( unsigned int ii = 0; ii < compareForces.size(); ii++ ){
(void) fprintf( getLog(), " ValidateOpenMMForces %u %s\n", ii, systemForceNameIndex[compareForces[ii]].c_str() );
}
(void) fflush( getLog() );
}
// get system copy and add forces to system
System* validationSystem = copySystemExcludingForces( system );
StringUIntMap forceNamesMap;
for( unsigned int ii = 0; ii < compareForces.size(); ii++ ){
const Force& forceToCopy = system.getForce( compareForces[ii] );
Force* force = copyForce( forceToCopy );
validationSystem->addForce( force );
forceNamesMap[systemForceNameIndex[compareForces[ii]]] = ii;
}
// include any missing dependencies (e.g, OBC force requires NB force for Cuda platform)
for( StringUIntMapI ii = forceNamesMap.begin(); ii != forceNamesMap.end(); ii++ ){
std::string forceName = (*ii).first;
StringVector dependencyVector;
getForceDependencies( forceName, dependencyVector );
for( unsigned int jj = 0; jj < dependencyVector.size(); jj++ ){
std::string dependentForceName = dependencyVector[jj];
StringUIntMapCI dependent = forceNamesMap.find( dependentForceName );
if( dependent == forceNamesMap.end() ){
forceNamesMap[dependentForceName] = 1;
int forceIndex = systemForceNameMap[dependentForceName];
const Force& forceToCopy = system.getForce( forceIndex );
validationSystem->addForce( copyForce( forceToCopy ) );
}
}
}
// create contexts
VerletIntegrator verletIntegrator( 0.001 );
Context* validationContext1 = new Context( *validationSystem, verletIntegrator, platform1);
Context* validationContext2 = new Context( *validationSystem, verletIntegrator, platform2);
// set positions
synchContexts( context, *validationContext1 );
synchContexts( context, *validationContext2 );
// diagnostics
if( 0 && getLog() ){
std::stringstream forceNames;
(void) fprintf( getLog(), " Validating system forces=%d\n", validationSystem->getNumForces() );
for( int ii = 0; ii < validationSystem->getNumForces(); ii++ ){
std::string forceName = getForceName( validationSystem->getForce( ii ) );
forceNames << forceName;
if( ii < (validationSystem->getNumForces()-1) ){
forceNames << "_";
} else {
forceNames << "Parameters.txt";
}
(void) fprintf( getLog(), " force %d %s\n", ii, forceName.c_str() );
}
writeParameterFile( *validationContext1, forceNames.str() );
(void) fflush( getLog() );
}
// calculate forces & build return result
ForceValidationResult* forceValidationResult = new ForceValidationResult( *validationContext1, *validationContext2, forceNamesMap );
delete validationContext1;
delete validationContext2;
delete validationSystem;
return forceValidationResult;
}
void ValidateOpenMMForces::checkForInconsistentForceEntries( std::vector<ForceValidationResult*>& forceValidationResults ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "ValidateOpenMMForces::checkForInconsistentForceEntries";
// ---------------------------------------------------------------------------------------
for( unsigned int ii = 0; ii < forceValidationResults.size(); ii++ ){
std::string forceName = forceValidationResults[ii]->getForceName();
double tolerance = getForceTolerance( forceName );
forceValidationResults[ii]->compareForces( tolerance );
}
}
int ValidateOpenMMForces::getTotalNumberOfInconsistentForceEntries( std::vector<ForceValidationResult*>& forceValidationResults ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "ValidateOpenMMForces::getTotalNumberOfInconsistentForceEntries";
// ---------------------------------------------------------------------------------------
int inconsistentEntries = 0;
for( unsigned int ii = 0; ii < forceValidationResults.size(); ii++ ){
inconsistentEntries += forceValidationResults[ii]->getNumberOfInconsistentForceEntries();
}
return inconsistentEntries;
}
int ValidateOpenMMForces::isExcludedForce( std::string forceName ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "ValidateOpenMMForces::isExcludedForce";
// ---------------------------------------------------------------------------------------
StringIntMapCI isPresent = _forcesToBeExcluded.find( forceName );
return isPresent != _forcesToBeExcluded.end() ? 1 : 0;
}
/*
* Get force tolerance for specified force
*
* @param forceName name of force
*
* @return force tolerance
*
* */
double ValidateOpenMMForces::getForceTolerance( const std::string& forceName ) const {
StringDoubleMapCI forceIsPresent = _forceTolerances.find( forceName );
if( forceIsPresent != _forceTolerances.end() ){
return (*forceIsPresent).second;
}
return _forceTolerance;
}
int ValidateOpenMMForces::getMaxErrorsToPrint() const {
return _maxErrorsToPrint;
}
void ValidateOpenMMForces::setMaxErrorsToPrint( int maxErrorsToPrint ){
_maxErrorsToPrint = maxErrorsToPrint;
}
/*
* Format line
*
* @param tab tab
* @param description description
* @param value value
*
* @return string containing contents
*
* */
std::string ValidateOpenMMForces::_getLine( const std::string& tab,
const std::string& description,
const std::string& value ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "GromacsOpenMMTest::_getLine";
static const unsigned int MAX_LINE_CHARS = 256;
char line[MAX_LINE_CHARS];
// ---------------------------------------------------------------------------------------
std::stringstream message;
memset( line, ' ', MAX_LINE_CHARS );
#ifdef WIN32
(void) sprintf_s( line, MAX_LINE_CHARS, "%s %-40s %s", tab.c_str(), description.c_str(), value.c_str() );
#else
(void) sprintf( line, "%s %-40s %s", tab.c_str(), description.c_str(), value.c_str() );
#endif
message << std::string( line ) << std::endl;
return message.str();
}
std::string ValidateOpenMMForces::getSummary( std::vector<ForceValidationResult*>& forceValidationResults ) const {
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "ValidateOpenMMForces::getSummary";
static const unsigned int MAX_LINE_CHARS = 256;
char value[MAX_LINE_CHARS];
// ---------------------------------------------------------------------------------------
#ifdef WIN32
#define LOCAL_SPRINTF0(a,b) sprintf_s( (a), MAX_LINE_CHARS, (b) );
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );
#define LOCAL_SPRINTF4(a,b,c,d) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d) );
#define LOCAL_SPRINTF5(a,b,c,d,e) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d), (e) );
#define LOCAL_SPRINTF6(a,b,c,d,e,f) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d), (e), (f) );
#define LOCAL_SPRINTF12(a,b,c,d,e,f,g,h,i,j,k,l) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d), (e), (f), (g), (h), (i), (j), (k), (l) );
#else
#define LOCAL_SPRINTF0(a,b) sprintf( (a), (b) );
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );
#define LOCAL_SPRINTF4(a,b,c,d) sprintf( (a), (b), (c), (d) );
#define LOCAL_SPRINTF5(a,b,c,d,e) sprintf( (a), (b), (c), (d), (e) );
#define LOCAL_SPRINTF6(a,b,c,d,e,f) sprintf( (a), (b), (c), (d), (e), (f) );
#define LOCAL_SPRINTF12(a,b,c,d,e,f,g,h,i,j,k,l) sprintf( (a), (b), (c), (d), (e), (f), (g), (h), (i), (j), (k), (l) );
#endif
unsigned int maxMissesToPrint = static_cast<unsigned int>(getMaxErrorsToPrint());
std::stringstream summary;
std::string tab = " ";
int useForces = 1;
for( unsigned int ii = 0; ii < forceValidationResults.size(); ii++ ){
// platforms
if( ii == 0 ){
std::string platform1 = forceValidationResults[ii]->getPlatformName( 0 );
std::string platform2 = forceValidationResults[ii]->getPlatformName( 1 );
(void) LOCAL_SPRINTF4( value, "%10s %10s\n", platform1.c_str(), platform2.c_str());
summary << _getLine( tab, "Platforms", value );
}
// force name
std::string forceName = forceValidationResults[ii]->getForceName();
double tolerance = getForceTolerance( forceName );
summary << _getLine( tab, "Force", forceName );
(void) LOCAL_SPRINTF( value, "%.3e ", tolerance );
summary << _getLine( tab, "Tolerance", value );
int index;
// delta
double delta = forceValidationResults[ii]->getMaxDeltaForceNorm( &index );
(void) LOCAL_SPRINTF4( value, "%.3e at index %6d", delta, index );
summary << _getLine( tab, "Max Delta", value );
// relative delta
delta = forceValidationResults[ii]->getMaxRelativeDeltaForceNorm( &index );
(void) LOCAL_SPRINTF4( value, "%.3e at index %6d", delta, index );
summary << _getLine( tab, "Max Relative Delta", value );
// PE delta
double pe1 = forceValidationResults[ii]->getPotentialEnergy( 0 );
double pe2 = forceValidationResults[ii]->getPotentialEnergy( 1 );
double sum = 0.5*(fabs( pe1 ) + fabs( pe2 ) );
delta = sum > 0.0 ? fabs( pe1 - pe2 )/sum : 0.0;
(void) LOCAL_SPRINTF5( value, "%10.4e PE[%14.6e %14.6e]", delta, pe1, pe2 );
summary << _getLine( tab, "Potential energies relative delta", value );
// misses
std::vector<int> inconsistentIndices;
forceValidationResults[ii]->getInconsistentForceIndexList( inconsistentIndices );
if( inconsistentIndices.size() > 0 ){
std::vector<Vec3> forces1 = forceValidationResults[ii]->getForces( 0 );
std::vector<Vec3> forces2 = forceValidationResults[ii]->getForces( 1 );
std::vector<double> forceNorms1 = forceValidationResults[ii]->getForceNorms( 0 );
std::vector<double> forceNorms2 = forceValidationResults[ii]->getForceNorms( 1 );
for( unsigned int kk = 0; kk < inconsistentIndices.size() && kk < maxMissesToPrint; kk++ ){
int jj = inconsistentIndices[kk];
if( isNanOrInfinity( forceNorms1[jj] ) || isNanOrInfinity( forceNorms2[jj] ) ){
(void) LOCAL_SPRINTF5( value, " nan at index %6d norms: [%12.5e %12.5e]", jj, forceNorms1[jj], forceNorms2[jj] );
summary << _getLine( tab, "Error", value );
} else {
double diff, sum;
if( useForces ){
double delta = std::sqrt( (forces1[jj][0]-forces2[jj][0])*(forces1[jj][0]-forces2[jj][0]) +
(forces1[jj][1]-forces2[jj][1])*(forces1[jj][1]-forces2[jj][1]) +
(forces1[jj][2]-forces2[jj][2])*(forces1[jj][2]-forces2[jj][2]) );
sum = 0.5*( fabs( forceNorms1[jj] ) + fabs( forceNorms2[jj] ) );
diff = delta > 0.0 ? delta/sum : 0.0;
} else {
sum = 0.5*( fabs( forceNorms1[jj] ) + fabs( forceNorms2[jj] ) );
diff = sum > 0.0 ? fabs( forceNorms1[jj] - forceNorms2[jj] )/sum : 0.0;
}
(void) LOCAL_SPRINTF12( value, "%11.5e at index %6d norms: [%11.5e %11.5e] forces: [%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]",
diff, jj, forceNorms1[jj], forceNorms2[jj],
forces1[jj][0], forces1[jj][1], forces1[jj][2],
forces2[jj][0], forces2[jj][1], forces2[jj][2] );
summary << _getLine( tab, "Error", value );
if( kk == (maxMissesToPrint-1) ){
value[0] = '\0';
summary << _getLine( tab, "No more errors will be printed: call ValidateOpenMMForces::setMaxErrorsToPrint() to modifiy this limit.", value );
}
}
}
}
(void) LOCAL_SPRINTF0( value, "\n" );
summary << _getLine( tab, " ", value );
if( inconsistentIndices.size() > 0 ){
(void) LOCAL_SPRINTF( value, "%u\n", static_cast<unsigned int>(inconsistentIndices.size()) );
summary << _getLine( tab, "Total errors", value );
}
if( forceValidationResults[ii]->nansDetected() ){
value[0] = '\0';
summary << _getLine( tab, "Nans detected.", value );
}
}
return summary.str();
}
olla/include/openmm/kernels.h
View file @ 167ae8a0
......@@ -48,7 +48,6 @@
#include "openmm/CustomManyParticleForce.h"
#include "openmm/CustomTorsionForce.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/GBVIForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/KernelImpl.h"
......@@ -667,35 +666,6 @@ public:
virtual void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) = 0;
};
/**
* This kernel is invoked by GBVIForce to calculate the forces acting on the system and the energy of the system.
*/
class CalcGBVIForceKernel : public KernelImpl {
public:
static std::string Name() {
return "CalcGBVIForce";
}
CalcGBVIForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the GBVIForce this kernel will be used for
* @param scaledRadii scaled radii
*/
virtual void initialize(const System& system, const GBVIForce& force, const std::vector<double>& scaledRadii) = 0;
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
};
/**
* This kernel is invoked by CustomGBForce to calculate the forces acting on the system and the energy of the system.
*/
......
openmmapi/include/OpenMM.h
View file @ 167ae8a0
......@@ -50,7 +50,6 @@
#include "openmm/CustomNonbondedForce.h"
#include "openmm/Force.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/GBVIForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/Integrator.h"
......
openmmapi/include/openmm/GBVIForce.h deleted 100644 → 0
View file @ b20c20fe
#ifndef OPENMM_GBVIFORCEFIELD_H_
#define OPENMM_GBVIFORCEFIELD_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* 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 "Force.h"
#include <vector>
#include "internal/windowsExport.h"
namespace OpenMM {
/**
* This class implements an implicit solvation force using the GB/VI model.
* <p>
* To use this class, create a GBVIForce object, then call addParticle() once for each particle in the
* System to define its parameters. The number of particles for which you define GB/VI parameters must
* be exactly equal to the number of particles in the System, or else an exception will be thrown when you
* try to create a Context. After a particle has been added, you can modify its force field parameters
* by calling setParticleParameters().
*
* @deprecated This class is not supported by most platforms, and will eventually be removed. You can implement the same force with CustomGBForce.
*/
class OPENMM_EXPORT GBVIForce : public Force {
public:
/**
* This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
*/
enum NonbondedMethod {
/**
* No cutoff is applied to nonbonded interactions. The full set of N^2 interactions is computed exactly.
* This necessarily means that periodic boundary conditions cannot be used. This is the default.
*/
NoCutoff = 0,
/**
* Interactions beyond the cutoff distance are ignored.
*/
CutoffNonPeriodic = 1,
/**
* Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
* each other particle. Interactions beyond the cutoff distance are ignored.
*/
CutoffPeriodic = 2,
};
/**
* This is an enumeration of the different methods that may be used for scaling of the Born radii.
*/
enum BornRadiusScalingMethod {
/**
* No scaling method is applied.
*/
NoScaling = 0,
/**
* Use quintic spline scaling function
*/
QuinticSpline = 1
};
/*
* Create a GBVIForce.
*/
GBVIForce();
/**
* Get the number of particles in the system.
*/
int getNumParticles() const {
return particles.size();
}
/**
* Add the GB/VI parameters for a particle. This should be called once for each particle
* in the System. When it is called for the i'th time, it specifies the parameters for the i'th particle.
*
* @param charge the charge of the particle, measured in units of the proton charge
* @param radius the GB/VI radius of the particle, measured in nm
* @param gamma the gamma parameter
* @return the index of the particle that was added
*/
int addParticle(double charge, double radius, double gamma);
/**
* Get the force field parameters for a particle.
*
* @param index the index of the particle for which to get parameters
* @param[out] charge the charge of the particle, measured in units of the proton charge
* @param[out] radius the GBSA radius of the particle, measured in nm
* @param[out] gamma the gamma parameter
*/
void getParticleParameters(int index, double& charge, double& radius, double& gamma) const;
/**
* Set the force field parameters for a particle.
*
* @param index the index of the particle for which to set parameters
* @param charge the charge of the particle, measured in units of the proton charge
* @param radius the GB/VI radius of the particle, measured in nm
* @param gamma the gamma parameter
*/
void setParticleParameters(int index, double charge, double radius, double gamma);
/**
* Add a bond
*
* @param particle1 the index of the first particle
* @param particle2 the index of the second particle
* @param distance the distance between the two particles, measured in nm
* @return the index of the bond that was added
*/
int addBond(int particle1, int particle2, double distance);
/**
* Get the parameters defining a bond
*
* @param index the index of the bond for which to get parameters
* @param[out] particle1 the index of the first particle involved in the bond
* @param[out] particle2 the index of the second particle involved in the bond
* @param[out] distance the distance between the two particles, measured in nm
*/
void getBondParameters(int index, int& particle1, int& particle2, double& distance) const;
/**
* Set 1-2 bonds
*
* @param index index of the bond for which to set parameters
* @param particle1 index of first atom in bond
* @param particle2 index of second atom in bond
* @param bondLength bond length, measured in nm
*/
void setBondParameters( int index, int particle1, int particle2, double bondLength);
/**
* Get number of bonds
*
* @return number of bonds
*/
int getNumBonds() const;
/**
* Get the dielectric constant for the solvent.
*/
double getSolventDielectric() const {
return solventDielectric;
}
/**
* Set the dielectric constant for the solvent.
*/
void setSolventDielectric(double dielectric) {
solventDielectric = dielectric;
}
/**
* Get the dielectric constant for the solute.
*/
double getSoluteDielectric() const {
return soluteDielectric;
}
/**
* Set the dielectric constant for the solute.
*/
void setSoluteDielectric(double dielectric) {
soluteDielectric = dielectric;
}
/**
* Get the method used for handling long range nonbonded interactions.
*/
NonbondedMethod getNonbondedMethod() const;
/**
* Set the method used for handling long range nonbonded interactions.
*/
void setNonbondedMethod(NonbondedMethod method);
/**
* Get the cutoff distance (in nm) being used for nonbonded interactions. If the NonbondedMethod in use
* is NoCutoff, this value will have no effect.
*
* @return the cutoff distance, measured in nm
*/
double getCutoffDistance() const;
/**
* Set the cutoff distance (in nm) being used for nonbonded interactions. If the NonbondedMethod in use
* is NoCutoff, this value will have no effect.
*
* @param distance the cutoff distance, measured in nm
*/
void setCutoffDistance(double distance);
/**
* Get Born radius scaling method
*/
BornRadiusScalingMethod getBornRadiusScalingMethod() const;
/**
* Set Born radius scaling method
*/
void setBornRadiusScalingMethod( BornRadiusScalingMethod method);
/**
* Get the lower limit factor used in the quintic spline scaling method (typically 0.5-0.8)
*/
double getQuinticLowerLimitFactor() const;
/**
* Set the lower limit factor used in the quintic spline scaling method (typically 0.5-0.8)
*/
void setQuinticLowerLimitFactor(double quinticLowerLimitFactor );
/**
* Get the upper limit used in the quintic spline scaling method, measured in nm (~5.0)
*/
double getQuinticUpperBornRadiusLimit() const;
/**
* Set the upper limit used in the quintic spline scaling method, measured in nm (~5.0)
*/
void setQuinticUpperBornRadiusLimit(double quinticUpperBornRadiusLimit);
/**
* Returns whether or not this force makes use of periodic boundary
* conditions.
*
* @returns true if force uses PBC and false otherwise
*/
bool usesPeriodicBoundaryConditions() const {
return nonbondedMethod == GBVIForce::CutoffPeriodic;
}
protected:
ForceImpl* createImpl() const;
private:
class ParticleInfo;
NonbondedMethod nonbondedMethod;
double cutoffDistance, solventDielectric, soluteDielectric;
BornRadiusScalingMethod scalingMethod;
double quinticLowerLimitFactor, quinticUpperBornRadiusLimit;
class BondInfo;
std::vector<ParticleInfo> particles;
std::vector<BondInfo> bonds;
};
/**
* This is an internal class used to record information about a particle.
* @private
*/
class GBVIForce::ParticleInfo {
public:
double charge, radius, gamma;
ParticleInfo() {
charge = radius = gamma = 0.0;
}
ParticleInfo(double charge, double radius, double gamma) :
charge(charge), radius(radius), gamma(gamma) {
}
};
/**
* This is an internal class used to record information about a bond.
* @private
*/
class GBVIForce::BondInfo {
public:
int particle1, particle2;
double bondLength;
BondInfo() {
bondLength = 0.0;
particle1 = -1;
particle2 = -1;
}
BondInfo(int atomIndex1, int atomIndex2, double bondLength) :
particle1(atomIndex1), particle2(atomIndex2), bondLength(bondLength) {
}
};
} // namespace OpenMM
#endif /*OPENMM_GBVIFORCEFIELD_H_*/
openmmapi/include/openmm/internal/GBVIForceImpl.h deleted 100644 → 0
View file @ b20c20fe
#ifndef OPENMM_GBVIFORCEFIELDIMPL_H_
#define OPENMM_GBVIFORCEFIELDIMPL_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* 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 "ForceImpl.h"
#include "openmm/GBVIForce.h"
#include "openmm/Kernel.h"
#include <string>
namespace OpenMM {
/**
* This is the internal implementation of GBVIForce.
*/
class GBVIForceImpl : public ForceImpl {
public:
GBVIForceImpl(const GBVIForce& owner);
void initialize(ContextImpl& context);
const GBVIForce& getOwner() const {
return owner;
}
// calculate scaled radii (Eq. 5 of Labute paper [JCC 29 1693-1698 2008])
void findScaledRadii( int numberOfParticles, const std::vector<std::vector<int> >& bondIndices,
const std::vector<double> & bondLengths, std::vector<double> & scaledRadii) const;
// if bond info not set, then use bond forces/constraints
int getBondsFromForces(ContextImpl& context);
void updateContextState(ContextImpl& context) {
// This force field doesn't update the state directly.
}
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
std::map<std::string, double> getDefaultParameters() {
return std::map<std::string, double>(); // This force field doesn't define any parameters.
}
std::vector<std::string> getKernelNames();
private:
const GBVIForce& owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_GBVIFORCEFIELDIMPL_H_*/
openmmapi/src/GBVIForce.cpp deleted 100644 → 0
View file @ b20c20fe
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* 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 "openmm/Force.h"
#include "openmm/OpenMMException.h"
#include "openmm/GBVIForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/internal/GBVIForceImpl.h"
#include <sstream>
using namespace OpenMM;
GBVIForce::GBVIForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), solventDielectric(78.3), soluteDielectric(1.0),
scalingMethod(QuinticSpline), quinticLowerLimitFactor(0.8), quinticUpperBornRadiusLimit(5.0) {
}
int GBVIForce::addParticle(double charge, double radius, double gamma) {
particles.push_back(ParticleInfo(charge, radius, gamma));
return particles.size()-1;
}
void GBVIForce::getParticleParameters(int index, double& charge, double& radius, double& gamma) const {
ASSERT_VALID_INDEX(index, particles);
charge = particles[index].charge;
radius = particles[index].radius;
gamma = particles[index].gamma;
}
void GBVIForce::setParticleParameters(int index, double charge, double radius, double gamma) {
ASSERT_VALID_INDEX(index, particles);
particles[index].charge = charge;
particles[index].radius = radius;
particles[index].gamma = gamma;
}
GBVIForce::NonbondedMethod GBVIForce::getNonbondedMethod() const {
return nonbondedMethod;
}
void GBVIForce::setNonbondedMethod(NonbondedMethod method) {
nonbondedMethod = method;
}
double GBVIForce::getCutoffDistance() const {
return cutoffDistance;
}
void GBVIForce::setCutoffDistance(double distance) {
cutoffDistance = distance;
}
GBVIForce::BornRadiusScalingMethod GBVIForce::getBornRadiusScalingMethod() const {
return scalingMethod;
}
void GBVIForce::setBornRadiusScalingMethod(BornRadiusScalingMethod method) {
scalingMethod = method;
}
double GBVIForce::getQuinticLowerLimitFactor() const {
return quinticLowerLimitFactor;
}
void GBVIForce::setQuinticLowerLimitFactor(double inputQuinticLowerLimitFactor ){
quinticLowerLimitFactor = inputQuinticLowerLimitFactor;
}
double GBVIForce::getQuinticUpperBornRadiusLimit() const {
return quinticUpperBornRadiusLimit;
}
void GBVIForce::setQuinticUpperBornRadiusLimit(double inputQuinticUpperBornRadiusLimit){
quinticUpperBornRadiusLimit = inputQuinticUpperBornRadiusLimit;
}
int GBVIForce::addBond(int particle1, int particle2, double bondLength) {
bonds.push_back(BondInfo(particle1, particle2, bondLength));
return bonds.size()-1;
}
void GBVIForce::setBondParameters( int index, int particle1, int particle2, double bondLength) {
ASSERT_VALID_INDEX(index, bonds);
bonds[index].particle1 = particle1;
bonds[index].particle2 = particle2;
bonds[index].bondLength = bondLength;
}
int GBVIForce::getNumBonds() const {
return (int) bonds.size();
}
void GBVIForce::getBondParameters(int index, int& bondIndex1, int& bondIndex2, double& bondLength) const {
ASSERT_VALID_INDEX(index, bonds);
bondIndex1 = bonds[index].particle1;
bondIndex2 = bonds[index].particle2;
bondLength = bonds[index].bondLength;
}
ForceImpl* GBVIForce::createImpl() const {
return new GBVIForceImpl(*this);
}
openmmapi/src/GBVIForceImpl.cpp deleted 100644 → 0
View file @ b20c20fe
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* 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 "openmm/internal/GBVIForceImpl.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/OpenMMException.h"
#include "openmm/kernels.h"
#include <vector>
#include <cmath>
#include <cstdio>
#include <sstream>
using namespace OpenMM;
using std::vector;
GBVIForceImpl::GBVIForceImpl(const GBVIForce& owner) : owner(owner) {
}
void GBVIForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcGBVIForceKernel::Name(), context);
if (owner.getNumParticles() != context.getSystem().getNumParticles())
throw OpenMMException("GBVIForce must have exactly as many particles as the System it belongs to.");
const System& system = context.getSystem();
int numberOfParticles = owner.getNumParticles();
int numberOfBonds = owner.getNumBonds();
// load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
// numberOfBonds < 1, indicating they were not set by the user
if( numberOfBonds < 1 && numberOfParticles > 1 ){
(void) fprintf( stderr, "Warning: no covalent bonds set for GB/VI force!\n" );
// getBondsFromForces( context );
// numberOfBonds = owner.getNumBonds();
}
std::vector< std::vector<int> > bondIndices;
bondIndices.resize( numberOfBonds );
std::vector<double> bondLengths;
bondLengths.resize( numberOfBonds );
for (int i = 0; i < numberOfBonds; i++) {
int particle1, particle2;
double bondLength;
owner.getBondParameters(i, particle1, particle2, bondLength);
if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
std::stringstream msg;
msg << "GBVISoftcoreForce: Illegal particle index: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
std::stringstream msg;
msg << "GBVISoftcoreForce: Illegal particle index: ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (bondLength < 0 ) {
std::stringstream msg;
msg << "GBVISoftcoreForce: negative bondlength: ";
msg << bondLength;
throw OpenMMException(msg.str());
}
bondIndices[i].push_back( particle1 );
bondIndices[i].push_back( particle2 );
bondLengths[i] = bondLength;
}
if (owner.getNonbondedMethod() == GBVIForce::CutoffPeriodic) {
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
double cutoff = owner.getCutoffDistance();
if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
throw OpenMMException("GBVIForce: The cutoff distance cannot be greater than half the periodic box size.");
}
vector<double> scaledRadii;
scaledRadii.resize(numberOfParticles);
findScaledRadii( numberOfParticles, bondIndices, bondLengths, scaledRadii);
kernel.getAs<CalcGBVIForceKernel>().initialize(context.getSystem(), owner, scaledRadii);
}
/*
int GBVIForceImpl::getBondsFromForces(ContextImpl& context) {
// load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
const System& system = context.getSystem();
for (int i = 0; i < system.getNumForces(); i++) {
if (dynamic_cast<const HarmonicBondForce*>(&system.getForce(i)) != NULL) {
const HarmonicBondForce& force = dynamic_cast<const HarmonicBondForce&>(system.getForce(i));
for (int j = 0; j < force.getNumBonds(); ++j) {
int particle1, particle2;
double length, k;
force.getBondParameters(j, particle1, particle2, length, k);
owner.addBond( particle1, particle2, length );
}
break;
}
}
// Also treat constrained distances as bonds if mass of one particle is < (2 + epsilon) (~2=deuterium)
for (int j = 0; j < system.getNumConstraints(); j++) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
double mass1 = system.getParticleMass( particle1 );
double mass2 = system.getParticleMass( particle2 );
if( mass1 < 2.1 || mass2 < 2.1 ){
owner.addBond( particle1, particle2, distance );
}
}
return 0;
}
*/
void GBVIForceImpl::findScaledRadii( int numberOfParticles, const std::vector<std::vector<int> >& bondIndices,
const std::vector<double> & bondLengths, std::vector<double> & scaledRadii) const {
// load 1-2 indicies for each atom
std::vector<std::vector<int> > bonded12(numberOfParticles);
for (int i = 0; i < (int) bondIndices.size(); ++i) {
bonded12[bondIndices[i][0]].push_back(i);
bonded12[bondIndices[i][1]].push_back(i);
}
int errors = 0;
// compute scaled radii (Eq. 5 of Labute paper [JCC 29 p. 1693-1698 2008])
for (int j = 0; j < (int) bonded12.size(); ++j){
double charge;
double gamma;
double radiusJ;
double scaledRadiusJ;
owner.getParticleParameters(j, charge, radiusJ, gamma);
if( bonded12[j].size() == 0 && numberOfParticles > 1 ){
(void) fprintf( stderr, "Warning GBVIForceImpl::findScaledRadii atom %d has no covalent bonds; using atomic radius=%.3f.\n", j, radiusJ );
scaledRadiusJ = radiusJ;
// errors++;
} else {
double rJ2 = radiusJ*radiusJ;
// loop over bonded neighbors of atom j, applying Eq. 5 in Labute
scaledRadiusJ = 0.0;
for (int i = 0; i < (int) bonded12[j].size(); ++i){
int index = bonded12[j][i];
int bondedAtomIndex = (j == bondIndices[index][0]) ? bondIndices[index][1] : bondIndices[index][0];
double radiusI;
owner.getParticleParameters(bondedAtomIndex, charge, radiusI, gamma);
double rI2 = radiusI*radiusI;
double a_ij = (radiusI - bondLengths[index]);
a_ij *= a_ij;
a_ij = (rJ2 - a_ij)/(2.0*bondLengths[index]);
double a_ji = radiusJ - bondLengths[index];
a_ji *= a_ji;
a_ji = (rI2 - a_ji)/(2.0*bondLengths[index]);
scaledRadiusJ += a_ij*a_ij*(3.0*radiusI - a_ij) + a_ji*a_ji*( 3.0*radiusJ - a_ji );
}
scaledRadiusJ = (radiusJ*radiusJ*radiusJ) - 0.125*scaledRadiusJ;
if( scaledRadiusJ > 0.0 ){
scaledRadiusJ = 0.95*pow( scaledRadiusJ, (1.0/3.0) );
} else {
scaledRadiusJ = 0.0;
}
}
scaledRadii[j] = scaledRadiusJ;
}
// abort if errors
if( errors ){
throw OpenMMException("GBVIForceImpl::findScaledRadii errors -- aborting");
}
}
double GBVIForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
if ((groups&(1<<owner.getForceGroup())) != 0)
return kernel.getAs<CalcGBVIForceKernel>().execute(context, includeForces, includeEnergy);
return 0.0;
}
std::vector<std::string> GBVIForceImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(CalcGBVIForceKernel::Name());
return names;
}
platforms/reference/include/GBVIParameters.h deleted 100644 → 0
View file @ b20c20fe
/* Portions copyright (c) 2006-2009 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 __GBVIParameters_H__
#define __GBVIParameters_H__
#include "RealVec.h"
#include <vector>
namespace OpenMM {
class GBVIParameters {
public:
/**
* This is an enumeration of the different methods that may be used for scaling of the Born radii.
*/
enum BornRadiusScalingMethod {
/**
* No scaling method is applied.
*/
NoScaling = 0,
/**
* Use quintic spline scaling function
*/
QuinticSpline = 1
};
private:
int _numberOfAtoms;
RealOpenMM _solventDielectric;
RealOpenMM _soluteDielectric;
RealOpenMM _electricConstant;
// parameter vectors
std::vector<RealOpenMM> _atomicRadii;
std::vector<RealOpenMM> _scaledRadii;
std::vector<RealOpenMM> _gammaParameters;
// cutoff and periodic boundary conditions
bool _cutoff;
bool _periodic;
OpenMM::RealVec _periodicBoxVectors[3];
RealOpenMM _cutoffDistance;
int _bornRadiusScalingMethod;
RealOpenMM _quinticLowerLimitFactor;
RealOpenMM _quinticUpperBornRadiusLimit;
public:
/**---------------------------------------------------------------------------------------
GBVIParameters constructor (Simbios)
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
GBVIParameters(int numberOfAtoms);
/**---------------------------------------------------------------------------------------
GBVIParameters destructor (Simbios)
--------------------------------------------------------------------------------------- */
~GBVIParameters();
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int getNumberOfAtoms() const;
/**---------------------------------------------------------------------------------------
Get electric constant
@return electric constant
--------------------------------------------------------------------------------------- */
RealOpenMM getElectricConstant() const;
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSolventDielectric() const;
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void setSolventDielectric(RealOpenMM solventDielectric);
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getSoluteDielectric() const;
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void setSoluteDielectric(RealOpenMM soluteDielectric);
/**---------------------------------------------------------------------------------------
Return scaled radii
@return array
--------------------------------------------------------------------------------------- */
const std::vector<RealOpenMM>& getScaledRadii() const;
/**---------------------------------------------------------------------------------------
Return scaled radii
@return array
--------------------------------------------------------------------------------------- */
void setScaledRadii(const std::vector<RealOpenMM>& scaledRadii);
/**---------------------------------------------------------------------------------------
Get AtomicRadii array w/ dielectric offset applied
@return array of atom volumes
--------------------------------------------------------------------------------------- */
const std::vector<RealOpenMM>& getAtomicRadii() const;
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
void setAtomicRadii(const std::vector<RealOpenMM>& atomicRadii);
/**---------------------------------------------------------------------------------------
Get GammaParameters array
@return array of gamma values
--------------------------------------------------------------------------------------- */
const std::vector<RealOpenMM>& getGammaParameters() const;
/**---------------------------------------------------------------------------------------
Set GammaParameters array
@param gammaParameters array of gamma parameters
--------------------------------------------------------------------------------------- */
void setGammaParameters(const std::vector<RealOpenMM>& gammaParameters);
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
--------------------------------------------------------------------------------------- */
void setUseCutoff(RealOpenMM distance);
/**---------------------------------------------------------------------------------------
Get whether to use a cutoff.
--------------------------------------------------------------------------------------- */
bool getUseCutoff();
/**---------------------------------------------------------------------------------------
Get the cutoff distance.
--------------------------------------------------------------------------------------- */
RealOpenMM getCutoffDistance();
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
Get whether to use periodic boundary conditions.
--------------------------------------------------------------------------------------- */
bool getPeriodic();
/**---------------------------------------------------------------------------------------
Get the periodic box dimension
--------------------------------------------------------------------------------------- */
const OpenMM::RealVec* getPeriodicBox();
/**---------------------------------------------------------------------------------------
Get tau prefactor
@return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM getTau() const;
/**---------------------------------------------------------------------------------------
Get bornRadiusScalingMethod
@return bornRadiusScalingMethod
--------------------------------------------------------------------------------------- */
int getBornRadiusScalingMethod() const;
/**---------------------------------------------------------------------------------------
Set bornRadiusScalingMethod
@param bornRadiusScalingMethod bornRadiusScalingMethod
--------------------------------------------------------------------------------------- */
void setBornRadiusScalingMethod(int bornRadiusScalingMethod);
/**---------------------------------------------------------------------------------------
Get quinticLowerLimitFactor
@return quinticLowerLimitFactor
--------------------------------------------------------------------------------------- */
RealOpenMM getQuinticLowerLimitFactor() const;
/**---------------------------------------------------------------------------------------
Set quinticLowerLimitFactor
@param quinticLowerLimitFactor quinticLowerLimitFactor
--------------------------------------------------------------------------------------- */
void setQuinticLowerLimitFactor(RealOpenMM quinticLowerLimitFactor);
/**---------------------------------------------------------------------------------------
Get quinticUpperBornRadiusLimit
@return quinticUpperBornRadiusLimit
--------------------------------------------------------------------------------------- */
RealOpenMM getQuinticUpperBornRadiusLimit() const;
/**---------------------------------------------------------------------------------------
Set quinticUpperBornRadiusLimit
@param quinticUpperBornRadiusLimit quinticUpperBornRadiusLimit
--------------------------------------------------------------------------------------- */
void setQuinticUpperBornRadiusLimit(RealOpenMM quinticUpperSplineLimit);
};
} // namespace OpenMM
#endif // __GBVIParameters_H__
platforms/reference/include/ReferenceGBVI.h deleted 100644 → 0
View file @ b20c20fe
/* 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 __ReferenceGBVI_H__
#define __ReferenceGBVI_H__
#include <vector>
#include "RealVec.h"
#include "GBVIParameters.h"
namespace OpenMM {
class ReferenceGBVI {
private:
// GB/VI parameters
GBVIParameters* _gbviParameters;
std::vector<RealOpenMM> _switchDeriviative;
public:
/**---------------------------------------------------------------------------------------
Constructor
@param implicitSolventParameters ImplicitSolventParameters reference
@return CpuImplicitSolvent object
--------------------------------------------------------------------------------------- */
ReferenceGBVI(GBVIParameters* gbviParameters);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceGBVI();
/**---------------------------------------------------------------------------------------
Return GBVIParameters
@return GBVIParameters
--------------------------------------------------------------------------------------- */
GBVIParameters* getGBVIParameters() const;
/**---------------------------------------------------------------------------------------
Set ImplicitSolventParameters
@param ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
void setGBVIParameters(GBVIParameters* gbviParameters);
/**---------------------------------------------------------------------------------------
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param gbviChain not used
--------------------------------------------------------------------------------------- */
void computeBornRadii(const std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<RealOpenMM>& bornRadii);
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
static RealOpenMM getVolume(RealOpenMM r, RealOpenMM R, RealOpenMM S);
/**---------------------------------------------------------------------------------------
Get L (analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM getL(RealOpenMM r, RealOpenMM x, RealOpenMM S);
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM dL_dr(RealOpenMM r, RealOpenMM x, RealOpenMM S);
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM dL_dx(RealOpenMM r, RealOpenMM x, RealOpenMM S);
/**---------------------------------------------------------------------------------------
Sgb function
@param t r*r*G_i*G_j
@return Sgb (p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM Sgb(RealOpenMM t);
/**---------------------------------------------------------------------------------------
Get GB/VI energy
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM computeBornEnergy(const std::vector<OpenMM::RealVec>& atomCoordinates, const std::vector<RealOpenMM>& partialCharges);
/**---------------------------------------------------------------------------------------
Get GB/VI forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces output forces
--------------------------------------------------------------------------------------- */
void computeBornForces(std::vector<OpenMM::RealVec>& atomCoordinates,
const std::vector<RealOpenMM>& partialCharges, std::vector<OpenMM::RealVec>& inputForces);
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
static double getVolumeD(double r, double R, double S);
/**---------------------------------------------------------------------------------------
Get L (analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static double getLD(double r, double x, double S);
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static double dL_drD(double r, double x, double S);
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static double dL_dxD(double r, double x, double S);
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
std::vector<RealOpenMM>& getSwitchDeriviative();
/**---------------------------------------------------------------------------------------
Compute quintic spline value and associated derviative
@param x value to compute spline at
@param rl lower cutoff value
@param ru upper cutoff value
@param outValue value of spline at x
@param outDerivative value of derivative of spline at x
--------------------------------------------------------------------------------------- */
void quinticSpline(RealOpenMM x, RealOpenMM rl, RealOpenMM ru,
RealOpenMM* outValue, RealOpenMM* outDerivative);
/**---------------------------------------------------------------------------------------
Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
and quintic splice switching function
@param atomicRadius3 atomic radius cubed
@param bornSum Born sum (volume integral)
@param gbviParameters Gbvi parameters (parameters used in spline
QuinticLowerLimitFactor & QuinticUpperBornRadiusLimit)
@param bornRadius output Born radius
@param switchDeriviative output switching function deriviative
--------------------------------------------------------------------------------------- */
void computeBornRadiiUsingQuinticSpline(RealOpenMM atomicRadius3, RealOpenMM bornSum,
GBVIParameters* gbviParameters,
RealOpenMM* bornRadius, RealOpenMM* switchDeriviative);
};
} // namespace OpenMM
#endif // __ReferenceGBVI_H__
platforms/reference/include/ReferenceKernels.h
View file @ 167ae8a0
......@@ -43,7 +43,6 @@
namespace OpenMM {
class ReferenceObc;
class ReferenceGBVI;
class ReferenceAndersenThermostat;
class ReferenceCustomCentroidBondIxn;
class ReferenceCustomCompoundBondIxn;
......@@ -683,37 +682,6 @@ private:
bool isPeriodic;
};
/**
* This kernel is invoked by GBVIForce to calculate the forces acting on the system.
*/
class ReferenceCalcGBVIForceKernel : public CalcGBVIForceKernel {
public:
ReferenceCalcGBVIForceKernel(std::string name, const Platform& platform) : CalcGBVIForceKernel(name, platform) {
}
~ReferenceCalcGBVIForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the GBVIForce this kernel will be used for
* @param scaled radii the scaled radii (Eq. 5 of Labute paper)
*/
void initialize(const System& system, const GBVIForce& force, const std::vector<double> & scaledRadii);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
ReferenceGBVI * gbvi;
std::vector<RealOpenMM> charges;
bool isPeriodic;
};
/**
* This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
*/
......
platforms/reference/src/ReferenceKernelFactory.cpp
View file @ 167ae8a0
......@@ -70,8 +70,6 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcCustomTorsionForceKernel(name, platform);
if (name == CalcGBSAOBCForceKernel::Name())
return new ReferenceCalcGBSAOBCForceKernel(name, platform);
if (name == CalcGBVIForceKernel::Name())
return new ReferenceCalcGBVIForceKernel(name, platform);
if (name == CalcCustomGBForceKernel::Name())
return new ReferenceCalcCustomGBForceKernel(name, platform);
if (name == CalcCustomExternalForceKernel::Name())
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 167ae8a0
......@@ -31,7 +31,6 @@
#include "ReferenceKernels.h"
#include "ReferenceObc.h"
#include "ReferenceGBVI.h"
#include "ReferenceAndersenThermostat.h"
#include "ReferenceAngleBondIxn.h"
#include "ReferenceBondForce.h"
......@@ -1224,69 +1223,6 @@ void ReferenceCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& conte
obcParameters->setScaledRadiusFactors(scaleFactors);
}
ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
if (gbvi) {
GBVIParameters * gBVIParameters = gbvi->getGBVIParameters();
delete gBVIParameters;
delete gbvi;
}
}
void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii) {
int numParticles = system.getNumParticles();
charges.resize(numParticles);
vector<RealOpenMM> atomicRadii(numParticles);
vector<RealOpenMM> scaledRadii(numParticles);
vector<RealOpenMM> gammas(numParticles);
for (int i = 0; i < numParticles; ++i) {
double charge, radius, gamma;
force.getParticleParameters(i, charge, radius, gamma);
charges[i] = static_cast<RealOpenMM>(charge);
atomicRadii[i] = static_cast<RealOpenMM>(radius);
gammas[i] = static_cast<RealOpenMM>(gamma);
scaledRadii[i] = static_cast<RealOpenMM>(inputScaledRadii[i]);
}
GBVIParameters * gBVIParameters = new GBVIParameters(numParticles);
gBVIParameters->setAtomicRadii(atomicRadii);
gBVIParameters->setGammaParameters(gammas);
gBVIParameters->setScaledRadii(scaledRadii);
gBVIParameters->setSolventDielectric(static_cast<RealOpenMM>(force.getSolventDielectric()));
gBVIParameters->setSoluteDielectric(static_cast<RealOpenMM>(force.getSoluteDielectric()));
gBVIParameters->setBornRadiusScalingMethod(force.getBornRadiusScalingMethod());
gBVIParameters->setQuinticUpperBornRadiusLimit(static_cast<RealOpenMM>(force.getQuinticUpperBornRadiusLimit()));
gBVIParameters->setQuinticLowerLimitFactor(static_cast<RealOpenMM>(force.getQuinticLowerLimitFactor()));
if (force.getNonbondedMethod() != GBVIForce::NoCutoff)
gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
isPeriodic = (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic);
gbvi = new ReferenceGBVI(gBVIParameters);
}
double ReferenceCalcGBVIForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
if (isPeriodic)
gbvi->getGBVIParameters()->setPeriodic(extractBoxVectors(context));
RealOpenMM energy;
if (includeForces) {
vector<RealVec>& forceData = extractForces(context);
gbvi->computeBornForces(posData, charges, forceData);
energy = 0.0;
}
if (includeEnergy) {
energy = gbvi->computeBornEnergy(posData, charges);
}
return static_cast<double>(energy);
}
ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
disposeRealArray(particleParamArray, numParticles);
if (neighborList != NULL)
......
platforms/reference/src/ReferencePlatform.cpp
View file @ 167ae8a0
......@@ -58,7 +58,6 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(CalcGBVIForceKernel::Name(), factory);
registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
......
platforms/reference/src/SimTKReference/GBVIParameters.cpp deleted 100644 → 0
View file @ b20c20fe
/* Portions copyright (c) 2006-2009 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 <string.h>
#include "openmm/OpenMMException.h"
#include "GBVIParameters.h"
using std::vector;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
GBVIParameters constructor
@param numberOfAtoms number of atoms
--------------------------------------------------------------------------------------- */
GBVIParameters::GBVIParameters(int numberOfAtoms) : _numberOfAtoms(numberOfAtoms),
_soluteDielectric(1.0),
_solventDielectric(78.3),
_electricConstant(-0.5*ONE_4PI_EPS0),
_cutoff(false),
_periodic(false),
_bornRadiusScalingMethod(0),
_quinticLowerLimitFactor(0.8),
_quinticUpperBornRadiusLimit(5.0) {
_atomicRadii.resize(numberOfAtoms);
_scaledRadii.resize(numberOfAtoms);
_gammaParameters.resize(numberOfAtoms);
}
/**---------------------------------------------------------------------------------------
GBVIParameters destructor
--------------------------------------------------------------------------------------- */
GBVIParameters::~GBVIParameters() {
}
/**---------------------------------------------------------------------------------------
Get number of atoms
@return number of atoms
--------------------------------------------------------------------------------------- */
int GBVIParameters::getNumberOfAtoms() const {
return _numberOfAtoms;
}
/**---------------------------------------------------------------------------------------
Get electric constant
@return electric constant
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getElectricConstant() const {
return _electricConstant;
}
/**---------------------------------------------------------------------------------------
Get solvent dielectric
@return solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getSolventDielectric() const {
return _solventDielectric;
}
/**---------------------------------------------------------------------------------------
Set solvent dielectric
@param solventDielectric solvent dielectric
--------------------------------------------------------------------------------------- */
void GBVIParameters::setSolventDielectric(RealOpenMM solventDielectric) {
_solventDielectric = solventDielectric;
}
/**---------------------------------------------------------------------------------------
Get solute dielectric
@return soluteDielectric
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getSoluteDielectric() const {
return _soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Set solute dielectric
@param soluteDielectric solute dielectric
--------------------------------------------------------------------------------------- */
void GBVIParameters::setSoluteDielectric(RealOpenMM soluteDielectric) {
_soluteDielectric = soluteDielectric;
}
/**---------------------------------------------------------------------------------------
Get AtomicRadii array
@return array of atomic radii
--------------------------------------------------------------------------------------- */
const vector<RealOpenMM>& GBVIParameters::getAtomicRadii() const {
return _atomicRadii;
}
/**---------------------------------------------------------------------------------------
Set AtomicRadii array
@param atomicRadii vector of atomic radii
--------------------------------------------------------------------------------------- */
void GBVIParameters::setAtomicRadii(const vector<RealOpenMM>& atomicRadii) {
if (atomicRadii.size() == _atomicRadii.size()) {
for (unsigned int ii = 0; ii < atomicRadii.size(); ii++) {
_atomicRadii[ii] = atomicRadii[ii];
}
} else {
std::stringstream msg;
msg << "GBVIParameters: input size for atomic radii does not agree w/ current size: input=";
msg << atomicRadii.size();
msg << " current size=" << _atomicRadii.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Return scaled radii
@return array
--------------------------------------------------------------------------------------- */
const vector<RealOpenMM>& GBVIParameters::getScaledRadii() const {
return _scaledRadii;
}
/**---------------------------------------------------------------------------------------
Set scaled radii
@param scaledRadii scaledRadii
--------------------------------------------------------------------------------------- */
void GBVIParameters::setScaledRadii(const vector<RealOpenMM>& scaledRadii) {
if (scaledRadii.size() == _scaledRadii.size()) {
for (unsigned int ii = 0; ii < scaledRadii.size(); ii++) {
_scaledRadii[ii] = scaledRadii[ii];
}
} else {
std::stringstream msg;
msg << "GBVIParameters: input size for scaled radii does not agree w/ current size: input=";
msg << scaledRadii.size();
msg << " current size=" << _scaledRadii.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Return gamma parameters
If not previously set, allocate space
@return array
--------------------------------------------------------------------------------------- */
const vector<RealOpenMM>& GBVIParameters::getGammaParameters() const {
return _gammaParameters;
}
/**---------------------------------------------------------------------------------------
Set gamma parameters
@param gammas gammas
--------------------------------------------------------------------------------------- */
void GBVIParameters::setGammaParameters(const vector<RealOpenMM>& gammas) {
if (gammas.size() == _gammaParameters.size()) {
for (unsigned int ii = 0; ii < gammas.size(); ii++) {
_gammaParameters[ii] = gammas[ii];
}
} else {
std::stringstream msg;
msg << "GBVIParameters: input size for gammas does not agree w/ current size: input=";
msg << gammas.size();
msg << " current size=" << _gammaParameters.size();
throw OpenMM::OpenMMException(msg.str());
}
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
--------------------------------------------------------------------------------------- */
void GBVIParameters::setUseCutoff(RealOpenMM distance) {
_cutoff = true;
_cutoffDistance = distance;
}
/**---------------------------------------------------------------------------------------
Get whether to use a cutoff.
--------------------------------------------------------------------------------------- */
bool GBVIParameters::getUseCutoff() {
return _cutoff;
}
/**---------------------------------------------------------------------------------------
Get the cutoff distance.
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getCutoffDistance() {
return _cutoffDistance;
}
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void GBVIParameters::setPeriodic(RealVec* vectors) {
assert(_cutoff);
assert(vectors[0][0] >= 2.0*_cutoffDistance);
assert(vectors[1][1] >= 2.0*_cutoffDistance);
assert(vectors[2][2] >= 2.0*_cutoffDistance);
_periodic = true;
_periodicBoxVectors[0] = vectors[0];
_periodicBoxVectors[1] = vectors[1];
_periodicBoxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Get whether to use periodic boundary conditions.
--------------------------------------------------------------------------------------- */
bool GBVIParameters::getPeriodic() {
return _periodic;
}
/**---------------------------------------------------------------------------------------
Get the periodic box dimension
--------------------------------------------------------------------------------------- */
const OpenMM::RealVec* GBVIParameters::getPeriodicBox() {
return _periodicBoxVectors;
}
/**---------------------------------------------------------------------------------------
Get tau prefactor
@return (1/e1 - 1/e0), where e1 = solute dielectric, e0 = solvent dielectric
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getTau() const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
RealOpenMM tau;
if (getSoluteDielectric() != zero && getSolventDielectric() != zero) {
tau = (one/getSoluteDielectric()) - (one/getSolventDielectric());
} else {
tau = zero;
}
return tau;
}
/**---------------------------------------------------------------------------------------
Get bornRadiusScalingMethod
@return bornRadiusScalingMethod
--------------------------------------------------------------------------------------- */
int GBVIParameters::getBornRadiusScalingMethod() const {
return _bornRadiusScalingMethod;
}
/**---------------------------------------------------------------------------------------
Set bornRadiusScalingMethod
@param bornRadiusScalingMethod bornRadiusScalingMethod
--------------------------------------------------------------------------------------- */
void GBVIParameters::setBornRadiusScalingMethod(int bornRadiusScalingMethod) {
_bornRadiusScalingMethod = bornRadiusScalingMethod;
}
/**---------------------------------------------------------------------------------------
Get quinticLowerLimitFactor
@return quinticLowerLimitFactor
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getQuinticLowerLimitFactor() const {
return _quinticLowerLimitFactor;
}
/**---------------------------------------------------------------------------------------
Set quinticLowerLimitFactor
@param quinticLowerLimitFactor quinticLowerLimitFactor
--------------------------------------------------------------------------------------- */
void GBVIParameters::setQuinticLowerLimitFactor(RealOpenMM quinticLowerLimitFactor) {
_quinticLowerLimitFactor = quinticLowerLimitFactor;
}
/**---------------------------------------------------------------------------------------
Get quinticUpperBornRadiusLimit
@return quinticUpperBornRadiusLimit
--------------------------------------------------------------------------------------- */
RealOpenMM GBVIParameters::getQuinticUpperBornRadiusLimit() const {
return _quinticUpperBornRadiusLimit;
}
/**---------------------------------------------------------------------------------------
Set quinticUpperBornRadiusLimit
@param quinticUpperBornRadiusLimit quinticUpperBornRadiusLimit
--------------------------------------------------------------------------------------- */
void GBVIParameters::setQuinticUpperBornRadiusLimit(RealOpenMM quinticUpperBornRadiusLimit) {
_quinticUpperBornRadiusLimit = quinticUpperBornRadiusLimit;
}
platforms/reference/src/SimTKReference/ReferenceGBVI.cpp deleted 100644 → 0
View file @ b20c20fe
/* Portions copyright (c) 2006-2009 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 <string.h>
#include <math.h>
#include <sstream>
#include <stdio.h>
#include "ReferenceForce.h"
#include "ReferenceGBVI.h"
using namespace std;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
ReferenceGBVI constructor
gbviParameters gbviParameters object
--------------------------------------------------------------------------------------- */
ReferenceGBVI::ReferenceGBVI(GBVIParameters* gbviParameters) : _gbviParameters(gbviParameters) {
_switchDeriviative.resize(gbviParameters->getNumberOfAtoms());
}
/**---------------------------------------------------------------------------------------
ReferenceGBVI destructor
--------------------------------------------------------------------------------------- */
ReferenceGBVI::~ReferenceGBVI() {
}
/**---------------------------------------------------------------------------------------
Get GBVIParameters reference
@return GBVIParameters reference
--------------------------------------------------------------------------------------- */
GBVIParameters* ReferenceGBVI::getGBVIParameters() const {
return _gbviParameters;
}
/**---------------------------------------------------------------------------------------
Set GBVIParameters reference
@param GBVIParameters reference
--------------------------------------------------------------------------------------- */
void ReferenceGBVI::setGBVIParameters(GBVIParameters* gbviParameters) {
_gbviParameters = gbviParameters;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcParameters->getNumberOfAtoms()
@return array
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& ReferenceGBVI::getSwitchDeriviative() {
return _switchDeriviative;
}
/**---------------------------------------------------------------------------------------
Compute quintic spline value and associated derviative
@param x value to compute spline at
@param rl lower cutoff value
@param ru upper cutoff value
@param outValue value of spline at x
@param outDerivative value of derivative of spline at x
--------------------------------------------------------------------------------------- */
void ReferenceGBVI::quinticSpline(RealOpenMM x, RealOpenMM rl, RealOpenMM ru,
RealOpenMM* outValue, RealOpenMM* outDerivative) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0);
static const RealOpenMM minusSix = static_cast<RealOpenMM>( -6.0);
static const RealOpenMM minusTen = static_cast<RealOpenMM>(-10.0);
static const RealOpenMM minusThirty = static_cast<RealOpenMM>(-30.0);
static const RealOpenMM fifteen = static_cast<RealOpenMM>( 15.0);
static const RealOpenMM sixty = static_cast<RealOpenMM>( 60.0);
// ---------------------------------------------------------------------------------------
RealOpenMM numerator = x - rl;
RealOpenMM denominator = ru - rl;
RealOpenMM ratio = numerator/denominator;
RealOpenMM ratio2 = ratio*ratio;
RealOpenMM ratio3 = ratio2*ratio;
*outValue = one + ratio3*(minusTen + fifteen*ratio + minusSix*ratio2);
*outDerivative = ratio2*(minusThirty + sixty*ratio + minusThirty*ratio2)/denominator;
}
/**---------------------------------------------------------------------------------------
Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
and quintic splice switching function
@param atomicRadius3 atomic radius cubed
@param bornSum Born sum (volume integral)
@param gbviParameters Gbvi parameters (parameters used in spline
QuinticLowerLimitFactor & QuinticUpperBornRadiusLimit)
@param bornRadius output Born radius
@param switchDeriviative output switching function deriviative
--------------------------------------------------------------------------------------- */
void ReferenceGBVI::computeBornRadiiUsingQuinticSpline(RealOpenMM atomicRadius3, RealOpenMM bornSum,
GBVIParameters* gbviParameters,
RealOpenMM* bornRadius, RealOpenMM* switchDeriviative) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0);
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0);
static const RealOpenMM minusOne = static_cast<RealOpenMM>(-1.0);
static const RealOpenMM minusThree = static_cast<RealOpenMM>(-3.0);
static const RealOpenMM oneEighth = static_cast<RealOpenMM>( 0.125);
static const RealOpenMM minusOneThird = static_cast<RealOpenMM>((-1.0/3.0));
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0);
// ---------------------------------------------------------------------------------------
// R = [ S(V)*(A - V) ]**(-1/3)
// S(V) = 1 V < L
// S(V) = qSpline + U/(A-V) L < V < A
// S(V) = U/(A-V) U < V
// dR/dr = (-1/3)*[ S(V)*(A - V) ]**(-4/3)*[ d{ S(V)*(A-V) }/dr
// d{ S(V)*(A-V) }/dr = (dV/dr)*[ (A-V)*dS/dV - S(V) ]
// (A - V)*dS/dV - S(V) = 0 - 1 V < L
// (A - V)*dS/dV - S(V) = (A-V)*d(qSpline) + (A-V)*U/(A-V)**2 - qSpline - U/(A-V)
// = (A-V)*d(qSpline) - qSpline L < V < A**(-3)
// (A - V)*dS/dV - S(V) = (A-V)*U*/(A-V)**2 - U/(A-V) = 0 U < V
RealOpenMM splineL = gbviParameters->getQuinticLowerLimitFactor()*atomicRadius3;
RealOpenMM sum;
if (bornSum > splineL) {
if (bornSum < atomicRadius3) {
RealOpenMM splineValue, splineDerivative;
quinticSpline(bornSum, splineL, atomicRadius3, &splineValue, &splineDerivative);
sum = (atomicRadius3 - bornSum)*splineValue + gbviParameters->getQuinticUpperBornRadiusLimit();
*switchDeriviative = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
} else {
sum = gbviParameters->getQuinticUpperBornRadiusLimit();
*switchDeriviative = zero;
}
} else {
sum = atomicRadius3 - bornSum;
*switchDeriviative = one;
}
*bornRadius = POW(sum, minusOneThird);
}
/**---------------------------------------------------------------------------------------
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param chain not used here
--------------------------------------------------------------------------------------- */
void ReferenceGBVI::computeBornRadii(const vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM minusThree = static_cast<RealOpenMM>(-3.0);
static const RealOpenMM oneEighth = static_cast<RealOpenMM>(0.125);
static const RealOpenMM minusOneThird = static_cast<RealOpenMM>((-1.0/3.0));
static const RealOpenMM three = static_cast<RealOpenMM>(3.0);
// ---------------------------------------------------------------------------------------
GBVIParameters* gbviParameters = getGBVIParameters();
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const vector<RealOpenMM>& atomicRadii = gbviParameters->getAtomicRadii();
const vector<RealOpenMM>& scaledRadii = gbviParameters->getScaledRadii();
vector<RealOpenMM>& switchDeriviatives = getSwitchDeriviative();
// ---------------------------------------------------------------------------------------
// calculate Born radii
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM sum = zero;
// sum over volumes
for (int atomJ = 0; atomJ < numberOfAtoms; atomJ++) {
if (atomJ != atomI) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
continue;
sum += ReferenceGBVI::getVolume(r, radiusI, scaledRadii[atomJ]);
}
}
RealOpenMM atomicRadius3 = POW(radiusI, minusThree);
if (_gbviParameters->getBornRadiusScalingMethod() != GBVIParameters::QuinticSpline) {
sum = atomicRadius3 - sum;
bornRadii[atomI] = POW(sum, minusOneThird);
switchDeriviatives[atomI] = one;
} else {
RealOpenMM bornRadius, switchDeriviative;
computeBornRadiiUsingQuinticSpline(atomicRadius3, sum, gbviParameters,
&bornRadius, &switchDeriviative);
bornRadii[atomI] = bornRadius;
switchDeriviatives[atomI] = switchDeriviative;
}
}
}
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceGBVI::getVolume(RealOpenMM r, RealOpenMM R, RealOpenMM S) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0);
static const RealOpenMM minusThree = static_cast<RealOpenMM>(-3.0);
RealOpenMM diff = (S - R);
if (FABS(diff) < r) {
RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
return (ReferenceGBVI::getL(r, (r + S), S) -
ReferenceGBVI::getL(r, lowerBound, S));
} else if (r <= diff) {
return ReferenceGBVI::getL(r, (r + S), S) -
ReferenceGBVI::getL(r, (r - S), S) +
POW(R, minusThree);
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
Get L (used in analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceGBVI::getL(RealOpenMM r, RealOpenMM x, RealOpenMM S) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM threeHalves = static_cast<RealOpenMM>(1.5);
static const RealOpenMM third = static_cast<RealOpenMM>((1.0/3.0));
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
static const RealOpenMM eighth = static_cast<RealOpenMM>(0.125);
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
return (threeHalves*xInv)*((xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv));
}
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceGBVI::dL_dr(RealOpenMM r, RealOpenMM x, RealOpenMM S) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM threeHalves = static_cast<RealOpenMM>(1.5);
static const RealOpenMM threeEights = static_cast<RealOpenMM>(0.375);
static const RealOpenMM third = static_cast<RealOpenMM>((1.0/3.0));
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
static const RealOpenMM eighth = static_cast<RealOpenMM>(0.125);
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM rInv2 = rInv*rInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
return ((-threeHalves*xInv2*rInv2)*(fourth + eighth*diff2*xInv2) + threeEights*xInv3*xInv);
}
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceGBVI::dL_dx(RealOpenMM r, RealOpenMM x, RealOpenMM S) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0);
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5);
static const RealOpenMM threeHalvesM = static_cast<RealOpenMM>(-1.5);
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0));
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff = (r + S)*(r - S);
return (threeHalvesM*xInv3)*((half*rInv) - xInv + (half*diff*xInv2*rInv));
}
/**---------------------------------------------------------------------------------------
Sgb function
@param t r*r*G_i*G_j
@return Sgb (top of p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceGBVI::Sgb(RealOpenMM t) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "ReferenceGBVI::Sgb";
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
// ---------------------------------------------------------------------------------------
return ((t != zero) ? one/SQRT((one + (fourth*EXP(-t))/t)) : zero);
}
/**---------------------------------------------------------------------------------------
Get GB/VI energy
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceGBVI::computeBornEnergy(const vector<RealVec>& atomCoordinates, const vector<RealOpenMM>& partialCharges) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM two = static_cast<RealOpenMM>(2.0);
static const RealOpenMM three = static_cast<RealOpenMM>(3.0);
static const RealOpenMM four = static_cast<RealOpenMM>(4.0);
static const RealOpenMM half = static_cast<RealOpenMM>(0.5);
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
static const RealOpenMM eighth = static_cast<RealOpenMM>(0.125);
// ---------------------------------------------------------------------------------------
const GBVIParameters* gbviParameters = getGBVIParameters();
const RealOpenMM preFactor = gbviParameters->getElectricConstant();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const vector<RealOpenMM>& atomicRadii = gbviParameters->getAtomicRadii();
const vector<RealOpenMM>& gammaParameters = gbviParameters->getGammaParameters();
// compute Born radii
vector<RealOpenMM> bornRadii(numberOfAtoms);
computeBornRadii(atomCoordinates, bornRadii);
// ---------------------------------------------------------------------------------------
// Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
// to minimze roundoff error sum cavityEnergy separately since in general much
// smaller than other contributions
RealOpenMM energy = zero;
RealOpenMM cavityEnergy = zero;
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
RealOpenMM partialChargeI = partialCharges[atomI];
// self-energy term
RealOpenMM atomIEnergy = half*partialChargeI/bornRadii[atomI];
// cavity term
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
cavityEnergy += gammaParameters[atomI]*ratio*ratio*ratio;
for (int atomJ = atomI + 1; atomJ < numberOfAtoms; atomJ++) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM t = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);
atomIEnergy += partialCharges[atomJ]*Sgb(t)/deltaR[ReferenceForce::RIndex];
}
energy += two*partialChargeI*atomIEnergy;
}
energy *= preFactor;
energy -= cavityEnergy;
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
return (conversion*energy);
}
/**---------------------------------------------------------------------------------------
Get GB/VI forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
void ReferenceGBVI::computeBornForces(std::vector<RealVec>& atomCoordinates, const vector<RealOpenMM>& partialCharges,
std::vector<OpenMM::RealVec>& inputForces) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM two = static_cast<RealOpenMM>(2.0);
static const RealOpenMM three = static_cast<RealOpenMM>(3.0);
static const RealOpenMM four = static_cast<RealOpenMM>(4.0);
static const RealOpenMM half = static_cast<RealOpenMM>(0.5);
static const RealOpenMM oneThird = static_cast<RealOpenMM>((1.0/3.0));
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
static const RealOpenMM eighth = static_cast<RealOpenMM>(0.125);
// ---------------------------------------------------------------------------------------
const GBVIParameters* gbviParameters = getGBVIParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const vector<RealOpenMM>& atomicRadii = gbviParameters->getAtomicRadii();
const vector<RealOpenMM>& gammaParameters = gbviParameters->getGammaParameters();
// ---------------------------------------------------------------------------------------
// constants
const RealOpenMM preFactor = two*gbviParameters->getElectricConstant();
// ---------------------------------------------------------------------------------------
// compute Born radii
vector<RealOpenMM> bornRadii(numberOfAtoms);
computeBornRadii(atomCoordinates, bornRadii);
// set energy/forces to zero
std::vector<OpenMM::RealVec> forces(numberOfAtoms);
for (int ii = 0; ii < numberOfAtoms; ii++) {
forces[ii][0] = zero;
forces[ii][1] = zero;
forces[ii][2] = zero;
}
vector<RealOpenMM> bornForces(numberOfAtoms, 0.0);
// ---------------------------------------------------------------------------------------
// first main loop
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
// partial of polar term wrt Born radius
// and (dGpol/dr)(dr/dx)
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for (int atomJ = atomI; atomJ < numberOfAtoms; atomJ++) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM expTerm = EXP(-D_ij);
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT(denominator2);
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*(one - fourth*expTerm)/denominator2;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*(one + D_ij)/denominator2;
if (atomI != atomJ) {
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
// ---------------------------------------------------------------------------------------
// second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
// dGpol/dBornRadius) = bornForces[]
// dBornRadius/dr = (1/3)*(bR**4)*(dV/dr)
const vector<RealOpenMM>& scaledRadii = gbviParameters->getScaledRadii();
const vector<RealOpenMM>& switchDeriviative = getSwitchDeriviative();
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
RealOpenMM R = atomicRadii[atomI];
// partial of cavity term wrt Born radius
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
bornForces[atomI] += (three*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
bornForces[atomI] *= switchDeriviative[atomI]*oneThird*b2*b2;
for (int atomJ = 0; atomJ < numberOfAtoms; atomJ++) {
if (atomJ != atomI) {
RealOpenMM deltaX = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
RealOpenMM deltaY = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
RealOpenMM deltaZ = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
deltaX = deltaR[ReferenceForce::XIndex];
deltaY = deltaR[ReferenceForce::YIndex];
deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = SQRT(r2);
RealOpenMM S = scaledRadii[atomJ];
RealOpenMM diff = (S - R);
RealOpenMM de = zero;
// find dRb/dr, where Rb is the Born radius
if (FABS(diff) < r) {
de = ReferenceGBVI::dL_dr(r, r+S, S) + ReferenceGBVI::dL_dx(r, r+S, S);
if (R > (r - S)) {
de -= ReferenceGBVI::dL_dr(r, R, S);
} else {
de -= (ReferenceGBVI::dL_dr(r, (r-S), S) + ReferenceGBVI::dL_dx(r, (r-S), S));
}
} else if (r < (S - R)) {
de = ReferenceGBVI::dL_dr(r, r+S, S) + ReferenceGBVI::dL_dx(r, r+S, S);
de -= (ReferenceGBVI::dL_dr(r, r-S, S) + ReferenceGBVI::dL_dx(r, r-S, S));
}
// de = (dG/dRb)(dRb/dr)
de *= bornForces[atomI]/r;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
}
}
// convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
inputForces[atomI][0] += conversion*forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1];
inputForces[atomI][2] += conversion*forces[atomI][2];
}
}
/**---------------------------------------------------------------------------------------
Use double precision
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
double ReferenceGBVI::getVolumeD(double r, double R, double S) {
// ---------------------------------------------------------------------------------------
static const double zero = 0.0;
static const double minusThree = -3.0;
double diff = (S - R);
if (fabs(diff) < r) {
double lowerBound = (R > (r - S)) ? R : (r - S);
return (ReferenceGBVI::getLD(r, (r + S), S) -
ReferenceGBVI::getLD(r, lowerBound, S));
} else if (r < diff) {
return ReferenceGBVI::getLD(r, (r + S), S) -
ReferenceGBVI::getLD(r, (r - S), S) +
pow(R, minusThree);
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
Use double precision
Get L (used in analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
double ReferenceGBVI::getLD(double r, double x, double S) {
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff2 = (r + S)*(r - S);
return (threeHalves*xInv)*((xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv));
}
/**---------------------------------------------------------------------------------------
Use double precision
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
double ReferenceGBVI::dL_drD(double r, double x, double S) {
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double threeEights = 0.375;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double rInv2 = rInv*rInv;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff2 = (r + S)*(r - S);
return ((-threeHalves*xInv2*rInv2)*(fourth + eighth*diff2*xInv2) + threeEights*xInv3*xInv);
}
/**---------------------------------------------------------------------------------------
Use double precision
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
double ReferenceGBVI::dL_dxD(double r, double x, double S) {
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double half = 0.5;
static const double threeHalvesM = -1.5;
static const double third = 1.0/3.0;
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff = (r + S)*(r - S);
return (threeHalvesM*xInv3)*((half*rInv) - xInv + (half*diff*xInv2*rInv));
}
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