diff --git a/modules/mol/mm/pymod/export_topology.cc b/modules/mol/mm/pymod/export_topology.cc
index c104fbf8c7de6fbe8f85fd6b259193c7bb2f8a58..2469e7610d024cde70f8db6d45a7cbae0567affc 100644
--- a/modules/mol/mm/pymod/export_topology.cc
+++ b/modules/mol/mm/pymod/export_topology.cc
@@ -388,6 +388,7 @@ void export_Topology()
.def("AddExclusion",&ost::mol::mm::Topology::AddExclusion,(arg("idx_one"),arg("idx_two")))
.def("AddPositionConstraint",&ost::mol::mm::Topology::AddPositionConstraint,(arg("idx")))
.def("ResetPositionConstraints",&ost::mol::mm::Topology::ResetPositionConstraints)
+ .def("SetDensity",&ost::mol::mm::Topology::SetDensity,(arg("density"),arg("resolution"),arg("scaling")=1.0))
.def("ResetExclusions",&ost::mol::mm::Topology::ResetExclusions)
.def("AddHarmonicPositionRestraint",&ost::mol::mm::Topology::AddHarmonicPositionRestraint,(arg("idx"),arg("ref_position"),arg("k"),arg("x_scale")=1.0,arg("y_scale")=1.0,arg("z_scale")=1.0))
.def("AddHarmonicDistanceRestraint",&ost::mol::mm::Topology::AddHarmonicDistanceRestraint,(arg("idx_one"),arg("idx_two"),arg("length"),arg("force_constant")))
@@ -438,6 +439,9 @@ void export_Topology()
.def("GetHarmonicDistanceRestraintParameters",&WrapGetHarmonicDistanceRestraintParam,(arg("interaction_idx")))
.def("GetFGMDHBondDonorParameters",&WrapGetFGMDHBondDonorParam,(arg("interaction_idx")))
.def("GetFGMDHBondAcceptorParameters",&WrapGetFGMDHBondAcceptorParam,(arg("interaction_idx")))
+ .def("GetDensity",&ost::mol::mm::Topology::GetDensity)
+ .def("GetDensityResolution",&ost::mol::mm::Topology::GetDensityResolution)
+ .def("GetDensityScaling",&ost::mol::mm::Topology::GetDensityResolution)
//setter functions for interaction parameters
.def("SetHarmonicBondParameters",&ost::mol::mm::Topology::SetHarmonicBondParameters,(arg("interaction_idx"),arg("bond_length"),arg("force_constant")))
diff --git a/modules/mol/mm/src/CMakeLists.txt b/modules/mol/mm/src/CMakeLists.txt
index f5d7ba62c87a26d480e10cd969fa0ac896f14c79..58cd9883c1518e6dacca772825999151e860b1d3 100644
--- a/modules/mol/mm/src/CMakeLists.txt
+++ b/modules/mol/mm/src/CMakeLists.txt
@@ -16,7 +16,11 @@ set(OST_MOL_MM_HEADERS
index.hh
topology.hh
steep.hh
-
+ openmm_plugins/density_plugin/openmmapi/density_force.hh
+ openmm_plugins/density_plugin/openmmapi/density_kernels.hh
+ openmm_plugins/density_plugin/openmmapi/density_force_impl.hh
+ openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.hh
+ openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.hh
)
set(OST_MOL_MM_SOURCES
@@ -34,6 +38,8 @@ set(OST_MOL_MM_SOURCES
topology_creator.cc
topology.cc
steep.cc
+ openmm_plugins/density_plugin/openmmapi/density_force.cc
+ openmm_plugins/density_plugin/openmmapi/density_force_impl.cc
)
# create settings.hh as configurational header, needed to set the plugins path
@@ -42,7 +48,7 @@ set(SETTINGS_HH_FILE "${CMAKE_CURRENT_SOURCE_DIR}/settings.hh")
configure_file(settings.hh.in ${SETTINGS_HH_FILE})
-set(MOL_MM_DEPS ost_mol ost_io)
+set(MOL_MM_DEPS ost_mol ost_io ost_img)
module(NAME mol_mm SOURCES ${OST_MOL_MM_SOURCES}
@@ -62,3 +68,41 @@ copy_if_different("${CMAKE_CURRENT_SOURCE_DIR}" "${STAGE_DIR}/share/openstructur
"CHARMM27.dat" "CHARMM_27_FORCEFIELD"
"ost_mol_mm")
install(FILES "CHARMM27.dat" DESTINATION "share/openstructure/forcefields/")
+
+#The plugins that get loaded by OpenMM have to be built seperately
+
+##########################
+#build the density plugin#
+##########################
+
+#define the sources (headers are already handled above)
+SET(DENSITY_PLUGIN_REFERENCE_SOURCES
+ openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.cc
+ openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.cc
+)
+
+# Set the library name
+SET(OPENMM_DENSITY_REFERENCE_LIBRARY_NAME density_plugin_reference)
+SET(SHARED_TARGET_DENSITY_PLUGIN ${OPENMM_DENSITY_REFERENCE_LIBRARY_NAME})
+
+ADD_LIBRARY(${SHARED_TARGET_DENSITY_PLUGIN} SHARED ${DENSITY_PLUGIN_REFERENCE_SOURCES})
+ADD_DEPENDENCIES(${SHARED_TARGET_DENSITY_PLUGIN} ost_mol_mm ost_img ost_mol_mm)
+
+SET_TARGET_PROPERTIES(${SHARED_TARGET_DENSITY_PLUGIN} PROPERTIES
+ LIBRARY_OUTPUT_DIRECTORY ${STAGE_DIR}/share/openstructure/openmm_plugins/
+ ARCHIVE_OUTPUT_DIRECTORY ${STAGE_DIR}/share/openstructure/openmm_plugins/
+ RUNTIME_OUTPUT_DIRECTORY ${STAGE_DIR}/share/openstructure/openmm_plugins/)
+
+foreach(_DEPENDENCY ${_MOL_MM_DEPS})
+ TARGET_LINK_LIBRARIES(${SHARED_TARGET_DENSITY_PLUGIN} ${_DEPENDENCY})
+endforeach()
+#also link against OpenMM and the mm module
+TARGET_LINK_LIBRARIES(${SHARED_TARGET_DENSITY_PLUGIN} ${OPEN_MM_LIBRARIES})
+TARGET_LINK_LIBRARIES(${SHARED_TARGET_DENSITY_PLUGIN} ost_mol_mm)
+
+file(MAKE_DIRECTORY ${STAGE_DIR}/share/openstructure/openmm_plugins/)
+INSTALL(TARGETS ${SHARED_TARGET_DENSITY_PLUGIN} LIBRARY DESTINATION "${STAGE_DIR}/share/openstructure/openmm_plugins/")
+
+
+
+
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force.cc b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8ff809d9a768edb8efd64acbde96fe2e51bcae31
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force.cc
@@ -0,0 +1,97 @@
+/* -------------------------------------------------------------------------- *
+ * 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) 2014 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 <ost/mol/mm/density_force.hh>
+#include <ost/mol/mm/density_force_impl.hh>
+#include <openmm/OpenMMException.h>
+#include <openmm/internal/AssertionUtilities.h>
+#include <string.h>
+
+namespace ost{ namespace mol{ namespace mm{
+
+OpenMM::ForceImpl* DensityForce::createImpl() const {
+ return new DensityForceImpl(*this);
+}
+
+void DensityForce::updateParametersInContext(OpenMM::Context& context) {
+ dynamic_cast<DensityForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
+}
+
+void DensityForce::addParticle(Real mass){
+ masses_.push_back(mass);
+ in_body_.push_back(false);
+}
+
+Real DensityForce::getParticleMass(int idx) const{
+ if(idx >= this->getNumParticles() || idx < 0){
+ throw OpenMM::OpenMMException("Invalid index provided!");
+ }
+ return masses_[idx];
+}
+
+bool DensityForce::isInBody(int idx) const{
+ if(idx >= this->getNumParticles() || idx < 0){
+ throw OpenMM::OpenMMException("Invalid index provided!");
+ }
+ return in_body_[idx];
+}
+
+void DensityForce::defineBody(const std::vector<int>& indices){
+ //first check, whether any of the indices is invalid or already part of a body
+ int num_particles = this->getNumParticles();
+ for(std::vector<int>::const_iterator i = indices.begin();
+ i != indices.end(); ++i){
+ if(*i >= num_particles || *i < 0){
+ throw OpenMM::OpenMMException("Invalid index provided!");
+ }
+ if(in_body_[*i]){
+ throw OpenMM::OpenMMException("Particle can be part of only one body!");
+ }
+ }
+
+ //let's change the state of the according particles
+ for(std::vector<int>::const_iterator i = indices.begin();
+ i != indices.end(); ++i){
+ in_body_[*i] = true;
+ }
+
+ //and finally define the body
+ bodies_.push_back(indices);
+}
+
+const std::vector<int>& DensityForce::getBody(int idx) const{
+ if(idx >= this->getNumBodies() || idx < 0){
+ throw OpenMM::OpenMMException("Invalid index provided!");
+ }
+ return bodies_[idx];
+}
+
+}}} //ns
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force.hh b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f932ff797b978cd54ebe8be72c03935f587f73af
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force.hh
@@ -0,0 +1,94 @@
+#ifndef OPENMM_DENSITYFORCE_H_
+#define OPENMM_DENSITYFORCE_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) 2014 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 <ost/img/image.hh>
+
+namespace ost { namespace mol{ namespace mm{
+
+
+class DensityForce : public OpenMM::Force {
+public:
+ /**
+ * Create an DensityForce. The density force only makes sense for simulation boxes
+ * with orthogonal base vectors.
+ */
+
+ //since the density is scaled in Angstrom, the resolution is also in Angstrom
+ DensityForce(const ost::img::ImageHandle& d, Real r, Real s): density_(d),
+ resolution_(r),
+ scaling_(s) { }
+
+ void addParticle(Real mass);
+
+ int getNumParticles() const { return masses_.size(); }
+
+ int getNumBodies() const { return bodies_.size(); }
+
+ Real getParticleMass(int idx) const;
+
+ bool isInBody(int idx) const;
+
+ void defineBody(const std::vector<int>& indices);
+
+ const std::vector<int>& getBody(int idx) const;
+
+ void setScaling(Real scaling) { scaling_ = scaling; }
+
+ Real getScaling() const { return scaling_; }
+
+ ost::img::ImageHandle getDensity() const { return density_; }
+
+ Real getResolution() const { return resolution_; }
+
+ void updateParametersInContext(OpenMM::Context& context);
+
+ bool usesPeriodicBoundaryConditions() const { return false; }
+
+protected:
+ OpenMM::ForceImpl* createImpl() const;
+private:
+
+ ost::img::ImageHandle density_;
+ std::vector<Real> masses_;
+ std::vector<bool> in_body_;
+ std::vector<std::vector<int> > bodies_;
+ Real resolution_;
+ Real scaling_;
+};
+
+}}} // ns
+
+#endif /*OPENMM_DENSITYFORCE_H_*/
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force_impl.cc b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force_impl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1d2b859a0eb006fe13c389a9818d2958bc69a801
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force_impl.cc
@@ -0,0 +1,71 @@
+/* -------------------------------------------------------------------------- *
+ * 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) 2014 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 <ost/mol/mm/density_force_impl.hh>
+#include <ost/mol/mm/density_kernels.hh>
+#include <openmm/OpenMMException.h>
+#include <openmm/internal/ContextImpl.h>
+#include <cmath>
+#include <map>
+#include <set>
+#include <sstream>
+
+namespace ost{ namespace mol{ namespace mm{
+
+DensityForceImpl::DensityForceImpl(const DensityForce& owner) : owner(owner) {
+}
+
+DensityForceImpl::~DensityForceImpl() {
+}
+
+void DensityForceImpl::initialize(OpenMM::ContextImpl& context) {
+ kernel = context.getPlatform().createKernel(CalcDensityForceKernel::Name(), context);
+ kernel.getAs<CalcDensityForceKernel>().initialize(context.getSystem(), owner);
+}
+
+double DensityForceImpl::calcForcesAndEnergy(OpenMM::ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+ if ((groups&(1<<owner.getForceGroup())) != 0)
+ return kernel.getAs<CalcDensityForceKernel>().execute(context, includeForces, includeEnergy);
+ return 0.0;
+}
+
+std::vector<std::string> DensityForceImpl::getKernelNames() {
+ std::vector<std::string> names;
+ names.push_back(CalcDensityForceKernel::Name());
+ return names;
+}
+
+void DensityForceImpl::updateParametersInContext(OpenMM::ContextImpl& context) {
+ kernel.getAs<CalcDensityForceKernel>().copyParametersToContext(context, owner);
+}
+
+}}} //ns
+
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force_impl.hh b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force_impl.hh
new file mode 100644
index 0000000000000000000000000000000000000000..850b2c5c23a33b48773bcc29b5001f1ad904353b
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_force_impl.hh
@@ -0,0 +1,74 @@
+#ifndef OPENMM_DENSITYFORCEIMPL_H_
+#define OPENMM_DENSITYFORCEIMPL_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) 2014 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 <ost/mol/mm/density_force.hh>
+#include <openmm/internal/ForceImpl.h>
+#include <openmm/Kernel.h>
+#include <utility>
+#include <set>
+#include <string>
+
+namespace ost{ namespace mol{ namespace mm{
+
+/**
+ * This is the internal implementation of DensityForce.
+ */
+
+class DensityForceImpl : public OpenMM::ForceImpl {
+public:
+ DensityForceImpl(const DensityForce& owner);
+ ~DensityForceImpl();
+ void initialize(OpenMM::ContextImpl& context);
+ const DensityForce& getOwner() const {
+ return owner;
+ }
+ void updateContextState(OpenMM::ContextImpl& context) {
+ // This force field doesn't update the state directly.
+ }
+ double calcForcesAndEnergy(OpenMM::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(OpenMM::ContextImpl& context);
+private:
+ const DensityForce& owner;
+ OpenMM::Kernel kernel;
+};
+
+}}} //ns
+
+#endif /*OPENMM_DENSITYFORCEIMPL_H_*/
+
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_kernels.hh b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_kernels.hh
new file mode 100644
index 0000000000000000000000000000000000000000..78e16cfa795519c82252b7edce168329ea404e97
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/openmmapi/density_kernels.hh
@@ -0,0 +1,80 @@
+#ifndef DENSITY_KERNELS_H_
+#define DENSITY_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) 2014 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 <ost/mol/mm/density_force.hh>
+#include <openmm/KernelImpl.h>
+#include <openmm/Platform.h>
+#include <openmm/System.h>
+#include <string>
+
+namespace ost{ namespace mol{ namespace mm{
+
+/**
+ * This kernel is invoked by DensityForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcDensityForceKernel : public OpenMM::KernelImpl {
+public:
+ static std::string Name() {
+ return "CalcDensityForce";
+ }
+ CalcDensityForceKernel(std::string name, const OpenMM::Platform& platform) : OpenMM::KernelImpl(name, platform) {
+ }
+ /**
+ * Initialize the kernel.
+ *
+ * @param system the System this kernel will be applied to
+ * @param force the DensityForce this kernel will be used for
+ */
+ virtual void initialize(const OpenMM::System& system, const DensityForce& 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(OpenMM::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 DensityForce to copy the parameters from
+ */
+ virtual void copyParametersToContext(OpenMM::ContextImpl& context, const DensityForce& force) = 0;
+};
+
+}}} //ns
+
+#endif /*DENSITY_KERNELS_H_*/
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.cc b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8df88868d90895b7b48b89e436799fd3bde97956
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.cc
@@ -0,0 +1,62 @@
+/* -------------------------------------------------------------------------- *
+ * OpenMMExample *
+ * -------------------------------------------------------------------------- *
+ * 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) 2014 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 <ost/mol/mm/reference_density_kernel_factory.hh>
+#include <ost/mol/mm/reference_density_kernels.hh>
+#include <openmm/reference/ReferencePlatform.h>
+#include <openmm/internal/ContextImpl.h>
+#include <openmm/OpenMMException.h>
+
+using namespace ost::mol::mm;
+
+extern "C" void registerPlatforms() {
+}
+
+extern "C" void registerKernelFactories() {
+ for (int i = 0; i < OpenMM::Platform::getNumPlatforms(); i++) {
+ OpenMM::Platform& platform = OpenMM::Platform::getPlatform(i);
+ if (dynamic_cast<OpenMM::ReferencePlatform*>(&platform) != NULL) {
+ ReferenceDensityKernelFactory* factory = new ReferenceDensityKernelFactory();
+ platform.registerKernelFactory(CalcDensityForceKernel::Name(), factory);
+ }
+ }
+}
+
+extern "C" void registerDensityReferenceKernelFactories() {
+ registerKernelFactories();
+}
+
+OpenMM::KernelImpl* ReferenceDensityKernelFactory::createKernelImpl(std::string name, const OpenMM::Platform& platform, OpenMM::ContextImpl& context) const {
+ OpenMM::ReferencePlatform::PlatformData& data = *static_cast<OpenMM::ReferencePlatform::PlatformData*>(context.getPlatformData());
+ if (name == CalcDensityForceKernel::Name())
+ return new ReferenceCalcDensityForceKernel(name, platform);
+ throw OpenMM::OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
+}
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.hh b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..14dfdf44c23060f65f12c3939815b8bd27664608
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernel_factory.hh
@@ -0,0 +1,53 @@
+#ifndef OPENMM_REFERENCEDensityKERNELFACTORY_H_
+#define OPENMM_REFERENCEDensityKERNELFACTORY_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) 2014 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>
+
+/**
+ * This KernelFactory creates kernels for the reference implementation of the Density plugin.
+ */
+
+namespace ost{ namespace mol{ namespace mm{
+
+class ReferenceDensityKernelFactory : public OpenMM::KernelFactory {
+public:
+ OpenMM::KernelImpl* createKernelImpl(std::string name,
+ const OpenMM::Platform& platform,
+ OpenMM::ContextImpl& context) const;
+};
+
+
+}}} //ns
+
+#endif /*OPENMM_REFERENCEDENSITYKERNELFACTORY_H_*/
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.cc b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.cc
new file mode 100644
index 0000000000000000000000000000000000000000..aeefb13871dd53b16c97501df1917fa31cf51e4b
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.cc
@@ -0,0 +1,342 @@
+/* -------------------------------------------------------------------------- *
+ * 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) 2014 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 <ost/mol/mm/reference_density_kernels.hh>
+#include <ost/mol/mm/density_force.hh>
+#include <openmm/OpenMMException.h>
+#include <openmm/internal/ContextImpl.h>
+#include <openmm/reference/RealVec.h>
+#include <openmm/reference/ReferencePlatform.h>
+
+using namespace OpenMM;
+using namespace std;
+
+namespace{
+
+static vector<RealVec>& extractPositions(ContextImpl& context) {
+ ReferencePlatform::PlatformData* data =
+ reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+ return *((vector<RealVec>*) data->positions);
+}
+
+static vector<RealVec>& extractForces(ContextImpl& context) {
+ ReferencePlatform::PlatformData* data =
+ reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+ return *((vector<RealVec>*) data->forces);
+}
+
+} //anon ns
+
+
+namespace ost{ namespace mol { namespace mm{
+
+ReferenceCalcDensityForceKernel::~ReferenceCalcDensityForceKernel(){
+ if(data_ != NULL) delete [] data_;
+}
+
+void ReferenceCalcDensityForceKernel::initialize(const System& system, const DensityForce& force) {
+
+ if(system.getNumParticles() != force.getNumParticles()){
+ std::string err = "DensityForce must contain same number of";
+ err += " particles as system!";
+ throw OpenMM::OpenMMException(err);
+ }
+
+ if(data_ != NULL) delete [] data_;
+
+ ost::img::ImageHandle density = force.getDensity();
+
+ x_sampling_ = density.GetSpatialSampling()[0] * 0.1;
+ y_sampling_ = density.GetSpatialSampling()[1] * 0.1;
+ z_sampling_ = density.GetSpatialSampling()[2] * 0.1;
+ one_over_x_sampling_ = Real(1.0) / x_sampling_;
+ one_over_y_sampling_ = Real(1.0) / y_sampling_;
+ one_over_z_sampling_ = Real(1.0) / z_sampling_;
+ half_x_sampling_ = Real(0.5) * x_sampling_;
+ half_y_sampling_ = Real(0.5) * y_sampling_;
+ half_z_sampling_ = Real(0.5) * z_sampling_;
+ x_origin_ = density.GetAbsoluteOrigin()[0] * 0.1;
+ y_origin_ = density.GetAbsoluteOrigin()[1] * 0.1;
+ z_origin_ = density.GetAbsoluteOrigin()[2] * 0.1;
+ x_extent_ = density.GetSize()[0];
+ y_extent_ = density.GetSize()[1];
+ z_extent_ = density.GetSize()[2];
+
+ data_ = new Real[x_extent_ * y_extent_ * z_extent_];
+ int current_idx = 0;
+ for(int i = 0; i < x_extent_; ++i){
+ for(int j = 0; j < y_extent_; ++j){
+ for(int k = 0; k < z_extent_; ++k){
+ data_[current_idx] = density.GetReal(ost::img::Point(i,j,k));
+ ++current_idx;
+ }
+ }
+ }
+
+ idx_helper_one_ = y_extent_ * z_extent_;
+ idx_helper_two_ = z_extent_;
+
+ //sigma for gaussian like representations of the atoms
+ resolution_ = force.getResolution() * 0.1;
+ s_ = resolution_ * 0.425;
+ one_s_ = Real(1.0) / s_;
+ square_one_s_ = Real(1.0) / (s_ * s_);
+ one_sqrt_2_pi_s_ = Real(1.0) / std::sqrt(Real(2.0) * Real(M_PI) * s_);
+ padding_dist_ = s_ * Real(2.0);
+ scaling_ = 10000 * force.getScaling();
+
+ //let's get body and particle info...
+
+ num_particles_ = system.getNumParticles();
+ num_bodies_ = force.getNumBodies();
+
+ particle_masses_.resize(num_particles_);
+
+ for(int i = 0; i < num_particles_; ++i){
+ particle_masses_[i] = force.getParticleMass(i);
+ }
+
+ bodies_.clear();
+ //let's set every particle not belonging to a body to its own body
+ for(int i = 0; i < num_particles_; ++i){
+ if(force.isInBody(i)) continue;
+ std::vector<int> body(1,i);
+ bodies_.push_back(body);
+ }
+
+ //let's fill the defined bodies!
+ for(int i = 0; i < num_bodies_; ++i){
+ const std::vector<int>& current_body = force.getBody(i);
+ bodies_.push_back(current_body);
+ }
+
+ //let's update the num_bodies_ !!!
+ num_bodies_ = bodies_.size();
+}
+
+double ReferenceCalcDensityForceKernel::execute(ContextImpl& context,
+ bool includeForces,
+ bool includeEnergy) {
+
+
+
+ vector<RealVec>& pos = extractPositions(context);
+ vector<RealVec>& force = extractForces(context);
+
+ //distance of a density position towards a particle
+ Real x_d;
+ Real y_d;
+ Real z_d;
+
+ //derivative of a full body
+ Real body_dx;
+ Real body_dy;
+ Real body_dz;
+
+ //return value for the full energy, will be updated when iterating over
+ //all bodies
+ Real overall_energy = 0.0;
+
+ //some helper variables
+ Real dd;
+ Real e_term;
+ Real density;
+ Real p_density;
+ std::vector<Real> p_dx;
+ std::vector<Real> p_dy;
+ std::vector<Real> p_dz;
+ std::vector<Real> p_m;
+ Real denominator;
+ Real nominator;
+ Real a, b;
+ Real temp_one, temp_two;
+ Real prefac;
+ Real min_x_pos, max_x_pos;
+ Real min_y_pos, max_y_pos;
+ Real min_z_pos, max_z_pos;
+ Real body_x_origin, body_y_origin, body_z_origin;
+ int body_x_extent, body_y_extent, body_z_extent;
+ int num_voxels;
+
+ for(int body_idx = 0; body_idx < num_bodies_; ++body_idx){
+
+ int body_size = bodies_[body_idx].size();
+
+ //set derivatives for every particle in each direction to zero
+ //and extract the required masses
+ p_dx.assign(body_size,0.0);
+ p_dy.assign(body_size,0.0);
+ p_dz.assign(body_size,0.0);
+ p_m.resize(body_size);
+
+ for(int p_idx = 0; p_idx < body_size; ++p_idx){
+ p_m[p_idx] = particle_masses_[bodies_[body_idx][p_idx]];
+ }
+
+ min_x_pos = std::numeric_limits<Real>::max();
+ max_x_pos = -std::numeric_limits<Real>::max();
+ min_y_pos = std::numeric_limits<Real>::max();
+ max_y_pos = -std::numeric_limits<Real>::max();
+ min_z_pos = std::numeric_limits<Real>::max();
+ max_z_pos = -std::numeric_limits<Real>::max();
+
+ //let's find out how big this thing is...
+ for(int p_idx = 0; p_idx < body_size; ++p_idx){
+ RealVec& p = pos[bodies_[body_idx][p_idx]];
+ min_x_pos = (p[0] < min_x_pos) ? p[0]: min_x_pos;
+ max_x_pos = (p[0] > max_x_pos) ? p[0]: max_x_pos;
+ min_y_pos = (p[1] < min_y_pos) ? p[1]: min_y_pos;
+ max_y_pos = (p[1] > max_y_pos) ? p[1]: max_y_pos;
+ min_z_pos = (p[2] < min_z_pos) ? p[2]: min_z_pos;
+ max_z_pos = (p[2] > max_z_pos) ? p[2]: max_z_pos;
+ }
+
+ //we add some padding
+ min_x_pos -= padding_dist_;
+ max_x_pos += padding_dist_;
+ min_y_pos -= padding_dist_;
+ max_y_pos += padding_dist_;
+ min_z_pos -= padding_dist_;
+ max_z_pos += padding_dist_;
+
+ body_x_extent = std::max(3,
+ static_cast<int>(std::ceil((max_x_pos - min_x_pos)*one_over_x_sampling_)));
+ body_y_extent = std::max(3,
+ static_cast<int>(std::ceil((max_y_pos - min_y_pos)*one_over_y_sampling_)));
+ body_z_extent = std::max(3,
+ static_cast<int>(std::ceil((max_z_pos - min_z_pos)*one_over_z_sampling_)));
+ num_voxels = body_x_extent * body_y_extent * body_z_extent;
+
+ Real overlap = min_x_pos + body_x_extent * x_sampling_ - max_x_pos;
+ body_x_origin = min_x_pos - 0.5 * overlap;
+ overlap = min_y_pos + body_y_extent * y_sampling_ - max_y_pos;
+ body_y_origin = min_y_pos - 0.5 * overlap;
+ overlap = min_z_pos + body_z_extent * z_sampling_ - max_z_pos;
+ body_z_origin = min_z_pos - 0.5 * overlap;
+
+
+ density_values_.resize(num_voxels);
+ body_values_.assign(num_voxels,0.0);
+
+ //let's fill the density values
+ int current_idx = 0;
+ for(int i = 0; i < body_x_extent; ++i){
+ for(int j = 0; j < body_y_extent; ++j){
+ for(int k = 0; k < body_z_extent; ++k){
+ density_values_[current_idx] =
+ this->GetIntpolValue(body_x_origin + i * x_sampling_ + half_x_sampling_,
+ body_y_origin + j * y_sampling_ + half_y_sampling_,
+ body_z_origin + k * z_sampling_ + half_z_sampling_);
+ ++current_idx;
+ }
+ }
+ }
+
+ for(int p_idx = 0; p_idx < body_size; ++p_idx){
+
+ //the actual position of the particle
+ RealVec& p = pos[bodies_[body_idx][p_idx]];
+ current_idx = 0;
+
+ for(int i = 0; i < body_x_extent; ++i){
+
+ x_d = p[0] - (body_x_origin + i * x_sampling_ + half_x_sampling_);
+
+ for(int j = 0; j < body_y_extent; ++j){
+
+ y_d = p[1] - (body_y_origin + j * y_sampling_ + half_y_sampling_);
+
+ for(int k = 0; k < body_z_extent; ++k){
+
+ z_d = p[2] - (body_z_origin + k * z_sampling_ + half_z_sampling_);
+
+ dd = x_d * x_d + y_d * y_d + z_d * z_d;
+
+ e_term = std::exp(Real(-0.5) * dd * square_one_s_);
+ density = density_values_[current_idx];
+ p_dx[p_idx] += (density * (x_d * one_s_) * e_term);
+ p_dy[p_idx] += (density * (y_d * one_s_) * e_term);
+ p_dz[p_idx] += (density * (z_d * one_s_) * e_term);
+ p_density = p_m[p_idx] * one_sqrt_2_pi_s_ * e_term;
+ body_values_[current_idx] += p_density;
+ ++current_idx;
+
+ }
+ }
+ }
+ }
+
+ nominator = 0.0;
+ denominator = 0.0;
+ temp_one = 0.0;
+ temp_two = 0.0;
+
+ for(int i = 0; i < num_voxels; ++i){
+ a = density_values_[i];
+ b = body_values_[i];
+ nominator += (a * b);
+ temp_one += (a * a);
+ temp_two += (b * b);
+ }
+
+ if(temp_one == 0.0 || temp_two == 0.0) continue;
+
+ denominator = std::sqrt(temp_one * temp_two);
+ prefac = scaling_ / denominator * one_sqrt_2_pi_s_;
+
+ body_dx = 0.0;
+ body_dy = 0.0;
+ body_dz = 0.0;
+
+ for(int p_idx = 0; p_idx < body_size; ++p_idx){
+ body_dx -= prefac * p_m[p_idx] * p_dx[p_idx];
+ body_dy -= prefac * p_m[p_idx] * p_dy[p_idx];
+ body_dz -= prefac * p_m[p_idx] * p_dz[p_idx];
+ }
+
+ //finally update the forces accordingly
+ for(int p_idx = 0; p_idx < body_size; ++p_idx){
+ force[bodies_[body_idx][p_idx]][0] += body_dx;
+ force[bodies_[body_idx][p_idx]][1] += body_dy;
+ force[bodies_[body_idx][p_idx]][2] += body_dz;
+ }
+
+ overall_energy -= (scaling_ * nominator/denominator);
+ }
+
+ return overall_energy;
+}
+
+void ReferenceCalcDensityForceKernel::copyParametersToContext(ContextImpl& context, const DensityForce& force) {
+ //ToDo
+}
+
+}}} //ns
diff --git a/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.hh b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7a06088ee572d7d91e606601d2fe9eb0b41725a2
--- /dev/null
+++ b/modules/mol/mm/src/openmm_plugins/density_plugin/platforms/reference/reference_density_kernels.hh
@@ -0,0 +1,164 @@
+#ifndef REFERENCE_DENSITY_KERNELS_H_
+#define REFERENCE_DENSITY_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) 2014 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 <ost/mol/mm/density_kernels.hh>
+#include <openmm/Platform.h>
+#include <vector>
+#include <set>
+#include <limits>
+
+namespace ost{ namespace mol{ namespace mm{
+
+/**
+ * This kernel is invoked by DensityForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcDensityForceKernel : public CalcDensityForceKernel {
+public:
+ ReferenceCalcDensityForceKernel(std::string name, const OpenMM::Platform& platform) : CalcDensityForceKernel(name, platform), data_(NULL) {
+ }
+
+ ~ReferenceCalcDensityForceKernel();
+
+ /**
+ * Initialize the kernel.
+ *
+ * @param system the System this kernel will be applied to
+ * @param force the DensityForce this kernel will be used for
+ */
+ void initialize(const OpenMM::System& system, const DensityForce& 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(OpenMM::ContextImpl& context, bool includeForces, bool includeEnergy);
+ /**
+ * Copy changed parameters over to a context.
+ *
+ * @param context the context to copy parameters to
+ * @param force the DensityForce to copy the parameters from
+ */
+ void copyParametersToContext(OpenMM::ContextImpl& context, const DensityForce& force);
+private:
+
+ inline Real GetValue(int a, int b, int c){
+ return data_[a * idx_helper_one_ + b * idx_helper_two_ + c];
+ }
+
+ inline Real GetIntpolValue(Real x, Real y, Real z){
+
+ //distances to origin
+ Real dx = x - x_origin_;
+ Real dy = y - y_origin_;
+ Real dz = z - z_origin_;
+
+ int x_bin = std::floor(dx * one_over_x_sampling_);
+ int y_bin = std::floor(dy * one_over_y_sampling_);
+ int z_bin = std::floor(dz * one_over_z_sampling_);
+
+ if(x_bin < 0 || x_bin >= (x_extent_-1) ||
+ y_bin < 0 || y_bin >= (y_extent_-1) ||
+ z_bin < 0 || z_bin >= (z_extent_-1)) return 0.0;
+
+ //distances in fraction of bin towards base bin
+ dx = (dx - x_bin * x_sampling_) * one_over_x_sampling_;
+ dy = (dy - y_bin * y_sampling_) * one_over_y_sampling_;
+ dz = (dz - z_bin * z_sampling_) * one_over_z_sampling_;
+
+ intpol_values_[0] = this->GetValue(x_bin,y_bin,z_bin);
+ intpol_values_[1] = this->GetValue(x_bin+1,y_bin,z_bin);
+ intpol_values_[2] = this->GetValue(x_bin,y_bin+1,z_bin);
+ intpol_values_[3] = this->GetValue(x_bin+1,y_bin+1,z_bin);
+ intpol_values_[4] = this->GetValue(x_bin,y_bin,z_bin+1);
+ intpol_values_[5] = this->GetValue(x_bin+1,y_bin,z_bin+1);
+ intpol_values_[6] = this->GetValue(x_bin,y_bin+1,z_bin+1);
+ intpol_values_[7] = this->GetValue(x_bin+1,y_bin+1,z_bin+1);
+
+ //do first bilinear interpolation
+ Real f11 = (1.0 - dx) * intpol_values_[0] + dx * intpol_values_[1];
+ Real f12 = (1.0 - dx) * intpol_values_[2] + dx * intpol_values_[3];
+ Real f21 = (1.0 - dx) * intpol_values_[4] + dx * intpol_values_[5];
+ Real f22 = (1.0 - dx) * intpol_values_[6] + dx * intpol_values_[7];
+
+ Real f1 = (1.0 - dy) * f11 + dy * f12;
+ Real f2 = (1.0 - dy) * f21 + dy * f22;
+
+ return (1.0 - dz) * f1 + dz * f2;
+ }
+
+ //this stuff is required to store and extract the internal density values
+ Real* data_;
+ Real intpol_values_[8];
+ int idx_helper_one_;
+ int idx_helper_two_;
+ Real x_sampling_;
+ Real y_sampling_;
+ Real z_sampling_;
+ Real one_over_x_sampling_;
+ Real one_over_y_sampling_;
+ Real one_over_z_sampling_;
+ Real half_x_sampling_;
+ Real half_y_sampling_;
+ Real half_z_sampling_;
+ Real x_origin_;
+ Real y_origin_;
+ Real z_origin_;
+ int x_extent_;
+ int y_extent_;
+ int z_extent_;
+
+ //this is stuff used in the actual force calculation
+ Real resolution_;
+ Real s_;
+ Real one_s_;
+ Real square_one_s_;
+ Real one_sqrt_2_pi_s_;
+ Real padding_dist_;
+ Real scaling_;
+ std::vector<Real> density_values_;
+ std::vector<Real> body_values_;
+
+ //this is stuff containing information of the simulated system
+ int num_bodies_;
+ int num_particles_;
+ std::vector<float> particle_masses_;
+ std::vector<std::vector<int> > bodies_;
+};
+
+}}} // ns
+
+#endif /*REFERENCE_DENSITY_KERNELS_H_*/
diff --git a/modules/mol/mm/src/settings.hh.in b/modules/mol/mm/src/settings.hh.in
index b3f5511a1d65f35e3f54bd8f4559549f3fdda4ba..50dd9b71abcbd4eb2e9f9229642f40a4843585a1 100644
--- a/modules/mol/mm/src/settings.hh.in
+++ b/modules/mol/mm/src/settings.hh.in
@@ -135,7 +135,7 @@ struct Settings{
use_dispersion_correction(true),
keep_ff_specific_naming(true),
openmm_plugin_directory("@OPEN_MM_PLUGIN_DIR@"),
- custom_plugin_directory("@OPEN_MM_PLUGIN_DIR@")
+ custom_plugin_directory("@STAGE_DIR@/share/openstructure/openmm_plugins")
{ }
diff --git a/modules/mol/mm/src/simulation.cc b/modules/mol/mm/src/simulation.cc
index 830b9bd5b6a1b14f92fdaf99248c95c279158dbb..d4e411becfbd9c25737c024cd5e3ca1793dd0f45 100644
--- a/modules/mol/mm/src/simulation.cc
+++ b/modules/mol/mm/src/simulation.cc
@@ -161,6 +161,7 @@ SimulationPtr Simulation::Load(const String& filename, SettingsPtr settings){
sim_ptr->integrator_ = settings->integrator;
OpenMM::Platform::loadPluginsFromDirectory (settings->openmm_plugin_directory);
+ OpenMM::Platform::loadPluginsFromDirectory (settings->custom_plugin_directory);
OpenMM::Platform* platform;
switch(settings->platform){
@@ -285,6 +286,7 @@ void Simulation::Init(const TopologyPtr top,
//to proceed in time, but first we have to load the proper platform
OpenMM::Platform::loadPluginsFromDirectory (settings->openmm_plugin_directory);
+ OpenMM::Platform::loadPluginsFromDirectory (settings->custom_plugin_directory);
OpenMM::Platform* platform;
std::map<String,String> context_properties;
diff --git a/modules/mol/mm/src/system_creator.cc b/modules/mol/mm/src/system_creator.cc
index 2912a5dcab88240f5d0c254bcc4659aa3ba922af..223f7ea9e540df0b56941f4cd25141622156297f 100644
--- a/modules/mol/mm/src/system_creator.cc
+++ b/modules/mol/mm/src/system_creator.cc
@@ -18,6 +18,7 @@
//------------------------------------------------------------------------------
#include <ost/mol/mm/system_creator.hh>
+#include <ost/mol/mm/density_force.hh>
#include <OpenMM.h>
namespace ost{ namespace mol{ namespace mm{
@@ -437,6 +438,27 @@ SystemPtr SystemCreator::Create(const TopologyPtr top,
settings->barostat_frequency);
sys->addForce(&barostat);
}
+
+ //add density if valid
+ ost::img::ImageHandle density = top->GetDensity();
+ if(density.IsValid()){
+ Real r = top->GetDensityResolution();
+ Real s = top->GetDensityScaling();
+ DensityForce* density_force_ptr = new DensityForce(density,r,s);
+ for(std::vector<Real>::iterator i = masses.begin();
+ i != masses.end(); ++i){
+ density_force_ptr->addParticle(*i);
+ }
+ const std::vector<std::vector<int> >& density_bodies = top->GetDensityBodies();
+ for(std::vector<std::vector<int> >::const_iterator i = density_bodies.begin();
+ i != density_bodies.end(); ++i){
+ density_force_ptr->defineBody(*i);
+ }
+ sys->addForce(density_force_ptr);
+ }
+
+
+
return sys;
}
}}} //ns
diff --git a/modules/mol/mm/src/topology.cc b/modules/mol/mm/src/topology.cc
index 448e25490f6b8486b60f44b1e34bb55750d2fbcc..819d2d500c12de758ec13436909c1645f01d52f3 100644
--- a/modules/mol/mm/src/topology.cc
+++ b/modules/mol/mm/src/topology.cc
@@ -255,6 +255,41 @@ void Topology::AddPositionConstraint(uint index){
position_constraints_.push_back(index);
}
+void Topology::SetDensity(const ost::img::ImageHandle& d, Real r, Real s){
+ density_ = d;
+ density_resolution_ = r;
+ density_scaling_ = s;
+ in_density_body_.assign(num_particles_,false);
+ density_bodies_.clear();
+}
+
+void Topology::DefineDensityBody(const std::vector<int>& indices){
+
+ if(!density_.IsValid()){
+ throw ost::Error("Must set valid density before you can define density bodies");
+ }
+
+ //first check, whether any of the indices is invalid or already part of a body
+ for(std::vector<int>::const_iterator i = indices.begin();
+ i != indices.end(); ++i){
+ if(*i >= static_cast<int>(num_particles_) || *i < 0){
+ throw ost::Error("Invalid index provided!");
+ }
+ if(in_density_body_[*i]){
+ throw ost::Error("Particle can be part of only one body!");
+ }
+ }
+
+ //let's change the state of the according particles
+ for(std::vector<int>::const_iterator i = indices.begin();
+ i != indices.end(); ++i){
+ in_density_body_[*i] = true;
+ }
+
+ //and finally define the body
+ density_bodies_.push_back(indices);
+}
+
uint Topology::AddHarmonicPositionRestraint(uint ind, const geom::Vec3& ref_position, Real k,
Real x_scale, Real y_scale, Real z_scale){
if(ind >= num_particles_){
diff --git a/modules/mol/mm/src/topology.hh b/modules/mol/mm/src/topology.hh
index e578e42001b8c7e56a9eef1a0aa0571b097ac8fc..486accf1eb55c6ba1fe0db5fea040a445d6bb182 100644
--- a/modules/mol/mm/src/topology.hh
+++ b/modules/mol/mm/src/topology.hh
@@ -38,6 +38,7 @@
#include <ost/mol/xcs_editor.hh>
#include <ost/mol/bond_handle.hh>
#include <ost/mol/residue_prop.hh>
+#include <ost/img/image_handle.hh>
#include <time.h>
@@ -127,6 +128,12 @@ public:
void AddPositionConstraint(uint index);
+ void SetDensity(const ost::img::ImageHandle& d, Real r, Real s = 1.0);
+
+ void DefineDensityBody(const std::vector<int>& indices);
+
+ int GetNumDensityBodies() const { return density_bodies_.size(); }
+
void ResetPositionConstraints() { position_constraints_.clear(); }
void ResetExclusions() { exclusions_.clear(); }
@@ -212,6 +219,12 @@ public:
void GetFGMDHBondAcceptorParameters(uint index, uint& index_one, uint& index_two) const;
+ ost::img::ImageHandle GetDensity() { return density_; }
+
+ Real GetDensityResolution() { return density_resolution_; }
+
+ Real GetDensityScaling() { return density_scaling_; }
+
void SetHarmonicBondParameters(uint index, const Real bond_length, const Real force_constant);
void SetHarmonicAngleParameters(uint index, const Real angle, const Real force_constant);
@@ -272,6 +285,9 @@ public:
const std::vector<Index<2> >& GetFGMDHBondAcceptors() const { return fgmd_hbond_acceptors_; }
+ const std::vector<std::vector<int> >& GetDensityBodies() const { return density_bodies_; }
+
+
std::vector<Real> GetSigmas() const { return sigmas_; }
std::vector<Real> GetEpsilons() const { return epsilons_; }
@@ -796,6 +812,12 @@ private:
std::set<Index<2> > added_distance_constraints_;
std::set<Index<2> > added_exclusions_;
+ //density map
+ ost::img::ImageHandle density_;
+ Real density_resolution_;
+ Real density_scaling_;
+ std::vector<std::vector<int> > density_bodies_;
+ std::vector<bool> in_density_body_;
};