Skip to content

GitLab

Commit 31d761d5 authored by Peter Eastman's avatar Peter Eastman
Browse files

Began implementing Drude plugin

parent 6fabca14
Hide whitespace changes
Inline Side-by-side
Showing with 2090 additions and 0 deletions
plugins/drude/CMakeLists.txt 0 → 100644
View file @ 31d761d5
#---------------------------------------------------
# OpenMM Drude Plugin
#
# Creates OpenMM Drude plugin library, base name=OpenMMDrude.
# Default libraries are shared & optimized. Variants
# are created for static (_static) and debug (_d).
#
# Windows:
# OpenMMDrude[_d].dll
# OpenMMDrude[_d].lib
# OpenMMDrude_static[_d].lib
# Unix:
# libOpenMMDrude[_d].so
# libOpenMMDrude_static[_d].a
#----------------------------------------------------
#INCLUDE(Dart)
# ----------------------------------------------------------------------------
SET(CREATE_SERIALIZABLE_OPENMM_DRUDE FALSE )
# ----------------------------------------------------------------------------
# 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_DRUDE_PLUGIN_SOURCE_SUBDIRS . openmmapi platforms/reference)
SET(OPENMM_DRUDE_LIBRARY_NAME OpenMMDrude)
SET(SHARED_DRUDE_TARGET ${OPENMM_DRUDE_LIBRARY_NAME})
SET(STATIC_DRUDE_TARGET ${OPENMM_DRUDE_LIBRARY_NAME}_static)
IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
SET(SHARED_DRUDE_SERIALIZABLE_TARGET ${OPENMM_DRUDE_LIBRARY_NAME}_serializable)
ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
# But on Unix or Cygwin we have to add the suffix manually
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_DRUDE_TARGET ${SHARED_DRUDE_TARGET}_d)
SET(STATIC_DRUDE_TARGET ${STATIC_DRUDE_TARGET}_d)
IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
SET(SHARED_DRUDE_SERIALIZABLE_TARGET ${SHARED_DRUDE_SERIALIZABLE_TARGET}_d)
ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
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_DRUDE_INCLUDE_DIRS) # start empty
FOREACH(subdir ${OPENMM_DRUDE_PLUGIN_SOURCE_SUBDIRS})
# append
SET(API_DRUDE_INCLUDE_DIRS ${API_DRUDE_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_DRUDE_REL_INCLUDE_FILES) # start these out empty
SET(API_DRUDE_ABS_INCLUDE_FILES)
FOREACH(dir ${API_DRUDE_INCLUDE_DIRS})
FILE(GLOB fullpaths ${dir}/*.h) # returns full pathnames
SET(API_DRUDE_ABS_INCLUDE_FILES ${API_DRUDE_ABS_INCLUDE_FILES} ${fullpaths})
FOREACH(pathname ${fullpaths})
GET_FILENAME_COMPONENT(filename ${pathname} NAME)
SET(API_DRUDE_REL_INCLUDE_FILES ${API_DRUDE_REL_INCLUDE_FILES} ${dir}/${filename})
ENDFOREACH(pathname)
ENDFOREACH(dir)
# collect up source files
SET(SOURCE_DRUDE_FILES) # empty
SET(SOURCE_DRUDE_INCLUDE_FILES)
FOREACH(subdir ${OPENMM_DRUDE_PLUGIN_SOURCE_SUBDIRS})
FILE(GLOB src_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.cpp ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*/*.cpp)
FILE(GLOB incl_files ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/src/*.h)
SET(SOURCE_DRUDE_FILES ${SOURCE_DRUDE_FILES} ${src_files}) #append
SET(SOURCE_DRUDE_INCLUDE_FILES ${SOURCE_DRUDE_INCLUDE_FILES} ${incl_files})
## Make sure we find these locally before looking in OpenMM/include if
## OpenMM was previously installed there.
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/${subdir}/include)
ENDFOREACH(subdir)
INCLUDE_DIRECTORIES(BEFORE ${OPENMM_DIR}/platforms/reference/src)
INCLUDE_DIRECTORIES(BEFORE ${OPENMM_DIR}/platforms/reference/src/SimTKReference)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/src/SimTKReference)
# ----------------------------------------------------------------------------
# If API_DRUDE wrappers are being generated, and add them to the build.
#IF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
# ADD_SUBDIRECTORY(wrappers)
# SET(SOURCE_DRUDE_FILES ${SOURCE_DRUDE_FILES} wrappers/DrudeOpenMMCWrapper.cpp wrappers/DrudeOpenMMFortranWrapper.cpp)
# SET_SOURCE_FILES_PROPERTIES(wrappers/DrudeOpenMMCWrapper.cpp wrappers/DrudeOpenMMFortranWrapper.cpp PROPERTIES GENERATED TRUE)
#ENDIF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
ADD_LIBRARY(${SHARED_DRUDE_TARGET} SHARED ${SOURCE_DRUDE_FILES} ${SOURCE_DRUDE_INCLUDE_FILES} ${API_DRUDE_ABS_INCLUDE_FILES})
SET_TARGET_PROPERTIES(${SHARED_DRUDE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_DRUDE_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY")
#IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
# ADD_LIBRARY(${SHARED_DRUDE_SERIALIZABLE_TARGET} SHARED ${SOURCE_DRUDE_FILES} ${SOURCE_DRUDE_INCLUDE_FILES} ${API_DRUDE_ABS_INCLUDE_FILES})
# SET_TARGET_PROPERTIES(${SHARED_DRUDE_SERIALIZABLE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_DRUDE_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_SERIALIZE")
# INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/../../serialization/include)
#ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
IF(OPENMM_BUILD_STATIC_LIB)
ADD_LIBRARY(${STATIC_DRUDE_TARGET} STATIC ${SOURCE_DRUDE_FILES} ${SOURCE_DRUDE_INCLUDE_FILES} ${API_DRUDE_ABS_INCLUDE_FILES})
SET_TARGET_PROPERTIES(${STATIC_DRUDE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_USE_STATIC_LIBRARIES -DLEPTON_BUILDING_STATIC_LIBRARY")
ENDIF(OPENMM_BUILD_STATIC_LIB)
#IF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
# ADD_DEPENDENCIES(${SHARED_DRUDE_TARGET} DrudeApiWrappers)
# IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
# ADD_DEPENDENCIES(${SHARED_DRUDE_SERIALIZABLE_TARGET} DrudeApiWrappers)
# ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
# IF(OPENMM_BUILD_STATIC_LIB)
# ADD_DEPENDENCIES(${STATIC_DRUDE_TARGET} DrudeApiWrappers)
# ENDIF(OPENMM_BUILD_STATIC_LIB)
#ENDIF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
# ----------------------------------------------------------------------------
# On Linux need to link to libdl
FIND_LIBRARY(DL_LIBRARY dl)
IF(DL_LIBRARY)
TARGET_LINK_LIBRARIES(${SHARED_DRUDE_TARGET} ${DL_LIBRARY})
IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
TARGET_LINK_LIBRARIES(${SHARED_DRUDE_SERIALIZABLE_TARGET} ${DL_LIBRARY})
ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
IF(OPENMM_BUILD_STATIC_LIB)
TARGET_LINK_LIBRARIES(${STATIC_DRUDE_TARGET} ${DL_LIBRARY})
ENDIF(OPENMM_BUILD_STATIC_LIB)
ENDIF(DL_LIBRARY)
SET( OpenMMLib OpenMM )
# But on Unix or Cygwin we have to add the suffix manually
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(OpenMMLib ${OpenMMLib}_d)
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
TARGET_LINK_LIBRARIES( ${SHARED_DRUDE_TARGET} ${SHARED_TARGET} )
IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
TARGET_LINK_LIBRARIES( ${SHARED_DRUDE_SERIALIZABLE_TARGET} ${SHARED_TARGET} )
ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
IF(OPENMM_BUILD_STATIC_LIB)
TARGET_LINK_LIBRARIES( ${STATIC_DRUDE_TARGET} ${STATIC_TARGET} )
ENDIF(OPENMM_BUILD_STATIC_LIB)
ADD_SUBDIRECTORY(platforms/reference/tests)
# Which hardware platforms to build
#SET(OPENMM_BUILD_DRUDE_PATH)
#SET(OPENMM_BUILD_DRUDE_CUDA_PATH)
#IF(OPENMM_BUILD_DRUDE_CUDA_LIB)
# ADD_SUBDIRECTORY(platforms/cuda)
# SET(OPENMM_BUILD_DRUDE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/platforms/cuda)
# SET(OPENMM_BUILD_DRUDE_CUDA_PATH ${CMAKE_CURRENT_SOURCE_DIR}/platforms/cuda)
# SET(OPENMM_DRUDE_CUDA_SOURCE_SUBDIRS . openmmapi olla platforms/cuda)
#ENDIF(OPENMM_BUILD_DRUDE_CUDA_LIB)
IF(OPENCL_FOUND)
SET(OPENMM_BUILD_DRUDE_OPENCL_LIB ON CACHE BOOL "Build Drude implementation for OpenCL")
ELSE(OPENCL_FOUND)
SET(OPENMM_BUILD_DRUDE_OPENCL_LIB OFF CACHE BOOL "Build Drude implementation for OpenCL")
ENDIF(OPENCL_FOUND)
IF(OPENMM_BUILD_DRUDE_OPENCL_LIB)
ADD_SUBDIRECTORY(platforms/opencl)
ENDIF(OPENMM_BUILD_DRUDE_OPENCL_LIB)
IF(CUDA_FOUND)
SET(OPENMM_BUILD_DRUDE_CUDA_LIB ON CACHE BOOL "Build Drude implementation for CUDA")
ELSE(CUDA_FOUND)
SET(OPENMM_BUILD_DRUDE_CUDA_LIB OFF CACHE BOOL "Build Drude implementation for CUDA")
ENDIF(CUDA_FOUND)
IF(OPENMM_BUILD_DRUDE_CUDA_LIB)
ADD_SUBDIRECTORY(platforms/cuda)
ENDIF(OPENMM_BUILD_DRUDE_CUDA_LIB)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_DRUDE_TARGET})
IF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_DRUDE_SERIALIZABLE_TARGET})
ENDIF( CREATE_SERIALIZABLE_OPENMM_DRUDE )
IF(OPENMM_BUILD_STATIC_LIB)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_DRUDE_TARGET})
ENDIF(OPENMM_BUILD_STATIC_LIB)
FILE(GLOB CORE_HEADERS include/*.h */include/*.h)
FILE(GLOB TOP_HEADERS include/openmm/*.h */include/openmm/*.h)
FILE(GLOB INTERNAL_HEADERS include/openmm/internal/*.h */include/openmm/internal/*.h )
INSTALL_FILES(/include FILES ${CORE_HEADERS})
INSTALL_FILES(/include/openmm FILES ${TOP_HEADERS})
INSTALL_FILES(/include/openmm/internal FILES ${INTERNAL_HEADERS})
#
# Testing
#
ENABLE_TESTING()
IF (EXECUTABLE_OUTPUT_PATH)
SET (TEST_PATH ${EXECUTABLE_OUTPUT_PATH})
ELSE (EXECUTABLE_OUTPUT_PATH)
SET (TEST_PATH .)
ENDIF (EXECUTABLE_OUTPUT_PATH)
#IF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
# ADD_SUBDIRECTORY(serialization)
#ENDIF (OPENMM_BUILD_SERIALIZATION_SUPPORT)
#INCLUDE(ApiDoxygen.cmake)
#ADD_SUBDIRECTORY(tests)
plugins/drude/openmmapi/include/openmm/DrudeForce.h 0 → 100644
View file @ 31d761d5
#ifndef OPENMM_DRUDEFORCE_H_
#define OPENMM_DRUDEFORCE_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/Context.h"
#include "openmm/Force.h"
#include <vector>
#include "internal/windowsExportDrude.h"
namespace OpenMM {
/**
* This class implements forces that are specific to Drude oscillators. There are two distinct forces
* it applies: an anisotropic harmonic force connecting each Drude particle to its parent particle; and
* a screened Coulomb interaction between specific pairs of dipoles. The latter is typically used between
* closely bonded particles whose Coulomb interaction would otherwise be fully excluded.
*
* To use this class, create a DrudeForce object, then call addParticle() once for each Drude particle in the
* System to define its parameters. After a particle has been added, you can modify its force field parameters
* by calling setParticleParameters(). This will have no effect on Contexts that already exist unless you
* call updateParametersInContext(). Likewise, call addScreenedPair() for each pair of dipoles (each dipole
* consisting of a Drude particle and its parent) that should be computed.
*/
class OPENMM_EXPORT_DRUDE DrudeForce : public Force {
public:
/**
* Create a DrudeForce.
*/
DrudeForce();
/**
* Get the number of particles for which force field parameters have been defined.
*/
int getNumParticles() const {
return particles.size();
}
/**
* Get the number of special interactions that should be calculated differently from other interactions.
*/
int getNumScreenedPairs() const {
return screenedPairs.size();
}
/**
* Add a Drude particle to which forces should be applied.
*
* @param particle the index within the System of the Drude particle
* @param particle1 the index within the System of the particle to which the Drude particle is attached
* @param particle2 the index within the System of the second particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso12 will be ignored.
* @param particle3 the index within the System of the third particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso34 will be ignored.
* @param particle4 the index within the System of the fourth particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso34 will be ignored.
* @param charge The charge on the Drude particle
* @param polarizability The isotropic polarizability
* @param aniso12 The scale factor for the polarizability along the direction defined by particle1 and particle2
* @param aniso34 The scale factor for the polarizability along the direction defined by particle3 and particle4
* @return the index of the particle that was added
*/
int addParticle(int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34);
/**
* Get the parameters for a Drude particle.
*
* @param index the index of the Drude particle for which to get parameters
* @param particle the index within the System of the Drude particle
* @param particle1 the index within the System of the particle to which the Drude particle is attached
* @param particle2 the index within the System of the second particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso12 will be ignored.
* @param particle3 the index within the System of the third particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso34 will be ignored.
* @param particle4 the index within the System of the fourth particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso34 will be ignored.
* @param charge The charge on the Drude particle
* @param polarizability The isotropic polarizability
* @param aniso12 The scale factor for the polarizability along the direction defined by particle1 and particle2
* @param aniso34 The scale factor for the polarizability along the direction defined by particle3 and particle4
*/
void getParticleParameters(int index, int& particle, int& particle1, int& particle2, int& particle3, int& particle4, double& charge, double& polarizability, double& aniso12, double& aniso34) const;
/**
* Set the parameters for a Drude particle.
*
* @param index the index of the Drude particle for which to set parameters
* @param particle the index within the System of the Drude particle
* @param particle1 the index within the System of the particle to which the Drude particle is attached
* @param particle2 the index within the System of the second particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso12 will be ignored.
* @param particle3 the index within the System of the third particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso34 will be ignored.
* @param particle4 the index within the System of the fourth particle used for defining anisotropic polarizability.
* This may be set to -1, in which case aniso34 will be ignored.
* @param charge The charge on the Drude particle
* @param polarizability The isotropic polarizability
* @param aniso12 The scale factor for the polarizability along the direction defined by particle1 and particle2
* @param aniso34 The scale factor for the polarizability along the direction defined by particle3 and particle4
*/
void setParticleParameters(int index, int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34);
/**
* Add an interaction to the list of screenedPairs.
*
* @param particle1 the index within this Force of the first particle involved in the interaction
* @param particle2 the index within this Force of the second particle involved in the interaction
* @param thole the Thole screening factor
* @return the index of the screenedPair that was added
*/
int addScreenedPair(int particle1, int particle2, double thole);
/**
* Get the force field parameters for screened pair.
*
* @param index the index of the pair for which to get parameters
* @param particle1 the index within this Force of the first particle involved in the interaction
* @param particle2 the index within this Force of the second particle involved in the interaction
* @param thole the Thole screening factor
*/
void getScreenedPairParameters(int index, int& particle1, int& particle2, double& thole) const;
/**
* Set the force field parameters for screened pair.
*
* @param index the index of the pair for which to get parameters
* @param particle1 the index within this Force of the first particle involved in the interaction
* @param particle2 the index within this Force of the second particle involved in the interaction
* @param thole the Thole screening factor
*/
void setScreenedPairParameters(int index, int particle1, int particle2, double thole);
/**
* Update the particle and screenedPair parameters in a Context to match those stored in this Force object. This method
* provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
* Simply call setParticleParameters() and setScreenedPairParameters() to modify this object's parameters, then call
* updateParametersInState() to copy them over to the Context.
*
* This method has several limitations. It can be used to modify the numeric parameters associated with a particle or
* screened pair (polarizability, thole, etc.), but not the identities of the particles they involve. It also cannot
* be used to add new particles or screenedPairs, only to change the parameters of existing ones.
*/
void updateParametersInContext(Context& context);
protected:
ForceImpl* createImpl() const;
private:
class ParticleInfo;
class ScreenedPairInfo;
std::vector<ParticleInfo> particles;
std::vector<ScreenedPairInfo> screenedPairs;
};
/**
* This is an internal class used to record information about a particle.
* @private
*/
class DrudeForce::ParticleInfo {
public:
int particle, particle1, particle2, particle3, particle4;
double charge, polarizability, aniso12, aniso34;
ParticleInfo() {
particle = particle1 = particle2 = particle3 = particle4 = -1;
charge = polarizability = aniso12 = aniso34 = 0.0;
}
ParticleInfo(int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34) :
particle(particle), particle1(particle1), particle2(particle2), particle3(particle3), particle4(particle4), charge(charge), polarizability(polarizability), aniso12(aniso12), aniso34(aniso34) {
}
};
/**
* This is an internal class used to record information about a screenedPair.
* @private
*/
class DrudeForce::ScreenedPairInfo {
public:
int particle1, particle2;
double thole;
ScreenedPairInfo() {
particle1 = particle2 = -1;
thole = 0.0;
}
ScreenedPairInfo(int particle1, int particle2, double thole) :
particle1(particle1), particle2(particle2), thole(thole) {
}
};
} // namespace OpenMM
#endif /*OPENMM_DRUDEFORCE_H_*/
plugins/drude/openmmapi/include/openmm/DrudeKernels.h 0 → 100644
View file @ 31d761d5
#ifndef DRUDE_KERNELS_H_
#define 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/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/Vec3.h"
#include <string>
#include <vector>
namespace OpenMM {
/**
* This kernel is invoked by DrudeForce to calculate the forces acting on the system and the energy of the system.
*/
class CalcDrudeForceKernel : public KernelImpl {
public:
static std::string Name() {
return "CalcDrudeForce";
}
CalcDrudeForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the DrudeForce this kernel will be used for
*/
virtual void initialize(const System& system, const DrudeForce& force) = 0;
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the DrudeForce to copy the parameters from
*/
virtual void copyParametersToContext(ContextImpl& context, const DrudeForce& force) = 0;
};
/**
* This kernel is invoked by DrudeLangevinIntegrator to take one time step.
*/
class IntegrateDrudeLangevinStepKernel : public KernelImpl {
public:
static std::string Name() {
return "IntegrateDrudeLangevinStep";
}
IntegrateDrudeLangevinStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* 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
*/
virtual void initialize(const System& system, const DrudeLangevinIntegrator& integrator, const DrudeForce& force) = 0;
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the DrudeLangevinIntegrator this kernel is being used for
*/
virtual void execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) = 0;
/**
* Compute the kinetic energy.
*/
virtual double computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) = 0;
};
} // namespace OpenMM
#endif /*DRUDE_KERNELS_H_*/
plugins/drude/openmmapi/include/openmm/DrudeLangevinIntegrator.h 0 → 100644
View file @ 31d761d5
#ifndef OPENMM_DRUDELANGEVININTEGRATOR_H_
#define OPENMM_DRUDELANGEVININTEGRATOR_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-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/Integrator.h"
#include "openmm/Kernel.h"
#include "openmm/internal/windowsExportDrude.h"
namespace OpenMM {
/**
* This Integrator simulates systems that include Drude particles. It applies two different Langevin
* thermostats to different parts of the system. The first is applied to ordinary particles (ones that
* are not part of a Drude particle pair), as well as to the center of mass of each Drude particle pair.
* A second thermostat, typically with a much lower temperature, is applied to the relative internal
* displacement of each pair.
*
* This Integrator requires the System to include a DrudeForce, which it uses to identify the Drude
* particles.
*/
class OPENMM_EXPORT_DRUDE DrudeLangevinIntegrator : public Integrator {
public:
/**
* Create a DrudeLangevinIntegrator.
*
* @param temperature the temperature of the main heat bath (in Kelvin)
* @param frictionCoeff the friction coefficient which couples the system to the main heat bath (in inverse picoseconds)
* @param drudeTemperature the temperature of the heat bath applied to internal coordinates of Drude particles (in Kelvin)
* @param drudeFrictionCoeff the friction coefficient which couples the system to the heat bath applied to internal coordinates of Drude particles (in inverse picoseconds)
* @param stepSize the step size with which to integrator the system (in picoseconds)
*/
DrudeLangevinIntegrator(double temperature, double frictionCoeff, double drudeTemperature, double drudeFrictionCoeff, double stepSize);
/**
* Get the temperature of the main heat bath (in Kelvin).
*
* @return the temperature of the heat bath, measured in Kelvin
*/
double getTemperature() const {
return temperature;
}
/**
* Set the temperature of the main heat bath (in Kelvin).
*
* @param temp the temperature of the heat bath, measured in Kelvin
*/
void setTemperature(double temp) {
temperature = temp;
}
/**
* Get the friction coefficient which determines how strongly the system is coupled to
* the main heat bath (in inverse ps).
*
* @return the friction coefficient, measured in 1/ps
*/
double getFriction() const {
return friction;
}
/**
* Set the friction coefficient which determines how strongly the system is coupled to
* the main heat bath (in inverse ps).
*
* @param coeff the friction coefficient, measured in 1/ps
*/
void setFriction(double coeff) {
friction = coeff;
}
/**
* Get the temperature of the heat bath applied to internal coordinates of Drude particles (in Kelvin).
*
* @return the temperature of the heat bath, measured in Kelvin
*/
double getDrudeTemperature() const {
return drudeTemperature;
}
/**
* Set the temperature of the heat bath applied to internal coordinates of Drude particles (in Kelvin).
*
* @param temp the temperature of the heat bath, measured in Kelvin
*/
void setDrudeTemperature(double temp) {
drudeTemperature = temp;
}
/**
* Get the friction coefficient which determines how strongly the internal coordinates of Drude particles
* are coupled to the heat bath (in inverse ps).
*
* @return the friction coefficient, measured in 1/ps
*/
double getDrudeFriction() const {
return drudeFriction;
}
/**
* Set the friction coefficient which determines how strongly the internal coordinates of Drude particles
* are coupled to the heat bath (in inverse ps).
*
* @param coeff the friction coefficient, measured in 1/ps
*/
void setDrudeFriction(double coeff) {
drudeFriction = coeff;
}
/**
* Get the random number seed. See setRandomNumberSeed() for details.
*/
int getRandomNumberSeed() const {
return randomNumberSeed;
}
/**
* Set the random number seed. The precise meaning of this parameter is undefined, and is left up
* to each Platform to interpret in an appropriate way. It is guaranteed that if two simulations
* are run with different random number seeds, the sequence of random forces will be different. On
* the other hand, no guarantees are made about the behavior of simulations that use the same seed.
* In particular, Platforms are permitted to use non-deterministic algorithms which produce different
* results on successive runs, even if those runs were initialized identically.
*/
void setRandomNumberSeed(int seed) {
randomNumberSeed = seed;
}
/**
* Advance a simulation through time by taking a series of time steps.
*
* @param steps the number of time steps to take
*/
void step(int steps);
protected:
/**
* This will be called by the Context when it is created. It informs the Integrator
* of what context it will be integrating, and gives it a chance to do any necessary initialization.
* It will also get called again if the application calls reinitialize() on the Context.
*/
void initialize(ContextImpl& context);
/**
* This will be called by the Context when it is destroyed to let the Integrator do any necessary
* cleanup. It will also get called again if the application calls reinitialize() on the Context.
*/
void cleanup();
/**
* When the user modifies the state, we need to mark that the forces need to be recalculated.
*/
void stateChanged(State::DataType changed);
/**
* Get the names of all Kernels used by this Integrator.
*/
std::vector<std::string> getKernelNames();
/**
* Compute the kinetic energy of the system at the current time.
*/
double computeKineticEnergy();
private:
double temperature, friction, drudeTemperature, drudeFriction;
int randomNumberSeed;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_DRUDELANGEVININTEGRATOR_H_*/
plugins/drude/openmmapi/include/openmm/internal/DrudeForceImpl.h 0 → 100644
View file @ 31d761d5
#ifndef OPENMM_DRUDEFORCEIMPL_H_
#define OPENMM_DRUDEFORCEIMPL_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/DrudeForce.h"
#include "openmm/internal/ForceImpl.h"
#include "openmm/Kernel.h"
#include <utility>
#include <set>
#include <string>
namespace OpenMM {
class System;
/**
* This is the internal implementation of DrudeForce.
*/
class OPENMM_EXPORT DrudeForceImpl : public ForceImpl {
public:
DrudeForceImpl(const DrudeForce& owner);
~DrudeForceImpl();
void initialize(ContextImpl& context);
const DrudeForce& getOwner() const {
return owner;
}
void updateContextState(ContextImpl& context) {
// This force field doesn't update the state directly.
}
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
std::map<std::string, double> getDefaultParameters() {
return std::map<std::string, double>(); // This force field doesn't define any parameters.
}
std::vector<std::string> getKernelNames();
void updateParametersInContext(ContextImpl& context);
private:
const DrudeForce& owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_DRUDEFORCEIMPL_H_*/
plugins/drude/openmmapi/include/openmm/internal/windowsExportDrude.h 0 → 100644
View file @ 31d761d5
#ifndef OPENMM_WINDOWSEXPORTDRUDE_H_
#define OPENMM_WINDOWSEXPORTDRUDE_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
* OPENMM_DRUDE_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_DRUDE_BUILDING_SHARED_LIBRARY)
#define OPENMM_EXPORT_DRUDE __declspec(dllexport)
#elif defined(OPENMM_DRUDE_BUILDING_STATIC_LIBRARY) || defined(OPENMM_DRUDE_USE_STATIC_LIBRARIES)
#define OPENMM_EXPORT_DRUDE
#else
#define OPENMM_EXPORT_DRUDE __declspec(dllimport) // i.e., a client of a shared library
#endif
#else
#define OPENMM_EXPORT_DRUDE // Linux, Mac
#endif
#endif // OPENMM_WINDOWSEXPORTDRUDE_H_
plugins/drude/openmmapi/src/DrudeForce.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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/Force.h"
#include "openmm/OpenMMException.h"
#include "openmm/DrudeForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/internal/DrudeForceImpl.h"
using namespace OpenMM;
using namespace std;
DrudeForce::DrudeForce() {
}
int DrudeForce::addParticle(int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34) {
particles.push_back(ParticleInfo(particle, particle1, particle2, particle3, particle4, charge, polarizability, aniso12, aniso34));
return particles.size()-1;
}
void DrudeForce::getParticleParameters(int index, int& particle, int& particle1, int& particle2, int& particle3, int& particle4, double& charge, double& polarizability, double& aniso12, double& aniso34) const {
ASSERT_VALID_INDEX(index, particles);
particle = particles[index].particle;
particle1 = particles[index].particle1;
particle2 = particles[index].particle2;
particle3 = particles[index].particle3;
particle4 = particles[index].particle4;
charge = particles[index].charge;
polarizability = particles[index].polarizability;
aniso12 = particles[index].aniso12;
aniso34 = particles[index].aniso34;
}
void DrudeForce::setParticleParameters(int index, int particle, int particle1, int particle2, int particle3, int particle4, double charge, double polarizability, double aniso12, double aniso34) {
ASSERT_VALID_INDEX(index, particles);
particles[index].particle = particle;
particles[index].particle1 = particle1;
particles[index].particle2 = particle2;
particles[index].particle3 = particle3;
particles[index].particle4 = particle4;
particles[index].charge = charge;
particles[index].polarizability = polarizability;
particles[index].aniso12 = aniso12;
particles[index].aniso34 = aniso34;
}
int DrudeForce::addScreenedPair(int particle1, int particle2, double thole) {
screenedPairs.push_back(ScreenedPairInfo(particle1, particle2, thole));
return screenedPairs.size()-1;
}
void DrudeForce::getScreenedPairParameters(int index, int& particle1, int& particle2, double& thole) const {
ASSERT_VALID_INDEX(index, screenedPairs);
particle1 = screenedPairs[index].particle1;
particle2 = screenedPairs[index].particle2;
thole = screenedPairs[index].thole;
}
void DrudeForce::setScreenedPairParameters(int index, int particle1, int particle2, double thole) {
ASSERT_VALID_INDEX(index, screenedPairs);
screenedPairs[index].particle1 = particle1;
screenedPairs[index].particle2 = particle2;
screenedPairs[index].thole = thole;
}
ForceImpl* DrudeForce::createImpl() const {
return new DrudeForceImpl(*this);
}
void DrudeForce::updateParametersInContext(Context& context) {
dynamic_cast<DrudeForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
}
plugins/drude/openmmapi/src/DrudeForceImpl.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/DrudeForceImpl.h"
#include "openmm/DrudeKernels.h"
#include <cmath>
#include <map>
#include <set>
#include <sstream>
using namespace OpenMM;
using namespace std;
DrudeForceImpl::DrudeForceImpl(const DrudeForce& owner) : owner(owner) {
}
DrudeForceImpl::~DrudeForceImpl() {
}
void DrudeForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcDrudeForceKernel::Name(), context);
const System& system = context.getSystem();
// Check for errors in the specification of particles.
set<int> usedParticles;
for (int i = 0; i < owner.getNumParticles(); i++) {
int particle, particle1, particle2, particle3, particle4;
double charge, k, k2, k3;
owner.getParticleParameters(i, particle, particle1, particle2, particle3, particle4, charge, k, k2, k3);
if (particle < 0 || particle >= system.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index: ";
msg << particle;
throw OpenMMException(msg.str());
}
if (particle1 < 0 || particle1 >= system.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < -1 || particle2 >= system.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index: ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (particle3 < -1 || particle3 >= system.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index: ";
msg << particle3;
throw OpenMMException(msg.str());
}
if (particle4 < -1 || particle4 >= system.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index: ";
msg << particle4;
throw OpenMMException(msg.str());
}
if (usedParticles.find(particle) != usedParticles.end()) {
stringstream msg;
msg << "DrudeForce: Particle index is used by two different Drude particles: ";
msg << particle;
throw OpenMMException(msg.str());
}
usedParticles.insert(particle);
if (usedParticles.find(particle1) != usedParticles.end()) {
stringstream msg;
msg << "DrudeForce: Particle index is used by two different Drude particles: ";
msg << particle1;
throw OpenMMException(msg.str());
}
usedParticles.insert(particle1);
}
// Check for errors in the specification of screened pairs.
vector<set<int> > screenedPairs(owner.getNumParticles());
for (int i = 0; i < owner.getNumScreenedPairs(); i++) {
int particle1, particle2;
double thole;
owner.getScreenedPairParameters(i, particle1, particle2, thole);
if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index for a screened pair: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
stringstream msg;
msg << "DrudeForce: Illegal particle index for a screened pair: ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (screenedPairs[particle1].count(particle2) > 0 || screenedPairs[particle2].count(particle1) > 0) {
stringstream msg;
msg << "DrudeForce: Multiple screened pairs are specified for particles ";
msg << particle1;
msg << " and ";
msg << particle2;
throw OpenMMException(msg.str());
}
screenedPairs[particle1].insert(particle2);
screenedPairs[particle2].insert(particle1);
}
kernel.getAs<CalcDrudeForceKernel>().initialize(context.getSystem(), owner);
}
double DrudeForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
return kernel.getAs<CalcDrudeForceKernel>().execute(context, includeForces, includeEnergy);
}
std::vector<std::string> DrudeForceImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(CalcDrudeForceKernel::Name());
return names;
}
void DrudeForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcDrudeForceKernel>().copyParametersToContext(context, owner);
}
plugins/drude/openmmapi/src/DrudeLangevinIntegrator.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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-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/DrudeLangevinIntegrator.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/DrudeKernels.h"
#include <cmath>
#include <ctime>
#include <string>
using namespace OpenMM;
using std::string;
using std::vector;
DrudeLangevinIntegrator::DrudeLangevinIntegrator(double temperature, double frictionCoeff, double drudeTemperature, double drudeFrictionCoeff, double stepSize) {
setTemperature(temperature);
setFriction(frictionCoeff);
setDrudeTemperature(drudeTemperature);
setDrudeFriction(drudeFrictionCoeff);
setStepSize(stepSize);
setConstraintTolerance(1e-5);
setRandomNumberSeed((int) time(NULL));
}
void DrudeLangevinIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
const DrudeForce* force = NULL;
const System& system = contextRef.getSystem();
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
if (force == NULL)
force = dynamic_cast<const DrudeForce*>(&system.getForce(i));
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (force == NULL)
throw OpenMMException("The System does not contain a DrudeForce");
context = &contextRef;
owner = &contextRef.getOwner();
kernel = context->getPlatform().createKernel(IntegrateDrudeLangevinStepKernel::Name(), contextRef);
kernel.getAs<IntegrateDrudeLangevinStepKernel>().initialize(contextRef.getSystem(), *this, *force);
}
void DrudeLangevinIntegrator::cleanup() {
kernel = Kernel();
}
void DrudeLangevinIntegrator::stateChanged(State::DataType changed) {
}
vector<string> DrudeLangevinIntegrator::getKernelNames() {
std::vector<std::string> names;
names.push_back(IntegrateDrudeLangevinStepKernel::Name());
return names;
}
double DrudeLangevinIntegrator::computeKineticEnergy() {
return kernel.getAs<IntegrateDrudeLangevinStepKernel>().computeKineticEnergy(*context, *this);
}
void DrudeLangevinIntegrator::step(int steps) {
for (int i = 0; i < steps; ++i) {
context->updateContextState();
context->calcForcesAndEnergy(true, false);
kernel.getAs<IntegrateDrudeLangevinStepKernel>().execute(*context, *this);
}
}
plugins/drude/platforms/reference/include/ReferenceDrudeKernelFactory.h 0 → 100644
View file @ 31d761d5
#ifndef OPENMM_REFERENCEDRUDEKERNELFACTORY_H_
#define OPENMM_REFERENCEDRUDEKERNELFACTORY_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 reference implementation of the Drude plugin.
*/
class ReferenceDrudeKernelFactory : public KernelFactory {
public:
KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCEDRUDEKERNELFACTORY_H_*/
plugins/drude/platforms/reference/src/ReferenceDrudeKernelFactory.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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 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 "ReferenceDrudeKernelFactory.h"
#include "ReferenceDrudeKernels.h"
#include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
#if defined(WIN32)
#include <windows.h>
extern "C" void initDrudeReferenceKernels();
BOOL WINAPI DllMain(HANDLE hModule, DWORD ul_reason_for_call, LPVOID lpReserved) {
if (ul_reason_for_call == DLL_PROCESS_ATTACH)
initDrudeReferenceKernels();
return TRUE;
}
#else
extern "C" void __attribute__((constructor)) initDrudeReferenceKernels();
#endif
extern "C" void initDrudeReferenceKernels() {
Platform& platform = Platform::getPlatformByName("Reference");
ReferenceDrudeKernelFactory* factory = new ReferenceDrudeKernelFactory();
platform.registerKernelFactory(CalcDrudeForceKernel::Name(), factory);
platform.registerKernelFactory(IntegrateDrudeLangevinStepKernel::Name(), factory);
}
KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
if (name == CalcDrudeForceKernel::Name())
return new ReferenceCalcDrudeForceKernel(name, platform);
if (name == IntegrateDrudeLangevinStepKernel::Name())
return new ReferenceIntegrateDrudeLangevinStepKernel(name, platform);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
}
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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) 2011-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 "ReferenceDrudeKernels.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include "SimTKReference/ReferenceCCMAAlgorithm.h"
#include <set>
using namespace OpenMM;
using namespace std;
static vector<RealVec>& extractPositions(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->positions);
}
static vector<RealVec>& extractVelocities(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->velocities);
}
static vector<RealVec>& extractForces(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->forces);
}
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
for (int i = 0; i < system.getNumForces(); i++) {
const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
if (force != NULL) {
for (int j = 0; j < force->getNumAngles(); j++) {
int atom1, atom2, atom3;
double angle, k;
force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
}
}
}
}
void ReferenceCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) {
// Initialize particle parameters.
int numParticles = force.getNumParticles();
particle.resize(numParticles);
particle1.resize(numParticles);
particle2.resize(numParticles);
particle3.resize(numParticles);
particle4.resize(numParticles);
charge.resize(numParticles);
polarizability.resize(numParticles);
aniso12.resize(numParticles);
aniso34.resize(numParticles);
for (int i = 0; i < numParticles; i++)
force.getParticleParameters(i, particle[i], particle1[i], particle2[i], particle3[i], particle4[i], charge[i], polarizability[i], aniso12[i], aniso34[i]);
// Initialize screened pair parameters.
int numPairs = force.getNumScreenedPairs();
pair1.resize(numPairs);
pair2.resize(numPairs);
pairThole.resize(numPairs);
for (int i = 0; i < numPairs; i++)
force.getScreenedPairParameters(i, pair1[i], pair2[i], pairThole[i]);
}
double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& pos = extractPositions(context);
vector<RealVec>& force = extractForces(context);
int numParticles = particle.size();
double energy = 0;
// Compute the interactions from the harmonic springs.
for (int i = 0; i < numParticles; i++) {
int p = particle[i];
int p1 = particle1[i];
int p2 = particle2[i];
int p3 = particle3[i];
int p4 = particle4[i];
RealOpenMM a1 = (p2 == -1 ? 1 : aniso12[i]);
RealOpenMM a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34[i]);
RealOpenMM a3 = 3-a1-a2;
RealOpenMM k3 = charge[i]*charge[i]/(polarizability[i]*a3);
RealOpenMM k1 = charge[i]*charge[i]/(polarizability[i]*a1) - k3;
RealOpenMM k2 = charge[i]*charge[i]/(polarizability[i]*a2) - k3;
// Compute the isotropic force.
RealVec delta = pos[p]-pos[p1];
RealOpenMM r2 = delta.dot(delta);
energy += 0.5*k3*r2;
force[p] -= delta*k3;
force[p1] += delta*k3;
// Compute the first anisotropic force.
if (p2 != -1) {
RealVec dir = pos[p1]-pos[p2];
RealOpenMM invDist = 1.0/sqrt(dir.dot(dir));
dir *= invDist;
RealOpenMM rprime = dir.dot(delta);
energy += 0.5*k1*rprime*rprime;
RealVec f1 = dir*(k1*rprime);
RealVec f2 = (delta-dir*rprime)*(k1*rprime*invDist);
force[p] -= f1;
force[p1] += f1-f2;
force[p2] += f2;
}
// Compute the second anisotropic force.
if (p3 != -1 && p4 != -1) {
RealVec dir = pos[p3]-pos[p4];
RealOpenMM invDist = 1.0/sqrt(dir.dot(dir));
dir *= invDist;
RealOpenMM rprime = dir.dot(delta);
energy += 0.5*k2*rprime*rprime;
RealVec f1 = dir*(k2*rprime);
RealVec f2 = (delta-dir*rprime)*(k2*rprime*invDist);
force[p] -= f1;
force[p1] += f1;
force[p3] -= f2;
force[p4] += f2;
}
}
// Compute the screened interaction between bonded dipoles.
int numPairs = pair1.size();
for (int i = 0; i < numPairs; i++) {
int dipole1 = pair1[i];
int dipole2 = pair2[i];
int dipole1Particles[] = {particle[dipole1], particle1[dipole1]};
int dipole2Particles[] = {particle[dipole2], particle1[dipole2]};
for (int j = 0; j < 2; j++)
for (int k = 0; k < 2; k++) {
int p1 = dipole1Particles[j];
int p2 = dipole2Particles[k];
RealOpenMM chargeProduct = charge[dipole1]*charge[dipole2]*(j == k ? 1 : -1);
RealVec delta = pos[p1]-pos[p2];
RealOpenMM r = sqrt(delta.dot(delta));
RealOpenMM u = r*pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0);
RealOpenMM screening = 1.0 - (1.0+0.5*u)*exp(-u);
energy += ONE_4PI_EPS0*chargeProduct*screening/r;
RealVec f = delta*(ONE_4PI_EPS0*chargeProduct*screening/(r*r*r));
force[p1] += f;
force[p2] -= f;
}
}
return energy;
}
void ReferenceCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context, const DrudeForce& force) {
if (force.getNumParticles() != particle.size())
throw OpenMMException("updateParametersInContext: The number of Drude particles has changed");
if (force.getNumScreenedPairs() != pair1.size())
throw OpenMMException("updateParametersInContext: The number of screened pairs has changed");
for (int i = 0; i < force.getNumParticles(); i++) {
int p, p1, p2, p3, p4;
force.getParticleParameters(i, p, p1, p2, p3, p4, charge[i], polarizability[i], aniso12[i], aniso34[i]);
if (p != particle[i] || p1 != particle1[i] || p2 != particle2[i] || p3 != particle3[i] || p4 != particle4[i])
throw OpenMMException("updateParametersInContext: A particle index has changed");
}
for (int i = 0; i < force.getNumScreenedPairs(); i++) {
int p1, p2;
force.getScreenedPairParameters(i, p1, p2, pairThole[i]);
if (p1 != pair1[i] || p2 != pair2[i])
throw OpenMMException("updateParametersInContext: A particle index for a screened pair has changed");
}
}
ReferenceIntegrateDrudeLangevinStepKernel::~ReferenceIntegrateDrudeLangevinStepKernel() {
if (constraints != NULL)
delete constraints;
}
void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system, const DrudeLangevinIntegrator& integrator, const DrudeForce& force) {
drudeForce = &force;
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
// Identify particle pairs and ordinary particles.
set<int> particles;
vector<RealOpenMM> particleMass;
for (int i = 0; i < system.getNumParticles(); i++) {
particles.insert(i);
double mass = system.getParticleMass(i);
particleMass.push_back(mass);
particleInvMass.push_back(mass == 0.0 ? 0.0 : 1.0/mass);
}
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);
pairParticles.push_back(make_pair(p, p1));
double m1 = system.getParticleMass(p);
double m2 = system.getParticleMass(p1);
pairInvTotalMass.push_back(1.0/(m1+m2));
pairInvReducedMass.push_back((m1+m2)/(m1*m2));
}
normalParticles.insert(normalParticles.begin(), particles.begin(), particles.end());
// Prepare constraints.
int numConstraints = system.getNumConstraints();
if (numConstraints > 0) {
vector<pair<int, int> > constraintIndices(numConstraints);
vector<RealOpenMM> constraintDistances(numConstraints);
for (int i = 0; i < numConstraints; ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
constraintIndices[i].first = particle1;
constraintIndices[i].second = particle2;
constraintDistances[i] = static_cast<RealOpenMM>(distance);
}
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, particleMass, angles, (RealOpenMM)integrator.getConstraintTolerance());
}
}
void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
// Update velocities of ordinary particles.
const RealOpenMM vscale = exp(-integrator.getStepSize()*integrator.getFriction());
const RealOpenMM fscale = (1-vscale)/integrator.getFriction();
const RealOpenMM kT = BOLTZ*integrator.getTemperature();
const RealOpenMM noisescale = sqrt(2*kT*integrator.getFriction())*sqrt(0.5*(1-vscale*vscale)/integrator.getFriction());
vector<RealVec>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context);
vector<RealVec>& force = extractForces(context);
for (int i = 0; i < (int) normalParticles.size(); i++) {
int index = normalParticles[i];
RealOpenMM invMass = particleInvMass[index];
if (invMass != 0.0) {
RealOpenMM sqrtInvMass = sqrt(invMass);
for (int j = 0; j < 3; j++)
vel[i][j] = vscale*vel[i][j] + fscale*invMass*force[i][j] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
// Update velocities of Drude particle pairs.
const RealOpenMM vscaleDrude = exp(-integrator.getStepSize()*integrator.getDrudeFriction());
const RealOpenMM fscaleDrude = (1-vscaleDrude)/integrator.getDrudeFriction();
const RealOpenMM kTDrude = BOLTZ*integrator.getDrudeTemperature();
const RealOpenMM noisescaleDrude = sqrt(2*kTDrude*integrator.getDrudeFriction())*sqrt(0.5*(1-vscaleDrude*vscaleDrude)/integrator.getDrudeFriction());
for (int i = 0; i < (int) pairParticles.size(); i++) {
int p1 = pairParticles[i].first;
int p2 = pairParticles[i].second;
RealOpenMM sqrtInvTotalMass = sqrt(pairInvTotalMass[i]);
RealOpenMM sqrtInvReducedMass = sqrt(pairInvReducedMass[i]);
RealVec cmVel = (vel[p2]+vel[p1])*0.5;
RealVec relVel = vel[p2]-vel[p1];
RealVec cmForce = force[p1]+force[p2];
RealVec relForce = force[p2]*(pairInvTotalMass[i]/particleInvMass[p1]) - force[p1]*(pairInvTotalMass[i]/particleInvMass[p2]);
for (int j = 0; j < 3; j++) {
cmVel[j] = vscale*cmVel[j] + fscale*pairInvTotalMass[i]*cmForce[j] + noisescale*sqrtInvTotalMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
relVel[j] = vscaleDrude*relVel[j] + fscaleDrude*pairInvReducedMass[i]*relForce[j] + noisescaleDrude*sqrtInvReducedMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
vel[p1] = cmVel-relVel*0.5;
vel[p2] = cmVel+relVel*0.5;
}
// Update the particle positions.
int numParticles = particleInvMass.size();
vector<RealVec> xPrime(numParticles);
RealOpenMM dt = integrator.getStepSize();
for (int i = 0; i < numParticles; i++)
if (particleInvMass[i] != 0.0)
xPrime[i] = pos[i]+vel[i]*dt;
// Apply constraints.
constraints->apply(numParticles, pos, xPrime, particleInvMass);
// Record the constrained positions and velocities.
RealOpenMM dtInv = 1.0/dt;
for (int i = 0; i < numParticles; i++) {
if (particleInvMass[i] != 0.0) {
vel[i] = (xPrime[i]-pos[i])*dtInv;
pos[i] = xPrime[i];
}
}
}
double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
const System& system = context.getSystem();
int numParticles = system.getNumParticles();
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context);
// Compute the shifted velocities.
vector<RealVec> shiftedVel(numParticles);
double timeShift = 0.5*integrator.getStepSize();
for (int i = 0; i < numParticles; ++i) {
if (particleInvMass[i] > 0)
shiftedVel[i] = velData[i]+forceData[i]*(timeShift/particleInvMass[i]);
else
shiftedVel[i] = velData[i];
}
// Apply constraints to them.
if (constraints != NULL) {
constraints->setTolerance(1e-4);
constraints->applyToVelocities(numParticles, posData, shiftedVel, particleInvMass);
}
// Compute the kinetic energy.
double energy = 0.0;
for (int i = 0; i < numParticles; ++i)
if (particleInvMass[i] > 0)
energy += (shiftedVel[i].dot(shiftedVel[i]))/particleInvMass[i];
return 0.5*energy;
}
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h 0 → 100644
View file @ 31d761d5
#ifndef REFERENCE_DRUDE_KERNELS_H_
#define REFERENCE_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 "ReferencePlatform.h"
#include "openmm/DrudeKernels.h"
#include "SimTKUtilities/RealVec.h"
#include <utility>
#include <vector>
class ReferenceConstraintAlgorithm;
namespace OpenMM {
/**
* This kernel is invoked by DrudeForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcDrudeForceKernel : public CalcDrudeForceKernel {
public:
ReferenceCalcDrudeForceKernel(std::string name, const Platform& platform) : CalcDrudeForceKernel(name, platform) {
}
/**
* 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:
std::vector<int> particle, particle1, particle2, particle3, particle4;
std::vector<double> charge, polarizability, aniso12, aniso34;
std::vector<int> pair1, pair2;
std::vector<double> pairThole;
};
/**
* This kernel is invoked by DrudeLangevinIntegrator to take one time step
*/
class ReferenceIntegrateDrudeLangevinStepKernel : public IntegrateDrudeLangevinStepKernel {
public:
ReferenceIntegrateDrudeLangevinStepKernel(std::string name, const Platform& platform) : IntegrateDrudeLangevinStepKernel(name, platform), constraints(NULL) {
}
~ReferenceIntegrateDrudeLangevinStepKernel();
/**
* 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:
const DrudeForce* drudeForce;
std::vector<int> normalParticles;
std::vector<std::pair<int, int> > pairParticles;
std::vector<double> particleInvMass;
std::vector<double> pairInvTotalMass;
std::vector<double> pairInvReducedMass;
ReferenceConstraintAlgorithm* constraints;
};
} // namespace OpenMM
#endif /*REFERENCE_DRUDE_KERNELS_H_*/
plugins/drude/platforms/reference/tests/CMakeLists.txt 0 → 100644
View file @ 31d761d5
#
# 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_OPENMM_DRUDE_TARGET OpenMMDrude)
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d)
SET(SHARED_OPENMM_DRUDE_TARGET ${SHARED_OPENMM_DRUDE_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_OPENMM_DRUDE_TARGET})
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
ENDFOREACH(TEST_PROG ${TEST_PROGS})
plugins/drude/platforms/reference/tests/TestReferenceDrudeForce.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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 Reference 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 "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
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("Reference");
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() {
try {
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/reference/tests/TestReferenceDrudeLangevinIntegrator.cpp 0 → 100644
View file @ 31d761d5
/* -------------------------------------------------------------------------- *
* 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 Reference 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/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
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, 10.0, temperatureDrude, 10.0, 0.003);
Platform& platform = Platform::getPlatformByName("Reference");
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);
}
int main() {
try {
testSinglePair();
}
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