Skip to content

GitLab

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

Created OpenCL implementation of DrudeForce and DrudeLangevinIntegrator

parent 09a8475d
Hide whitespace changes
Inline Side-by-side
Showing with 1364 additions and 1 deletion
CMakeLists.txt
View file @ f8c98f32
...@@ -414,6 +414,15 @@ IF(OPENMM_BUILD_RPMD_PLUGIN) ...@@ -414,6 +414,15 @@ IF(OPENMM_BUILD_RPMD_PLUGIN)
ADD_SUBDIRECTORY(plugins/rpmd) ADD_SUBDIRECTORY(plugins/rpmd)
ENDIF(OPENMM_BUILD_RPMD_PLUGIN) ENDIF(OPENMM_BUILD_RPMD_PLUGIN)
# Drude plugin
SET(OPENMM_BUILD_DRUDE_PLUGIN ON CACHE BOOL "Build Drude plugin")
SET(OPENMM_BUILD_DRUDE_PATH)
IF(OPENMM_BUILD_DRUDE_PLUGIN)
SET(OPENMM_BUILD_DRUDE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/plugins/drude)
ADD_SUBDIRECTORY(plugins/drude)
ENDIF(OPENMM_BUILD_DRUDE_PLUGIN)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_TARGET}) INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_TARGET})
IF(OPENMM_BUILD_STATIC_LIB) IF(OPENMM_BUILD_STATIC_LIB)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_TARGET}) INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_TARGET})
...@@ -445,6 +454,9 @@ IF(OPENMM_BUILD_PYTHON_WRAPPERS) ...@@ -445,6 +454,9 @@ IF(OPENMM_BUILD_PYTHON_WRAPPERS)
IF(NOT OPENMM_BUILD_RPMD_PLUGIN) IF(NOT OPENMM_BUILD_RPMD_PLUGIN)
MESSAGE(SEND_ERROR "The Python wrappers require that the RPMD plugin be built.") MESSAGE(SEND_ERROR "The Python wrappers require that the RPMD plugin be built.")
ENDIF(NOT OPENMM_BUILD_RPMD_PLUGIN) ENDIF(NOT OPENMM_BUILD_RPMD_PLUGIN)
IF(NOT OPENMM_BUILD_DRUDE_PLUGIN)
MESSAGE(SEND_ERROR "The Python wrappers require that the Drude plugin be built.")
ENDIF(NOT OPENMM_BUILD_DRUDE_PLUGIN)
ADD_SUBDIRECTORY(wrappers/python) ADD_SUBDIRECTORY(wrappers/python)
ENDIF(OPENMM_BUILD_PYTHON_WRAPPERS) ENDIF(OPENMM_BUILD_PYTHON_WRAPPERS)
......
plugins/drude/platforms/cuda/src/CudaDrudeKernels.cpp
View file @ f8c98f32
...@@ -107,7 +107,7 @@ CudaCalcDrudeForceKernel::~CudaCalcDrudeForceKernel() { ...@@ -107,7 +107,7 @@ CudaCalcDrudeForceKernel::~CudaCalcDrudeForceKernel() {
if (pairParams != NULL) if (pairParams != NULL)
delete pairParams; delete pairParams;
} }
#include <iostream>
void CudaCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) { void CudaCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) {
cu.setAsCurrent(); cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size(); int numContexts = cu.getPlatformData().contexts.size();
......
plugins/drude/platforms/opencl/CMakeLists.txt 0 → 100644
View file @ f8c98f32
#---------------------------------------------------
# OpenMM OpenCL Drude Integrator
#
# Creates OpenMM library, base name=OpenMMDrudeOpenCL.
# Default libraries are shared & optimized. Variants
# are created for debug (_d).
#
# Windows:
# OpenMMDrudeOpenCL[_d].dll
# OpenMMDrudeOpenCL[_d].lib
# Unix:
# libOpenMMDrudeOpenCL[_d].so
#----------------------------------------------------
IF (APPLE)
SET (CMAKE_OSX_DEPLOYMENT_TARGET "10.6")
SET (CMAKE_OSX_SYSROOT "/Developer/SDKs/MacOSX10.6.sdk")
ENDIF (APPLE)
# The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS.
SET(OPENMM_SOURCE_SUBDIRS .)
# 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(OPENMMDRUDEOPENCL_LIBRARY_NAME OpenMMDrudeOpenCL)
SET(SHARED_TARGET ${OPENMMDRUDEOPENCL_LIBRARY_NAME})
# 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)
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_INCLUDE_DIRS) # start empty
FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
# append
SET(API_INCLUDE_DIRS ${API_INCLUDE_DIRS}
${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include
${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include/internal)
ENDFOREACH(subdir)
# We'll need both *relative* path names, starting with their API_INCLUDE_DIRS,
# and absolute pathnames.
SET(API_REL_INCLUDE_FILES) # start these out empty
SET(API_ABS_INCLUDE_FILES)
FOREACH(dir ${API_INCLUDE_DIRS})
FILE(GLOB fullpaths ${dir}/*.h) # returns full pathnames
SET(API_ABS_INCLUDE_FILES ${API_ABS_INCLUDE_FILES} ${fullpaths})
FOREACH(pathname ${fullpaths})
GET_FILENAME_COMPONENT(filename ${pathname} NAME)
SET(API_REL_INCLUDE_FILES ${API_REL_INCLUDE_FILES} ${dir}/${filename})
ENDFOREACH(pathname)
ENDFOREACH(dir)
# collect up source files
SET(SOURCE_FILES) # empty
SET(SOURCE_INCLUDE_FILES)
FOREACH(subdir ${OPENMM_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})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
ENDFOREACH(subdir)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/opencl/include)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/opencl/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_BINARY_DIR}/platforms/opencl/src)
# Set variables needed for encoding kernel sources into a C++ class
SET(CL_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
SET(CL_SOURCE_CLASS OpenCLDrudeKernelSources)
SET(CL_KERNELS_CPP ${CMAKE_CURRENT_BINARY_DIR}/src/${CL_SOURCE_CLASS}.cpp)
SET(CL_KERNELS_H ${CMAKE_CURRENT_BINARY_DIR}/src/${CL_SOURCE_CLASS}.h)
SET(SOURCE_FILES ${SOURCE_FILES} ${CL_KERNELS_CPP} ${CL_KERNELS_H})
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}/src)
# Create the library
INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR})
FILE(GLOB OPENCL_KERNELS ${CL_SOURCE_DIR}/kernels/*.cl)
ADD_CUSTOM_COMMAND(OUTPUT ${CL_KERNELS_CPP} ${CL_KERNELS_H}
COMMAND ${CMAKE_COMMAND}
ARGS -D CL_SOURCE_DIR=${CL_SOURCE_DIR} -D CL_KERNELS_CPP=${CL_KERNELS_CPP} -D CL_KERNELS_H=${CL_KERNELS_H} -D CL_SOURCE_CLASS=${CL_SOURCE_CLASS} -P ${CMAKE_SOURCE_DIR}/platforms/opencl/EncodeCLFiles.cmake
DEPENDS ${OPENCL_KERNELS}
)
SET_SOURCE_FILES_PROPERTIES(${CL_KERNELS_CPP} ${CL_KERNELS_H} PROPERTIES GENERATED TRUE)
ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME}_d)
ELSE (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME})
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
TARGET_LINK_LIBRARIES(${SHARED_TARGET} ${MAIN_OPENMM_LIB} ${OPENCL_LIBRARIES} ${PTHREADS_LIB})
TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}OpenCL_d optimized ${OPENMM_LIBRARY_NAME}OpenCL)
TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${SHARED_DRUDE_TARGET} optimized ${SHARED_DRUDE_TARGET})
SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY")
INSTALL(TARGETS ${SHARED_TARGET} DESTINATION ${CMAKE_INSTALL_PREFIX}/lib/plugins)
# Ensure that links to the main OpenCL library will be resolved.
IF (APPLE)
IF (CMAKE_BUILD_TYPE MATCHES Debug)
SET(OPENCL_LIBRARY libOpenMMOpenCL_d.dylib)
ELSE (CMAKE_BUILD_TYPE MATCHES Debug)
SET(OPENCL_LIBRARY libOpenMMOpenCL.dylib)
ENDIF (CMAKE_BUILD_TYPE MATCHES Debug)
INSTALL(CODE "EXECUTE_PROCESS(COMMAND install_name_tool -change ${OPENCL_LIBRARY} @loader_path/${OPENCL_LIBRARY} ${CMAKE_INSTALL_PREFIX}/lib/plugins/lib${SHARED_TARGET}.dylib)")
ENDIF (APPLE)
if(OPENMM_BUILD_OPENCL_TESTS)
SUBDIRS (tests)
endif(OPENMM_BUILD_OPENCL_TESTS)
plugins/drude/platforms/opencl/include/OpenCLDrudeKernelFactory.h 0 → 100644
View file @ f8c98f32
#ifndef OPENMM_OPENCLDRUDEKERNELFACTORY_H_
#define OPENMM_OPENCLDRUDEKERNELFACTORY_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/KernelFactory.h"
namespace OpenMM {
/**
* This KernelFactory creates kernels for the OpenCL implementation of the Drude plugin.
*/
class OpenCLDrudeKernelFactory : public KernelFactory {
public:
KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
};
} // namespace OpenMM
#endif /*OPENMM_OPENCLDRUDEKERNELFACTORY_H_*/
plugins/drude/platforms/opencl/src/OpenCLDrudeKernelFactory.cpp 0 → 100644
View file @ f8c98f32
/* -------------------------------------------------------------------------- *
* OpenMMAmoeba *
* -------------------------------------------------------------------------- *
* 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) 2011-2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* 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 <exception>
#include "OpenCLDrudeKernelFactory.h"
#include "OpenCLDrudeKernels.h"
#include "openmm/internal/windowsExport.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
extern "C" void registerPlatforms() {
}
extern "C" void registerKernelFactories() {
try {
Platform& platform = Platform::getPlatformByName("OpenCL");
OpenCLDrudeKernelFactory* factory = new OpenCLDrudeKernelFactory();
platform.registerKernelFactory(CalcDrudeForceKernel::Name(), factory);
platform.registerKernelFactory(IntegrateDrudeLangevinStepKernel::Name(), factory);
}
catch (std::exception ex) {
// Ignore
}
}
extern "C" OPENMM_EXPORT void registerDrudeOpenCLKernelFactories() {
try {
Platform::getPlatformByName("OpenCL");
}
catch (...) {
Platform::registerPlatform(new OpenCLPlatform());
}
registerKernelFactories();
}
KernelImpl* OpenCLDrudeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
OpenCLContext& cl = *static_cast<OpenCLPlatform::PlatformData*>(context.getPlatformData())->contexts[0];
if (name == CalcDrudeForceKernel::Name())
return new OpenCLCalcDrudeForceKernel(name, platform, cl);
if (name == IntegrateDrudeLangevinStepKernel::Name())
return new OpenCLIntegrateDrudeLangevinStepKernel(name, platform, cl);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
}
plugins/drude/platforms/opencl/src/OpenCLDrudeKernelSources.cpp.in 0 → 100644
View file @ f8c98f32
/* -------------------------------------------------------------------------- *
* 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) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* 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 "OpenCLDrudeKernelSources.h"
using namespace OpenMM;
using namespace std;
plugins/drude/platforms/opencl/src/OpenCLDrudeKernelSources.h.in 0 → 100644
View file @ f8c98f32
#ifndef OPENMM_OPENCLDRUDEKERNELSOURCES_H_
#define OPENMM_OPENCLDRUDEKERNELSOURCES_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* 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 <string>
namespace OpenMM {
/**
* This class is a central holding place for the source code of OpenCL kernels.
* The CMake build script inserts declarations into it based on the .cu files in the
* kernels subfolder.
*/
class OpenCLDrudeKernelSources {
public:
@CL_FILE_DECLARATIONS@
};
} // namespace OpenMM
#endif /*OPENMM_OPENCLDRUDEKERNELSOURCES_H_*/
plugins/drude/platforms/opencl/src/OpenCLDrudeKernels.cpp 0 → 100644
View file @ f8c98f32
/* -------------------------------------------------------------------------- *
* 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) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "OpenCLDrudeKernels.h"
#include "OpenCLDrudeKernelSources.h"
#include "openmm/internal/ContextImpl.h"
#include "OpenCLBondedUtilities.h"
#include "OpenCLForceInfo.h"
#include "OpenCLIntegrationUtilities.h"
#include "OpenCLKernelSources.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <set>
using namespace OpenMM;
using namespace std;
class OpenCLDrudeForceInfo : public OpenCLForceInfo {
public:
OpenCLDrudeForceInfo(const DrudeForce& force) : OpenCLForceInfo(0), force(force) {
}
int getNumParticleGroups() {
return force.getNumParticles()+force.getNumScreenedPairs();
}
void getParticlesInGroup(int index, vector<int>& particles) {
particles.clear();
if (index < force.getNumParticles()) {
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(index, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
particles.push_back(p);
particles.push_back(p1);
if (p2 != -1)
particles.push_back(p2);
if (p3 != -1)
particles.push_back(p3);
if (p4 != -1)
particles.push_back(p4);
}
else {
int drude1, drude2;
double thole;
force.getScreenedPairParameters(index-force.getNumParticles(), drude1, drude2, thole);
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(drude1, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
particles.push_back(p);
particles.push_back(p1);
force.getParticleParameters(drude2, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
particles.push_back(p);
particles.push_back(p1);
}
}
bool areGroupsIdentical(int group1, int group2) {
if (group1 < force.getNumParticles() && group2 < force.getNumParticles()) {
int p, p1, p2, p3, p4;
double charge1, polarizability1, aniso12_1, aniso34_1;
double charge2, polarizability2, aniso12_2, aniso34_2;
force.getParticleParameters(group1, p, p1, p2, p3, p4, charge1, polarizability1, aniso12_1, aniso34_1);
force.getParticleParameters(group2, p, p1, p2, p3, p4, charge2, polarizability2, aniso12_2, aniso34_2);
return (charge1 == charge2 && polarizability1 == polarizability2 && aniso12_1 == aniso12_2 && aniso34_1 == aniso34_2);
}
if (group1 >= force.getNumParticles() && group2 >= force.getNumParticles()) {
int drude1, drude2;
double thole1, thole2;
force.getScreenedPairParameters(group1-force.getNumParticles(), drude1, drude2, thole1);
force.getScreenedPairParameters(group1-force.getNumParticles(), drude1, drude2, thole2);
return (thole1 == thole2);
}
return false;
}
private:
const DrudeForce& force;
};
OpenCLCalcDrudeForceKernel::~OpenCLCalcDrudeForceKernel() {
if (particleParams != NULL)
delete particleParams;
if (pairParams != NULL)
delete pairParams;
}
void OpenCLCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) {
int numContexts = cl.getPlatformData().contexts.size();
int startParticleIndex = cl.getContextIndex()*force.getNumParticles()/numContexts;
int endParticleIndex = (cl.getContextIndex()+1)*force.getNumParticles()/numContexts;
int numParticles = endParticleIndex-startParticleIndex;
if (numParticles > 0) {
// Create the harmonic interaction .
vector<vector<int> > atoms(numParticles, vector<int>(5));
particleParams = OpenCLArray::create<mm_float4>(cl, numParticles, "drudeParticleParams");
vector<mm_float4> paramVector(numParticles);
for (int i = 0; i < numParticles; i++) {
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(startParticleIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], atoms[i][4], charge, polarizability, aniso12, aniso34);
double a1 = (atoms[i][2] == -1 ? 1 : aniso12);
double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34);
double a3 = 3-a1-a2;
double k3 = charge*charge/(polarizability*a3);
double k1 = charge*charge/(polarizability*a1) - k3;
double k2 = charge*charge/(polarizability*a2) - k3;
if (atoms[i][2] == -1) {
atoms[i][2] = 0;
k1 = 0;
}
if (atoms[i][3] == -1 || atoms[i][4] == -1) {
atoms[i][3] = 0;
atoms[i][4] = 0;
k2 = 0;
}
paramVector[i] = mm_float4((float) k1, (float) k2, (float) k3, 0.0f);
}
particleParams->upload(paramVector);
map<string, string> replacements;
replacements["PARAMS"] = cl.getBondedUtilities().addArgument(particleParams->getDeviceBuffer(), "float4");
cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLDrudeKernelSources::drudeParticleForce, replacements), force.getForceGroup());
}
int startPairIndex = cl.getContextIndex()*force.getNumScreenedPairs()/numContexts;
int endPairIndex = (cl.getContextIndex()+1)*force.getNumScreenedPairs()/numContexts;
int numPairs = endPairIndex-startPairIndex;
if (numPairs > 0) {
// Create the screened interaction between dipole pairs.
vector<vector<int> > atoms(numPairs, vector<int>(4));
pairParams = OpenCLArray::create<mm_float2>(cl, numPairs, "drudePairParams");
vector<mm_float2> paramVector(numPairs);
for (int i = 0; i < numPairs; i++) {
int drude1, drude2;
double thole;
force.getScreenedPairParameters(startPairIndex+i, drude1, drude2, thole);
int p2, p3, p4;
double charge1, charge2, polarizability1, polarizability2, aniso12, aniso34;
force.getParticleParameters(drude1, atoms[i][0], atoms[i][1], p2, p3, p4, charge1, polarizability1, aniso12, aniso34);
force.getParticleParameters(drude2, atoms[i][2], atoms[i][3], p2, p3, p4, charge2, polarizability2, aniso12, aniso34);
double screeningScale = thole/pow(polarizability1*polarizability2, 1.0/6.0);
double energyScale = ONE_4PI_EPS0*charge1*charge2;
paramVector[i] = mm_float2((float) screeningScale, (float) energyScale);
}
pairParams->upload(paramVector);
map<string, string> replacements;
replacements["PARAMS"] = cl.getBondedUtilities().addArgument(pairParams->getDeviceBuffer(), "float2");
cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLDrudeKernelSources::drudePairForce, replacements), force.getForceGroup());
}
cl.addForce(new OpenCLDrudeForceInfo(force));
}
double OpenCLCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return 0.0;
}
void OpenCLCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) {
}
OpenCLIntegrateDrudeLangevinStepKernel::~OpenCLIntegrateDrudeLangevinStepKernel() {
if (normalParticles != NULL)
delete normalParticles;
if (pairParticles != NULL)
delete pairParticles;
}
void OpenCLIntegrateDrudeLangevinStepKernel::initialize(const System& system, const DrudeLangevinIntegrator& integrator, const DrudeForce& force) {
cl.getPlatformData().initializeContexts(system);
cl.getIntegrationUtilities().initRandomNumberGenerator((unsigned int) integrator.getRandomNumberSeed());
// Identify particle pairs and ordinary particles.
set<int> particles;
vector<int> normalParticleVec;
vector<mm_int2> pairParticleVec;
for (int i = 0; i < system.getNumParticles(); i++)
particles.insert(i);
for (int i = 0; i < force.getNumParticles(); i++) {
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
force.getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
particles.erase(p);
particles.erase(p1);
pairParticleVec.push_back(mm_int2(p, p1));
}
normalParticleVec.insert(normalParticleVec.begin(), particles.begin(), particles.end());
normalParticles = OpenCLArray::create<int>(cl, max((int) normalParticleVec.size(), 1), "drudeNormalParticles");
pairParticles = OpenCLArray::create<cl_int2>(cl, max((int) pairParticleVec.size(), 1), "drudePairParticles");
if (normalParticleVec.size() > 0)
normalParticles->upload(normalParticleVec);
if (pairParticleVec.size() > 0)
pairParticles->upload(pairParticleVec);
// Create kernels.
map<string, string> defines;
defines["NUM_ATOMS"] = cl.intToString(cl.getNumAtoms());
defines["PADDED_NUM_ATOMS"] = cl.intToString(cl.getPaddedNumAtoms());
defines["NUM_NORMAL_PARTICLES"] = cl.intToString(normalParticleVec.size());
defines["NUM_PAIRS"] = cl.intToString(pairParticleVec.size());
map<string, string> replacements;
cl::Program program = cl.createProgram(OpenCLDrudeKernelSources::drudeLangevin, defines, "");
kernel1 = cl::Kernel(program, "integrateDrudeLangevinPart1");
kernel2 = cl::Kernel(program, "integrateDrudeLangevinPart2");
prevStepSize = -1.0;
}
void OpenCLIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
int numAtoms = cl.getNumAtoms();
if (!hasInitializedKernels) {
hasInitializedKernels = true;
kernel1.setArg<cl::Buffer>(0, cl.getVelm().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(1, cl.getForce().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(2, integration.getPosDelta().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(3, normalParticles->getDeviceBuffer());
kernel1.setArg<cl::Buffer>(4, pairParticles->getDeviceBuffer());
kernel1.setArg<cl::Buffer>(5, integration.getStepSize().getDeviceBuffer());
kernel1.setArg<cl::Buffer>(12, integration.getRandom().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(0, cl.getPosq().getDeviceBuffer());
if (cl.getUseMixedPrecision())
kernel2.setArg<cl::Buffer>(1, cl.getPosqCorrection().getDeviceBuffer());
else
kernel2.setArg<void*>(1, NULL);
kernel2.setArg<cl::Buffer>(2, integration.getPosDelta().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(3, cl.getVelm().getDeviceBuffer());
kernel2.setArg<cl::Buffer>(4, integration.getStepSize().getDeviceBuffer());
}
// Compute integrator coefficients.
double stepSize = integrator.getStepSize();
double vscale = exp(-stepSize*integrator.getFriction());
double fscale = (1-vscale)/integrator.getFriction();
double noisescale = sqrt(2*BOLTZ*integrator.getTemperature()*integrator.getFriction())*sqrt(0.5*(1-vscale*vscale)/integrator.getFriction());
double vscaleDrude = exp(-stepSize*integrator.getDrudeFriction());
double fscaleDrude = (1-vscaleDrude)/integrator.getDrudeFriction();
double noisescaleDrude = sqrt(2*BOLTZ*integrator.getDrudeTemperature()*integrator.getDrudeFriction())*sqrt(0.5*(1-vscaleDrude*vscaleDrude)/integrator.getDrudeFriction());
if (stepSize != prevStepSize) {
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
mm_double2 ss = mm_double2(0, stepSize);
integration.getStepSize().upload(&ss);
}
else {
mm_float2 ss = mm_float2(0, (float) stepSize);
integration.getStepSize().upload(&ss);
}
prevStepSize = stepSize;
}
if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {
kernel1.setArg<cl_double>(6, vscale);
kernel1.setArg<cl_double>(7, fscale);
kernel1.setArg<cl_double>(8, noisescale);
kernel1.setArg<cl_double>(9, vscaleDrude);
kernel1.setArg<cl_double>(10, fscaleDrude);
kernel1.setArg<cl_double>(11, noisescaleDrude);
}
else {
kernel1.setArg<cl_float>(6, (cl_float) vscale);
kernel1.setArg<cl_float>(7, (cl_float) fscale);
kernel1.setArg<cl_float>(8, (cl_float) noisescale);
kernel1.setArg<cl_float>(9, (cl_float) vscaleDrude);
kernel1.setArg<cl_float>(10, (cl_float) fscaleDrude);
kernel1.setArg<cl_float>(11, (cl_float) noisescaleDrude);
}
// Call the first integration kernel.
kernel1.setArg<cl_uint>(13, integration.prepareRandomNumbers(normalParticles->getSize()+2*pairParticles->getSize()));
cl.executeKernel(kernel1, numAtoms);
// Apply constraints.
integration.applyConstraints(integrator.getConstraintTolerance());
// Call the second integration kernel.
cl.executeKernel(kernel2, numAtoms);
integration.computeVirtualSites();
// Update the time and step count.
cl.setTime(cl.getTime()+stepSize);
cl.setStepCount(cl.getStepCount()+1);
}
double OpenCLIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
return cl.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
}
plugins/drude/platforms/opencl/src/OpenCLDrudeKernels.h 0 → 100644
View file @ f8c98f32
#ifndef OPENCL_DRUDE_KERNELS_H_
#define OPENCL_DRUDE_KERNELS_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/DrudeKernels.h"
#include "OpenCLContext.h"
#include "OpenCLArray.h"
namespace OpenMM {
/**
* This kernel is invoked by DrudeForce to calculate the forces acting on the system and the energy of the system.
*/
class OpenCLCalcDrudeForceKernel : public CalcDrudeForceKernel {
public:
OpenCLCalcDrudeForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) :
CalcDrudeForceKernel(name, platform), cl(cl), particleParams(NULL), pairParams(NULL) {
}
~OpenCLCalcDrudeForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the DrudeForce this kernel will be used for
*/
void initialize(const System& system, const DrudeForce& 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);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the DrudeForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const DrudeForce& force);
private:
OpenCLContext& cl;
OpenCLArray* particleParams;
OpenCLArray* pairParams;
};
/**
* This kernel is invoked by DrudeLangevinIntegrator to take one time step
*/
class OpenCLIntegrateDrudeLangevinStepKernel : public IntegrateDrudeLangevinStepKernel {
public:
OpenCLIntegrateDrudeLangevinStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) :
IntegrateDrudeLangevinStepKernel(name, platform), cl(cl), hasInitializedKernels(false), normalParticles(NULL), pairParticles(NULL) {
}
~OpenCLIntegrateDrudeLangevinStepKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the DrudeLangevinIntegrator this kernel will be used for
* @param force the DrudeForce to get particle parameters from
*/
void initialize(const System& system, const DrudeLangevinIntegrator& integrator, const DrudeForce& force);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the DrudeLangevinIntegrator this kernel is being used for
*/
void execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the DrudeLangevinIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator);
private:
OpenCLContext& cl;
bool hasInitializedKernels;
double prevStepSize;
OpenCLArray* normalParticles;
OpenCLArray* pairParticles;
cl::Kernel kernel1, kernel2;
};
} // namespace OpenMM
#endif /*OPENCL_DRUDE_KERNELS_H_*/
plugins/drude/platforms/opencl/src/kernels/drudeLangevin.cl 0 → 100644
View file @ f8c98f32
/**
* Perform the first step of Langevin integration.
*/
__kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __global const real4* restrict force, __global mixed4* restrict posDelta,
__global const int* restrict normalParticles, __global const int2* restrict pairParticles, __global const mixed2* restrict dt, mixed vscale, mixed fscale,
mixed noisescale, mixed vscaleDrude, mixed fscaleDrude, mixed noisescaleDrude, __global const float4* restrict random, unsigned int randomIndex) {
mixed stepSize = dt[0].y;
// Update normal particles.
for (int i = get_global_id(0); i < NUM_NORMAL_PARTICLES; i += get_global_size(0)) {
int index = normalParticles[i];
mixed4 velocity = velm[index];
if (velocity.w != 0) {
mixed sqrtInvMass = SQRT(velocity.w);
float4 rand = random[randomIndex+index];
real4 f = force[index];
velocity.xyz = vscale*velocity.xyz + fscale*velocity.w*f.xyz + noisescale*sqrtInvMass*rand.xyz;
velm[index] = velocity;
posDelta[index] = (mixed4) (stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0);
}
}
// Update Drude particle pairs.
randomIndex += NUM_NORMAL_PARTICLES;
for (int i = get_global_id(0); i < NUM_PAIRS; i += get_global_size(0)) {
int2 particles = pairParticles[i];
mixed4 velocity1 = velm[particles.x];
mixed4 velocity2 = velm[particles.y];
mixed mass1 = RECIP(velocity1.w);
mixed mass2 = RECIP(velocity2.w);
mixed invTotalMass = RECIP(mass1+mass2);
mixed invReducedMass = (mass1+mass2)*velocity1.w*velocity2.w;
mixed mass1fract = invTotalMass*mass1;
mixed mass2fract = invTotalMass*mass2;
mixed sqrtInvTotalMass = SQRT(invTotalMass);
mixed sqrtInvReducedMass = SQRT(invReducedMass);
mixed4 cmVel = velocity1*mass1fract+velocity2*mass2fract;
mixed4 relVel = velocity2-velocity1;
mixed4 force1 = force[particles.x];
mixed4 force2 = force[particles.y];
mixed4 cmForce = force1+force2;
mixed4 relForce = force2*mass1fract - force1*mass2fract;
float4 rand1 = random[randomIndex+2*i];
float4 rand2 = random[randomIndex+2*i+1];
cmVel.xyz = vscale*cmVel.xyz + fscale*invTotalMass*cmForce.xyz + noisescale*sqrtInvTotalMass*rand1.xyz;
relVel.xyz = vscaleDrude*relVel.xyz + fscaleDrude*invReducedMass*relForce.xyz + noisescaleDrude*sqrtInvReducedMass*rand2.xyz;
velocity1.xyz = cmVel.xyz-relVel.xyz*mass2fract;
velocity2.xyz = cmVel.xyz+relVel.xyz*mass1fract;
velm[particles.x] = velocity1;
velm[particles.y] = velocity2;
posDelta[particles.x] = (mixed4) (stepSize*velocity1.x, stepSize*velocity1.y, stepSize*velocity1.z, 0);
posDelta[particles.y] = (mixed4) (stepSize*velocity2.x, stepSize*velocity2.y, stepSize*velocity2.z, 0);
}
}
/**
* Perform the second step of Langevin integration.
*/
__kernel void integrateDrudeLangevinPart2(__global real4* restrict posq, __global real4* restrict posqCorrection, __global const mixed4* restrict posDelta, __global mixed4* restrict velm, __global const mixed2* restrict dt) {
#ifdef SUPPORTS_DOUBLE_PRECISION
double invStepSize = 1.0/dt[0].y;
#else
float invStepSize = 1.0f/dt[0].y;
#endif
int index = get_global_id(0);
while (index < NUM_ATOMS) {
mixed4 vel = velm[index];
if (vel.w != 0.0) {
#ifdef USE_MIXED_PRECISION
real4 pos1 = posq[index];
real4 pos2 = posqCorrection[index];
mixed4 pos = (mixed4) (pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
#else
real4 pos = posq[index];
#endif
mixed4 delta = posDelta[index];
pos.xyz += delta.xyz;
#ifdef SUPPORTS_DOUBLE_PRECISION
vel.xyz = convert_mixed4(invStepSize*convert_double4(delta)).xyz;
#else
vel.xyz = invStepSize*delta.xyz;
#endif
#ifdef USE_MIXED_PRECISION
posq[index] = convert_real4(pos);
posqCorrection[index] = (real4) (pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
#else
posq[index] = pos;
#endif
velm[index] = vel;
}
index += get_global_size(0);
}
}
plugins/drude/platforms/opencl/src/kernels/drudePairForce.cl 0 → 100644
View file @ f8c98f32
float2 drudeParams = PARAMS[index];
real4 force1 = 0;
real4 force2 = 0;
real4 force3 = 0;
real4 force4 = 0;
// First pair.
real4 delta = (real4) (pos1.xyz-pos3.xyz, 0);
real rInv = RSQRT(dot(delta, delta));
real r = RECIP(rInv);
real u = drudeParams.x*r;
real screening = 1-(1+0.5f*u)*EXP(-u);
real pairEnergy = drudeParams.y*screening*rInv;
energy += pairEnergy;
real4 f = delta*(pairEnergy*rInv*rInv);
force1 += f;
force3 -= f;
// Second pair.
delta = (real4) (pos1.xyz-pos4.xyz, 0);
rInv = RSQRT(dot(delta, delta));
r = RECIP(rInv);
u = drudeParams.x*r;
screening = 1-(1+0.5f*u)*EXP(-u);
pairEnergy = -drudeParams.y*screening*rInv;
energy += pairEnergy;
f = delta*(pairEnergy*rInv*rInv);
force1 += f;
force4 -= f;
// Third pair.
delta = (real4) (pos2.xyz-pos3.xyz, 0);
rInv = RSQRT(dot(delta, delta));
r = RECIP(rInv);
u = drudeParams.x*r;
screening = 1-(1+0.5f*u)*EXP(-u);
pairEnergy = -drudeParams.y*screening*rInv;
energy += pairEnergy;
f = delta*(pairEnergy*rInv*rInv);
force2 += f;
force3 -= f;
// Fourth pair.
delta = (real4) (pos2.xyz-pos4.xyz, 0);
rInv = RSQRT(dot(delta, delta));
r = RECIP(rInv);
u = drudeParams.x*r;
screening = 1-(1+0.5f*u)*EXP(-u);
pairEnergy = drudeParams.y*screening*rInv;
energy += pairEnergy;
f = delta*(pairEnergy*rInv*rInv);
force2 += f;
force4 -= f;
plugins/drude/platforms/opencl/src/kernels/drudeParticleForce.cl 0 → 100644
View file @ f8c98f32
real4 delta = (real4) (pos1.xyz-pos2.xyz, 0);
real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
float4 drudeParams = PARAMS[index];
float k1 = drudeParams.x;
float k2 = drudeParams.y;
float k3 = drudeParams.z;
// Compute the isotropic force.
energy += 0.5f*k3*r2;
real4 force1 = -delta*k3;
real4 force2 = delta*k3;
real4 force3 = 0;
real4 force4 = 0;
real4 force5 = 0;
// Compute the first anisotropic force.
if (k1 != 0) {
real4 dir = (real4) (pos2.xyz-pos3.xyz, 0);
real invDist = RSQRT(dot(dir, dir));
dir *= invDist;
real rprime = dot(dir, delta);
energy += 0.5f*k1*rprime*rprime;
real4 f1 = dir*(k1*rprime);
real4 f2 = (delta-dir*rprime)*(k1*rprime*invDist);
force1 -= f1;
force2 += f1-f2;
force3 += f2;
}
// Compute the second anisotropic force.
if (k2 != 0) {
real4 dir = (real4) (pos3.xyz-pos4.xyz, 0);
real invDist = RSQRT(dot(dir, dir));
dir *= invDist;
real rprime = dot(dir, delta);
energy += 0.5f*k2*rprime*rprime;
real4 f1 = dir*(k2*rprime);
real4 f2 = (delta-dir*rprime)*(k2*rprime*invDist);
force1 -= f1;
force2 += f1;
force4 -= f2;
force5 += f2;
}
plugins/drude/platforms/opencl/tests/CMakeLists.txt 0 → 100644
View file @ f8c98f32
#
# Testing
#
ENABLE_TESTING()
INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR})
# 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_DRUDE_TARGET} ${SHARED_TARGET})
ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
IF (OPENMM_BUILD_OPENCL_DOUBLE_PRECISION_TESTS)
ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed)
ADD_TEST(${TEST_ROOT}Double ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} double)
ENDIF(OPENMM_BUILD_OPENCL_DOUBLE_PRECISION_TESTS)
ENDFOREACH(TEST_PROG ${TEST_PROGS})
plugins/drude/platforms/opencl/tests/TestOpenCLDrudeForce.cpp 0 → 100644
View file @ f8c98f32
/* -------------------------------------------------------------------------- *
* 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) 2013 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 OpenCL implementation of DrudeForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/DrudeForce.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
extern "C" void registerDrudeOpenCLKernelFactories();
void validateForce(System& system, vector<Vec3>& positions, double expectedEnergy) {
// Given a System containing a Drude force, check that its energy has the expected value.
VerletIntegrator integ(1.0);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Energy | State::Forces);
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
// Try moving each particle along each axis, and see if the energy changes by the correct amount.
double offset = 1e-3;
for (int i = 0; i < system.getNumParticles(); i++)
for (int j = 0; j < 3; j++) {
vector<Vec3> offsetPos = positions;
offsetPos[i][j] = positions[i][j]-offset;
context.setPositions(offsetPos);
double e1 = context.getState(State::Energy | State::Forces).getPotentialEnergy();
offsetPos[i][j] = positions[i][j]+offset;
context.setPositions(offsetPos);
double e2 = context.getState(State::Energy | State::Forces).getPotentialEnergy();
ASSERT_EQUAL_TOL(state.getForces()[i][j], (e1-e2)/(2*offset), 1e-3);
}
}
void testSingleParticle() {
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
system.addForce(drude);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(2, 0, 0);
validateForce(system, positions, 0.5*k*3*3);
}
void testAnisotropicParticle() {
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
const double a1 = 0.8;
const double a2 = 1.1;
const double k1 = k/a1;
const double k2 = k/a2;
const double k3 = k/(3-a1-a2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, 2, 3, 4, charge, alpha, a1, a2);
system.addForce(drude);
vector<Vec3> positions(5);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0.1, -0.5, 0.8);
positions[2] = Vec3(0, 2, 0);
positions[3] = Vec3(1, 2, 0);
positions[4] = Vec3(1, 2, 3);
validateForce(system, positions, 0.5*k1*0.5*0.5 + 0.5*k2*0.8*0.8 + 0.5*k3*0.1*0.1);
}
double computeScreening(double r, double thole, double alpha1, double alpha2) {
double u = r*thole/pow(alpha1*alpha2, 1.0/6.0);
return 1.0-(1.0+u/2)*exp(-u);
}
void testThole() {
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
const double thole = 2.5;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
drude->addParticle(3, 2, -1, -1, -1, charge, alpha, 1, 1);
drude->addScreenedPair(0, 1, thole);
system.addForce(drude);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, -0.5, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 0, 0.3);
double energySpring1 = 0.5*k*0.5*0.5;
double energySpring2 = 0.5*k*0.3*0.3;
double energyDipole = 0.0;
double q[] = {-charge, charge, -charge, charge};
for (int i = 0; i < 2; i++)
for (int j = 2; j < 4; j++) {
Vec3 delta = positions[i]-positions[j];
double r = sqrt(delta.dot(delta));
energyDipole += ONE_4PI_EPS0*q[i]*q[j]*computeScreening(r, thole, alpha, alpha)/r;
}
validateForce(system, positions, energySpring1+energySpring2+energyDipole);
}
int main(int argc, char* argv[]) {
try {
registerDrudeOpenCLKernelFactories();
if (argc > 1)
Platform::getPlatformByName("OpenCL").setPropertyDefaultValue("OpenCLPrecision", std::string(argv[1]));
testSingleParticle();
testAnisotropicParticle();
testThole();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
}
plugins/drude/platforms/opencl/tests/TestOpenCLDrudeLangevinIntegrator.cpp 0 → 100644
View file @ f8c98f32
/* -------------------------------------------------------------------------- *
* 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) 2013 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 OpenCL implementation of DrudeLangevinIntegrator.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/VirtualSite.h"
#include "openmm/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
extern "C" OPENMM_EXPORT void registerDrudeOpenCLKernelFactories();
void testSinglePair() {
const double temperature = 300.0;
const double temperatureDrude = 10.0;
const double k = 1.5;
const double charge = 0.1;
const double alpha = charge*charge/k;
const double mass1 = 1.0;
const double mass2 = 0.1;
const double totalMass = mass1+mass2;
const double reducedMass = (mass1*mass2)/(mass1+mass2);
System system;
system.addParticle(mass1);
system.addParticle(mass2);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
system.addForce(drude);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 0, 0);
DrudeLangevinIntegrator integ(temperature, 20.0, temperatureDrude, 20.0, 0.003);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
context.setPositions(positions);
// Equilibrate.
integ.step(1000);
// Compute the internal and center of mass temperatures.
double keCM = 0, keInternal = 0;
int numSteps = 10000;
for (int i = 0; i < numSteps; i++) {
integ.step(10);
State state = context.getState(State::Velocities);
const vector<Vec3>& vel = state.getVelocities();
Vec3 velCM = vel[0]*(mass1/totalMass) + vel[1]*(mass2/totalMass);
keCM += 0.5*totalMass*velCM.dot(velCM);
Vec3 velInternal = vel[0]-vel[1];
keInternal += 0.5*reducedMass*velInternal.dot(velInternal);
}
ASSERT_USUALLY_EQUAL_TOL(3*0.5*BOLTZ*temperature, keCM/numSteps, 0.1);
ASSERT_USUALLY_EQUAL_TOL(3*0.5*BOLTZ*temperatureDrude, keInternal/numSteps, 0.01);
}
void testWater() {
// Create a box of SWM4-NDP water molecules. This involves constraints, virtual sites,
// and Drude particles.
const int gridSize = 4;
const int numMolecules = gridSize*gridSize*gridSize;
const double spacing = 0.6;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
const double temperatureDrude = 10.0;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
DrudeForce* drude = new DrudeForce();
system.addForce(nonbonded);
system.addForce(drude);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(1.0);
for (int i = 0; i < numMolecules; i++) {
int startIndex = system.getNumParticles();
system.addParticle(15.6); // O
system.addParticle(0.4); // D
system.addParticle(1.0); // H1
system.addParticle(1.0); // H2
system.addParticle(0.0); // M
nonbonded->addParticle(1.71636, 0.318395, 0.21094*4.184);
nonbonded->addParticle(-1.71636, 1, 0);
nonbonded->addParticle(0.55733, 1, 0);
nonbonded->addParticle(0.55733, 1, 0);
nonbonded->addParticle(-1.11466, 1, 0);
for (int j = 0; j < 5; j++)
for (int k = 0; k < j; k++)
nonbonded->addException(startIndex+j, startIndex+k, 0, 1, 0);
system.addConstraint(startIndex, startIndex+2, 0.09572);
system.addConstraint(startIndex, startIndex+3, 0.09572);
system.addConstraint(startIndex+2, startIndex+3, 0.15139);
system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
}
vector<Vec3> positions;
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos(i*spacing, j*spacing, k*spacing);
positions.push_back(pos);
positions.push_back(pos);
positions.push_back(pos+Vec3(0.09572, 0, 0));
positions.push_back(pos+Vec3(-0.023999, 0.092663, 0));
positions.push_back(pos);
}
// Simulate it and check the temperature.
DrudeLangevinIntegrator integ(temperature, 50.0, temperatureDrude, 50.0, 0.0005);
Platform& platform = Platform::getPlatformByName("OpenCL");
Context context(system, integ, platform);
context.setPositions(positions);
context.applyConstraints(1e-5);
// Equilibrate.
integ.step(1000);
// Compute the internal and center of mass temperatures.
double ke = 0;
int numSteps = 10000;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
ke += context.getState(State::Energy).getKineticEnergy();
}
ke /= numSteps;
int numStandardDof = 3*3*numMolecules-system.getNumConstraints();
int numDrudeDof = 3*numMolecules;
int numDof = numStandardDof+numDrudeDof;
double expectedTemp = (numStandardDof*temperature+numDrudeDof*temperatureDrude)/numDof;
ASSERT_USUALLY_EQUAL_TOL(expectedTemp, ke/(0.5*numDof*BOLTZ), 0.02);
}
int main(int argc, char* argv[]) {
try {
registerDrudeOpenCLKernelFactories();
if (argc > 1)
Platform::getPlatformByName("OpenCL").setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testSinglePair();
testWater();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment