Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
tsoc
openmm
Commits
4868314a
"vscode:/vscode.git/clone" did not exist on "3c16ceff2cc1964149da9f0f3a794d365efb0695"
Commit
4868314a
authored
Sep 08, 2010
by
Mark Friedrichs
Browse files
Initial Amnoeba reference code
parent
49173b75
Changes
16
Hide whitespace changes
Inline
Side-by-side
Showing
16 changed files
with
2923 additions
and
0 deletions
+2923
-0
plugins/amoeba/platforms/reference/CMakeLists.txt
plugins/amoeba/platforms/reference/CMakeLists.txt
+138
-0
plugins/amoeba/platforms/reference/include/AmoebaReferenceKernelFactory.h
...latforms/reference/include/AmoebaReferenceKernelFactory.h
+45
-0
plugins/amoeba/platforms/reference/include/windowsExportAmoebaReference.h
...latforms/reference/include/windowsExportAmoebaReference.h
+41
-0
plugins/amoeba/platforms/reference/sharedTarget/CMakeLists.txt
...ns/amoeba/platforms/reference/sharedTarget/CMakeLists.txt
+76
-0
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp
.../platforms/reference/src/AmoebaReferenceKernelFactory.cpp
+103
-0
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
...amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
+774
-0
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
...s/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
+399
-0
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.cpp
...rms/reference/src/SimTKReference/AmoebaReferenceForce.cpp
+150
-0
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.h
...forms/reference/src/SimTKReference/AmoebaReferenceForce.h
+121
-0
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.cpp
.../src/SimTKReference/AmoebaReferenceHarmonicAngleForce.cpp
+183
-0
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.h
...ce/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.h
+105
-0
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.cpp
...e/src/SimTKReference/AmoebaReferenceHarmonicBondForce.cpp
+89
-0
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.h
...nce/src/SimTKReference/AmoebaReferenceHarmonicBondForce.h
+80
-0
plugins/amoeba/platforms/reference/tests/CMakeLists.txt
plugins/amoeba/platforms/reference/tests/CMakeLists.txt
+81
-0
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaHarmonicAngleForce.cpp
...reference/tests/TestReferenceAmoebaHarmonicAngleForce.cpp
+314
-0
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaHarmonicBondForce.cpp
.../reference/tests/TestReferenceAmoebaHarmonicBondForce.cpp
+224
-0
No files found.
plugins/amoeba/platforms/reference/CMakeLists.txt
0 → 100644
View file @
4868314a
#---------------------------------------------------
# OpenMMAmoeba REFERENCE Platform
#
# Creates OpenMM library, base name=OpenMMAmoebaReference.
# Default libraries are shared & optimized. Variants
# are created for static (_static) and debug (_d).
#
# Windows:
# OpenMMAmoebaReference[_d].dll
# OpenMMAmoebaReference[_d].lib
# OpenMMAmoebaReference_static[_d].lib
# Unix:
# libOpenMMAmoebaReference[_d].so
# libOpenMMAmoebaReference_static[_d].a
#----------------------------------------------------
# ----------------------------------------------------------------------------
# logging
SET
(
LOG TRUE
)
IF
(
LOG
)
SET
(
LOG_FILE
"CMakeLog.txt"
)
FILE
(
WRITE
${
LOG_FILE
}
"In plugins/amoeba/platforms/reference
\n
"
)
ENDIF
(
LOG
)
IF
(
LOG
)
MACRO
(
LOG_DIR LOG_FILE DIR_LIST
)
FILE
(
APPEND
${
LOG_FILE
}
"
\n
${
DIR_LIST
}
\n
"
)
FOREACH
(
currentFile
${
ARGN
}
)
FILE
(
APPEND
${
LOG_FILE
}
"
${
currentFile
}
\n
"
)
ENDFOREACH
(
currentFile
)
ENDMACRO
(
LOG_DIR
)
ENDIF
(
LOG
)
MESSAGE
(
"YYY Reference tests"
)
# ----------------------------------------------------------------------------
SET
(
DO_TESTS TRUE
)
# 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_AMOEBA_SOURCE_SUBDIRS .
)
# Collect up information about the version of the OpenMM library we're building
# and make it available to the code so it can be built into the binaries.
SET
(
OPENMM_REFERENCE_LIBRARY_NAME OpenMMAmoebaReference
)
SET
(
SHARED_TARGET
${
OPENMM_REFERENCE_LIBRARY_NAME
}
)
SET
(
STATIC_TARGET
${
OPENMM_REFERENCE_LIBRARY_NAME
}
_static
)
# Ensure that debug libraries have "_d" appended to their names.
# CMake gets this right on Windows automatically with this definition.
IF
(
${
CMAKE_GENERATOR
}
MATCHES
"Visual Studio"
)
SET
(
CMAKE_DEBUG_POSTFIX
"_d"
CACHE INTERNAL
""
FORCE
)
ENDIF
(
${
CMAKE_GENERATOR
}
MATCHES
"Visual Studio"
)
# But on Unix or Cygwin we have to add the suffix manually
IF
(
UNIX AND CMAKE_BUILD_TYPE MATCHES Debug
)
SET
(
SHARED_TARGET
${
SHARED_TARGET
}
_d
)
SET
(
STATIC_TARGET
${
STATIC_TARGET
}
_d
)
ENDIF
(
UNIX AND CMAKE_BUILD_TYPE MATCHES Debug
)
# These are all the places to search for header files which are
# to be part of the API.
SET
(
API_AMOEBA_INCLUDE_DIRS
)
# start empty
FOREACH
(
subdir
${
OPENMM_AMOEBA_SOURCE_SUBDIRS
}
)
# append
SET
(
API_AMOEBA_INCLUDE_DIRS
${
API_AMOEBA_INCLUDE_DIRS
}
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/include
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/include/internal
)
ENDFOREACH
(
subdir
)
## ----------------------------------------------------------------------------
IF
(
LOG
)
LOG_DIR
(
${
LOG_FILE
}
"OPENMM_SOURCE_SUBDIRS"
${
OPENMM_SOURCE_SUBDIRS
}
)
LOG_DIR
(
${
LOG_FILE
}
"OPENMM_AMOEBA_SOURCE_SUBDIRS"
${
OPENMM_AMOEBA_SOURCE_SUBDIRS
}
)
LOG_DIR
(
${
LOG_FILE
}
"API_AMOEBA_INCLUDE_DIRS"
${
API_AMOEBA_INCLUDE_DIRS
}
)
LOG_DIR
(
${
LOG_FILE
}
"CMAKE_CURRENT_SOURCE_DIR"
${
CMAKE_CURRENT_SOURCE_DIR
}
)
ENDIF
(
LOG
)
## ----------------------------------------------------------------------------
# We'll need both *relative* path names, starting with their API_AMOEBA_INCLUDE_DIRS,
# and absolute pathnames.
SET
(
API_AMOEBA_REL_INCLUDE_FILES
)
# start these out empty
SET
(
API_AMOEBA_ABS_INCLUDE_FILES
)
FOREACH
(
dir
${
API_AMOEBA_INCLUDE_DIRS
}
)
FILE
(
GLOB fullpaths
${
dir
}
/*.h
)
# returns full pathnames
SET
(
API_AMOEBA_ABS_INCLUDE_FILES
${
API_AMOEBA_ABS_INCLUDE_FILES
}
${
fullpaths
}
)
FOREACH
(
pathname
${
fullpaths
}
)
GET_FILENAME_COMPONENT
(
filename
${
pathname
}
NAME
)
SET
(
API_AMOEBA_REL_INCLUDE_FILES
${
API_AMOEBA_REL_INCLUDE_FILES
}
${
dir
}
/
${
filename
}
)
ENDFOREACH
(
pathname
)
ENDFOREACH
(
dir
)
## ----------------------------------------------------------------------------
IF
(
LOG
)
LOG_DIR
(
${
LOG_FILE
}
"API_AMOEBA_REL_INCLUDE_FILES"
${
API_AMOEBA_REL_INCLUDE_FILES
}
)
LOG_DIR
(
${
LOG_FILE
}
"OPENMM_DIR"
${
OPENMM_DIR
}
)
ENDIF
(
LOG
)
## ----------------------------------------------------------------------------
# collect up source files
SET
(
SOURCE_FILES
)
# empty
SET
(
SOURCE_INCLUDE_FILES
)
FOREACH
(
subdir
${
OPENMM_AMOEBA_SOURCE_SUBDIRS
}
)
FILE
(
GLOB_RECURSE src_files
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/src/*.cpp
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/src/*.c
)
FILE
(
GLOB incl_files
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/src/*.h
)
SET
(
SOURCE_FILES
${
SOURCE_FILES
}
${
src_files
}
)
#append
SET
(
SOURCE_INCLUDE_FILES
${
SOURCE_INCLUDE_FILES
}
${
incl_files
}
)
IF
(
LOG
)
LOG_DIR
(
${
LOG_FILE
}
"Adding include dir: "
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/include
)
ENDIF
(
LOG
)
INCLUDE_DIRECTORIES
(
${
CMAKE_CURRENT_SOURCE_DIR
}
/
${
subdir
}
/include
)
ENDFOREACH
(
subdir
)
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/platforms/reference/include
)
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/platforms/reference/src
)
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/platforms/reference/src/SimTKReference
)
IF
(
LOG
)
LOG_DIR
(
${
LOG_FILE
}
"CMAKE_CURRENT_SOURCE_DIR"
${
CMAKE_CURRENT_SOURCE_DIR
}
)
FILE
(
APPEND
${
LOG_FILE
}
"CMAKE_CURRENT_SOURCE_DIR=
${
CMAKE_CURRENT_SOURCE_DIR
}
\n
"
)
LOG_DIR
(
${
LOG_FILE
}
"SOURCE_FILES"
${
SOURCE_FILES
}
)
LOG_DIR
(
${
LOG_FILE
}
"SOURCE_INCLUDE_FILES"
${
SOURCE_INCLUDE_FILES
}
)
ENDIF
(
LOG
)
INCLUDE_DIRECTORIES
(
${
CMAKE_CURRENT_SOURCE_DIR
}
/src
)
MESSAGE
(
"YYY Reference tests"
)
# SUBDIRS (sharedTarget staticTarget)
SUBDIRS
(
sharedTarget
)
SUBDIRS
(
tests
)
plugins/amoeba/platforms/reference/include/AmoebaReferenceKernelFactory.h
0 → 100644
View file @
4868314a
#ifndef AMOEBA_OPENMM_REFERERENCE_KERNEL_FACTORY_H_
#define AMOEBA_OPENMM_REFERERENCE_KERNEL_FACTORY_H_
/* -------------------------------------------------------------------------- *
* AmoebaOpenMM *
* -------------------------------------------------------------------------- *
* 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: *
* 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/KernelFactory.h"
namespace
OpenMM
{
/**
* This KernelFactory creates all kernels for AmoebaReferencePlatform.
*/
class
AmoebaReferenceKernelFactory
:
public
KernelFactory
{
public:
KernelImpl
*
createKernelImpl
(
std
::
string
name
,
const
Platform
&
platform
,
ContextImpl
&
context
)
const
;
};
}
// namespace OpenMM
#endif
/*AMOEBA_OPENMM_REFERERENCE_KERNEL_FACTORY_H_*/
plugins/amoeba/platforms/reference/include/windowsExportAmoebaReference.h
0 → 100644
View file @
4868314a
#ifndef OPENMM_WINDOWS_EXPORT_REFERENCE_H_
#define OPENMM_WINDOWS_EXPORT_REFERENCE_H_
/*
* Shared libraries are messy in Visual Studio. We have to distinguish three
* cases:
* (1) this header is being used to build the OpenMM shared library
* (dllexport)
* (2) this header is being used by a *client* of the OpenMM shared
* library (dllimport)
* (3) we are building the OpenMM static library, or the client is
* being compiled with the expectation of linking with the
* OpenMM static library (nothing special needed)
* In the CMake script for building this library, we define one of the symbols
* OpenMMCUDA_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_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 OpenMM 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)
// Keep MS VC++ quiet about lack of dll export of private members.
#pragma warning(disable:4251)
#if defined(OPENMM_AMOEBA_REFERENCE_BUILDING_SHARED_LIBRARY)
#define OPENMM_AMOEBA_REFERENCE_EXPORT __declspec(dllexport)
#elif defined(OPENMM_AMOEBA_REFERENCE_BUILDING_STATIC_LIBRARY) || defined(OPENMM_AMOEBA_REFERENCE_USE_STATIC_LIBRARIES)
#define OPENMM_AMOEBA_REFERENCE_EXPORT
#else
#define OPENMM_AMOEBA_REFERENCE_EXPORT __declspec(dllimport) // i.e., a client of a shared library
#endif
#else
#define OPENMM_AMOEBA_REFERENCE_EXPORT // Linux, Mac
#endif
#endif // OPENMM_WINDOWS_EXPORT_REFERENCE_H_
plugins/amoeba/platforms/reference/sharedTarget/CMakeLists.txt
0 → 100644
View file @
4868314a
#
# Include REFERENCE related files.
#
# ----------------------------------------------------------------------------
# logging
SET
(
LOG FALSE
)
IF
(
LOG
)
SET
(
LOG_FILE
"CMakeLog.txt"
)
FILE
(
WRITE
${
LOG_FILE
}
"In amoeba/platforms/reference/sharedTarget Cmake
\n
"
)
# FILE( APPEND ${LOG_FILE} "BROOK_LIB_PATH=${BROOK_LIB_PATH}\n")
ENDIF
(
LOG
)
IF
(
LOG
)
MACRO
(
LOG_DIR LOG_FILE DIR_LIST
)
FILE
(
APPEND
${
LOG_FILE
}
"
\n
${
DIR_LIST
}
\n
"
)
FOREACH
(
currentFile
${
ARGN
}
)
FILE
(
APPEND
${
LOG_FILE
}
"
${
currentFile
}
\n
"
)
ENDFOREACH
(
currentFile
)
ENDMACRO
(
LOG_DIR
)
ENDIF
(
LOG
)
# ----------------------------------------------------------------------------
SET
(
OPENMM_BUILD_AMOEBA_PATH
${
CMAKE_SOURCE_DIR
}
/plugins/amoeba
)
# ----------------------------------------------------------------------------
IF
(
LOG
)
LOG_DIR
(
${
LOG_FILE
}
"Pre OPENMM_SOURCE_SUBDIRS"
${
OPENMM_SOURCE_SUBDIRS
}
)
LOG_DIR
(
${
LOG_FILE
}
"Pre OPENMM_AMOEBA_SOURCE_SUBDIRS "
${
OPENMM_AMOEBA_SOURCE_SUBDIRS
}
)
LOG_DIR
(
${
LOG_FILE
}
"Pre SOURCE_FILES"
${
SOURCE_FILES
}
)
ENDIF
(
LOG
)
## ----------------------------------------------------------------------------
INCLUDE_DIRECTORIES
(
${
REFERENCE_INCLUDE
}
)
LINK_DIRECTORIES
(
${
REFERENCE_TARGET_LINK
}
)
FOREACH
(
subdir
${
OPENMM_AMOEBA_SOURCE_SUBDIRS
}
)
FILE
(
GLOB src_files
${
OPENMM_BUILD_AMOEBA_PATH
}
/platforms/reference/
${
subdir
}
/src/*.cu
${
OPENMM_BUILD_AMOEBA_PATH
}
/platforms/reference/src/*/*.cu
)
SET
(
SOURCE_FILES
${
SOURCE_FILES
}
${
src_files
}
)
INCLUDE_DIRECTORIES
(
BEFORE
${
OPENMM_BUILD_AMOEBA_PATH
}
/platforms/reference/../
${
subdir
}
/include
)
ENDFOREACH
(
subdir
)
# ----------------------------------------------------------------------------
IF
(
LOG
)
LOG_DIR
(
${
LOG_FILE
}
"OPENMM_BUILD_AMOEBA_PATH"
${
OPENMM_BUILD_AMOEBA_PATH
}
)
FILE
(
APPEND
${
LOG_FILE
}
"OPENMM_BUILD_AMOEBA_PATH=
${
OPENMM_BUILD_AMOEBA_PATH
}
\n
"
)
LOG_DIR
(
${
LOG_FILE
}
"OPENMM_SOURCE_SUBDIRS"
${
OPENMM_SOURCE_SUBDIRS
}
)
LOG_DIR
(
${
LOG_FILE
}
"CMAKE_SOURCE_DIR"
${
CMAKE_SOURCE_DIR
}
)
LOG_DIR
(
${
LOG_FILE
}
"REFERENCE_INCLUDE"
${
REFERENCE_INCLUDE
}
)
LOG_DIR
(
${
LOG_FILE
}
"REFERENCE_TARGET_LINK"
${
REFERENCE_TARGET_LINK
}
)
LOG_DIR
(
${
LOG_FILE
}
"SHARED_TARGET"
${
SHARED_TARGET
}
)
LOG_DIR
(
${
LOG_FILE
}
"OPENMM_DIR"
${
OPENMM_DIR
}
)
LOG_DIR
(
${
LOG_FILE
}
"SOURCE_FILES"
${
SOURCE_FILES
}
)
ENDIF
(
LOG
)
## ----------------------------------------------------------------------------
# REFERENCE_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/jama/include)
INCLUDE_DIRECTORIES
(
${
OPENMM_BUILD_AMOEBA_PATH
}
/platforms/reference/../src
${
OPENMM_BUILD_AMOEBA_PATH
}
/platforms/reference/include
${
OPENMM_DIR
}
/platforms/reference/src
${
OPENMM_DIR
}
/platforms/reference/include
${
OPENMM_DIR
}
/platforms/reference/src/kernels
${
OPENMM_DIR
}
/openmmapi/include
)
ADD_LIBRARY
(
${
SHARED_TARGET
}
SHARED
${
SOURCE_FILES
}
${
SOURCE_INCLUDE_FILES
}
${
API_ABS_INCLUDE_FILES
}
)
TARGET_LINK_LIBRARIES
(
${
SHARED_TARGET
}
debug
${
OPENMM_LIBRARY_NAME
}
_d optimized
${
OPENMM_LIBRARY_NAME
}
)
TARGET_LINK_LIBRARIES
(
${
SHARED_TARGET
}
debug
${
OPENMM_LIBRARY_NAME
}
Cuda_d optimized
${
OPENMM_LIBRARY_NAME
}
Cuda
)
TARGET_LINK_LIBRARIES
(
${
SHARED_TARGET
}
debug
${
OPENMM_AMOEBA_LIBRARY_NAME
}
_d optimized
${
OPENMM_AMOEBA_LIBRARY_NAME
}
)
SET_TARGET_PROPERTIES
(
${
SHARED_TARGET
}
PROPERTIES COMPILE_FLAGS
"-DOPENMMREFERENCE_BUILDING_SHARED_LIBRARY -DOPENMMREFERENCEAMOEBA_BUILDING_SHARED_LIBRARY"
)
INSTALL_TARGETS
(
/lib/plugins RUNTIME_DIRECTORY /lib/plugins
${
SHARED_TARGET
}
)
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp
0 → 100644
View file @
4868314a
/* -------------------------------------------------------------------------- *
* AmoebaOpenMM *
* -------------------------------------------------------------------------- *
* 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: *
* 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 "AmoebaReferenceKernelFactory.h"
#include "AmoebaReferenceKernels.h"
#include "ReferencePlatform.h"
#include "windowsExportAmoebaReference.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/OpenMMException.h"
using
namespace
OpenMM
;
extern
"C"
void
OPENMM_AMOEBA_REFERENCE_EXPORT
registerKernelFactories
()
{
fprintf
(
stderr
,
"In registerKernelFactories AmoebaReferenceKernelFactory
\n
"
);
fflush
(
stderr
);
for
(
int
ii
=
0
;
ii
<
Platform
::
getNumPlatforms
();
ii
++
){
Platform
&
platform
=
Platform
::
getPlatform
(
ii
);
if
(
platform
.
getName
()
==
"Reference"
){
AmoebaReferenceKernelFactory
*
factory
=
new
AmoebaReferenceKernelFactory
();
platform
.
registerKernelFactory
(
CalcAmoebaHarmonicBondForceKernel
::
Name
(),
factory
);
platform
.
registerKernelFactory
(
CalcAmoebaHarmonicAngleForceKernel
::
Name
(),
factory
);
/*
platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaWcaDispersionForceKernel::Name(), factory);
*/
}
}
}
KernelImpl
*
AmoebaReferenceKernelFactory
::
createKernelImpl
(
std
::
string
name
,
const
Platform
&
platform
,
ContextImpl
&
context
)
const
{
ReferencePlatform
::
PlatformData
&
referencePlatformData
=
*
static_cast
<
ReferencePlatform
::
PlatformData
*>
(
context
.
getPlatformData
());
// create AmoebaReferenceData object if contextToAmoebaDataMap does not contain
// key equal to current context
if
(
name
==
CalcAmoebaHarmonicBondForceKernel
::
Name
())
return
new
ReferenceCalcAmoebaHarmonicBondForceKernel
(
name
,
platform
,
context
.
getSystem
());
if
(
name
==
CalcAmoebaHarmonicAngleForceKernel
::
Name
())
return
new
ReferenceCalcAmoebaHarmonicAngleForceKernel
(
name
,
platform
,
context
.
getSystem
());
/*
if (name == CalcAmoebaHarmonicInPlaneAngleForceKernel::Name())
return new ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaTorsionForceKernel::Name())
return new ReferenceCalcAmoebaTorsionForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaPiTorsionForceKernel::Name())
return new ReferenceCalcAmoebaPiTorsionForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaStretchBendForceKernel::Name())
return new ReferenceCalcAmoebaStretchBendForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaOutOfPlaneBendForceKernel::Name())
return new ReferenceCalcAmoebaOutOfPlaneBendForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaTorsionTorsionForceKernel::Name())
return new ReferenceCalcAmoebaTorsionTorsionForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaMultipoleForceKernel::Name())
return new ReferenceCalcAmoebaMultipoleForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaGeneralizedKirkwoodForceKernel::Name())
return new ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaVdwForceKernel::Name())
return new ReferenceCalcAmoebaVdwForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaWcaDispersionForceKernel::Name())
return new ReferenceCalcAmoebaWcaDispersionForceKernel(name, platform, context.getSystem());
*/
throw
OpenMMException
((
std
::
string
(
"Tried to create kernel with illegal kernel name '"
)
+
name
+
"'"
).
c_str
());
}
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
0 → 100644
View file @
4868314a
/* -------------------------------------------------------------------------- *
* AmoebaOpenMM *
* -------------------------------------------------------------------------- *
* 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: *
* 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 "AmoebaReferenceKernels.h"
#include "AmoebaReferenceHarmonicBondForce.h"
#include "AmoebaReferenceHarmonicAngleForce.h"
#include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h"
//#include "internal/AmoebaMultipoleForceImpl.h"
//#include "internal/AmoebaWcaDispersionForceImpl.h"
#include <cmath>
#ifdef _MSC_VER
#include <windows.h>
#endif
using
namespace
OpenMM
;
using
namespace
std
;
static
RealOpenMM
**
extractPositions
(
ContextImpl
&
context
)
{
ReferencePlatform
::
PlatformData
*
data
=
reinterpret_cast
<
ReferencePlatform
::
PlatformData
*>
(
context
.
getPlatformData
());
return
(
RealOpenMM
**
)
data
->
positions
;
}
static
RealOpenMM
**
extractVelocities
(
ContextImpl
&
context
)
{
ReferencePlatform
::
PlatformData
*
data
=
reinterpret_cast
<
ReferencePlatform
::
PlatformData
*>
(
context
.
getPlatformData
());
return
(
RealOpenMM
**
)
data
->
velocities
;
}
static
RealOpenMM
**
extractForces
(
ContextImpl
&
context
)
{
ReferencePlatform
::
PlatformData
*
data
=
reinterpret_cast
<
ReferencePlatform
::
PlatformData
*>
(
context
.
getPlatformData
());
return
(
RealOpenMM
**
)
data
->
forces
;
}
static
RealOpenMM
*
extractBoxSize
(
ContextImpl
&
context
)
{
ReferencePlatform
::
PlatformData
*
data
=
reinterpret_cast
<
ReferencePlatform
::
PlatformData
*>
(
context
.
getPlatformData
());
return
(
RealOpenMM
*
)
data
->
periodicBoxSize
;
}
// ***************************************************************************
ReferenceCalcAmoebaHarmonicBondForceKernel
::
ReferenceCalcAmoebaHarmonicBondForceKernel
(
std
::
string
name
,
const
Platform
&
platform
,
System
&
system
)
:
CalcAmoebaHarmonicBondForceKernel
(
name
,
platform
),
system
(
system
)
{
}
ReferenceCalcAmoebaHarmonicBondForceKernel
::~
ReferenceCalcAmoebaHarmonicBondForceKernel
()
{
}
void
ReferenceCalcAmoebaHarmonicBondForceKernel
::
initialize
(
const
System
&
system
,
const
AmoebaHarmonicBondForce
&
force
)
{
numBonds
=
force
.
getNumBonds
();
for
(
int
ii
=
0
;
ii
<
numBonds
;
ii
++
)
{
int
particle1Index
,
particle2Index
;
double
lengthValue
,
kValue
;
force
.
getBondParameters
(
ii
,
particle1Index
,
particle2Index
,
lengthValue
,
kValue
);
particle1
.
push_back
(
particle1Index
);
particle2
.
push_back
(
particle2Index
);
length
.
push_back
(
static_cast
<
RealOpenMM
>
(
lengthValue
)
);
kQuadratic
.
push_back
(
static_cast
<
RealOpenMM
>
(
kValue
)
);
}
globalHarmonicBondCubic
=
static_cast
<
RealOpenMM
>
(
force
.
getAmoebaGlobalHarmonicBondCubic
());
globalHarmonicBondQuartic
=
static_cast
<
RealOpenMM
>
(
force
.
getAmoebaGlobalHarmonicBondQuartic
());
}
double
ReferenceCalcAmoebaHarmonicBondForceKernel
::
execute
(
ContextImpl
&
context
,
bool
includeForces
,
bool
includeEnergy
)
{
RealOpenMM
**
posData
=
extractPositions
(
context
);
RealOpenMM
**
forceData
=
extractForces
(
context
);
RealOpenMM
energy
=
0.0
;
for
(
int
ii
=
0
;
ii
<
numBonds
;
ii
++
)
{
int
particle1Index
=
particle1
[
ii
];
int
particle2Index
=
particle2
[
ii
];
RealOpenMM
bondLength
=
length
[
ii
];
RealOpenMM
bondK
=
kQuadratic
[
ii
];
RealOpenMM
*
forces
[
2
];
forces
[
0
]
=
forceData
[
particle1Index
];
forces
[
1
]
=
forceData
[
particle2Index
];
energy
+=
AmoebaReferenceHarmonicBondForce
::
calculateForceAndEnergy
(
posData
[
particle1Index
],
posData
[
particle2Index
],
bondLength
,
bondK
,
globalHarmonicBondCubic
,
globalHarmonicBondQuartic
,
forces
);
}
return
static_cast
<
double
>
(
energy
);
}
ReferenceCalcAmoebaHarmonicAngleForceKernel
::
ReferenceCalcAmoebaHarmonicAngleForceKernel
(
std
::
string
name
,
const
Platform
&
platform
,
System
&
system
)
:
CalcAmoebaHarmonicAngleForceKernel
(
name
,
platform
),
system
(
system
)
{
}
ReferenceCalcAmoebaHarmonicAngleForceKernel
::~
ReferenceCalcAmoebaHarmonicAngleForceKernel
()
{
}
void
ReferenceCalcAmoebaHarmonicAngleForceKernel
::
initialize
(
const
System
&
system
,
const
AmoebaHarmonicAngleForce
&
force
)
{
numAngles
=
force
.
getNumAngles
();
for
(
int
ii
=
0
;
ii
<
numAngles
;
ii
++
)
{
int
particle1Index
,
particle2Index
,
particle3Index
;
double
angleValue
,
k
;
force
.
getAngleParameters
(
ii
,
particle1Index
,
particle2Index
,
particle3Index
,
angleValue
,
k
);
particle1
.
push_back
(
particle1Index
);
particle2
.
push_back
(
particle2Index
);
particle3
.
push_back
(
particle3Index
);
angle
.
push_back
(
static_cast
<
RealOpenMM
>
(
angleValue
)
);
kQuadratic
.
push_back
(
static_cast
<
RealOpenMM
>
(
k
)
);
}
globalHarmonicAngleCubic
=
static_cast
<
RealOpenMM
>
(
force
.
getAmoebaGlobalHarmonicAngleCubic
());
globalHarmonicAngleQuartic
=
static_cast
<
RealOpenMM
>
(
force
.
getAmoebaGlobalHarmonicAngleQuartic
());
globalHarmonicAnglePentic
=
static_cast
<
RealOpenMM
>
(
force
.
getAmoebaGlobalHarmonicAnglePentic
());
globalHarmonicAngleSextic
=
static_cast
<
RealOpenMM
>
(
force
.
getAmoebaGlobalHarmonicAngleSextic
());
}
double
ReferenceCalcAmoebaHarmonicAngleForceKernel
::
execute
(
ContextImpl
&
context
,
bool
includeForces
,
bool
includeEnergy
)
{
RealOpenMM
**
posData
=
extractPositions
(
context
);
RealOpenMM
**
forceData
=
extractForces
(
context
);
RealOpenMM
energy
=
0.0
;
for
(
unsigned
int
ii
=
0
;
ii
<
numAngles
;
ii
++
)
{
int
particle1Index
=
particle1
[
ii
];
int
particle2Index
=
particle2
[
ii
];
int
particle3Index
=
particle3
[
ii
];
RealOpenMM
idealAngle
=
angle
[
ii
];
RealOpenMM
angleK
=
kQuadratic
[
ii
];
RealOpenMM
*
forces
[
3
];
forces
[
0
]
=
forceData
[
particle1Index
];
forces
[
1
]
=
forceData
[
particle2Index
];
forces
[
2
]
=
forceData
[
particle3Index
];
energy
+=
AmoebaReferenceHarmonicAngleForce
::
calculateForceAndEnergy
(
posData
[
particle1Index
],
posData
[
particle2Index
],
posData
[
particle3Index
],
idealAngle
,
angleK
,
globalHarmonicAngleCubic
,
globalHarmonicAngleQuartic
,
globalHarmonicAnglePentic
,
globalHarmonicAngleSextic
,
forces
);
}
return
static_cast
<
double
>
(
energy
);
}
//ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform), system(system) {
//}
//
//ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::~ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel() {
//}
//
//void ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& force) {
//
// numAngles = force.getNumAngles();
//
// std::vector<int> particle1(numAngles);
// std::vector<int> particle2(numAngles);
// std::vector<int> particle3(numAngles);
// std::vector<int> particle4(numAngles);
// std::vector<RealOpenMM> angle(numAngles);
// std::vector<RealOpenMM> k(numAngles);
//
// for (int i = 0; i < numAngles; i++) {
// double angleValue, kQuadratic;
// force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], angleValue, kQuadratic);
// //angle[i] = static_cast<RealOpenMM>( (angleValue*RadiansToDegrees) );
// angle[i] = static_cast<RealOpenMM>( angleValue );
// k[i] = static_cast<RealOpenMM>( kQuadratic );
// }
///*
// gpuSetAmoebaInPlaneAngleParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, angle, k,
// static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAngleCubic()),
// static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAngleQuartic()),
// static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAnglePentic()),
// static_cast<RealOpenMM>( force.getAmoebaGlobalHarmonicInPlaneAngleSextic() ) );
//*/
//
//}
//
//double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// if( data.getAmoebaLocalForcesKernel() == this ){
// computeAmoebaLocalForces( data );
// }
// return 0.0;
//}
//
//ReferenceCalcAmoebaTorsionForceKernel::ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaTorsionForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaTorsionForceKernel::~ReferenceCalcAmoebaTorsionForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaTorsionForceKernel::initialize(const System& system, const AmoebaTorsionForce& force) {
//
// data.setAmoebaLocalForcesKernel( this );
// numTorsions = force.getNumTorsions();
// std::vector<int> particle1(numTorsions);
// std::vector<int> particle2(numTorsions);
// std::vector<int> particle3(numTorsions);
// std::vector<int> particle4(numTorsions);
//
// std::vector< std::vector<RealOpenMM> > torsionParameters1(numTorsions);
// std::vector< std::vector<RealOpenMM> > torsionParameters2(numTorsions);
// std::vector< std::vector<RealOpenMM> > torsionParameters3(numTorsions);
//
// for (int i = 0; i < numTorsions; i++) {
//
// std::vector<double> torsionParameter1;
// std::vector<double> torsionParameter2;
// std::vector<double> torsionParameter3;
//
// std::vector<RealOpenMM> torsionParameters1F(3);
// std::vector<RealOpenMM> torsionParameters2F(3);
// std::vector<RealOpenMM> torsionParameters3F(3);
//
// force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], torsionParameter1, torsionParameter2, torsionParameter3 );
// for ( unsigned int jj = 0; jj < torsionParameter1.size(); jj++) {
// torsionParameters1F[jj] = static_cast<RealOpenMM>(torsionParameter1[jj]);
// torsionParameters2F[jj] = static_cast<RealOpenMM>(torsionParameter2[jj]);
// torsionParameters3F[jj] = static_cast<RealOpenMM>(torsionParameter3[jj]);
// }
// torsionParameters1[i] = torsionParameters1F;
// torsionParameters2[i] = torsionParameters2F;
// torsionParameters3[i] = torsionParameters3F;
// }
// gpuSetAmoebaTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, torsionParameters1, torsionParameters2, torsionParameters3 );
//
//}
//
//double ReferenceCalcAmoebaTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// if( data.getAmoebaLocalForcesKernel() == this ){
// computeAmoebaLocalForces( data );
// }
// return 0.0;
//}
//
//ReferenceCalcAmoebaPiTorsionForceKernel::ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaPiTorsionForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaPiTorsionForceKernel::~ReferenceCalcAmoebaPiTorsionForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaPiTorsionForceKernel::initialize(const System& system, const AmoebaPiTorsionForce& force) {
//
// data.setAmoebaLocalForcesKernel( this );
// numPiTorsions = force.getNumPiTorsions();
//
// std::vector<int> particle1(numPiTorsions);
// std::vector<int> particle2(numPiTorsions);
// std::vector<int> particle3(numPiTorsions);
// std::vector<int> particle4(numPiTorsions);
// std::vector<int> particle5(numPiTorsions);
// std::vector<int> particle6(numPiTorsions);
//
// std::vector<RealOpenMM> torsionKParameters(numPiTorsions);
//
// for (int i = 0; i < numPiTorsions; i++) {
//
// double torsionKParameter;
//
// force.getPiTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], particle5[i], particle6[i], torsionKParameter);
// torsionKParameters[i] = static_cast<RealOpenMM>(torsionKParameter);
// }
// gpuSetAmoebaPiTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, particle5, particle6, torsionKParameters);
//}
//
//double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// if( data.getAmoebaLocalForcesKernel() == this ){
// computeAmoebaLocalForces( data );
// }
// return 0.0;
//}
//
//ReferenceCalcAmoebaStretchBendForceKernel::ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaStretchBendForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaStretchBendForceKernel::~ReferenceCalcAmoebaStretchBendForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system, const AmoebaStretchBendForce& force) {
//
// data.setAmoebaLocalForcesKernel( this );
// numStretchBends = force.getNumStretchBends();
//
// std::vector<int> particle1(numStretchBends);
// std::vector<int> particle2(numStretchBends);
// std::vector<int> particle3(numStretchBends);
// std::vector<RealOpenMM> lengthABParameters(numStretchBends);
// std::vector<RealOpenMM> lengthCBParameters(numStretchBends);
// std::vector<RealOpenMM> angleParameters(numStretchBends);
// std::vector<RealOpenMM> kParameters(numStretchBends);
//
// for (int i = 0; i < numStretchBends; i++) {
//
// double lengthAB, lengthCB, angle, k;
//
// force.getStretchBendParameters(i, particle1[i], particle2[i], particle3[i], lengthAB, lengthCB, angle, k);
// lengthABParameters[i] = static_cast<RealOpenMM>(lengthAB);
// lengthCBParameters[i] = static_cast<RealOpenMM>(lengthCB);
// angleParameters[i] = static_cast<RealOpenMM>(angle);
// kParameters[i] = static_cast<RealOpenMM>(k);
// }
// gpuSetAmoebaStretchBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, lengthABParameters, lengthCBParameters, angleParameters, kParameters);
//
//}
//
//double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// if( data.getAmoebaLocalForcesKernel() == this ){
// computeAmoebaLocalForces( data );
// }
// return 0.0;
//}
//ReferenceCalcAmoebaOutOfPlaneBendForceKernel::ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaOutOfPlaneBendForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaOutOfPlaneBendForceKernel::~ReferenceCalcAmoebaOutOfPlaneBendForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& system, const AmoebaOutOfPlaneBendForce& force) {
//
// data.setAmoebaLocalForcesKernel( this );
// numOutOfPlaneBends = force.getNumOutOfPlaneBends();
//
// std::vector<int> particle1(numOutOfPlaneBends);
// std::vector<int> particle2(numOutOfPlaneBends);
// std::vector<int> particle3(numOutOfPlaneBends);
// std::vector<int> particle4(numOutOfPlaneBends);
// std::vector<RealOpenMM> kParameters(numOutOfPlaneBends);
//
// for (int i = 0; i < numOutOfPlaneBends; i++) {
//
// double k;
//
// force.getOutOfPlaneBendParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], k);
// kParameters[i] = static_cast<RealOpenMM>(k);
// }
// gpuSetAmoebaOutOfPlaneBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, kParameters,
// static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendCubic()),
// static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendQuartic()),
// static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendPentic()),
// static_cast<RealOpenMM>( force.getAmoebaGlobalOutOfPlaneBendSextic() ) );
//
//}
//
//double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// if( data.getAmoebaLocalForcesKernel() == this ){
// computeAmoebaLocalForces( data );
// }
// return 0.0;
//}
//
//ReferenceCalcAmoebaTorsionTorsionForceKernel::ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaTorsionTorsionForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaTorsionTorsionForceKernel::~ReferenceCalcAmoebaTorsionTorsionForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, const AmoebaTorsionTorsionForce& force) {
//
// data.setAmoebaLocalForcesKernel( this );
// numTorsionTorsions = force.getNumTorsionTorsions();
//
// // torsion-torsion parameters
//
// std::vector<int> particle1(numTorsionTorsions);
// std::vector<int> particle2(numTorsionTorsions);
// std::vector<int> particle3(numTorsionTorsions);
// std::vector<int> particle4(numTorsionTorsions);
// std::vector<int> particle5(numTorsionTorsions);
// std::vector<int> chiralCheckAtomIndex(numTorsionTorsions);
// std::vector<int> gridIndices(numTorsionTorsions);
//
// for (int i = 0; i < numTorsionTorsions; i++) {
// force.getTorsionTorsionParameters(i, particle1[i], particle2[i], particle3[i],
// particle4[i], particle5[i],
// chiralCheckAtomIndex[i], gridIndices[i]);
// }
// gpuSetAmoebaTorsionTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, particle5, chiralCheckAtomIndex, gridIndices );
//
// // torsion-torsion grids
//
// numTorsionTorsionGrids = force.getNumTorsionTorsionGrids();
// std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > > RealOpenMMGrids;
//
// RealOpenMMGrids.resize(numTorsionTorsionGrids);
// for (int i = 0; i < numTorsionTorsionGrids; i++) {
//
// TorsionTorsionGrid grid;
// force.getTorsionTorsionGrid(i, grid );
//
// RealOpenMMGrids[i].resize( grid.size() );
// for (unsigned int ii = 0; ii < grid.size(); ii++) {
//
// RealOpenMMGrids[i][ii].resize( grid[ii].size() );
// for (unsigned int jj = 0; jj < grid[ii].size(); jj++) {
//
// RealOpenMMGrids[i][ii][jj].resize( grid[ii][jj].size() );
// for (unsigned int kk = 0; kk < grid[ii][kk].size(); kk++) {
// RealOpenMMGrids[i][ii][jj][kk] = static_cast<RealOpenMM>(grid[ii][jj][kk]);
// }
// }
// }
// }
// gpuSetAmoebaTorsionTorsionGrids(data.getAmoebaGpu(), RealOpenMMGrids );
//
//}
//
//double CudaCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// if( data.getAmoebaLocalForcesKernel() == this ){
// computeAmoebaLocalForces( data );
// }
// return 0.0;
//}
//
///* -------------------------------------------------------------------------- *
// * AmoebaMultipole *
// * -------------------------------------------------------------------------- */
//
//static void computeAmoebaMultipoleForce( AmoebaCudaData& data ) {
//
// amoebaGpuContext gpu = data.getAmoebaGpu();
// data.initializeGpu();
//
// if( 0 && data.getLog() ){
// (void) fprintf( data.getLog(), "computeAmoebaMultipoleForce\n" );
// (void) fflush( data.getLog());
// }
//
// // calculate Born radii
//
// if( data.getHasAmoebaGeneralizedKirkwood() ){
// kCalculateObcGbsaBornSum(gpu->gpuContext);
// kReduceObcGbsaBornSum(gpu->gpuContext);
// }
//
// // multipoles
//
// kCalculateAmoebaMultipoleForces(gpu, data.getHasAmoebaGeneralizedKirkwood() );
//
////kClearForces(gpu->gpuContext);
////kClearEnergy(gpu->gpuContext);
////(void) fprintf( data.getLog(), "computeAmoebaMultipoleForce clearing forces/energy after kCalculateAmoebaMultipoleForces()\n" );
//
// // GK
//
// if( data.getHasAmoebaGeneralizedKirkwood() ){
// kCalculateAmoebaKirkwood(gpu);
// }
//
// if( 0 && data.getLog() ){
// (void) fprintf( data.getLog(), "completed computeAmoebaMultipoleForce\n" );
// (void) fflush( data.getLog());
// }
//}
//
//CudaCalcAmoebaMultipoleForceKernel::CudaCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, AmoebaCudaData& data, System& system) :
// CalcAmoebaMultipoleForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//CudaCalcAmoebaMultipoleForceKernel::~CudaCalcAmoebaMultipoleForceKernel() {
// data.decrementKernelCount();
//}
//
//void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const AmoebaMultipoleForce& force) {
//
// numMultipoles = force.getNumMultipoles();
//
// data.setHasAmoebaMultipole( true );
//
// std::vector<RealOpenMM> charges(numMultipoles);
// std::vector<RealOpenMM> dipoles(3*numMultipoles);
// std::vector<RealOpenMM> quadrupoles(9*numMultipoles);
// std::vector<RealOpenMM> tholes(numMultipoles);
// std::vector<RealOpenMM> dampingFactors(numMultipoles);
// std::vector<RealOpenMM> polarity(numMultipoles);
// std::vector<int> axisTypes(numMultipoles);
// std::vector<int> multipoleAtomId1s(numMultipoles);
// std::vector<int> multipoleAtomId2s(numMultipoles);
// std::vector< std::vector< std::vector<int> > > multipoleAtomCovalentInfo(numMultipoles);
// std::vector<int> minCovalentIndices(numMultipoles);
// std::vector<int> minCovalentPolarizationIndices(numMultipoles);
//
// RealOpenMM scalingDistanceCutoff = static_cast<RealOpenMM>(force.getScalingDistanceCutoff());
//
// std::vector<AmoebaMultipoleForce::CovalentType> covalentList;
// covalentList.push_back( AmoebaMultipoleForce::Covalent12 );
// covalentList.push_back( AmoebaMultipoleForce::Covalent13 );
// covalentList.push_back( AmoebaMultipoleForce::Covalent14 );
// covalentList.push_back( AmoebaMultipoleForce::Covalent15 );
//
// std::vector<AmoebaMultipoleForce::CovalentType> polarizationCovalentList;
// polarizationCovalentList.push_back( AmoebaMultipoleForce::PolarizationCovalent11 );
// polarizationCovalentList.push_back( AmoebaMultipoleForce::PolarizationCovalent12 );
// polarizationCovalentList.push_back( AmoebaMultipoleForce::PolarizationCovalent13 );
// polarizationCovalentList.push_back( AmoebaMultipoleForce::PolarizationCovalent14 );
//
// std::vector<int> covalentDegree;
// AmoebaMultipoleForceImpl::getCovalentDegree( force, covalentDegree );
// int dipoleIndex = 0;
// int quadrupoleIndex = 0;
// int maxCovalentRange = 0;
// double totalCharge = 0.0;
// for (int i = 0; i < numMultipoles; i++) {
//
// // multipoles
//
// int axisType, multipoleAtomId1, multipoleAtomId2;
// double charge, tholeD, dampingFactorD, polarityD;
// std::vector<double> dipolesD;
// std::vector<double> quadrupolesD;
// force.getMultipoleParameters(i, charge, dipolesD, quadrupolesD, axisType, multipoleAtomId1, multipoleAtomId2,
// tholeD, dampingFactorD, polarityD );
//
// totalCharge += charge;
// axisTypes[i] = axisType;
// multipoleAtomId1s[i] = multipoleAtomId1;
// multipoleAtomId2s[i] = multipoleAtomId2;
//
// charges[i] = static_cast<RealOpenMM>(charge);
// tholes[i] = static_cast<RealOpenMM>(tholeD);
// dampingFactors[i] = static_cast<RealOpenMM>(dampingFactorD);
// polarity[i] = static_cast<RealOpenMM>(polarityD);
//
// dipoles[dipoleIndex++] = static_cast<RealOpenMM>(dipolesD[0]);
// dipoles[dipoleIndex++] = static_cast<RealOpenMM>(dipolesD[1]);
// dipoles[dipoleIndex++] = static_cast<RealOpenMM>(dipolesD[2]);
//
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[0]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[1]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[2]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[3]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[4]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[5]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[6]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[7]);
// quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[8]);
//
// // covalent info
//
// std::vector< std::vector<int> > covalentLists;
// force.getCovalentMaps(i, covalentLists );
// multipoleAtomCovalentInfo[i] = covalentLists;
//
// int minCovalentIndex, maxCovalentIndex;
// AmoebaMultipoleForceImpl::getCovalentRange( force, i, covalentList, &minCovalentIndex, &maxCovalentIndex );
// minCovalentIndices[i] = minCovalentIndex;
// if( maxCovalentRange < (maxCovalentIndex - minCovalentIndex) ){
// maxCovalentRange = maxCovalentIndex - minCovalentIndex;
// }
//
// AmoebaMultipoleForceImpl::getCovalentRange( force, i, polarizationCovalentList, &minCovalentIndex, &maxCovalentIndex );
// minCovalentPolarizationIndices[i] = minCovalentIndex;
// if( maxCovalentRange < (maxCovalentIndex - minCovalentIndex) ){
// maxCovalentRange = maxCovalentIndex - minCovalentIndex;
// }
// }
//
// int iterativeMethod = static_cast<int>(force.getMutualInducedIterationMethod());
// if( iterativeMethod != 0 ){
// throw OpenMMException("Iterative method for mutual induced dipoles not recognized.\n");
// }
//
// int nonbondedMethod = static_cast<int>(force.getNonbondedMethod());
// if( nonbondedMethod != 0 && nonbondedMethod != 1 ){
// throw OpenMMException("AmoebaMultipoleForce nonbonded method not recognized.\n");
// }
//
// gpuSetAmoebaMultipoleParameters(data.getAmoebaGpu(), charges, dipoles, quadrupoles, axisTypes, multipoleAtomId1s, multipoleAtomId2s,
// tholes, scalingDistanceCutoff, dampingFactors, polarity,
// multipoleAtomCovalentInfo, covalentDegree, minCovalentIndices, minCovalentPolarizationIndices, (maxCovalentRange+2),
// static_cast<int>(force.getMutualInducedIterationMethod()),
// force.getMutualInducedMaxIterations(),
// static_cast<RealOpenMM>( force.getMutualInducedTargetEpsilon()),
// nonbondedMethod,
// static_cast<RealOpenMM>( force.getCutoffDistance()),
// static_cast<RealOpenMM>( force.getAEwald()),
// static_cast<RealOpenMM>( force.getElectricConstant()) );
//
//}
//
//double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// computeAmoebaMultipoleForce( data );
// return 0.0;
//}
//
///* -------------------------------------------------------------------------- *
// * AmoebaGeneralizedKirkwood *
// * -------------------------------------------------------------------------- */
//
//ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaGeneralizedKirkwoodForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force) {
//
// data.setHasAmoebaGeneralizedKirkwood( true );
//
// int numParticles = system.getNumParticles();
//
// std::vector<RealOpenMM> radius(numParticles);
// std::vector<RealOpenMM> scale(numParticles);
// std::vector<RealOpenMM> charge(numParticles);
//
// for( int ii = 0; ii < numParticles; ii++ ){
// double particleCharge, particleRadius, scalingFactor;
// force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor);
// radius[ii] = static_cast<RealOpenMM>( particleRadius );
// scale[ii] = static_cast<RealOpenMM>( scalingFactor );
// charge[ii] = static_cast<RealOpenMM>( particleCharge );
// }
// gpuSetAmoebaObcParameters( data.getAmoebaGpu(), static_cast<RealOpenMM>(force.getSoluteDielectric() ),
// static_cast<RealOpenMM>( force.getSolventDielectric() ),
// static_cast<RealOpenMM>( force.getDielectricOffset() ), radius, scale, charge,
// force.getIncludeCavityTerm(),
// static_cast<RealOpenMM>( force.getProbeRadius() ),
// static_cast<RealOpenMM>( force.getSurfaceAreaFactor() ) );
//}
//
//double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// // handled in computeAmoebaMultipoleForce()
// return 0.0;
//}
//
//static void computeAmoebaVdwForce( AmoebaReferenceData& data ) {
//
// amoebaGpuContext gpu = data.getAmoebaGpu();
// data.initializeGpu();
//
// // Vdw14_7F
//
// kCalculateAmoebaVdw14_7Forces(gpu);
//}
//
//ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaVdwForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
//
// // per-particle parameters
//
// int numParticles = system.getNumParticles();
//
// std::vector<int> indexIVs(numParticles);
// std::vector<int> indexClasses(numParticles);
// std::vector< std::vector<int> > allExclusions(numParticles);
// std::vector<RealOpenMM> sigmas(numParticles);
// std::vector<RealOpenMM> epsilons(numParticles);
// std::vector<RealOpenMM> reductions(numParticles);
// for( int ii = 0; ii < numParticles; ii++ ){
//
// int indexIV, indexClass;
// double sigma, epsilon, reduction;
// std::vector<int> exclusions;
//
// force.getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
// force.getParticleExclusions( ii, exclusions );
// for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){
// allExclusions[ii].push_back( exclusions[jj] );
// }
//
// indexIVs[ii] = indexIV;
// indexClasses[ii] = indexClass;
// sigmas[ii] = static_cast<RealOpenMM>( sigma );
// epsilons[ii] = static_cast<RealOpenMM>( epsilon );
// reductions[ii] = static_cast<RealOpenMM>( reduction );
// }
//
// gpuSetAmoebaVdwParameters( data.getAmoebaGpu(), indexIVs, indexClasses, sigmas, epsilons, reductions,
// force.getSigmaCombiningRule(), force.getEpsilonCombiningRule(),
// allExclusions );
//}
//
//double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// computeAmoebaVdwForce( data );
// return 0.0;
//}
//
///* -------------------------------------------------------------------------- *
// * AmoebaWcaDispersion *
// * -------------------------------------------------------------------------- */
//
//static void computeAmoebaWcaDispersionForce( AmoebaReferenceData& data ) {
//
// data.initializeGpu();
// if( 0 && data.getLog() ){
// (void) fprintf( data.getLog(), "Calling computeAmoebaWcaDispersionForce " ); (void) fflush( data.getLog() );
// }
//
// kCalculateAmoebaWcaDispersionForces( data.getAmoebaGpu() );
//
// if( 0 && data.getLog() ){
// (void) fprintf( data.getLog(), " -- completed\n" ); (void) fflush( data.getLog() );
// }
//}
//
//ReferenceCalcAmoebaWcaDispersionForceKernel::ReferenceCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaWcaDispersionForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaWcaDispersionForceKernel::~ReferenceCalcAmoebaWcaDispersionForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaWcaDispersionForceKernel::initialize(const System& system, const AmoebaWcaDispersionForce& force) {
//
// // per-particle parameters
//
// int numParticles = system.getNumParticles();
// std::vector<RealOpenMM> radii(numParticles);
// std::vector<RealOpenMM> epsilons(numParticles);
// for( int ii = 0; ii < numParticles; ii++ ){
//
// double radius, epsilon;
// force.getParticleParameters( ii, radius, epsilon );
//
// radii[ii] = static_cast<RealOpenMM>( radius );
// epsilons[ii] = static_cast<RealOpenMM>( epsilon );
// }
// RealOpenMM totalMaximumDispersionEnergy = static_cast<RealOpenMM>( AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy( force ) );
// gpuSetAmoebaWcaDispersionParameters( data.getAmoebaGpu(), radii, epsilons, totalMaximumDispersionEnergy,
// static_cast<RealOpenMM>( force.getEpso( )),
// static_cast<RealOpenMM>( force.getEpsh( )),
// static_cast<RealOpenMM>( force.getRmino( )),
// static_cast<RealOpenMM>( force.getRminh( )),
// static_cast<RealOpenMM>( force.getAwater( )),
// static_cast<RealOpenMM>( force.getShctd( )),
// static_cast<RealOpenMM>( force.getDispoff( ) ) );
//}
//
//double ReferenceCalcAmoebaWcaDispersionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// computeAmoebaWcaDispersionForce( data );
// return 0.0;
//}
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
0 → 100644
View file @
4868314a
#ifndef AMOEBA_OPENMM_REFERENCE_KERNELS_H_
#define AMOEBA_OPENMM_REFERENCE_KERNELS_H_
/* -------------------------------------------------------------------------- *
* AmoebaOpenMM *
* -------------------------------------------------------------------------- *
* 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: *
* 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 "amoebaKernels.h"
#include "openmm/System.h"
#include "SimTKUtilities/SimTKOpenMMRealType.h"
namespace
OpenMM
{
/**
* This kernel is invoked by AmoebaHarmonicBondForce to calculate the forces acting on the system and the energy of the system.
*/
class
ReferenceCalcAmoebaHarmonicBondForceKernel
:
public
CalcAmoebaHarmonicBondForceKernel
{
public:
ReferenceCalcAmoebaHarmonicBondForceKernel
(
std
::
string
name
,
const
Platform
&
platform
,
System
&
system
);
~
ReferenceCalcAmoebaHarmonicBondForceKernel
();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the AmoebaHarmonicBondForce this kernel will be used for
*/
void
initialize
(
const
System
&
system
,
const
AmoebaHarmonicBondForce
&
force
);
/**
* 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:
int
numBonds
;
std
::
vector
<
int
>
particle1
;
std
::
vector
<
int
>
particle2
;
std
::
vector
<
RealOpenMM
>
length
;
std
::
vector
<
RealOpenMM
>
kQuadratic
;
RealOpenMM
globalHarmonicBondCubic
;
RealOpenMM
globalHarmonicBondQuartic
;
System
&
system
;
};
/**
* This kernel is invoked by AmoebaHarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
*/
class
ReferenceCalcAmoebaHarmonicAngleForceKernel
:
public
CalcAmoebaHarmonicAngleForceKernel
{
public:
ReferenceCalcAmoebaHarmonicAngleForceKernel
(
std
::
string
name
,
const
Platform
&
platform
,
System
&
system
);
~
ReferenceCalcAmoebaHarmonicAngleForceKernel
();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the AmoebaHarmonicAngleForce this kernel will be used for
*/
void
initialize
(
const
System
&
system
,
const
AmoebaHarmonicAngleForce
&
force
);
/**
* 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:
int
numAngles
;
std
::
vector
<
int
>
particle1
;
std
::
vector
<
int
>
particle2
;
std
::
vector
<
int
>
particle3
;
std
::
vector
<
RealOpenMM
>
angle
;
std
::
vector
<
RealOpenMM
>
kQuadratic
;
RealOpenMM
globalHarmonicAngleCubic
;
RealOpenMM
globalHarmonicAngleQuartic
;
RealOpenMM
globalHarmonicAnglePentic
;
RealOpenMM
globalHarmonicAngleSextic
;
System
&
system
;
};
// /**
// * This kernel is invoked by AmoebaHarmonicInPlaneAngleForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel : public CalcAmoebaHarmonicInPlaneAngleForceKernel {
// public:
// ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaHarmonicInPlaneAngleForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaHarmonicInPlaneAngleForce& force);
// /**
// * 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:
// int numAngles;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaTorsionForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaTorsionForceKernel : public CalcAmoebaTorsionForceKernel {
// public:
// ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaTorsionForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaTorsionForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaTorsionForce& force);
// /**
// * 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:
// int numTorsions;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel {
// public:
// ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaPiTorsionForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaPiTorsionForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaPiTorsionForce& force);
// /**
// * 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:
// int numPiTorsions;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
// public:
// ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaStretchBendForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaStretchBendForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaStretchBendForce& force);
// /**
// * 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:
// int numStretchBends;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaOutOfPlaneBendForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
// public:
// ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaOutOfPlaneBendForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaOutOfPlaneBendForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaOutOfPlaneBendForce& force);
// /**
// * 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:
// int numOutOfPlaneBends;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel {
// public:
// ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaTorsionTorsionForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaTorsionTorsionForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaTorsionTorsionForce& force);
// /**
// * 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:
// int numTorsionTorsions;
// int numTorsionTorsionGrids;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaMultipoleForceKernel : public CalcAmoebaMultipoleForceKernel {
// public:
// ReferenceCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaMultipoleForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaMultipoleForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaMultipoleForce& force);
// /**
// * 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:
// int numMultipoles;
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
// public:
// ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaMultipoleForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force);
// /**
// * 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:
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
// public:
// ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaVdwForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaMultipoleForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaVdwForce& force);
// /**
// * 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:
// AmoebaReferenceData& data;
// System& system;
// };
//
// /**
// * This kernel is invoked to calculate the WCA dispersion forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaWcaDispersionForceKernel : public CalcAmoebaWcaDispersionForceKernel {
// public:
// ReferenceCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, AmoebaReferenceData& data, System& system);
// ~ReferenceCalcAmoebaWcaDispersionForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaMultipoleForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaWcaDispersionForce& force);
// /**
// * 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:
// AmoebaReferenceData& data;
// System& system;
// };
}
// namespace OpenMM
#endif
/*AMOEBA_OPENMM_REFERENCE_KERNELS_H*/
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.cpp
0 → 100644
View file @
4868314a
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "AmoebaReferenceForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Load delta of two vectors
@param xVector first vector
@param yVector second vector
@param deltaR output vector: y - x
--------------------------------------------------------------------------------------- */
void
AmoebaReferenceForce
::
loadDeltaR
(
const
RealOpenMM
*
xVector
,
const
RealOpenMM
*
yVector
,
std
::
vector
<
RealOpenMM
>&
deltaR
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::loadDeltaR";
// ---------------------------------------------------------------------------------------
deltaR
.
resize
(
0
);
deltaR
.
push_back
(
yVector
[
0
]
-
xVector
[
0
]
);
deltaR
.
push_back
(
yVector
[
1
]
-
xVector
[
1
]
);
deltaR
.
push_back
(
yVector
[
2
]
-
xVector
[
2
]
);
}
/**---------------------------------------------------------------------------------------
Calculate norm squared of 3d vector
@param inputVector vector whose norm squared is to be computed
@return norm squared
--------------------------------------------------------------------------------------- */
RealOpenMM
AmoebaReferenceForce
::
getNormSquared3
(
const
std
::
vector
<
RealOpenMM
>&
inputVector
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getNorm3";
// ---------------------------------------------------------------------------------------
// get 3 norm
return
(
inputVector
[
0
]
*
inputVector
[
0
]
+
inputVector
[
1
]
*
inputVector
[
1
]
+
inputVector
[
2
]
*
inputVector
[
2
]
);
}
/**---------------------------------------------------------------------------------------
Calculate norm of 3d vector
@param inputVector vector whose norm is to be computed
@return norm
--------------------------------------------------------------------------------------- */
RealOpenMM
AmoebaReferenceForce
::
getNorm3
(
const
std
::
vector
<
RealOpenMM
>&
inputVector
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getNorm3";
// ---------------------------------------------------------------------------------------
// get 3 norm
return
SQRT
(
inputVector
[
0
]
*
inputVector
[
0
]
+
inputVector
[
1
]
*
inputVector
[
1
]
+
inputVector
[
2
]
*
inputVector
[
2
]
);
}
/**---------------------------------------------------------------------------------------
Calculate dot product of 3d vectors
@param xVector first vector
@param yVector second vector
@return dot product
--------------------------------------------------------------------------------------- */
RealOpenMM
AmoebaReferenceForce
::
getDotProduct3
(
const
std
::
vector
<
RealOpenMM
>&
xVector
,
const
std
::
vector
<
RealOpenMM
>&
yVector
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
// ---------------------------------------------------------------------------------------
// get dot product
return
xVector
[
0
]
*
yVector
[
0
]
+
xVector
[
1
]
*
yVector
[
1
]
+
xVector
[
2
]
*
yVector
[
2
];
}
/**---------------------------------------------------------------------------------------
Calculate z = x X y
@param xVector input vector
@param yVector input vector
@param zVector output vector: z = x X y
--------------------------------------------------------------------------------------- */
void
AmoebaReferenceForce
::
getCrossProduct
(
const
std
::
vector
<
RealOpenMM
>&
xVector
,
const
std
::
vector
<
RealOpenMM
>&
yVector
,
std
::
vector
<
RealOpenMM
>&
zVector
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getCrossProduct";
// ---------------------------------------------------------------------------------------
zVector
[
0
]
=
xVector
[
1
]
*
yVector
[
2
]
-
xVector
[
2
]
*
yVector
[
1
];
zVector
[
1
]
=
xVector
[
2
]
*
yVector
[
0
]
-
xVector
[
0
]
*
yVector
[
2
];
zVector
[
2
]
=
xVector
[
0
]
*
yVector
[
1
]
-
xVector
[
1
]
*
yVector
[
0
];
return
;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.h
0 → 100644
View file @
4868314a
/* 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 __AmoebaReferenceHarmonicBondForce_H__
#define __AmoebaReferenceHarmonicBondForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
class
AmoebaReferenceForce
{
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceForce
(
);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~
AmoebaReferenceForce
(
);
/**---------------------------------------------------------------------------------------
Load delta of two vectors
@param xVector first vector
@param yVector second vector
@param deltaR output vector: y - x
--------------------------------------------------------------------------------------- */
static
void
loadDeltaR
(
const
RealOpenMM
*
xVector
,
const
RealOpenMM
*
yVector
,
std
::
vector
<
RealOpenMM
>&
deltaR
);
/**---------------------------------------------------------------------------------------
Calculate norm squared of 3d vector
@param inputVector vector whose norm squared is to be computed
@return norm squared
--------------------------------------------------------------------------------------- */
static
RealOpenMM
getNormSquared3
(
const
std
::
vector
<
RealOpenMM
>&
inputVector
);
/**---------------------------------------------------------------------------------------
Calculate norm of 3d vector
@param inputVector vector whose norm is to be computed
@return norm
--------------------------------------------------------------------------------------- */
static
RealOpenMM
getNorm3
(
const
std
::
vector
<
RealOpenMM
>&
inputVector
);
/**---------------------------------------------------------------------------------------
Calculate dot product of 3d vectors
@param xVector first vector
@param yVector second vector
@return dot product
--------------------------------------------------------------------------------------- */
static
RealOpenMM
getDotProduct3
(
const
std
::
vector
<
RealOpenMM
>&
xVector
,
const
std
::
vector
<
RealOpenMM
>&
yVector
);
/**---------------------------------------------------------------------------------------
Calculate z = x X y
@param xVector input vector
@param yVector input vector
@param zVector output vector: z = x X y
--------------------------------------------------------------------------------------- */
static
void
getCrossProduct
(
const
std
::
vector
<
RealOpenMM
>&
xVector
,
const
std
::
vector
<
RealOpenMM
>&
yVector
,
std
::
vector
<
RealOpenMM
>&
zVector
);
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceHarmonicBondForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.cpp
0 → 100644
View file @
4868314a
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceHarmonicAngleForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Get dEdT and energy prefactor given cosine of angle :: the calculation for different
the force types is identical
@param cosine cosine of angle
@param idealAngle ideal angle
@param angleK angle k (quadratic prefactor)
@param angleCubic cubic prefactor
@param angleQuartic quartiic prefactor
@param anglePentic pentic prefactor
@param angleSextic sextic prefactor
@param dEdR dEdR
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM
AmoebaReferenceHarmonicAngleForce
::
getPrefactorsGivenAngleCosine
(
RealOpenMM
cosine
,
RealOpenMM
idealAngle
,
RealOpenMM
angleK
,
RealOpenMM
angleCubic
,
RealOpenMM
angleQuartic
,
RealOpenMM
anglePentic
,
RealOpenMM
angleSextic
,
RealOpenMM
*
dEdR
){
// ---------------------------------------------------------------------------------------
static
const
RealOpenMM
zero
=
0.0
;
static
const
RealOpenMM
one
=
1.0
;
static
const
RealOpenMM
two
=
2.0
;
static
const
RealOpenMM
three
=
3.0
;
static
const
RealOpenMM
four
=
4.0
;
static
const
RealOpenMM
five
=
5.0
;
static
const
RealOpenMM
six
=
6.0
;
// static const std::string methodName = "AmoebaReferenceHarmonicAngleForce::getPrefactorsGivenAngleCosine";
// ---------------------------------------------------------------------------------------
RealOpenMM
angle
;
if
(
cosine
>=
one
){
angle
=
zero
;
}
else
if
(
cosine
<=
-
one
){
angle
=
RADIAN
*
PI_M
;
}
else
{
angle
=
RADIAN
*
ACOS
(
cosine
);
}
RealOpenMM
deltaIdeal
=
angle
-
idealAngle
;
RealOpenMM
deltaIdeal2
=
deltaIdeal
*
deltaIdeal
;
RealOpenMM
deltaIdeal3
=
deltaIdeal
*
deltaIdeal2
;
RealOpenMM
deltaIdeal4
=
deltaIdeal2
*
deltaIdeal2
;
*
dEdR
=
(
two
+
three
*
angleCubic
*
deltaIdeal
+
four
*
angleQuartic
*
deltaIdeal2
+
five
*
anglePentic
*
deltaIdeal3
+
six
*
angleSextic
*
deltaIdeal4
);
*
dEdR
*=
RADIAN
*
angleK
*
deltaIdeal
;
RealOpenMM
energy
=
1.0
f
+
angleCubic
*
deltaIdeal
+
angleQuartic
*
deltaIdeal2
+
anglePentic
*
deltaIdeal3
+
angleSextic
*
deltaIdeal4
;
energy
*=
angleK
*
deltaIdeal2
;
return
energy
;
}
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param angleLength angle
@param angleK quadratic angle force
@param angleCubic cubic angle force parameter
@param angleQuartic quartic angle force parameter
@param anglePentic pentic angle force parameter
@param angleSextic sextic angle force parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM
AmoebaReferenceHarmonicAngleForce
::
calculateForceAndEnergy
(
RealOpenMM
*
positionAtomA
,
RealOpenMM
*
positionAtomB
,
RealOpenMM
*
positionAtomC
,
RealOpenMM
angle
,
RealOpenMM
angleK
,
RealOpenMM
angleCubic
,
RealOpenMM
angleQuartic
,
RealOpenMM
anglePentic
,
RealOpenMM
angleSextic
,
RealOpenMM
**
forces
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceHarmonicAngleForce::calculateHarmonicForce";
static
const
RealOpenMM
zero
=
0.0
;
static
const
RealOpenMM
one
=
1.0
;
static
const
RealOpenMM
onePt5
=
1.5
;
static
const
RealOpenMM
two
=
2.0
;
// ---------------------------------------------------------------------------------------
std
::
vector
<
RealOpenMM
>
deltaR
[
2
];
AmoebaReferenceForce
::
loadDeltaR
(
positionAtomA
,
positionAtomB
,
deltaR
[
0
]
);
RealOpenMM
rAB2
=
AmoebaReferenceForce
::
getNormSquared3
(
deltaR
[
0
]
);
RealOpenMM
rAB
=
SQRT
(
rAB2
);
AmoebaReferenceForce
::
loadDeltaR
(
positionAtomC
,
positionAtomB
,
deltaR
[
1
]
);
RealOpenMM
rCB2
=
AmoebaReferenceForce
::
getNormSquared3
(
deltaR
[
1
]
);
RealOpenMM
rCB
=
SQRT
(
rCB2
);
if
(
rAB
<=
zero
||
rCB
<=
zero
){
return
zero
;
}
std
::
vector
<
RealOpenMM
>
pVector
(
3
);
AmoebaReferenceForce
::
getCrossProduct
(
deltaR
[
0
],
deltaR
[
1
],
pVector
);
RealOpenMM
rp
=
AmoebaReferenceForce
::
getNorm3
(
pVector
);
if
(
rp
<
1.0e-06
){
rp
=
1.0e-06
;
}
RealOpenMM
dot
=
AmoebaReferenceForce
::
getDotProduct3
(
deltaR
[
0
],
deltaR
[
1
]
);
RealOpenMM
cosine
=
dot
/
(
rAB
*
rCB
);
RealOpenMM
dEdR
;
RealOpenMM
energy
=
getPrefactorsGivenAngleCosine
(
cosine
,
angle
,
angleK
,
angleCubic
,
angleQuartic
,
anglePentic
,
angleSextic
,
&
dEdR
);
RealOpenMM
termA
=
dEdR
/
(
rAB2
*
rp
);
RealOpenMM
termC
=
-
dEdR
/
(
rCB2
*
rp
);
std
::
vector
<
RealOpenMM
>
deltaCrossP
[
3
];
deltaCrossP
[
0
].
resize
(
3
);
deltaCrossP
[
1
].
resize
(
3
);
deltaCrossP
[
2
].
resize
(
3
);
AmoebaReferenceForce
::
getCrossProduct
(
deltaR
[
0
],
pVector
,
deltaCrossP
[
0
]
);
AmoebaReferenceForce
::
getCrossProduct
(
deltaR
[
1
],
pVector
,
deltaCrossP
[
2
]
);
for
(
unsigned
int
ii
=
0
;
ii
<
3
;
ii
++
){
deltaCrossP
[
0
][
ii
]
*=
termA
;
deltaCrossP
[
2
][
ii
]
*=
termC
;
deltaCrossP
[
1
][
ii
]
=
-
1.0
f
*
(
deltaCrossP
[
0
][
ii
]
+
deltaCrossP
[
2
][
ii
]);
}
// accumulate forces
for
(
int
jj
=
0
;
jj
<
3
;
jj
++
){
forces
[
jj
][
0
]
+=
deltaCrossP
[
jj
][
0
];
forces
[
jj
][
1
]
+=
deltaCrossP
[
jj
][
1
];
forces
[
jj
][
2
]
+=
deltaCrossP
[
jj
][
2
];
}
return
energy
;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.h
0 → 100644
View file @
4868314a
/* 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 __AmoebaReferenceHarmonicAngleForce_H__
#define __AmoebaReferenceHarmonicAngleForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
// ---------------------------------------------------------------------------------------
class
AmoebaReferenceHarmonicAngleForce
{
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceHarmonicAngleForce
(
);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~
AmoebaReferenceHarmonicAngleForce
(
);
/**---------------------------------------------------------------------------------------
Get dEdT and energy prefactor given cosine of angle :: the calculation for different
the force types is identical
@param cosine cosine of angle
@param idealAngle ideal angle
@param angleK angle k (quadratic prefactor)
@param angleCubic cubic prefactor
@param angleQuartic quartiic prefactor
@param anglePentic pentic prefactor
@param angleSextic sextic prefactor
@param dEdR dEdR
@return energy
--------------------------------------------------------------------------------------- */
static
RealOpenMM
getPrefactorsGivenAngleCosine
(
RealOpenMM
cosine
,
RealOpenMM
idealAngle
,
RealOpenMM
angleK
,
RealOpenMM
angleCubic
,
RealOpenMM
angleQuartic
,
RealOpenMM
anglePentic
,
RealOpenMM
angleSextic
,
RealOpenMM
*
dEdR
);
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param angleLength angle
@param angleK quadratic angle force
@param angleCubic cubic angle force parameter
@param angleQuartic quartic angle force parameter
@param anglePentic pentic angle force parameter
@param angleSextic sextic angle force parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
static
RealOpenMM
calculateForceAndEnergy
(
RealOpenMM
*
positionAtomA
,
RealOpenMM
*
positionAtomB
,
RealOpenMM
*
positionAtomC
,
RealOpenMM
angle
,
RealOpenMM
angleK
,
RealOpenMM
angleCubic
,
RealOpenMM
angleQuartic
,
RealOpenMM
anglePentic
,
RealOpenMM
angleSextic
,
RealOpenMM
**
forces
);
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceHarmonicAngleForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.cpp
0 → 100644
View file @
4868314a
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "AmoebaReferenceHarmonicBondForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic bond ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param bondLength bond length
@param bondK bond force
@param bondCubic cubic bond force parameter
@param bondQuartic quartic bond force parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM
AmoebaReferenceHarmonicBondForce
::
calculateForceAndEnergy
(
RealOpenMM
*
positionAtomA
,
RealOpenMM
*
positionAtomB
,
RealOpenMM
bondLength
,
RealOpenMM
bondK
,
RealOpenMM
bondCubic
,
RealOpenMM
bondQuartic
,
RealOpenMM
**
forces
){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceHarmonicBondForce::calculateHarmonicForce";
static
const
RealOpenMM
zero
=
0.0
;
static
const
RealOpenMM
one
=
1.0
;
static
const
RealOpenMM
onePt5
=
1.5
;
static
const
RealOpenMM
two
=
2.0
;
// ---------------------------------------------------------------------------------------
// get deltaR, R2, and R between 2 atoms
std
::
vector
<
RealOpenMM
>
deltaR
;
deltaR
.
push_back
(
positionAtomB
[
0
]
-
positionAtomA
[
0
]
);
deltaR
.
push_back
(
positionAtomB
[
1
]
-
positionAtomA
[
1
]
);
deltaR
.
push_back
(
positionAtomB
[
2
]
-
positionAtomA
[
2
]
);
RealOpenMM
r
=
SQRT
(
deltaR
[
0
]
*
deltaR
[
0
]
+
deltaR
[
1
]
*
deltaR
[
1
]
+
deltaR
[
2
]
*
deltaR
[
2
]
);
// deltaIdeal = r - r_0
RealOpenMM
deltaIdeal
=
r
-
bondLength
;
RealOpenMM
deltaIdeal2
=
deltaIdeal
*
deltaIdeal
;
RealOpenMM
dEdR
=
(
one
+
onePt5
*
bondCubic
*
deltaIdeal
+
two
*
bondQuartic
*
deltaIdeal2
);
dEdR
*=
two
*
bondK
*
deltaIdeal
;
dEdR
=
r
>
zero
?
(
dEdR
/
r
)
:
zero
;
forces
[
0
][
0
]
+=
dEdR
*
deltaR
[
0
];
forces
[
0
][
1
]
+=
dEdR
*
deltaR
[
1
];
forces
[
0
][
2
]
+=
dEdR
*
deltaR
[
2
];
forces
[
1
][
0
]
-=
dEdR
*
deltaR
[
0
];
forces
[
1
][
1
]
-=
dEdR
*
deltaR
[
1
];
forces
[
1
][
2
]
-=
dEdR
*
deltaR
[
2
];
RealOpenMM
energy
=
bondK
*
deltaIdeal2
*
(
one
+
bondCubic
*
deltaIdeal
+
bondQuartic
*
deltaIdeal2
);
return
energy
;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.h
0 → 100644
View file @
4868314a
/* 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 __AmoebaReferenceHarmonicBondForce_H__
#define __AmoebaReferenceHarmonicBondForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
// ---------------------------------------------------------------------------------------
class
AmoebaReferenceHarmonicBondForce
{
private:
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceHarmonicBondForce
(
){};
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~
AmoebaReferenceHarmonicBondForce
(
);
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic bond ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param bondLength bond length
@param bondK bond force
@param bondCubic cubic bond force parameter
@param bondQuartic quartic bond force parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
static
RealOpenMM
calculateForceAndEnergy
(
RealOpenMM
*
positionAtomA
,
RealOpenMM
*
positionAtomB
,
RealOpenMM
bondLength
,
RealOpenMM
bondK
,
RealOpenMM
bondCubic
,
RealOpenMM
bondQuartic
,
RealOpenMM
**
forces
);
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceHarmonicBondForce___
plugins/amoeba/platforms/reference/tests/CMakeLists.txt
0 → 100644
View file @
4868314a
#
# Testing
#
ENABLE_TESTING
()
#INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/platforms/reference/include
)
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/openmmapi/include/openmm
)
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/platforms/reference/src
)
INCLUDE_DIRECTORIES
(
${
OPENMM_DIR
}
/platforms/reference/src/kernels
)
#SET(SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET "AmoebaTinkerParameterFile" )
#SET(AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES "AmoebaTinkerParameterFile.cpp" )
#SET(AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES "AmoebaTinkerParameterFile.h" )
Set
(
SHARED_OPENMM__AMOEBA_TARGET OpenMMAmoeba
)
Set
(
STATIC_OPENMM_TARGET OpenMMAmoeba_static
)
Set
(
SHARED_CUDA_TARGET OpenMMAmoebaCuda
)
Set
(
STATIC_CUDA_TARGET OpenMMCuda_static OpenMMAmoebaCuda_static
)
#ADD_LIBRARY(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} SHARED ${AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES} ${AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES} )
#SET_TARGET_PROPERTIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_VALIDATE_BUILDING_SHARED_LIBRARY")
#TARGET_LINK_LIBRARIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME})
#TARGET_LINK_LIBRARIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} debug ${SHARED_OPENMM__AMOEBA_TARGET}_d optimized ${SHARED_OPENMM__AMOEBA_TARGET})
#TARGET_LINK_LIBRARIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} debug ${SHARED_CUDA_TARGET}_d optimized ${SHARED_CUDA_TARGET})
IF
(
UNIX AND CMAKE_BUILD_TYPE MATCHES Debug
)
SET
(
SHARED_CUDA_TARGET
${
SHARED_CUDA_TARGET
}
_d
)
SET
(
SHARED_OPENMM__AMOEBA_TARGET
${
SHARED_OPENMM__AMOEBA_TARGET
}
_d
)
#SET(STATIC_CUDA_TARGET ${STATIC_CUDA_TARGET}_d)
#Set(STATIC_OPENMM_TARGET ${STATIC_OPENMM_TARGET}_d)
ENDIF
(
UNIX AND CMAKE_BUILD_TYPE MATCHES Debug
)
#LINK_DIRECTORIES
# Automatically create tests using files named "Test*.cpp"
FILE
(
GLOB TEST_PROGS
"*Test*.cpp"
)
FOREACH
(
TEST_PROG
${
TEST_PROGS
}
)
GET_FILENAME_COMPONENT
(
TEST_ROOT
${
TEST_PROG
}
NAME_WE
)
# Link with shared library
ADD_EXECUTABLE
(
${
TEST_ROOT
}
${
TEST_PROG
}
)
#TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET} ${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET})
TARGET_LINK_LIBRARIES
(
${
TEST_ROOT
}
${
SHARED_TARGET
}
${
SHARED_OPENMM_TARGET
}
${
SHARED_CUDA_TARGET
}
)
ADD_TEST
(
${
TEST_ROOT
}
${
EXECUTABLE_OUTPUT_PATH
}
/
${
TEST_ROOT
}
)
# Link with static library
# SET(TEST_STATIC ${TEST_ROOT}Static)
# CUDA_ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
# SET_TARGET_PROPERTIES(${TEST_STATIC}
# PROPERTIES
# COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
# )
# TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET} ${STATIC_OPENMM_TARGET} ${STATIC_CUDA_TARGET})
# ADD_TEST(${TEST_STATIC} ${EXECUTABLE_OUTPUT_PATH}/${TEST_STATIC})
ENDFOREACH
(
TEST_PROG
${
TEST_PROGS
}
)
# TestCudaUsingParameterFile customized w/ command-line argument (input file name used in test)
#ADD_EXECUTABLE(TestAmoebaCudaUsingParameterFile TstAmoebaCudaUsingParameterFile.cpp)
#TARGET_LINK_LIBRARIES(TestAmoebaCudaUsingParameterFile ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET})
#ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/lambdaSdObcParameters.txt")
#ADD_TEST(TestCudaUsingParameterFile "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFile" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfNoPbcParameters.txt")
#
#SET(TEST_ROOT TestCudaUsingParameterFile)
#SET(TEST_PROG TstCudaUsingParameterFile.cpp)
#SET(TEST_STATIC ${TEST_ROOT}Static)
#SET(INCLUDE_CUDA_STATIC 1)
#IF(INCLUDE_CUDA_STATIC)
# ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
# SET_TARGET_PROPERTIES(${TEST_STATIC}
# PROPERTIES
# COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
# )
# TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET} ${STATIC_BROOK_TARGET})
# ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/lambdaSdObcParameters.txt")
# ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfNoPbcParameters.txt")
# ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfPbcParameters.txt" " +checkEnergyForceConsistent -checkForces" )
#ENDIF(INCLUDE_CUDA_STATIC)
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaHarmonicAngleForce.cpp
0 → 100644
View file @
4868314a
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
/**
* This tests the Cuda implementation of CudaAmoebaHarmonicAngleForce.
*/
#include "../../../tests/AssertionUtilities.h"
//#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "AmoebaOpenMM.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
using
namespace
OpenMM
;
const
double
TOL
=
1e-5
;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
#define RADIAN_TO_DEGREE 57.29577951308
#define DEGREE_TO_RADIAN 0.01745329252
#define RADIAN_INVERSE 0.01745329252
/* ---------------------------------------------------------------------------------------
Compute cross product of two 3-vectors and place in 3rd vector
vectorZ = vectorX x vectorY
@param vectorX x-vector
@param vectorY y-vector
@param vectorZ z-vector
@return vector is vectorZ
--------------------------------------------------------------------------------------- */
static
void
crossProductVector3
(
double
*
vectorX
,
double
*
vectorY
,
double
*
vectorZ
){
vectorZ
[
0
]
=
vectorX
[
1
]
*
vectorY
[
2
]
-
vectorX
[
2
]
*
vectorY
[
1
];
vectorZ
[
1
]
=
vectorX
[
2
]
*
vectorY
[
0
]
-
vectorX
[
0
]
*
vectorY
[
2
];
vectorZ
[
2
]
=
vectorX
[
0
]
*
vectorY
[
1
]
-
vectorX
[
1
]
*
vectorY
[
0
];
return
;
}
static
void
getPrefactorsGivenAngleCosine
(
double
cosine
,
double
idealAngle
,
double
quadraticK
,
double
cubicK
,
double
quarticK
,
double
penticK
,
double
sexticK
,
double
*
dEdR
,
double
*
energyTerm
,
FILE
*
log
)
{
double
angle
;
if
(
cosine
>=
1.0
){
angle
=
0.0
f
;
}
else
if
(
cosine
<=
-
1.0
){
angle
=
RADIAN
*
PI_M
;
}
else
{
angle
=
RADIAN
*
acos
(
cosine
);
}
if
(
log
){
(
void
)
fprintf
(
log
,
"getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e
\n
"
,
cosine
,
angle
,
idealAngle
);
(
void
)
fflush
(
log
);
}
double
deltaIdeal
=
angle
-
idealAngle
;
double
deltaIdeal2
=
deltaIdeal
*
deltaIdeal
;
double
deltaIdeal3
=
deltaIdeal
*
deltaIdeal2
;
double
deltaIdeal4
=
deltaIdeal2
*
deltaIdeal2
;
// deltaIdeal = r - r_0
*
dEdR
=
(
2.0
+
3.0
*
cubicK
*
deltaIdeal
+
4.0
*
quarticK
*
deltaIdeal2
+
5.0
*
penticK
*
deltaIdeal3
+
6.0
*
sexticK
*
deltaIdeal4
);
*
dEdR
*=
RADIAN
*
quadraticK
*
deltaIdeal
;
*
energyTerm
=
1.0
f
+
cubicK
*
deltaIdeal
+
quarticK
*
deltaIdeal2
+
penticK
*
deltaIdeal3
+
sexticK
*
deltaIdeal4
;
*
energyTerm
*=
quadraticK
*
deltaIdeal2
;
return
;
}
static
void
computeAmoebaHarmonicAngleForce
(
int
bondIndex
,
std
::
vector
<
Vec3
>&
positions
,
AmoebaHarmonicAngleForce
&
amoebaHarmonicAngleForce
,
std
::
vector
<
Vec3
>&
forces
,
double
*
energy
,
FILE
*
log
)
{
int
particle1
,
particle2
,
particle3
;
double
idealAngle
;
double
quadraticK
;
amoebaHarmonicAngleForce
.
getAngleParameters
(
bondIndex
,
particle1
,
particle2
,
particle3
,
idealAngle
,
quadraticK
);
double
cubicK
=
amoebaHarmonicAngleForce
.
getAmoebaGlobalHarmonicAngleCubic
();
double
quarticK
=
amoebaHarmonicAngleForce
.
getAmoebaGlobalHarmonicAngleQuartic
();
double
penticK
=
amoebaHarmonicAngleForce
.
getAmoebaGlobalHarmonicAnglePentic
();
double
sexticK
=
amoebaHarmonicAngleForce
.
getAmoebaGlobalHarmonicAngleSextic
();
if
(
log
){
(
void
)
fprintf
(
log
,
"computeAmoebaHarmonicAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]
\n
"
,
bondIndex
,
particle1
,
particle2
,
particle3
,
idealAngle
,
quadraticK
,
cubicK
,
quarticK
,
penticK
,
sexticK
);
(
void
)
fflush
(
log
);
}
double
deltaR
[
2
][
3
];
double
r2_0
=
0.0
;
double
r2_1
=
0.0
;
for
(
int
ii
=
0
;
ii
<
3
;
ii
++
){
deltaR
[
0
][
ii
]
=
positions
[
particle1
][
ii
]
-
positions
[
particle2
][
ii
];
r2_0
+=
deltaR
[
0
][
ii
]
*
deltaR
[
0
][
ii
];
deltaR
[
1
][
ii
]
=
positions
[
particle3
][
ii
]
-
positions
[
particle2
][
ii
];
r2_1
+=
deltaR
[
1
][
ii
]
*
deltaR
[
1
][
ii
];
}
double
pVector
[
3
];
crossProductVector3
(
deltaR
[
0
],
deltaR
[
1
],
pVector
);
double
rp
=
sqrt
(
pVector
[
0
]
*
pVector
[
0
]
+
pVector
[
1
]
*
pVector
[
1
]
+
pVector
[
2
]
*
pVector
[
2
]
);
if
(
rp
<
1.0e-06
){
rp
=
1.0e-06
;
}
double
dot
=
deltaR
[
0
][
0
]
*
deltaR
[
1
][
0
]
+
deltaR
[
0
][
1
]
*
deltaR
[
1
][
1
]
+
deltaR
[
0
][
2
]
*
deltaR
[
1
][
2
];
double
cosine
=
dot
/
sqrt
(
r2_0
*
r2_1
);
if
(
log
){
(
void
)
fprintf
(
log
,
"dot=%10.3e r2_0=%10.3e r2_1=%10.3e
\n
"
,
dot
,
r2_0
,
r2_1
);
(
void
)
fflush
(
log
);
}
double
dEdR
;
double
energyTerm
;
getPrefactorsGivenAngleCosine
(
cosine
,
idealAngle
,
quadraticK
,
cubicK
,
quarticK
,
penticK
,
sexticK
,
&
dEdR
,
&
energyTerm
,
log
);
double
termA
=
-
dEdR
/
(
r2_0
*
rp
);
double
termC
=
dEdR
/
(
r2_1
*
rp
);
double
deltaCrossP
[
3
][
3
];
crossProductVector3
(
deltaR
[
0
],
pVector
,
deltaCrossP
[
0
]
);
crossProductVector3
(
deltaR
[
1
],
pVector
,
deltaCrossP
[
2
]
);
for
(
int
ii
=
0
;
ii
<
3
;
ii
++
){
deltaCrossP
[
0
][
ii
]
*=
termA
;
deltaCrossP
[
2
][
ii
]
*=
termC
;
deltaCrossP
[
1
][
ii
]
=
-
1.0
*
(
deltaCrossP
[
0
][
ii
]
+
deltaCrossP
[
2
][
ii
]);
}
forces
[
particle1
][
0
]
+=
deltaCrossP
[
0
][
0
];
forces
[
particle1
][
1
]
+=
deltaCrossP
[
0
][
1
];
forces
[
particle1
][
2
]
+=
deltaCrossP
[
0
][
2
];
forces
[
particle2
][
0
]
+=
deltaCrossP
[
1
][
0
];
forces
[
particle2
][
1
]
+=
deltaCrossP
[
1
][
1
];
forces
[
particle2
][
2
]
+=
deltaCrossP
[
1
][
2
];
forces
[
particle3
][
0
]
+=
deltaCrossP
[
2
][
0
];
forces
[
particle3
][
1
]
+=
deltaCrossP
[
2
][
1
];
forces
[
particle3
][
2
]
+=
deltaCrossP
[
2
][
2
];
*
energy
+=
energyTerm
;
}
static
void
computeAmoebaHarmonicAngleForces
(
Context
&
context
,
AmoebaHarmonicAngleForce
&
amoebaHarmonicAngleForce
,
std
::
vector
<
Vec3
>&
expectedForces
,
double
*
expectedEnergy
,
FILE
*
log
)
{
// get positions and zero forces
State
state
=
context
.
getState
(
State
::
Positions
);
std
::
vector
<
Vec3
>
positions
=
state
.
getPositions
();
expectedForces
.
resize
(
positions
.
size
()
);
for
(
unsigned
int
ii
=
0
;
ii
<
expectedForces
.
size
();
ii
++
){
expectedForces
[
ii
][
0
]
=
expectedForces
[
ii
][
1
]
=
expectedForces
[
ii
][
2
]
=
0.0
;
}
// calculates forces/energy
*
expectedEnergy
=
0.0
;
for
(
int
ii
=
0
;
ii
<
amoebaHarmonicAngleForce
.
getNumAngles
();
ii
++
){
computeAmoebaHarmonicAngleForce
(
ii
,
positions
,
amoebaHarmonicAngleForce
,
expectedForces
,
expectedEnergy
,
log
);
}
if
(
log
){
(
void
)
fprintf
(
log
,
"computeAmoebaHarmonicAngleForces: expected energy=%14.7e
\n
"
,
*
expectedEnergy
);
for
(
unsigned
int
ii
=
0
;
ii
<
positions
.
size
();
ii
++
){
(
void
)
fprintf
(
log
,
"%6u [%14.7e %14.7e %14.7e]
\n
"
,
ii
,
expectedForces
[
ii
][
0
],
expectedForces
[
ii
][
1
],
expectedForces
[
ii
][
2
]
);
}
(
void
)
fflush
(
log
);
}
return
;
}
void
compareWithExpectedForceAndEnergy
(
Context
&
context
,
AmoebaHarmonicAngleForce
&
amoebaHarmonicAngleForce
,
double
tolerance
,
const
std
::
string
&
idString
,
FILE
*
log
)
{
std
::
vector
<
Vec3
>
expectedForces
;
double
expectedEnergy
;
computeAmoebaHarmonicAngleForces
(
context
,
amoebaHarmonicAngleForce
,
expectedForces
,
&
expectedEnergy
,
log
);
State
state
=
context
.
getState
(
State
::
Forces
|
State
::
Energy
);
const
std
::
vector
<
Vec3
>
forces
=
state
.
getForces
();
if
(
log
){
(
void
)
fprintf
(
log
,
"computeAmoebaHarmonicAngleForces: expected energy=%14.7e %14.7e
\n
"
,
expectedEnergy
,
state
.
getPotentialEnergy
()
);
for
(
unsigned
int
ii
=
0
;
ii
<
forces
.
size
();
ii
++
){
(
void
)
fprintf
(
log
,
"%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]
\n
"
,
ii
,
expectedForces
[
ii
][
0
],
expectedForces
[
ii
][
1
],
expectedForces
[
ii
][
2
],
forces
[
ii
][
0
],
forces
[
ii
][
1
],
forces
[
ii
][
2
]
);
}
(
void
)
fflush
(
log
);
}
for
(
unsigned
int
ii
=
0
;
ii
<
forces
.
size
();
ii
++
){
ASSERT_EQUAL_VEC
(
expectedForces
[
ii
],
forces
[
ii
],
tolerance
);
}
ASSERT_EQUAL_TOL
(
expectedEnergy
,
state
.
getPotentialEnergy
(),
tolerance
);
}
void
testOneAngle
(
FILE
*
log
)
{
System
system
;
int
numberOfParticles
=
3
;
for
(
int
ii
=
0
;
ii
<
numberOfParticles
;
ii
++
){
system
.
addParticle
(
1.0
);
}
LangevinIntegrator
integrator
(
0.0
,
0.1
,
0.01
);
AmoebaHarmonicAngleForce
*
amoebaHarmonicAngleForce
=
new
AmoebaHarmonicAngleForce
();
double
angle
=
100.0
;
double
quadraticK
=
1.0
;
double
cubicK
=
1.0e-01
;
double
quarticK
=
1.0e-02
;
double
penticK
=
1.0e-03
;
double
sexticK
=
1.0e-04
;
amoebaHarmonicAngleForce
->
addAngle
(
0
,
1
,
2
,
angle
,
quadraticK
);
amoebaHarmonicAngleForce
->
setAmoebaGlobalHarmonicAngleCubic
(
cubicK
);
amoebaHarmonicAngleForce
->
setAmoebaGlobalHarmonicAngleQuartic
(
quarticK
);
amoebaHarmonicAngleForce
->
setAmoebaGlobalHarmonicAnglePentic
(
penticK
);
amoebaHarmonicAngleForce
->
setAmoebaGlobalHarmonicAngleSextic
(
sexticK
);
system
.
addForce
(
amoebaHarmonicAngleForce
);
Context
context
(
system
,
integrator
,
Platform
::
getPlatformByName
(
"Reference"
));
std
::
vector
<
Vec3
>
positions
(
numberOfParticles
);
positions
[
0
]
=
Vec3
(
0
,
1
,
0
);
positions
[
1
]
=
Vec3
(
0
,
0
,
0
);
positions
[
2
]
=
Vec3
(
0
,
0
,
1
);
context
.
setPositions
(
positions
);
compareWithExpectedForceAndEnergy
(
context
,
*
amoebaHarmonicAngleForce
,
TOL
,
"testOneAngle"
,
log
);
}
int
main
(
int
numberOfArguments
,
char
*
argv
[]
)
{
try
{
std
::
cout
<<
"TestCudaAmoebaHarmonicAngleForce running test..."
<<
std
::
endl
;
Platform
::
loadPluginsFromDirectory
(
Platform
::
getDefaultPluginsDirectory
()
);
//FILE* log = fopen( "AmoebaHarmonicAngleForce.log", "w" );;
//FILE* log = NULL;
FILE
*
log
=
stderr
;
testOneAngle
(
log
);
if
(
log
&&
log
!=
stderr
)
(
void
)
fclose
(
log
);
}
catch
(
const
std
::
exception
&
e
)
{
std
::
cout
<<
"exception: "
<<
e
.
what
()
<<
std
::
endl
;
std
::
cout
<<
"FAIL - ERROR. Test failed."
<<
std
::
endl
;
return
1
;
}
std
::
cout
<<
"PASS - Test succeeded."
<<
std
::
endl
;
return
0
;
}
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaHarmonicBondForce.cpp
0 → 100644
View file @
4868314a
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of HarmonicBondForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "AmoebaOpenMM.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
using
namespace
OpenMM
;
const
double
TOL
=
1e-5
;
static
void
computeAmoebaHarmonicBondForce
(
int
bondIndex
,
std
::
vector
<
Vec3
>&
positions
,
AmoebaHarmonicBondForce
&
amoebaHarmonicBondForce
,
std
::
vector
<
Vec3
>&
forces
,
double
*
energy
)
{
int
particle1
,
particle2
;
double
bondLength
;
double
quadraticK
;
double
cubicK
=
amoebaHarmonicBondForce
.
getAmoebaGlobalHarmonicBondCubic
();
double
quarticK
=
amoebaHarmonicBondForce
.
getAmoebaGlobalHarmonicBondQuartic
();
amoebaHarmonicBondForce
.
getBondParameters
(
bondIndex
,
particle1
,
particle2
,
bondLength
,
quadraticK
);
double
deltaR
[
3
];
double
r2
=
0.0
;
for
(
int
ii
=
0
;
ii
<
3
;
ii
++
){
deltaR
[
ii
]
=
positions
[
particle2
][
ii
]
-
positions
[
particle1
][
ii
];
r2
+=
deltaR
[
ii
]
*
deltaR
[
ii
];
}
double
r
=
sqrt
(
r2
);
double
bondDelta
=
(
r
-
bondLength
);
double
bondDelta2
=
bondDelta
*
bondDelta
;
double
dEdR
=
1.0
+
1.5
*
cubicK
*
bondDelta
+
2.0
*
quarticK
*
bondDelta2
;
dEdR
*=
(
r
>
0.0
)
?
(
2.0
*
quadraticK
*
bondDelta
)
/
r
:
0.0
;
forces
[
particle1
][
0
]
+=
dEdR
*
deltaR
[
0
];
forces
[
particle1
][
1
]
+=
dEdR
*
deltaR
[
1
];
forces
[
particle1
][
2
]
+=
dEdR
*
deltaR
[
2
];
forces
[
particle2
][
0
]
-=
dEdR
*
deltaR
[
0
];
forces
[
particle2
][
1
]
-=
dEdR
*
deltaR
[
1
];
forces
[
particle2
][
2
]
-=
dEdR
*
deltaR
[
2
];
*
energy
+=
(
1.0
f
+
cubicK
*
bondDelta
+
quarticK
*
bondDelta2
)
*
quadraticK
*
bondDelta2
;
}
static
void
computeAmoebaHarmonicBondForces
(
Context
&
context
,
AmoebaHarmonicBondForce
&
amoebaHarmonicBondForce
,
std
::
vector
<
Vec3
>&
expectedForces
,
double
*
expectedEnergy
,
FILE
*
log
)
{
// get positions and zero forces
State
state
=
context
.
getState
(
State
::
Positions
);
std
::
vector
<
Vec3
>
positions
=
state
.
getPositions
();
expectedForces
.
resize
(
positions
.
size
()
);
for
(
unsigned
int
ii
=
0
;
ii
<
expectedForces
.
size
();
ii
++
){
expectedForces
[
ii
][
0
]
=
expectedForces
[
ii
][
1
]
=
expectedForces
[
ii
][
2
]
=
0.0
;
}
// calculates forces/energy
*
expectedEnergy
=
0.0
;
for
(
int
ii
=
0
;
ii
<
amoebaHarmonicBondForce
.
getNumBonds
();
ii
++
){
computeAmoebaHarmonicBondForce
(
ii
,
positions
,
amoebaHarmonicBondForce
,
expectedForces
,
expectedEnergy
);
}
if
(
log
){
(
void
)
fprintf
(
log
,
"computeAmoebaHarmonicBondForces: expected energy=%15.7e
\n
"
,
*
expectedEnergy
);
for
(
unsigned
int
ii
=
0
;
ii
<
positions
.
size
();
ii
++
){
(
void
)
fprintf
(
log
,
"%6u [%15.7e %15.7e %15.7e]
\n
"
,
ii
,
expectedForces
[
ii
][
0
],
expectedForces
[
ii
][
1
],
expectedForces
[
ii
][
2
]
);
}
(
void
)
fflush
(
log
);
}
return
;
}
void
compareWithExpectedForceAndEnergy
(
Context
&
context
,
AmoebaHarmonicBondForce
&
amoebaHarmonicBondForce
,
double
tolerance
,
const
std
::
string
&
idString
,
FILE
*
log
)
{
std
::
vector
<
Vec3
>
expectedForces
;
double
expectedEnergy
;
computeAmoebaHarmonicBondForces
(
context
,
amoebaHarmonicBondForce
,
expectedForces
,
&
expectedEnergy
,
NULL
);
State
state
=
context
.
getState
(
State
::
Forces
|
State
::
Energy
);
const
std
::
vector
<
Vec3
>
forces
=
state
.
getForces
();
if
(
log
){
(
void
)
fprintf
(
log
,
"computeAmoebaHarmonicBondForces: expected energy=%15.7e %15.7e
\n
"
,
expectedEnergy
,
state
.
getPotentialEnergy
()
);
for
(
unsigned
int
ii
=
0
;
ii
<
forces
.
size
();
ii
++
){
(
void
)
fprintf
(
log
,
"%6u [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e]
\n
"
,
ii
,
expectedForces
[
ii
][
0
],
expectedForces
[
ii
][
1
],
expectedForces
[
ii
][
2
],
forces
[
ii
][
0
],
forces
[
ii
][
1
],
forces
[
ii
][
2
]
);
}
(
void
)
fflush
(
log
);
}
for
(
unsigned
int
ii
=
0
;
ii
<
forces
.
size
();
ii
++
){
ASSERT_EQUAL_VEC
(
expectedForces
[
ii
],
forces
[
ii
],
tolerance
);
}
ASSERT_EQUAL_TOL
(
expectedEnergy
,
state
.
getPotentialEnergy
(),
tolerance
);
}
void
testOneBond
(
FILE
*
log
)
{
System
system
;
system
.
addParticle
(
1.0
);
system
.
addParticle
(
1.0
);
LangevinIntegrator
integrator
(
0.0
,
0.1
,
0.01
);
AmoebaHarmonicBondForce
*
amoebaHarmonicBondForce
=
new
AmoebaHarmonicBondForce
();
double
bondLength
=
1.5
;
double
quadraticK
=
1.0
;
double
cubicK
=
2.0
;
double
quarticicK
=
3.0
;
amoebaHarmonicBondForce
->
setAmoebaGlobalHarmonicBondCubic
(
cubicK
);
amoebaHarmonicBondForce
->
setAmoebaGlobalHarmonicBondQuartic
(
quarticicK
);
amoebaHarmonicBondForce
->
addBond
(
0
,
1
,
bondLength
,
quadraticK
);
system
.
addForce
(
amoebaHarmonicBondForce
);
Context
context
(
system
,
integrator
,
Platform
::
getPlatformByName
(
"Reference"
));
std
::
vector
<
Vec3
>
positions
(
2
);
positions
[
0
]
=
Vec3
(
0
,
1
,
0
);
positions
[
1
]
=
Vec3
(
0
,
0
,
0
);
context
.
setPositions
(
positions
);
compareWithExpectedForceAndEnergy
(
context
,
*
amoebaHarmonicBondForce
,
TOL
,
"testOneBond"
,
log
);
}
void
testTwoBond
(
FILE
*
log
)
{
System
system
;
system
.
addParticle
(
1.0
);
system
.
addParticle
(
1.0
);
system
.
addParticle
(
1.0
);
LangevinIntegrator
integrator
(
0.0
,
0.1
,
0.01
);
AmoebaHarmonicBondForce
*
amoebaHarmonicBondForce
=
new
AmoebaHarmonicBondForce
();
double
bondLength
=
1.5
;
double
quadraticK
=
1.0
;
double
cubicK
=
2.0
;
double
quarticicK
=
3.0
;
amoebaHarmonicBondForce
->
setAmoebaGlobalHarmonicBondCubic
(
cubicK
);
amoebaHarmonicBondForce
->
setAmoebaGlobalHarmonicBondQuartic
(
quarticicK
);
amoebaHarmonicBondForce
->
addBond
(
0
,
1
,
bondLength
,
quadraticK
);
amoebaHarmonicBondForce
->
addBond
(
1
,
2
,
bondLength
,
quadraticK
);
system
.
addForce
(
amoebaHarmonicBondForce
);
Context
context
(
system
,
integrator
,
Platform
::
getPlatformByName
(
"Reference"
));
std
::
vector
<
Vec3
>
positions
(
3
);
positions
[
0
]
=
Vec3
(
0
,
1
,
0
);
positions
[
1
]
=
Vec3
(
0
,
0
,
0
);
positions
[
2
]
=
Vec3
(
1
,
0
,
1
);
context
.
setPositions
(
positions
);
compareWithExpectedForceAndEnergy
(
context
,
*
amoebaHarmonicBondForce
,
TOL
,
"testTwoBond"
,
log
);
}
int
main
(
int
numberOfArguments
,
char
*
argv
[]
)
{
try
{
std
::
cout
<<
"TestReferenceAmoebaHarmonicBondForce running test..."
<<
std
::
endl
;
Platform
::
loadPluginsFromDirectory
(
Platform
::
getDefaultPluginsDirectory
()
);
//FILE* log = NULL;
FILE
*
log
=
stderr
;
//testOneBond( log );
testTwoBond
(
log
);
if
(
log
&&
log
!=
stderr
)
(
void
)
fclose
(
log
);
}
catch
(
const
std
::
exception
&
e
)
{
std
::
cout
<<
"exception: "
<<
e
.
what
()
<<
std
::
endl
;
std
::
cout
<<
"FAIL - ERROR. Test failed."
<<
std
::
endl
;
return
1
;
}
std
::
cout
<<
"PASS - Test succeeded."
<<
std
::
endl
;
return
0
;
}
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment
