Commit 18a837c0 authored by Peter Eastman's avatar Peter Eastman
Browse files

Initial implementation of RPMD plugin

parent 839ab51e
This diff is collapsed.
...@@ -32,10 +32,34 @@ extern "C" { ...@@ -32,10 +32,34 @@ extern "C" {
#endif #endif
typedef struct { class t_complex {
RealOpenMM re; public:
RealOpenMM im; RealOpenMM re;
} t_complex; RealOpenMM im;
t_complex() : re(0.0), im(0.0) {
}
t_complex(RealOpenMM re, RealOpenMM im) : re(re), im(im) {
}
t_complex(const t_complex& c) : re(c.re), im(c.im) {
}
t_complex operator*(RealOpenMM r) {
return t_complex(re*r, im*r);
}
t_complex operator+(const t_complex& c) const {
return t_complex(re+c.re, im+c.im);
}
t_complex operator-(const t_complex& c) const {
return t_complex(re-c.re, im-c.im);
}
t_complex operator+=(const t_complex& c) {
re += c.re;
im += c.im;
}
t_complex operator-=(const t_complex& c) {
re -= c.re;
im -= c.im;
}
};
/*! \brief Datatype for FFT setup /*! \brief Datatype for FFT setup
......
#---------------------------------------------------
# OpenMM RPMD Plugin
#
# Creates OpenMM RPMD plugin library, base name=OpenMMRPMD.
# Default libraries are shared & optimized. Variants
# are created for static (_static) and debug (_d).
#
# Windows:
# OpenMMRPMD[_d].dll
# OpenMMRPMD[_d].lib
# OpenMMRPMD_static[_d].lib
# Unix:
# libOpenMMRPMD[_d].so
# libOpenMMRPMD_static[_d].a
#----------------------------------------------------
#INCLUDE(Dart)
# ----------------------------------------------------------------------------
SET(CREATE_SERIALIZABLE_OPENMM_RPMD 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_RPMD_PLUGIN_SOURCE_SUBDIRS . openmmapi platforms/reference)
SET(OPENMM_RPMD_LIBRARY_NAME OpenMMRPMD)
SET(SHARED_RPMD_TARGET ${OPENMM_RPMD_LIBRARY_NAME})
SET(STATIC_RPMD_TARGET ${OPENMM_RPMD_LIBRARY_NAME}_static)
IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
SET(SHARED_RPMD_SERIALIZABLE_TARGET ${OPENMM_RPMD_LIBRARY_NAME}_serializable)
ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
# But on Unix or Cygwin we have to add the suffix manually
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_RPMD_TARGET ${SHARED_RPMD_TARGET}_d)
SET(STATIC_RPMD_TARGET ${STATIC_RPMD_TARGET}_d)
IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
SET(SHARED_RPMD_SERIALIZABLE_TARGET ${SHARED_RPMD_SERIALIZABLE_TARGET}_d)
ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
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_RPMD_INCLUDE_DIRS) # start empty
FOREACH(subdir ${OPENMM_RPMD_PLUGIN_SOURCE_SUBDIRS})
# append
SET(API_RPMD_INCLUDE_DIRS ${API_RPMD_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_RPMD_REL_INCLUDE_FILES) # start these out empty
SET(API_RPMD_ABS_INCLUDE_FILES)
FOREACH(dir ${API_RPMD_INCLUDE_DIRS})
FILE(GLOB fullpaths ${dir}/*.h) # returns full pathnames
SET(API_RPMD_ABS_INCLUDE_FILES ${API_RPMD_ABS_INCLUDE_FILES} ${fullpaths})
FOREACH(pathname ${fullpaths})
GET_FILENAME_COMPONENT(filename ${pathname} NAME)
SET(API_RPMD_REL_INCLUDE_FILES ${API_RPMD_REL_INCLUDE_FILES} ${dir}/${filename})
ENDFOREACH(pathname)
ENDFOREACH(dir)
# collect up source files
SET(SOURCE_RPMD_FILES) # empty
SET(SOURCE_RPMD_INCLUDE_FILES)
FOREACH(subdir ${OPENMM_RPMD_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_RPMD_FILES ${SOURCE_RPMD_FILES} ${src_files}) #append
SET(SOURCE_RPMD_INCLUDE_FILES ${SOURCE_RPMD_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_RPMD wrappers are being generated, and add them to the build.
#IF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
# ADD_SUBDIRECTORY(wrappers)
# SET(SOURCE_RPMD_FILES ${SOURCE_RPMD_FILES} wrappers/RPMDOpenMMCWrapper.cpp wrappers/RPMDOpenMMFortranWrapper.cpp)
# SET_SOURCE_FILES_PROPERTIES(wrappers/RPMDOpenMMCWrapper.cpp wrappers/RPMDOpenMMFortranWrapper.cpp PROPERTIES GENERATED TRUE)
#ENDIF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
ADD_LIBRARY(${SHARED_RPMD_TARGET} SHARED ${SOURCE_RPMD_FILES} ${SOURCE_RPMD_INCLUDE_FILES} ${API_RPMD_ABS_INCLUDE_FILES})
SET_TARGET_PROPERTIES(${SHARED_RPMD_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY")
#IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
# ADD_LIBRARY(${SHARED_RPMD_SERIALIZABLE_TARGET} SHARED ${SOURCE_RPMD_FILES} ${SOURCE_RPMD_INCLUDE_FILES} ${API_RPMD_ABS_INCLUDE_FILES})
# SET_TARGET_PROPERTIES(${SHARED_RPMD_SERIALIZABLE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_SERIALIZE")
# INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/../../serialization/include)
#ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
IF(OPENMM_BUILD_STATIC_LIB)
ADD_LIBRARY(${STATIC_RPMD_TARGET} STATIC ${SOURCE_RPMD_FILES} ${SOURCE_RPMD_INCLUDE_FILES} ${API_RPMD_ABS_INCLUDE_FILES})
SET_TARGET_PROPERTIES(${STATIC_RPMD_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_RPMD_TARGET} RPMDApiWrappers)
# IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
# ADD_DEPENDENCIES(${SHARED_RPMD_SERIALIZABLE_TARGET} RPMDApiWrappers)
# ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
# IF(OPENMM_BUILD_STATIC_LIB)
# ADD_DEPENDENCIES(${STATIC_RPMD_TARGET} RPMDApiWrappers)
# 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_RPMD_TARGET} ${DL_LIBRARY})
IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
TARGET_LINK_LIBRARIES(${SHARED_RPMD_SERIALIZABLE_TARGET} ${DL_LIBRARY})
ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
IF(OPENMM_BUILD_STATIC_LIB)
TARGET_LINK_LIBRARIES(${STATIC_RPMD_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_RPMD_TARGET} ${SHARED_TARGET} )
IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
TARGET_LINK_LIBRARIES( ${SHARED_RPMD_SERIALIZABLE_TARGET} ${SHARED_TARGET} )
ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
IF(OPENMM_BUILD_STATIC_LIB)
TARGET_LINK_LIBRARIES( ${STATIC_RPMD_TARGET} ${STATIC_TARGET} )
ENDIF(OPENMM_BUILD_STATIC_LIB)
#ADD_SUBDIRECTORY(platforms/reference/tests)
# Which hardware platforms to build
IF(CUDA_FOUND)
SET(OPENMM_BUILD_RPMD_CUDA_LIB ON CACHE BOOL "Build OpenMMRPMDCuda library for Nvidia GPUs")
ELSE(CUDA_FOUND)
SET(OPENMM_BUILD_RPMD_CUDA_LIB OFF CACHE BOOL "Build OpenMMRPMDCuda library for Nvidia GPUs")
ENDIF(CUDA_FOUND)
#SET(OPENMM_BUILD_RPMD_PATH)
#SET(OPENMM_BUILD_RPMD_CUDA_PATH)
#IF(OPENMM_BUILD_RPMD_CUDA_LIB)
# ADD_SUBDIRECTORY(platforms/cuda)
# SET(OPENMM_BUILD_RPMD_PATH ${CMAKE_CURRENT_SOURCE_DIR}/platforms/cuda)
# SET(OPENMM_BUILD_RPMD_CUDA_PATH ${CMAKE_CURRENT_SOURCE_DIR}/platforms/cuda)
# SET(OPENMM_RPMD_CUDA_SOURCE_SUBDIRS . openmmapi olla platforms/cuda)
#ENDIF(OPENMM_BUILD_RPMD_CUDA_LIB)
#SET(OPENMM_BUILD_OPENCL_LIB OFF CACHE BOOL "Build OpenMMOpenCL library for Nvidia GPUs")
#IF(OPENMM_BUILD_OPENCL_LIB)
# ADD_SUBDIRECTORY(platforms/opencl)
#ENDIF(OPENMM_BUILD_OPENCL_LIB)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_RPMD_TARGET})
IF( CREATE_SERIALIZABLE_OPENMM_RPMD )
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_RPMD_SERIALIZABLE_TARGET})
ENDIF( CREATE_SERIALIZABLE_OPENMM_RPMD )
IF(OPENMM_BUILD_STATIC_LIB)
INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_RPMD_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)
#ADD_SUBDIRECTORY(examples)
#ifndef OPENMM_RPMDINTEGRATOR_H_
#define OPENMM_RPMDINTEGRATOR_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-2011 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/State.h"
#include "openmm/Vec3.h"
#include "openmm/internal/windowsExport.h"
namespace OpenMM {
/**
* This is an Integrator which simulates a System using RPMD dynamics.
*/
class OPENMM_EXPORT RPMDIntegrator : public Integrator {
public:
/**
* Create a RPMDIntegrator.
*
* @param numCopies the number of copies of the system that should be simulated
* @param temperature the temperature of the heat bath (in Kelvin)
* @param frictionCoeff the friction coefficient which couples the system to the heat bath (in inverse picoseconds)
* @param stepSize the step size with which to integrator the system (in picoseconds)
*/
RPMDIntegrator(int numCopies, double temperature, double frictionCoeff, double stepSize);
/**
* Get the number of copies of the system being simulated.
*/
int getNumCopies() const {
return numCopies;
}
/**
* Get the temperature of the heat bath (in Kelvin).
*
* @return the temperature of the heat bath, measured in Kelvin
*/
double getTemperature() const {
return temperature;
}
/**
* Set the temperature of the 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 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 heat bath (in inverse ps).
*
* @param coeff the friction coefficient, measured in 1/ps
*/
void setFriction(double coeff) {
friction = 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;
}
/**
* Set the positions of all particles in one copy of the system.
*
* @param copy the index of the copy for which to set positions
* @param positions the positions of all particles in the system
*/
void setPositions(int copy, const std::vector<Vec3>& positions);
/**
* Get the velocities of all particles in one copy of the system.
*
* @param copy the index of the copy for which to set velocities
* @param velocities the velocities of all particles in the system
*/
void setVelocities(int copy, const std::vector<Vec3>& velocities);
/**
* Get a State object recording the current state information about one copy of the system.
*
* @param copy the index of the copy for which to retrieve state information
* @param types the set of data types which should be stored in the State object. This
* should be a union of DataType values, e.g. (State::Positions | State::Velocities).
*/
State getState(int copy, int types);
/**
* 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);
/**
* Get the names of all Kernels used by this Integrator.
*/
std::vector<std::string> getKernelNames();
private:
double temperature, friction;
int numCopies, randomNumberSeed;
ContextImpl* context;
Context* owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_RPMDINTEGRATOR_H_*/
#ifndef RPMD_KERNELS_H_
#define RPMD_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) 2011 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/RPMDIntegrator.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/Vec3.h"
#include <string>
#include <vector>
namespace OpenMM {
/**
* This kernel is invoked by RPMDIntegrator to take one time step, and to get and
* set the state of system copies.
*/
class IntegrateRPMDStepKernel : public KernelImpl {
public:
static std::string Name() {
return "IntegrateRPMDStep";
}
IntegrateRPMDStepKernel(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 RPMDIntegrator this kernel will be used for
*/
virtual void initialize(const System& system, const RPMDIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the RPMDIntegrator this kernel is being used for
*/
virtual void execute(ContextImpl& context, const RPMDIntegrator& integrator) = 0;
/**
* Get the positions of all particles in one copy of the system.
*/
virtual void setPositions(int copy, const std::vector<Vec3>& positions) = 0;
/**
* Get the velocities of all particles in one copy of the system.
*/
virtual void setVelocities(int copy, const std::vector<Vec3>& velocities) = 0;
/**
* Copy positions and velocities for one copy into the context.
*/
virtual void copyToContext(int copy, ContextImpl& context) const = 0;
};
} // namespace OpenMM
#endif /*RPMD_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) 2008-2011 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/RPMDIntegrator.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/RPMDKernels.h"
#include <ctime>
#include <string>
using namespace OpenMM;
using std::string;
using std::vector;
RPMDIntegrator::RPMDIntegrator(int numCopies, double temperature, double frictionCoeff, double stepSize) : owner(NULL), numCopies(numCopies) {
setTemperature(temperature);
setFriction(frictionCoeff);
setStepSize(stepSize);
setConstraintTolerance(1e-4);
setRandomNumberSeed((int) time(NULL));
}
void RPMDIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
context = &contextRef;
owner = &contextRef.getOwner();
kernel = context->getPlatform().createKernel(IntegrateRPMDStepKernel::Name(), contextRef);
dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
}
vector<string> RPMDIntegrator::getKernelNames() {
std::vector<std::string> names;
names.push_back(IntegrateRPMDStepKernel::Name());
return names;
}
void RPMDIntegrator::setPositions(int copy, const vector<Vec3>& positions) {
dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).setPositions(copy, positions);
}
void RPMDIntegrator::setVelocities(int copy, const vector<Vec3>& velocities) {
dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).setVelocities(copy, velocities);
}
State RPMDIntegrator::getState(int copy, int types) {
dynamic_cast<const IntegrateRPMDStepKernel&>(kernel.getImpl()).copyToContext(copy, *context);
return context->getOwner().getState(types);
}
void RPMDIntegrator::step(int steps) {
for (int i = 0; i < steps; ++i) {
context->updateContextState();
context->calcForcesAndEnergy(true, false);
dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).execute(*context, *this);
}
}
/* -------------------------------------------------------------------------- *
* 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 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 "ReferenceRpmdKernels.h"
#include "openmm/internal/ContextImpl.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
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);
}
ReferenceIntegrateRPMDStepKernel::~ReferenceIntegrateRPMDStepKernel() {
if (fft != NULL)
fftpack_destroy(fft);
}
void ReferenceIntegrateRPMDStepKernel::initialize(const System& system, const RPMDIntegrator& integrator) {
int numCopies = integrator.getNumCopies();
int numParticles = system.getNumParticles();
positions.resize(numCopies);
velocities.resize(numCopies);
forces.resize(numCopies);
for (int i = 0; i < numCopies; i++) {
positions[i].resize(numParticles);
velocities[i].resize(numParticles);
forces[i].resize(numParticles);
}
fftpack_init_1d(&fft, numParticles);
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
}
void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator) {
int numCopies = positions.size();
int numParticles = positions[0].size();
vector<RealVec>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context);
vector<RealVec>& f = extractForces(context);
if (!hasSetPosition) {
// Initialize the positions from the context.
for (int i = 0; i < numCopies; i++)
for (int j = 0; j < numParticles; j++)
positions[i][j] = pos[j];
hasSetPosition = true;
}
if (!hasSetVelocity) {
// Initialize the velocities from the context.
for (int i = 0; i < numCopies; i++)
for (int j = 0; j < numParticles; j++)
velocities[i][j] = vel[j];
hasSetVelocity = true;
}
// Loop over copies and compute the force on each one.
for (int i = 0; i < numCopies; i++) {
pos = positions[i];
vel = velocities[i];
context.calcForcesAndEnergy(true, false);
forces[i] = f;
}
// Update velocities.
const RealOpenMM dt = integrator.getStepSize();
const System& system = context.getSystem();
for (int i = 0; i < numCopies; i++)
for (int j = 0; j < numParticles; j++)
velocities[i][j] += forces[i][j]*(dt/system.getParticleMass(j));
// Evolve the free ring polymer by transforming to the frequency domain.
vector<t_complex> p(numCopies);
vector<t_complex> q(numCopies);
const RealOpenMM scale = 1.0/sqrt(numCopies);
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/1000;
const RealOpenMM nkT = numCopies*BOLTZ*integrator.getTemperature();
const RealOpenMM twown = 2.0*nkT/hbar;
for (int particle = 0; particle < numParticles; particle++) {
for (int component = 0; component < 3; component++) {
for (int k = 0; k < numCopies; k++) {
q[k] = t_complex(scale*positions[k][particle][component], 0.0);
p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
}
fftpack_exec_1d(fft, FFTPACK_FORWARD, &q[0], &q[0]);
fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
q[0] += p[0]*(dt/system.getParticleMass(particle));
for (int k = 1; k < numCopies; k++) {
const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
const RealOpenMM wt = wk*dt;
const RealOpenMM sinwt2 = sin(wt/2);
const RealOpenMM wm = wk*system.getParticleMass(particle);
const t_complex pprime = p[k] - q[k]*2.0*wm*sinwt2; // Advance momentum from t-dt/2 to t+dt/2
q[k] = p[k]*(2.0*sinwt2/wm) + q[k]*(1.0-4.0*sinwt2*sinwt2); // Advance position from t to t+dt
p[k] = pprime;
}
fftpack_exec_1d(fft, FFTPACK_BACKWARD, &q[0], &q[0]);
fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
for (int k = 0; k < numCopies; k++) {
positions[k][particle][component] = scale*q[k].re;
velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
}
}
}
// Apply the PILE-L thermostat.
const RealOpenMM c1_0 = exp(-dt*integrator.getFriction());
const RealOpenMM c2_0 = sqrt(1.0-c1_0*c1_0);
for (int particle = 0; particle < numParticles; particle++) {
const RealOpenMM localc3 = c2_0*sqrt(system.getParticleMass(particle));
for (int component = 0; component < 3; component++) {
for (int k = 0; k < numCopies; k++)
p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
// Apply a local Langevin thermostat to the centroid mode.
p[0].re = p[0].re*c1_0 + localc3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
// Use critical damping white noise for the remaining modes.
for (int k = 1; k < numCopies/2; k++) {
const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
const RealOpenMM c1 = exp(-2.0*wk*dt);
const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
const RealOpenMM c3 = c2*sqrt(system.getParticleMass(particle)*nkT);
RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
p[k] = p[k]*c1 + t_complex(rand1, rand2);
if (k < numCopies-k)
p[numCopies-k] = p[numCopies-k]*c1 + t_complex(rand1, -rand2);
}
fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
for (int k = 0; k < numCopies; k++)
velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
}
}
}
void ReferenceIntegrateRPMDStepKernel::setPositions(int copy, const vector<Vec3>& pos) {
int numParticles = positions[copy].size();
for (int i = 0; i < numParticles; i++)
positions[copy][i] = pos[i];
hasSetPosition = true;
}
void ReferenceIntegrateRPMDStepKernel::setVelocities(int copy, const vector<Vec3>& vel) {
int numParticles = velocities[copy].size();
for (int i = 0; i < numParticles; i++)
velocities[copy][i] = vel[i];
hasSetVelocity = true;
}
void ReferenceIntegrateRPMDStepKernel::copyToContext(int copy, ContextImpl& context) const {
extractPositions(context) = positions[copy];
extractVelocities(context) = velocities[copy];
}
#ifndef REFERENCE_RPMD_KERNELS_H_
#define REFERENCE_RPMD_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) 2011 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/RpmdKernels.h"
#include "SimTKUtilities/RealVec.h"
#include "SimTKReference/fftpack.h"
namespace OpenMM {
/**
* This kernel is invoked by RPMDIntegrator to take one time step, and to get and
* set the state of system copies.
*/
class ReferenceIntegrateRPMDStepKernel : public IntegrateRPMDStepKernel {
public:
ReferenceIntegrateRPMDStepKernel(std::string name, const Platform& platform) :
IntegrateRPMDStepKernel(name, platform), fft(NULL), hasSetPosition(false), hasSetVelocity(false) {
}
~ReferenceIntegrateRPMDStepKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the RPMDIntegrator this kernel will be used for
*/
void initialize(const System& system, const RPMDIntegrator& integrator);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the RPMDIntegrator this kernel is being used for
*/
void execute(ContextImpl& context, const RPMDIntegrator& integrator);
/**
* Get the positions of all particles in one copy of the system.
*/
void setPositions(int copy, const std::vector<Vec3>& positions);
/**
* Get the velocities of all particles in one copy of the system.
*/
void setVelocities(int copy, const std::vector<Vec3>& velocities);
/**
* Copy positions and velocities for one copy into the context.
*/
void copyToContext(int copy, ContextImpl& context) const;
private:
std::vector<std::vector<RealVec> > positions;
std::vector<std::vector<RealVec> > velocities;
std::vector<std::vector<RealVec> > forces;
fftpack* fft;
bool hasSetPosition, hasSetVelocity;
};
} // namespace OpenMM
#endif /*REFERENCE_RPMD_KERNELS_H_*/
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