#include "gmxapi/md.h"
#include "gmxapi/session.h"
#include "gmxapi/status.h"
-#include "gmxapi/system.h"
#include "system_impl.h"
#include "workflow.h"
{
GMX_ASSERT(impl_, "Constructor requires valid implementation object.");
}
+System::Impl* System::get() const
+{
+ return impl_.get();
+}
System::~System() = default;
//! \endcond
return system;
}
+std::shared_ptr<Workflow> getWork(const System::Impl& system)
+{
+ return system.workflow_;
+}
+
System::Impl::Impl(std::unique_ptr<gmxapi::Workflow> workflow) noexcept :
workflow_(std::move(workflow)), spec_(std::make_shared<MDWorkSpec>())
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*/
std::shared_ptr<Session> launch(const std::shared_ptr<Context>& context);
-private:
//! Description of simulation work.
std::shared_ptr<Workflow> workflow_;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
* 2. A new Session is created using the ContextImpl and the runner
*
* Then, for each module available through getSpec()->getModules(),
- * the session and module are passed to gmxapi::setSessionRestraint().
+ * the session and module are passed to gmxapi::addSessionRestraint().
* 1. A gmx::IRestraintPotential is retrieved from the module.
* 2. A unique, named SessionResources is created for the module and attached to the SessionImpl.
* 1. The module is added as a signaller to the session SignalManager
*/
std::shared_ptr<Session> launch(const std::shared_ptr<Context>& context);
+ [[nodiscard]] Impl* get() const;
+
private:
/*!
* \brief Opaque pointer to implementation.
*/
System fromTprFile(const std::string& filename);
+class Workflow;
+std::shared_ptr<Workflow> getWork(const System::Impl& system);
+
} // end namespace gmxapi
#endif // include guard
)
# Ensure that "make tests" builds all nblib tests so the top-level
# "make check" will work.
-add_dependencies(tests nblib-tests)
+if (BUILD_TESTING)
+ add_dependencies(tests nblib-tests)
-# this allows all nblib tests to be run with "make check-nblib"
-add_custom_target(check-nblib
- COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -R NbLib
- COMMENT "Running nblib tests"
- USES_TERMINAL VERBATIM)
-add_dependencies(check-nblib nblib-tests)
+ # this allows all nblib tests to be run with "make check-nblib"
+ add_custom_target(check-nblib
+ COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -R NbLib
+ COMMENT "Running nblib tests"
+ USES_TERMINAL VERBATIM)
+ add_dependencies(check-nblib nblib-tests)
+endif()
set(NBLIB_MAJOR 0)
set(NBLIB_MINOR 1)
# DIRECTORIES_TO_CHECKSUM - List of directories to hash
# VERSION_CMAKEIN - path to an input template file
# VERSION_OUT - path to the output file
+# HAVE_FUNCTIONING_PYTHON - Whether the python used can do the checking.
#
# The following a possible additional definitions
# PYTHON_EXECUTABLE - Needed to run checking stage of current tree
# version of GROMACS is in use.
set(RELEASE_CHECKSUM_FILE "${PROJECT_SOURCE_DIR}/src/reference_checksum")
if(NOT VERSION_STRING_OF_FORK OR "${VERSION_STRING_OF_FORK}" STREQUAL "")
- if(NOT EXISTS ${RELEASE_CHECKSUM_FILE})
- set(GMX_VERSION_STRING_FULL "${GMX_VERSION_STRING_FULL}-UNCHECKED")
- set(GMX_RELEASE_SOURCE_FILE_CHECKSUM "NoChecksumFile")
- set(GMX_CURRENT_SOURCE_FILE_CHECKSUM "NoChecksumFile")
- message(WARNING "Could not valdiate the GROMACS source due to missing reference checksum file.")
- elseif(NOT PYTHON_EXECUTABLE)
- set(GMX_VERSION_STRING_FULL "${GMX_VERSION_STRING_FULL}-UNCHECKED")
- set(GMX_RELEASE_SOURCE_FILE_CHECKSUM "NoPythonAvailable")
- set(GMX_CURRENT_SOURCE_FILE_CHECKSUM "NoPythonAvailable")
- message(STATUS "Could not calculate checksum of source files without Python")
- else()
+ if(EXISTS ${RELEASE_CHECKSUM_FILE} AND HAVE_FULL_FUNCTIONING_PYTHON)
file(READ ${RELEASE_CHECKSUM_FILE} GMX_RELEASE_SOURCE_FILE_CHECKSUM)
string(STRIP ${GMX_RELEASE_SOURCE_FILE_CHECKSUM} GMX_RELEASE_SOURCE_FILE_CHECKSUM)
set(CHECKSUM_RESULT_FILE "${CMAKE_CURRENT_BINARY_DIR}/computed_checksum")
set(GMX_VERSION_STRING_FULL "${GMX_VERSION_STRING_FULL}-MODIFIED")
message(STATUS "The source code for this GROMACS installation is different from the officially released version.")
endif()
+ elseif(HAVE_FULL_FUNCTIONING_PYTHON)
+ set(GMX_VERSION_STRING_FULL "${GMX_VERSION_STRING_FULL}-UNCHECKED")
+ set(GMX_RELEASE_SOURCE_FILE_CHECKSUM "NoChecksumFile")
+ set(GMX_CURRENT_SOURCE_FILE_CHECKSUM "NoChecksumFile")
+ message(WARNING "Could not valdiate the GROMACS source due to missing reference checksum file.")
+ else()
+ set(GMX_VERSION_STRING_FULL "${GMX_VERSION_STRING_FULL}-UNCHECKED")
+ set(GMX_RELEASE_SOURCE_FILE_CHECKSUM "NoPythonAvailable")
+ set(GMX_CURRENT_SOURCE_FILE_CHECKSUM "NoPythonAvailable")
+ message(STATUS "Could not calculate checksum of source files without functioning Python")
endif()
endif()
else()
endif()
include(gmxCustomCommandUtilities)
-
+include(FindPythonModule)
# The first two are also for use outside this file, encapsulating the details
# of how to use the generated VersionInfo.cmake.
set(VERSION_INFO_CMAKE_FILE ${PROJECT_BINARY_DIR}/VersionInfo.cmake)
# step 1, so they get regenerated only when the static version info
# changes.
+# Check if we have all necessary python modules available
+if (Python3_Interpreter_FOUND)
+ set(HAVE_FULL_FUNCTIONING_PYTHON Python3_Interpreter_FOUND)
+ foreach(module argparse hashlib hmac os stat re) # add further modules if necessary
+ find_python_module(${module} QUIET)
+ string(TOUPPER ${module} module_upper)
+ if(NOT PYTHONMODULE_${module_upper})
+ message(STATUS
+ "Python module ${module} not found - disabling checksum validation")
+ unset(HAVE_FULL_FUNCTIONING_PYTHON)
+ endif()
+ endforeach()
+endif()
+
# Configure information known at this time into a partially filled
# version info file.
set(VERSION_INFO_CMAKEIN_FILE_PARTIAL
OUTPUT ${VERSION_INFO_CMAKE_FILE}
COMMAND ${CMAKE_COMMAND}
-D PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}
+ -D HAVE_FULL_FUNCTIONING_PYTHON=${HAVE_FULL_FUNCTIONING_PYTHON}
-D PROJECT_VERSION=${GMX_VERSION_STRING}
-D PROJECT_SOURCE_DIR=${PROJECT_SOURCE_DIR}
-D DIRECTORIES_TO_CHECKSUM=${DIRECTORIES_TO_CHECKSUM_STRING}
# not been tampered with.
# Note: The RUN_ALWAYS here is to regenerate the hash file only, it does not
# mean that the target is run in all builds
-if (Python3_Interpreter_FOUND)
+if (HAVE_FULL_FUNCTIONING_PYTHON)
# We need the full path to the directories after passing it through
set(FULL_PATH_DIRECTORIES "")
foreach(DIR ${SET_OF_DIRECTORIES_TO_CHECKSUM})
arithmetic, by paying attention to summation order or reorganization of
mathematical expressions. The default configuration uses mixed
precision, but it is easy to turn on double precision by adding the
-option ``-DGMX\_DOUBLE=on`` to ``cmake``. Double precision will be 20 to 100%
+option ``-DGMX_DOUBLE=on`` to ``cmake``. Double precision will be 20 to 100%
slower than mixed precision depending on the architecture you are
running on. Double precision will use somewhat more memory and run
input, energy and full-precision trajectory files will be almost twice
.. _refBernetti2020:
-:sup:`184` Bernetti, M. and Bussi, G.,
+:sup:`184` Bernetti, M. and Bussi G.,
"Pressure control using stochastic cell rescaling", *J. Chem. Phys.*, **153**, 114107 (2020).
+.. raw:: html
+
+ </div>
+
+ <div id="ref-Lidmar2012">
+
+.. _refLidmar2012:
+
+:sup:`185` Lidmar J.,
+"Improving the efficiency of extended ensemble simulations: The accelerated weight histogram method", *Phys. Rev. E*, **85**, 0256708 (2012).
+
+.. raw:: html
+
+ </div>
+
+ <div id="ref-Lindahl2018">
+
+.. _refLindahl2018:
+
+:sup:`186` Lindahl V., Lidmar J. and Hess B.,
+"Riemann metric approach to optimal sampling of multidimensional free-energy landscapes", *Phys. Rev. E*, **98**, 023312 (2018).
+
+.. raw:: html
+
+ </div>
+
+ <div id="ref-LundBorg2021">
+
+.. _refLundborg2021:
+
+:sup:`187` Lundborg M., Lidmar J. and Hess B.,
+"The accelerated weight histogram method for alchemical free energy calculations", *J. Chem. Phys.*, **154**, 204103 (2021).
+
.. raw:: html
</div>
Adaptive biasing with AWH
-------------------------
-The accelerated weight histogram method
+The accelerated weight histogram method :ref:`185 <refLidmar2012>`
:ref:`137 <reflindahl2014accelerated>` calculates the PMF along a reaction coordinate by adding
an adaptively determined biasing potential. AWH flattens free energy
barriers along the reaction coordinate by applying a history-dependent
toward :math:`\rho(\lambda)`.
It is also possible to directly control the :math:`\lambda` state
-of, e.g., alchemical free energy perturbations. In that case there is no harmonic
+of, e.g., alchemical free energy perturbations :ref:`187 <reflundborg2021>`. In that case there is no harmonic
potential and :math:`\lambda` changes in discrete steps along the reaction coordinate
depending on the biased free energy difference between the :math:`\lambda` states.
N.b., it is not yet possible to use AWH in combination with perturbed masses or
is the deviation of the force. The factors :math:`\omega(\lambda|x(t))`,
see :eq:`Eq %s <eqawhomega>`, reweight the samples.
:math:`\eta_{\mu\nu}(\lambda)` is a friction
-tensor \ :ref:`144 <refsivak2012thermodynamic>`. Its matrix elements are inversely proportional to local
+tensor :ref:`186 <reflindahl2018>` and :ref:`144 <refsivak2012thermodynamic>`. Its matrix elements are inversely proportional to local
diffusion coefficients. A measure of sampling (in)efficiency at each
:math:`\lambda` is given by
molecules.
- **AWH code:** currently acts on coordinates provided by the pull
- code.
+ code or the free-energy lambda parameter.
- **free-energy code with harmonic bonds or constraints:** between
single atoms.
.. |NM4| replace:: (kJ mol\ :math:`^{-1}`\ nm\ :math:`^{-4}`
.. |DKJ| replace:: :math:`D` (kJ mol\ :math:`^{-1}`
.. |BETA| replace:: :math:`\beta` (nm\ :math:`^{-1}`
-.. |C23| replace:: :math:`C_{i=2,3}` (kJ mol\ :math:`^{-1}\ nm\ :math:`^{-i}`
+.. |C23| replace:: :math:`C_{i=2,3}` (kJ mol\ :math:`^{-1}`\ nm\ :math:`^{-i}`
.. |BMM| replace:: :math:`b_m`
.. |GE0| replace:: :math:`\geq 0`
.. |KO| replace:: :math:`k`
Fixes where mdrun could behave incorrectly
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Fixed gmxapi MD plugin binding
+""""""""""""""""""""""""""""""
+
+Molecular Dynamics extension code was not properly handled when added to a
+simulation through the gmxapi Python interface.
+This meant that restraint potentials would silently fail to be applied with
+gmxapi versions >= 0.1.
+Updates have been applied internally to gmxapi.
+Third party code should not need to be updated, but developers will
+note an additional "null restraint" in
+https://gitlab.com/gromacs/gromacs/-/tree/master/python_packaging/sample_restraint
+(for illustration and testing purposes).
+
+:issue:`4078`
+
Fixes for ``gmx`` tools
^^^^^^^^^^^^^^^^^^^^^^^
+Fix gmx nmr -viol option
+""""""""""""""""""""""""
+
+The tool would previously fail with a cryptic error.
+Also enforces that this option is exclusive with other analysis modes.
+
+:issue:`4060`
+
+Fixed gmx dipoles -quad option
+""""""""""""""""""""""""""""""
+
+The tool now reports correct values.
+
+:issue:`4080`
+
+Make sure gmx convert-tpr -until works
+""""""""""""""""""""""""""""""""""""""
+
+This got broken during reworking the internals of the tool and didn't
+calculate the number of remaining steps correctly.
+
+:issue:`4056`
+
Fixes that affect portability
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Check that necessary python modules are available
+"""""""""""""""""""""""""""""""""""""""""""""""""
+
+The source code validation could otherwise fail a build with cryptic errors.
+
+:issue:`3985`
+
+Ensure that NB-LIB and gmxapi can be build even without tests enabled
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Could otherwise lead to cryptic build errors.
+
Miscellaneous
^^^^^^^^^^^^^
+Removed performance loss in the mdrun domain decomposition
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+With 16 or more so-called PP MPI ranks, the domain decomposition
+repartitioning could incur large performance overheads due to a sub-optimally
+sized hash table. This has now been fixed.
+
+:issue:`4054`
cmake_minimum_required(VERSION 3.16.3)
# If you are using this repository as a template, you should probably change the
# project name and adopt your own versioning scheme.
-project(sample_restraint VERSION 0.0.8)
+project(sample_restraint VERSION 0.2.1)
find_package(PythonInterp)
# Defines targets for the C++ restraints implemented here. These CMake targets are used by the
# unit tests and by the Python module target defined in ../pythonmodule/CMakeLists.txt
-# Create a shared object library for our restrained ensemble plugin.
-add_library(gmxapi_extension_ensemblepotential STATIC
- ensemblepotential.h
- ensemblepotential.cpp
+# These targets depend both on Gromacs::libgromacs and on Gromacs::gmxapi.
+# These dependencies are likely to evolve.
+# In particular, Gromacs::libgromacs has been deprecated for some time.
+# https://gitlab.com/gromacs/gromacs/-/issues?label_name%5B%5D=API+restructuring
+
+add_library(gmxapi_extension_resources STATIC
sessionresources.cpp)
+set_target_properties(gmxapi_extension_resources PROPERTIES POSITION_INDEPENDENT_CODE ON)
+target_include_directories(gmxapi_extension_resources PUBLIC
+ ${CMAKE_CURRENT_SOURCE_DIR})
+target_link_libraries(gmxapi_extension_resources PUBLIC
+ Gromacs::libgromacs
+ Gromacs::gmxapi)
+
+# Create an archive library for our restrained ensemble plugin.
+add_library(gmxapi_extension_ensemblepotential STATIC
+ ensemblepotential.cpp)
set_target_properties(gmxapi_extension_ensemblepotential PROPERTIES POSITION_INDEPENDENT_CODE ON)
target_include_directories(gmxapi_extension_ensemblepotential PUBLIC
# If building with setuptools, CMake will not be performing the install
set_target_properties(gmxapi_extension_ensemblepotential PROPERTIES BUILD_WITH_INSTALL_RPATH TRUE)
-target_link_libraries(gmxapi_extension_ensemblepotential PUBLIC Gromacs::libgromacs Gromacs::gmxapi)
+target_link_libraries(gmxapi_extension_ensemblepotential PRIVATE
+ gmxapi_extension_resources
+ Gromacs::libgromacs
+ Gromacs::gmxapi)
+
+# Create an archive library for a test plugin.
+add_library(gmxapi_extension_test STATIC
+ nullpotential.cpp
+ )
+set_target_properties(gmxapi_extension_test PROPERTIES POSITION_INDEPENDENT_CODE ON)
+
+target_include_directories(gmxapi_extension_test PUBLIC
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ )
+set_target_properties(gmxapi_extension_test PROPERTIES SKIP_BUILD_RPATH FALSE)
+set_target_properties(gmxapi_extension_test PROPERTIES BUILD_WITH_INSTALL_RPATH TRUE)
+target_link_libraries(gmxapi_extension_test PRIVATE
+ gmxapi_extension_resources)
+target_link_libraries(gmxapi_extension_test PUBLIC
+ Gromacs::libgromacs
+ Gromacs::gmxapi)
* This file currently contains boilerplate that will not be necessary in future gmxapi releases, as
* well as additional code used in implementing the restrained ensemble example workflow.
*
- * A simpler restraint potential would only update the calculate() function. If a callback function is
- * not needed or desired, remove the callback() code from this file and from ensemblepotential.h
+ * A simpler restraint potential would only update the calculate() function. If a callback function
+ * is not needed or desired, remove the callback() code from this file and from ensemblepotential.h
*
* \author M. Eric Irrgang <ericirrgang@gmail.com>
*/
*/
class BlurToGrid
{
- public:
- /*!
- * \brief Construct the blurring functor.
- *
- * \param low The coordinate value of the first grid point.
- * \param gridSpacing Distance between grid points.
- * \param sigma Gaussian parameter for blurring inputs onto the grid.
- */
- BlurToGrid(double low,
- double gridSpacing,
- double sigma) :
- low_{low},
- binWidth_{gridSpacing},
- sigma_{sigma}
- {
- };
-
- /*!
- * \brief Callable for the functor.
- *
- * \param samples A list of values to be blurred onto the grid.
- * \param grid Pointer to the container into which to accumulate a blurred histogram of samples.
- *
- * Example:
- *
- * # Acquire 3 samples to be discretized with blurring.
- * std::vector<double> someData = {3.7, 8.1, 4.2};
- *
- * # Create an empty grid to store magnitudes for points 0.5, 1.0, ..., 10.0.
- * std::vector<double> histogram(20, 0.);
- *
- * # Specify the above grid and a Gaussian parameter of 0.8.
- * auto blur = BlurToGrid(0.5, 0.5, 0.8);
- *
- * # Collect the density grid for the samples.
- * blur(someData, &histogram);
- *
- */
- void operator()(const std::vector<double>& samples,
- std::vector<double>* grid)
+public:
+ /*!
+ * \brief Construct the blurring functor.
+ *
+ * \param low The coordinate value of the first grid point.
+ * \param gridSpacing Distance between grid points.
+ * \param sigma Gaussian parameter for blurring inputs onto the grid.
+ */
+ BlurToGrid(double low, double gridSpacing, double sigma) :
+ low_{ low }, binWidth_{ gridSpacing }, sigma_{ sigma } {};
+
+ /*!
+ * \brief Callable for the functor.
+ *
+ * \param samples A list of values to be blurred onto the grid.
+ * \param grid Pointer to the container into which to accumulate a blurred histogram of samples.
+ *
+ * Example:
+ *
+ * # Acquire 3 samples to be discretized with blurring.
+ * std::vector<double> someData = {3.7, 8.1, 4.2};
+ *
+ * # Create an empty grid to store magnitudes for points 0.5, 1.0, ..., 10.0.
+ * std::vector<double> histogram(20, 0.);
+ *
+ * # Specify the above grid and a Gaussian parameter of 0.8.
+ * auto blur = BlurToGrid(0.5, 0.5, 0.8);
+ *
+ * # Collect the density grid for the samples.
+ * blur(someData, &histogram);
+ *
+ */
+ void operator()(const std::vector<double>& samples, std::vector<double>* grid)
+ {
+ const auto nbins = grid->size();
+ const double& dx{ binWidth_ };
+ const auto num_samples = samples.size();
+
+ const double denominator = 1.0 / (2 * sigma_ * sigma_);
+ const double normalization = 1.0 / (num_samples * sqrt(2.0 * M_PI * sigma_ * sigma_));
+ // We aren't doing any filtering of values too far away to contribute meaningfully, which
+ // is admittedly wasteful for large sigma...
+ for (size_t i = 0; i < nbins; ++i)
{
- const auto nbins = grid->size();
- const double& dx{binWidth_};
- const auto num_samples = samples.size();
-
- const double denominator = 1.0 / (2 * sigma_ * sigma_);
- const double normalization = 1.0 / (num_samples * sqrt(2.0 * M_PI * sigma_ * sigma_));
- // We aren't doing any filtering of values too far away to contribute meaningfully, which
- // is admittedly wasteful for large sigma...
- for (size_t i = 0;i < nbins;++i)
+ double bin_value{ 0 };
+ const double bin_x{ low_ + i * dx };
+ for (const auto distance : samples)
{
- double bin_value{0};
- const double bin_x{low_ + i * dx};
- for (const auto distance : samples)
- {
- const double relative_distance{bin_x - distance};
- const auto numerator = -relative_distance * relative_distance;
- bin_value += normalization * exp(numerator * denominator);
- }
- grid->at(i) = bin_value;
+ const double relative_distance{ bin_x - distance };
+ const auto numerator = -relative_distance * relative_distance;
+ bin_value += normalization * exp(numerator * denominator);
}
- };
+ grid->at(i) = bin_value;
+ }
+ };
- private:
- /// Minimum value of bin zero
- const double low_;
+private:
+ /// Minimum value of bin zero
+ const double low_;
- /// Size of each bin
- const double binWidth_;
+ /// Size of each bin
+ const double binWidth_;
- /// Smoothing factor
- const double sigma_;
+ /// Smoothing factor
+ const double sigma_;
};
-EnsemblePotential::EnsemblePotential(size_t nbins,
- double binWidth,
- double minDist,
- double maxDist,
- PairHist experimental,
- unsigned int nSamples,
- double samplePeriod,
- unsigned int nWindows,
- double k,
- double sigma) :
- nBins_{nbins},
- binWidth_{binWidth},
- minDist_{minDist},
- maxDist_{maxDist},
- histogram_(nbins,
- 0),
- experimental_{std::move(experimental)},
- nSamples_{nSamples},
- currentSample_{0},
- samplePeriod_{samplePeriod},
+EnsemblePotential::EnsemblePotential(size_t nbins,
+ double binWidth,
+ double minDist,
+ double maxDist,
+ PairHist experimental,
+ unsigned int nSamples,
+ double samplePeriod,
+ unsigned int nWindows,
+ double k,
+ double sigma) :
+ nBins_{ nbins },
+ binWidth_{ binWidth },
+ minDist_{ minDist },
+ maxDist_{ maxDist },
+ histogram_(nbins, 0),
+ experimental_{ std::move(experimental) },
+ nSamples_{ nSamples },
+ currentSample_{ 0 },
+ samplePeriod_{ samplePeriod },
// In actuality, we have nsamples at (samplePeriod - dt), but we don't have access to dt.
- nextSampleTime_{samplePeriod},
+ nextSampleTime_{ samplePeriod },
distanceSamples_(nSamples),
- nWindows_{nWindows},
- currentWindow_{0},
- windowStartTime_{0},
- nextWindowUpdateTime_{nSamples * samplePeriod},
+ nWindows_{ nWindows },
+ currentWindow_{ 0 },
+ windowStartTime_{ 0 },
+ nextWindowUpdateTime_{ nSamples * samplePeriod },
windows_{},
- k_{k},
- sigma_{sigma}
-{}
+ k_{ k },
+ sigma_{ sigma }
+{
+}
EnsemblePotential::EnsemblePotential(const input_param_type& params) :
EnsemblePotential(params.nBins,
- params.binWidth,
- params.minDist,
- params.maxDist,
- params.experimental,
- params.nSamples,
- params.samplePeriod,
- params.nWindows,
- params.k,
- params.sigma)
+ params.binWidth,
+ params.minDist,
+ params.maxDist,
+ params.experimental,
+ params.nSamples,
+ params.samplePeriod,
+ params.nWindows,
+ params.k,
+ params.sigma)
{
}
// a parallelized simulation).
//
//
-void EnsemblePotential::callback(gmx::Vector v,
- gmx::Vector v0,
- double t,
- const Resources& resources)
+void EnsemblePotential::callback(gmx::Vector v, gmx::Vector v0, double t, const Resources& resources)
{
- const auto rdiff = v - v0;
- const auto Rsquared = dot(rdiff,
- rdiff);
- const auto R = sqrt(Rsquared);
+ const auto rdiff = v - v0;
+ const auto Rsquared = dot(rdiff, rdiff);
+ const auto R = sqrt(Rsquared);
// Store historical data every sample_period steps
if (t >= nextSampleTime_)
if (t >= nextWindowUpdateTime_)
{
// Get next histogram array, recycling old one if available.
- std::unique_ptr<Matrix<double>> new_window = std::make_unique<Matrix<double>>(1,
- nBins_);
+ std::unique_ptr<Matrix<double>> new_window = std::make_unique<Matrix<double>>(1, nBins_);
std::unique_ptr<Matrix<double>> temp_window;
if (windows_.size() == nWindows_)
{
}
else
{
- auto new_temp_window = std::make_unique<Matrix<double>>(1,
- nBins_);
+ auto new_temp_window = std::make_unique<Matrix<double>>(1, nBins_);
assert(new_temp_window);
temp_window.swap(new_temp_window);
}
// Reduce sampled data for this restraint in this simulation, applying a Gaussian blur to fill a grid.
- auto blur = BlurToGrid(0.0,
- binWidth_,
- sigma_);
+ auto blur = BlurToGrid(0.0, binWidth_, sigma_);
assert(new_window != nullptr);
assert(distanceSamples_.size() == nSamples_);
assert(currentSample_ == nSamples_);
- blur(distanceSamples_,
- new_window->vector());
- // We can just do the blur locally since there aren't many bins. Bundling these operations for
- // all restraints could give us a chance at some parallelism. We should at least use some
- // threading if we can.
+ blur(distanceSamples_, new_window->vector());
+ // We can just do the blur locally since there aren't many bins. Bundling these operations
+ // for all restraints could give us a chance at some parallelism. We should at least use
+ // some threading if we can.
// We request a handle each time before using resources to make error handling easier if there is a failure in
// one of the ensemble member processes and to give more freedom to how resources are managed from step to step.
// Get global reduction (sum) and checkpoint.
assert(temp_window != nullptr);
// Todo: in reduce function, give us a mean instead of a sum.
- ensemble.reduce(*new_window,
- temp_window.get());
+ ensemble.reduce(*new_window, temp_window.get());
// Update window list with smoothed data.
windows_.emplace_back(std::move(new_window));
}
for (const auto& window : windows_)
{
- for (size_t i = 0;i < window->cols();++i)
+ for (size_t i = 0; i < window->cols(); ++i)
{
histogram_.at(i) += (window->vector()->at(i) - experimental_.at(i)) / windows_.size();
}
// with the same number of MD steps in each interval, and the interval will effectively lose digits as the
// simulation progresses, so _update_period should be cleanly representable in binary. When we extract this
// to a facility, we can look for a part of the code with access to the current timestep.
- windowStartTime_ = t;
+ windowStartTime_ = t;
nextWindowUpdateTime_ = nSamples_ * samplePeriod_ + windowStartTime_;
++currentWindow_; // This is currently never used. I'm not sure it will be, either...
// Reset sample times.
nextSampleTime_ = t + samplePeriod_;
};
-
}
// HERE is the function that does the calculation of the restraint force.
//
//
-gmx::PotentialPointData EnsemblePotential::calculate(gmx::Vector v,
- gmx::Vector v0,
- double /* t */)
+gmx::PotentialPointData EnsemblePotential::calculate(gmx::Vector v, gmx::Vector v0, double /* t */)
{
// This is not the vector from v to v0. It is the position of a site
// at v, relative to the origin v0. This is a potentially confusing convention...
- const auto rdiff = v - v0;
- const auto Rsquared = dot(rdiff,
- rdiff);
- const auto R = sqrt(Rsquared);
+ const auto rdiff = v - v0;
+ const auto Rsquared = dot(rdiff, rdiff);
+ const auto R = sqrt(Rsquared);
// Compute output
gmx::PotentialPointData output;
// Energy not needed right now.
-// output.energy = 0;
+ // output.energy = 0;
if (R != 0) // Direction of force is ill-defined when v == v0
{
- double f{0};
+ double f{ 0 };
if (R > maxDist_)
{
}
else
{
- double f_scal{0};
+ double f_scal{ 0 };
- const size_t numBins = histogram_.size();
- double normConst = sqrt(2 * M_PI) * sigma_ * sigma_ * sigma_;
+ const size_t numBins = histogram_.size();
+ double normConst = sqrt(2 * M_PI) * sigma_ * sigma_ * sigma_;
- for (size_t n = 0;n < numBins;n++)
+ for (size_t n = 0; n < numBins; n++)
{
- const double x{n * binWidth_ - R};
- const double argExp{-0.5 * x * x / (sigma_ * sigma_)};
+ const double x{ n * binWidth_ - R };
+ const double argExp{ -0.5 * x * x / (sigma_ * sigma_) };
f_scal += histogram_.at(n) * exp(argExp) * x / normConst;
}
f = -k_ * f_scal;
}
const auto magnitude = f / norm(rdiff);
- output.force = rdiff * static_cast<decltype(rdiff[0])>(magnitude);
+ output.force = rdiff * static_cast<decltype(rdiff[0])>(magnitude);
}
return output;
}
-std::unique_ptr<ensemble_input_param_type>
-makeEnsembleParams(size_t nbins,
- double binWidth,
- double minDist,
- double maxDist,
- const std::vector<double>& experimental,
- unsigned int nSamples,
- double samplePeriod,
- unsigned int nWindows,
- double k,
- double sigma)
+std::unique_ptr<ensemble_input_param_type> makeEnsembleParams(size_t nbins,
+ double binWidth,
+ double minDist,
+ double maxDist,
+ const std::vector<double>& experimental,
+ unsigned int nSamples,
+ double samplePeriod,
+ unsigned int nWindows,
+ double k,
+ double sigma)
{
using std::make_unique;
- auto params = make_unique<ensemble_input_param_type>();
- params->nBins = nbins;
- params->binWidth = binWidth;
- params->minDist = minDist;
- params->maxDist = maxDist;
+ auto params = make_unique<ensemble_input_param_type>();
+ params->nBins = nbins;
+ params->binWidth = binWidth;
+ params->minDist = minDist;
+ params->maxDist = maxDist;
params->experimental = experimental;
- params->nSamples = nSamples;
+ params->nSamples = nSamples;
params->samplePeriod = samplePeriod;
- params->nWindows = nWindows;
- params->k = k;
- params->sigma = sigma;
+ params->nWindows = nWindows;
+ params->k = k;
+ params->sigma = sigma;
return params;
};
-// Important: Explicitly instantiate a definition for the templated class declared in ensemblepotential.h.
-// Failing to do this will cause a linker error.
-template
-class ::plugin::RestraintModule<EnsembleRestraint>;
+// Important: Explicitly instantiate a definition for the templated class declared in
+// ensemblepotential.h. Failing to do this will cause a linker error.
+template class ::plugin::RestraintModule<EnsembleRestraint>;
} // end namespace plugin
struct ensemble_input_param_type
{
/// distance histogram parameters
- size_t nBins{0};
- double binWidth{0.};
+ size_t nBins{ 0 };
+ double binWidth{ 0. };
/// Flat-bottom potential boundaries.
- double minDist{0};
- double maxDist{0};
+ double minDist{ 0 };
+ double maxDist{ 0 };
/// Experimental reference distribution.
PairHist experimental{};
/// Number of samples to store during each window.
- unsigned int nSamples{0};
- double samplePeriod{0};
+ unsigned int nSamples{ 0 };
+ double samplePeriod{ 0 };
/// Number of windows to use for smoothing histogram updates.
- unsigned int nWindows{0};
+ unsigned int nWindows{ 0 };
/// Harmonic force coefficient
- double k{0};
+ double k{ 0 };
/// Smoothing factor: width of Gaussian interpolation for histogram
- double sigma{0};
-
+ double sigma{ 0 };
};
// \todo We should be able to automate a lot of the parameter setting stuff
// The statically compiled fast parameter structure would be generated with a recursive variadic template
// the way a tuple is. ref https://eli.thegreenplace.net/2014/variadic-templates-in-c/
-std::unique_ptr<ensemble_input_param_type>
-makeEnsembleParams(size_t nbins,
- double binWidth,
- double minDist,
- double maxDist,
- const std::vector<double>& experimental,
- unsigned int nSamples,
- double samplePeriod,
- unsigned int nWindows,
- double k,
- double sigma);
+std::unique_ptr<ensemble_input_param_type> makeEnsembleParams(size_t nbins,
+ double binWidth,
+ double minDist,
+ double maxDist,
+ const std::vector<double>& experimental,
+ unsigned int nSamples,
+ double samplePeriod,
+ unsigned int nWindows,
+ double k,
+ double sigma);
/*!
* \brief a residue-pair bias calculator for use in restrained-ensemble simulations.
* Applies a force between two sites according to the difference between an experimentally observed
* site pair distance distribution and the distance distribution observed earlier in the simulation
* trajectory. The sampled distribution is averaged from the previous `nwindows` histograms from all
- * ensemble members. Each window contains a histogram populated with `nsamples` distances recorded at
- * `sample_period` step intervals.
+ * ensemble members. Each window contains a histogram populated with `nsamples` distances recorded
+ * at `sample_period` step intervals.
*
* \internal
- * During a the window_update_period steps of a window, the potential applied is a harmonic function of
- * the difference between the sampled and experimental histograms. At the beginning of the window, this
- * difference is found and a Gaussian blur is applied.
+ * During a the window_update_period steps of a window, the potential applied is a harmonic function
+ * of the difference between the sampled and experimental histograms. At the beginning of the
+ * window, this difference is found and a Gaussian blur is applied.
*/
class EnsemblePotential
{
- public:
- using input_param_type = ensemble_input_param_type;
-
- /* No default constructor. Parameters must be provided. */
- EnsemblePotential() = delete;
-
- /*!
- * \brief Constructor called by the wrapper code to produce a new instance.
- *
- * This constructor is called once per simulation per GROMACS process. Note that until
- * gmxapi 0.0.8 there is only one instance per simulation in a thread-MPI simulation.
- *
- * \param params
- */
- explicit EnsemblePotential(const input_param_type& params);
-
- /*!
- * \brief Deprecated constructor taking a parameter list.
- *
- * \param nbins
- * \param binWidth
- * \param minDist
- * \param maxDist
- * \param experimental
- * \param nSamples
- * \param samplePeriod
- * \param nWindows
- * \param k
- * \param sigma
- */
- EnsemblePotential(size_t nbins,
- double binWidth,
- double minDist,
- double maxDist,
- PairHist experimental,
- unsigned int nSamples,
- double samplePeriod,
- unsigned int nWindows,
- double k,
- double sigma);
-
- /*!
- * \brief Evaluates the pair restraint potential.
- *
- * In parallel simulations, the gmxapi framework does not make guarantees about where or
- * how many times this function is called. It should be simple and stateless; it should not
- * update class member data (see ``ensemblepotential.cpp``. For a more controlled API hook
- * and to manage state in the object, use ``callback()``.
- *
- * \param v position of the site for which force is being calculated.
- * \param v0 reference site (other member of the pair).
- * \param t current simulation time (ps).
- * \return container for force and potential energy data.
- */
- // Implementation note for the future: If dispatching this virtual function is not fast
- // enough, the compiler may be able to better optimize a free
- // function that receives the current restraint as an argument.
- gmx::PotentialPointData calculate(gmx::Vector v,
- gmx::Vector v0,
- gmx_unused double t);
-
- /*!
- * \brief An update function to be called on the simulation master rank/thread periodically by the Restraint framework.
- *
- * Defining this function in a plugin potential is optional. If the function is defined,
- * the restraint framework calls this function (on the first rank only in a parallel simulation) before calling calculate().
- *
- * The callback may use resources provided by the Session in the callback to perform updates
- * to the local or global state of an ensemble of simulations. Future gmxapi releases will
- * include additional optimizations, allowing call-back frequency to be expressed, and more
- * general Session resources, as well as more flexible call signatures.
- */
- void callback(gmx::Vector v,
- gmx::Vector v0,
- double t,
- const Resources& resources);
-
- private:
- /// Width of bins (distance) in histogram
- size_t nBins_;
- double binWidth_;
-
- /// Flat-bottom potential boundaries.
- double minDist_;
- double maxDist_;
- /// Smoothed historic distribution for this restraint. An element of the array of restraints in this simulation.
- // Was `hij` in earlier code.
- PairHist histogram_;
- PairHist experimental_;
-
- /// Number of samples to store during each window.
- unsigned int nSamples_;
- unsigned int currentSample_;
- double samplePeriod_;
- double nextSampleTime_;
- /// Accumulated list of samples during a new window.
- std::vector<double> distanceSamples_;
-
- /// Number of windows to use for smoothing histogram updates.
- size_t nWindows_;
- size_t currentWindow_;
- double windowStartTime_;
- double nextWindowUpdateTime_;
- /// The history of nwindows histograms for this restraint.
- std::vector<std::unique_ptr<plugin::Matrix<double>>> windows_;
-
- /// Harmonic force coefficient
- double k_;
- /// Smoothing factor: width of Gaussian interpolation for histogram
- double sigma_;
+public:
+ using input_param_type = ensemble_input_param_type;
+
+ /* No default constructor. Parameters must be provided. */
+ EnsemblePotential() = delete;
+
+ /*!
+ * \brief Constructor called by the wrapper code to produce a new instance.
+ *
+ * This constructor is called once per simulation per GROMACS process. Note that until
+ * gmxapi 0.0.8 there is only one instance per simulation in a thread-MPI simulation.
+ *
+ * \param params
+ */
+ explicit EnsemblePotential(const input_param_type& params);
+
+ /*!
+ * \brief Deprecated constructor taking a parameter list.
+ *
+ * \param nbins
+ * \param binWidth
+ * \param minDist
+ * \param maxDist
+ * \param experimental
+ * \param nSamples
+ * \param samplePeriod
+ * \param nWindows
+ * \param k
+ * \param sigma
+ */
+ EnsemblePotential(size_t nbins,
+ double binWidth,
+ double minDist,
+ double maxDist,
+ PairHist experimental,
+ unsigned int nSamples,
+ double samplePeriod,
+ unsigned int nWindows,
+ double k,
+ double sigma);
+
+ /*!
+ * \brief Evaluates the pair restraint potential.
+ *
+ * In parallel simulations, the gmxapi framework does not make guarantees about where or
+ * how many times this function is called. It should be simple and stateless; it should not
+ * update class member data (see ``ensemblepotential.cpp``. For a more controlled API hook
+ * and to manage state in the object, use ``callback()``.
+ *
+ * \param v position of the site for which force is being calculated.
+ * \param v0 reference site (other member of the pair).
+ * \param t current simulation time (ps).
+ * \return container for force and potential energy data.
+ */
+ // Implementation note for the future: If dispatching this virtual function is not fast
+ // enough, the compiler may be able to better optimize a free
+ // function that receives the current restraint as an argument.
+ gmx::PotentialPointData calculate(gmx::Vector v, gmx::Vector v0, gmx_unused double t);
+
+ /*!
+ * \brief An update function to be called on the simulation master rank/thread periodically by the Restraint framework.
+ *
+ * Defining this function in a plugin potential is optional. If the function is defined,
+ * the restraint framework calls this function (on the first rank only in a parallel simulation) before calling calculate().
+ *
+ * The callback may use resources provided by the Session in the callback to perform updates
+ * to the local or global state of an ensemble of simulations. Future gmxapi releases will
+ * include additional optimizations, allowing call-back frequency to be expressed, and more
+ * general Session resources, as well as more flexible call signatures.
+ */
+ void callback(gmx::Vector v, gmx::Vector v0, double t, const Resources& resources);
+
+private:
+ /// Width of bins (distance) in histogram
+ size_t nBins_;
+ double binWidth_;
+
+ /// Flat-bottom potential boundaries.
+ double minDist_;
+ double maxDist_;
+ /// Smoothed historic distribution for this restraint. An element of the array of restraints in this simulation.
+ // Was `hij` in earlier code.
+ PairHist histogram_;
+ PairHist experimental_;
+
+ /// Number of samples to store during each window.
+ unsigned int nSamples_;
+ unsigned int currentSample_;
+ double samplePeriod_;
+ double nextSampleTime_;
+ /// Accumulated list of samples during a new window.
+ std::vector<double> distanceSamples_;
+
+ /// Number of windows to use for smoothing histogram updates.
+ size_t nWindows_;
+ size_t currentWindow_;
+ double windowStartTime_;
+ double nextWindowUpdateTime_;
+ /// The history of nwindows histograms for this restraint.
+ std::vector<std::unique_ptr<plugin::Matrix<double>>> windows_;
+
+ /// Harmonic force coefficient
+ double k_;
+ /// Smoothing factor: width of Gaussian interpolation for histogram
+ double sigma_;
};
/*!
*/
class EnsembleRestraint : public ::gmx::IRestraintPotential, private EnsemblePotential
{
- public:
- using EnsemblePotential::input_param_type;
-
- EnsembleRestraint(std::vector<int> sites,
- const input_param_type& params,
- std::shared_ptr<Resources> resources
- ) :
- EnsemblePotential(params),
- sites_{std::move(sites)},
- resources_{std::move(resources)}
- {}
-
- ~EnsembleRestraint() override = default;
-
- /*!
- * \brief Implement required interface of gmx::IRestraintPotential
- *
- * \return list of configured site indices.
- *
- * \todo remove to template header
- * \todo abstraction of site references
- */
- std::vector<int> sites() const override
- {
- return sites_;
- }
-
- /*!
- * \brief Implement the interface gmx::IRestraintPotential
- *
- * Dispatch to calculate() method.
- *
- * \param r1 coordinate of first site
- * \param r2 reference coordinate (second site)
- * \param t simulation time
- * \return calculated force and energy
- *
- * \todo remove to template header.
- */
- gmx::PotentialPointData evaluate(gmx::Vector r1,
- gmx::Vector r2,
- double t) override
- {
- return calculate(r1,
- r2,
- t);
- };
-
- /*!
- * \brief An update function to be called on the simulation master rank/thread periodically by the Restraint framework.
- *
- * Implements optional override of gmx::IRestraintPotential::update
- *
- * This boilerplate will disappear into the Restraint template in an upcoming gmxapi release.
- */
- void update(gmx::Vector v,
- gmx::Vector v0,
- double t) override
- {
- // Todo: use a callback period to mostly bypass this and avoid excessive mutex locking.
- callback(v,
- v0,
- t,
- *resources_);
- };
-
- /*!
- * \brief Implement the binding protocol that allows access to Session resources.
- *
- * The client receives a non-owning pointer to the session and cannot extent the life of the session. In
- * the future we can use a more formal handle mechanism.
- *
- * \param session pointer to the current session
- */
- void bindSession(gmxapi::SessionResources* session) override
- {
- resources_->setSession(session);
- }
-
- void setResources(std::unique_ptr<Resources>&& resources)
- {
- resources_ = std::move(resources);
- }
-
- private:
- std::vector<int> sites_;
-// double callbackPeriod_;
-// double nextCallback_;
- std::shared_ptr<Resources> resources_;
+public:
+ using EnsemblePotential::input_param_type;
+
+ EnsembleRestraint(std::vector<int> sites,
+ const input_param_type& params,
+ std::shared_ptr<Resources> resources) :
+ EnsemblePotential(params), sites_{ std::move(sites) }, resources_{ std::move(resources) }
+ {
+ }
+
+ ~EnsembleRestraint() override = default;
+
+ /*!
+ * \brief Implement required interface of gmx::IRestraintPotential
+ *
+ * \return list of configured site indices.
+ *
+ * \todo remove to template header
+ * \todo abstraction of site references
+ */
+ std::vector<int> sites() const override { return sites_; }
+
+ /*!
+ * \brief Implement the interface gmx::IRestraintPotential
+ *
+ * Dispatch to calculate() method.
+ *
+ * \param r1 coordinate of first site
+ * \param r2 reference coordinate (second site)
+ * \param t simulation time
+ * \return calculated force and energy
+ *
+ * \todo remove to template header.
+ */
+ gmx::PotentialPointData evaluate(gmx::Vector r1, gmx::Vector r2, double t) override
+ {
+ return calculate(r1, r2, t);
+ };
+
+ /*!
+ * \brief An update function to be called on the simulation master rank/thread periodically by the Restraint framework.
+ *
+ * Implements optional override of gmx::IRestraintPotential::update
+ *
+ * This boilerplate will disappear into the Restraint template in an upcoming gmxapi release.
+ */
+ void update(gmx::Vector v, gmx::Vector v0, double t) override
+ {
+ // Todo: use a callback period to mostly bypass this and avoid excessive mutex locking.
+ callback(v, v0, t, *resources_);
+ };
+
+ /*!
+ * \brief Implement the binding protocol that allows access to Session resources.
+ *
+ * The client receives a non-owning pointer to the session and cannot extent the life of the
+ * session. In the future we can use a more formal handle mechanism.
+ *
+ * \param session pointer to the current session
+ */
+ void bindSession(gmxapi::SessionResources* session) override
+ {
+ resources_->setSession(session);
+ }
+
+ void setResources(std::unique_ptr<Resources>&& resources) { resources_ = std::move(resources); }
+
+private:
+ std::vector<int> sites_;
+ // double callbackPeriod_;
+ // double nextCallback_;
+ std::shared_ptr<Resources> resources_;
};
// Important: Just declare the template instantiation here for client code.
// We will explicitly instantiate a definition in the .cpp file where the input_param_type is defined.
-extern template
-class RestraintModule<EnsembleRestraint>;
+extern template class RestraintModule<EnsembleRestraint>;
} // end namespace plugin
-#endif //HARMONICRESTRAINT_ENSEMBLEPOTENTIAL_H
+#endif // HARMONICRESTRAINT_ENSEMBLEPOTENTIAL_H
--- /dev/null
+/*! \file
+ * \brief Code to implement the potential declared in nullpotential.h
+ *
+ * This restraint is implemented in a transitional style. We are moving in the direction of
+ * callback based data flow. There is also a preference amongst the GROMACS developers for
+ * stateless objects or free functions. State can be provided by library-managed facilities
+ * rather than stored in long-lived objects.
+ *
+ * The IRestraintPotential framework is due for revision in conjunction with ongoing evolution of
+ * the gmx::MDModules interactions. Until then, we try to use a forward looking design.
+ *
+ * \author M. Eric Irrgang <ericirrgang@gmail.com>
+ */
+
+#include "nullpotential.h"
+
+#include <utility>
+
+namespace plugin
+{
+
+null_input_param_type makeNullParams(std::vector<int>&& sites)
+{
+ return null_input_param_type{ std::move(sites), 0 };
+}
+
+std::vector<int> sites(const null_input_param_type& input)
+{
+ return input.sites_;
+}
+
+gmx::PotentialPointData evaluate(
+ gmx::Vector /*r1*/,
+ gmx::Vector /*r2*/,
+ double /*t*/,
+ null_input_param_type* input)
+{
+ ++input->count_;
+ return gmx::PotentialPointData();
+}
+
+int count(const null_input_param_type& input)
+{
+ return input.count_;
+}
+
+gmx::PotentialPointData NullRestraint::evaluate(gmx::Vector r1, gmx::Vector r2, double t)
+{
+ return ::plugin::evaluate(r1, r2, t, &data_);
+}
+
+std::vector<int> NullRestraint::sites() const
+{
+ return ::plugin::sites(data_);
+}
+
+NullRestraint::NullRestraint(std::vector<int> sites,
+ const NullRestraint::input_param_type& params,
+ std::shared_ptr<Resources> /*resources*/) :
+ data_{ std::move(sites), params.count_ }
+{
+}
+
+// Important: Explicitly instantiate a definition for the templated class declared in
+// ensemblepotential.h. Failing to do this will cause a linker error.
+template class ::plugin::RestraintModule<NullRestraint>;
+
+} // end namespace plugin
--- /dev/null
+/*! \file
+ * \brief Provide a minimal pluggable restraint potential for illustration and testing purposes.
+ *
+ * \author M. Eric Irrgang <ericirrgang@gmail.com>
+ */
+
+#ifndef GMXAPI_EXTENSION_NULLPOTENTIAL_H
+#define GMXAPI_EXTENSION_NULLPOTENTIAL_H
+
+
+#include <array>
+#include <memory>
+#include <mutex>
+#include <vector>
+
+#include "gmxapi/gromacsfwd.h"
+#include "gmxapi/session.h"
+#include "gmxapi/md/mdmodule.h"
+
+#include "gromacs/restraint/restraintpotential.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/real.h"
+
+#include "sessionresources.h"
+
+
+// Ultimately, the shared object library for the Python module should not export any symbols,
+// but we use a generic namespace here for tidiness.
+namespace plugin
+{
+
+struct null_input_param_type
+{
+ std::vector<int> sites_;
+ int count_;
+};
+
+//! Creation function for NullRestraint input and state data.
+null_input_param_type makeNullParams(std::vector<int>&& sites);
+
+//! Support the IRestraintPotential protocol.
+std::vector<int> sites(const null_input_param_type& input);
+
+//! Implement the NullRestraint force-provider.
+gmx::PotentialPointData evaluate(gmx::Vector r1, gmx::Vector r2, double t, null_input_param_type* input);
+
+//! Get the number of times evaluate has been called.
+int count(const null_input_param_type& input);
+
+class NullRestraint : public ::gmx::IRestraintPotential
+{
+public:
+ using input_param_type = null_input_param_type;
+
+ NullRestraint(std::vector<int> sites,
+ const input_param_type& params,
+ std::shared_ptr<Resources> resources);
+
+ ~NullRestraint() override = default;
+ [[nodiscard]] std::vector<int> sites() const override;
+ gmx::PotentialPointData evaluate(gmx::Vector r1, gmx::Vector r2, double t) override;
+
+ input_param_type data_;
+};
+
+// Important: Just declare the template instantiation here for client code.
+// We will explicitly instantiate a definition in the .cpp file where the input_param_type is defined.
+extern template class RestraintModule<NullRestraint>;
+
+} // end namespace plugin
+
+#endif // GMXAPI_EXTENSION_NULLPOTENTIAL_H
{
// Explicit instantiation.
-template
-class ::plugin::Matrix<double>;
+template class ::plugin::Matrix<double>;
-void ResourcesHandle::reduce(const Matrix<double>& send,
- Matrix<double>* receive) const
+void ResourcesHandle::reduce(const Matrix<double>& send, Matrix<double>* receive) const
{
assert(reduce_);
if (*reduce_)
{
- (*reduce_)(send,
- receive);
+ (*reduce_)(send, receive);
}
else
{
void ResourcesHandle::stop()
{
assert(session_);
- auto signaller = gmxapi::getMdrunnerSignal(session_,
- gmxapi::md::signals::STOP);
+ auto signaller = gmxapi::getMdrunnerSignal(session_, gmxapi::md::signals::STOP);
// Should probably check that the function object has been initialized...
signaller();
if (!bool(reduce_))
{
- throw gmxapi::ProtocolError("reduce operation functor is not set, which should not happen...");
+ throw gmxapi::ProtocolError(
+ "reduce operation functor is not set, which should not happen...");
}
handle.reduce_ = &reduce_;
if (!session_)
{
- throw gmxapi::ProtocolError("Resources::getHandle() must not be called before setSession() has been called.");
+ throw gmxapi::ProtocolError(
+ "Resources::getHandle() must not be called before setSession() has been called.");
}
handle.session_ = session_;
{
if (!session)
{
- throw gmxapi::ProtocolError("Resources::setSession received a null SessionResources pointer.");
+ throw gmxapi::ProtocolError(
+ "Resources::setSession received a null SessionResources pointer.");
}
session_ = session;
}
-} // end namespace myplugin
-
+} // namespace plugin
namespace plugin
{
-// Stop-gap for cross-language data exchange pending SharedData implementation and inclusion of Eigen.
-// Adapted from pybind docs.
+// Stop-gap for cross-language data exchange pending SharedData implementation and inclusion of
+// Eigen. Adapted from pybind docs.
template<class T>
class Matrix
{
- public:
- Matrix(size_t rows,
- size_t cols) :
- rows_(rows),
- cols_(cols),
- data_(rows_ * cols_,
- 0)
- {
- }
+public:
+ Matrix(size_t rows, size_t cols) : rows_(rows), cols_(cols), data_(rows_ * cols_, 0) {}
- explicit Matrix(std::vector<T>&& captured_data) :
- rows_{1},
- cols_{captured_data.size()},
- data_{std::move(captured_data)}
- {
- }
+ explicit Matrix(std::vector<T>&& captured_data) :
+ rows_{ 1 }, cols_{ captured_data.size() }, data_{ std::move(captured_data) }
+ {
+ }
- std::vector<T>* vector()
- { return &data_; }
+ std::vector<T>* vector() { return &data_; }
- T* data()
- { return data_.data(); };
+ T* data() { return data_.data(); };
- size_t rows() const
- { return rows_; }
+ [[nodiscard]] size_t rows() const { return rows_; }
- size_t cols() const
- { return cols_; }
+ [[nodiscard]] size_t cols() const { return cols_; }
- private:
- size_t rows_;
- size_t cols_;
- std::vector<T> data_;
+private:
+ size_t rows_;
+ size_t cols_;
+ std::vector<T> data_;
};
// Defer implicit instantiation to ensemblepotential.cpp
-extern template
-class Matrix<double>;
+extern template class Matrix<double>;
/*!
* \brief An active handle to ensemble resources provided by the Context.
*/
class ResourcesHandle
{
- public:
- /*!
- * \brief Ensemble reduce.
- *
- * \param send Matrices to be summed across the ensemble using Context resources.
- * \param receive destination of reduced data instead of updating internal Matrix.
- */
- void reduce(const Matrix<double>& send,
- Matrix<double>* receive) const;
-
- /*!
- * \brief Issue a stop condition event.
- *
- * Can be called on any or all ranks. Sets a condition that will cause the current simulation to shut down
- * after the current step.
- */
- void stop();
-
- // to be abstracted and hidden...
- const std::function<void(const Matrix<double>&,
- Matrix<double>*)>* reduce_;
-
- gmxapi::SessionResources* session_;
+public:
+ /*!
+ * \brief Ensemble reduce.
+ *
+ * \param send Matrices to be summed across the ensemble using Context resources.
+ * \param receive destination of reduced data instead of updating internal Matrix.
+ */
+ void reduce(const Matrix<double>& send, Matrix<double>* receive) const;
+
+ /*!
+ * \brief Issue a stop condition event.
+ *
+ * Can be called on any or all ranks. Sets a condition that will cause the current simulation to
+ * shut down after the current step.
+ */
+ void stop();
+
+ // to be abstracted and hidden...
+ const std::function<void(const Matrix<double>&, Matrix<double>*)>* reduce_;
+
+ gmxapi::SessionResources* session_;
};
/*!
*/
class Resources
{
- public:
- /*!
- * \brief Create a new resources object.
- *
- * This constructor is called by the framework during Session launch to provide the plugin
- * potential with external resources.
- *
- * \note If getHandle() is going to be used, setSession() must be called first.
- *
- * \param reduce ownership of a function object providing ensemble averaging of a 2D matrix.
- */
- explicit Resources(std::function<void(const Matrix<double>&,
- Matrix<double>*)>&& reduce) :
- reduce_(reduce),
- session_(nullptr)
- {};
-
- /*!
- * \brief Grant the caller an active handle for the currently executing block of code.
- *
- * Objects should not keep resource handles open for longer than a single block of code.
- * calculate() and callback() functions get a handle to the resources for the current time step
- * by calling getHandle().
- *
- * \note setSession() must be called before this function can be used.
- * This clumsy protocol requires other infrastructure before it can be
- * cleaned up for gmxapi 0.1
- *
- * \return resource handle
- *
- * In this release, the only facility provided by the resources is a function object for
- * the ensemble averaging function provided by the Context.
- */
- ResourcesHandle getHandle() const;
-
- /*!
- * \brief Acquires a pointer to a Session managing these resources.
- *
- * \param session non-owning pointer to Session resources.
- */
- void setSession(gmxapi::SessionResources* session);
-
- private:
- //! bound function object to provide ensemble reduce facility.
- std::function<void(const Matrix<double>&,
- Matrix<double>*)> reduce_;
-
- // Raw pointer to the session in which these resources live.
- gmxapi::SessionResources* session_;
+public:
+ /*!
+ * \brief Create a new resources object.
+ *
+ * This constructor is called by the framework during Session launch to provide the plugin
+ * potential with external resources.
+ *
+ * \note If getHandle() is going to be used, setSession() must be called first.
+ *
+ * \param reduce ownership of a function object providing ensemble averaging of a 2D matrix.
+ */
+ explicit Resources(std::function<void(const Matrix<double>&, Matrix<double>*)>&& reduce) :
+ reduce_(reduce), session_(nullptr){};
+
+ /*!
+ * \brief Grant the caller an active handle for the currently executing block of code.
+ *
+ * Objects should not keep resource handles open for longer than a single block of code.
+ * calculate() and callback() functions get a handle to the resources for the current time step
+ * by calling getHandle().
+ *
+ * \note setSession() must be called before this function can be used.
+ * This clumsy protocol requires other infrastructure before it can be
+ * cleaned up for gmxapi 0.1
+ *
+ * \return resource handle
+ *
+ * In this release, the only facility provided by the resources is a function object for
+ * the ensemble averaging function provided by the Context.
+ */
+ [[nodiscard]] ResourcesHandle getHandle() const;
+
+ /*!
+ * \brief Acquires a pointer to a Session managing these resources.
+ *
+ * \param session non-owning pointer to Session resources.
+ */
+ void setSession(gmxapi::SessionResources* session);
+
+private:
+ //! bound function object to provide ensemble reduce facility.
+ std::function<void(const Matrix<double>&, Matrix<double>*)> reduce_;
+
+ // Raw pointer to the session in which these resources live.
+ gmxapi::SessionResources* session_;
};
/*!
*
* \tparam R a class implementing the gmx::IRestraintPotential interface.
*
- * The template type parameter should define a ``input_param_type`` member type.
+ * The template type parameter is a class that defines a ``input_param_type`` member type
+ * and which implements a constructor with the signature
+ * ``R(std::vector<int>, const R::input_param_type&, std::shared_ptr<Resources>)``.
*
* \todo move this to a template header in gmxapi */
template<class R>
class RestraintModule : public gmxapi::MDModule // consider names
{
- public:
- using param_t = typename R::input_param_type;
-
- /*!
- * \brief Construct a named restraint module.
- *
- * Objects of this type are created during Session launch, so this code really doesn't belong
- * here. The Director / Builder for the restraint uses a generic interface to pass standard
- * parameters for pair restraints: a list of sites, a (custom) parameters structure, and
- * resources provided by the Session.
- *
- * \param name
- * \param sites
- * \param params
- * \param resources
- */
- RestraintModule(std::string name,
- std::vector<int> sites,
- const typename R::input_param_type& params,
- std::shared_ptr<Resources> resources) :
- sites_{std::move(sites)},
- params_{params},
- resources_{std::move(resources)},
- name_{std::move(name)}
- {
+public:
+ using param_t = typename R::input_param_type;
+
+ /*!
+ * \brief Construct a named restraint module.
+ *
+ * Objects of this type are created during Session launch, so this code really doesn't belong
+ * here. The Director / Builder for the restraint uses a generic interface to pass standard
+ * parameters for pair restraints: a list of sites, a (custom) parameters structure, and
+ * resources provided by the Session.
+ *
+ * \param name
+ * \param sites
+ * \param params
+ * \param resources
+ */
+ RestraintModule(std::string name,
+ std::vector<int> sites,
+ const typename R::input_param_type& params,
+ std::shared_ptr<Resources> resources) :
+ sites_{ std::move(sites) },
+ params_{ params },
+ resources_{ std::move(resources) },
+ name_{ std::move(name) } {
};
- ~RestraintModule() override = default;
-
- /*!
- * \brief Implement gmxapi::MDModule interface to get module name.
- *
- * name is provided during the building stage.
- * \return
- */
- // \todo make member function const
- const char* name() const override
+ ~RestraintModule() override = default;
+
+ /*!
+ * \brief Implement gmxapi::MDModule interface to get module name.
+ *
+ * name is provided during the building stage.
+ * \return
+ */
+ // \todo make member function const
+ [[nodiscard]] const char* name() const override { return name_.c_str(); }
+
+ /*!
+ * \brief Implement gmxapi::MDModule interface to create a restraint for libgromacs.
+ *
+ * \return (Possibly shared) Ownership of a restraint instance
+ *
+ * Creates the restraint instance if it does not already exist. Only creates one restraint
+ * instance in the lifetime of the RestraintModule.
+ *
+ * Note this interface is not stable but requires other GROMACS and gmxapi infrastructure
+ * to mature before it is clear whether we will be creating a new instance or sharing ownership
+ * of the object. A future version may use a std::unique_ptr.
+ */
+ std::shared_ptr<gmx::IRestraintPotential> getRestraint() override
+ {
+ std::lock_guard<std::mutex> lock(restraintInstantiation_);
+ if (!restraint_)
{
- return name_.c_str();
+ restraint_ = std::make_shared<R>(sites_, params_, resources_);
}
+ return restraint_;
+ }
- /*!
- * \brief Implement gmxapi::MDModule interface to create a restraint for libgromacs.
- *
- * \return (Possibly shared) Ownership of a restraint instance
- *
- * Creates the restraint instance if it does not already exist. Only creates one restraint
- * instance in the lifetime of the RestraintModule.
- *
- * Note this interface is not stable but requires other GROMACS and gmxapi infrastructure
- * to mature before it is clear whether we will be creating a new instance or sharing ownership
- * of the object. A future version may use a std::unique_ptr.
- */
- std::shared_ptr<gmx::IRestraintPotential> getRestraint() override
- {
- std::lock_guard<std::mutex> lock(restraintInstantiation_);
- if (!restraint_)
- {
- restraint_ = std::make_shared<R>(sites_,
- params_,
- resources_);
- }
- return restraint_;
- }
-
- private:
- std::vector<int> sites_;
- param_t params_;
+private:
+ std::vector<int> sites_;
+ param_t params_;
- // Need to figure out if this is copyable or who owns it.
- std::shared_ptr<Resources> resources_;
+ // Need to figure out if this is copyable or who owns it.
+ std::shared_ptr<Resources> resources_;
- const std::string name_;
- std::shared_ptr<R> restraint_{nullptr};
- std::mutex restraintInstantiation_;
+ const std::string name_;
+ std::shared_ptr<R> restraint_{ nullptr };
+ std::mutex restraintInstantiation_;
};
/*!
* \brief Filehandle management helper class.
*
- * Use the RAII pattern to make sure that a (newly) constructed object has an open filehandle and that
- * a the filehandle for a destructed object is closed. Closing a file is not guaranteed to be error-free,
- * so the programmer should explicitly call close() and check for errors (see the C library docs
- * for fclose()).
+ * Use the RAII pattern to make sure that a (newly) constructed object has an open filehandle and
+ * that a the filehandle for a destructed object is closed. Closing a file is not guaranteed to be
+ * error-free, so the programmer should explicitly call close() and check for errors (see the C
+ * library docs for fclose()).
*
* RAIIFile makes sure that fclose() is called exactly once, whether client code issues close()
* or not.
*/
class RAIIFile
{
- public:
-
- /*!
- * \brief Open a file in the chosen access mode.
- *
- * \param filename Name of file to be opened.
- * \param mode access mode as described for the fopen C library call.
- */
- RAIIFile(const char* filename,
- const char* mode) :
- fh_{fopen(filename,
- mode)}
- {}
-
- /*!
- * \brief Open a file for writing.
- *
- * \param filename Name of file to be opened.
- *
- * File is opened in mode "w", which truncates data if the file already exists.
- * For other file access modes, use RAIIFile(const char* filename, const char* mode)
- */
- explicit RAIIFile(const char* filename) :
- RAIIFile(filename,
- "w")
- {}
-
- /*!
- * \brief Explicitly close the associated filehandle.
- *
- * It is good practice to explicitly close the file at a known point in the client code, though
- * it is not strictly necessary. If the filehandle is still open when the RAIIFile object is
- * destroyed, the fclose will be called then.
- *
- * Calling close() additional times on the same RAIIFile object is fine and has no effect in
- * single-threaded code. However, the destructor and close() routines are not thread-safe, so
- * the client code should make sure that close() is not called at the same time by multiple threads.
- * Standard reference-counting constructs, like std::shared_ptr, can be used to make sure the
- * object destructor is called exactly once if it needs to be shared.
- *
- * Refer to documentation on fclose() on checking for and interpreting `errno`.
- */
- void close()
+public:
+ /*!
+ * \brief Open a file in the chosen access mode.
+ *
+ * \param filename Name of file to be opened.
+ * \param mode access mode as described for the fopen C library call.
+ */
+ RAIIFile(const char* filename, const char* mode) : fh_{ fopen(filename, mode) } {}
+
+ /*!
+ * \brief Open a file for writing.
+ *
+ * \param filename Name of file to be opened.
+ *
+ * File is opened in mode "w", which truncates data if the file already exists.
+ * For other file access modes, use RAIIFile(const char* filename, const char* mode)
+ */
+ explicit RAIIFile(const char* filename) : RAIIFile(filename, "w") {}
+
+ /*!
+ * \brief Explicitly close the associated filehandle.
+ *
+ * It is good practice to explicitly close the file at a known point in the client code, though
+ * it is not strictly necessary. If the filehandle is still open when the RAIIFile object is
+ * destroyed, the fclose will be called then.
+ *
+ * Calling close() additional times on the same RAIIFile object is fine and has no effect in
+ * single-threaded code. However, the destructor and close() routines are not thread-safe, so
+ * the client code should make sure that close() is not called at the same time by multiple
+ * threads. Standard reference-counting constructs, like std::shared_ptr, can be used to make
+ * sure the object destructor is called exactly once if it needs to be shared.
+ *
+ * Refer to documentation on fclose() on checking for and interpreting `errno`.
+ */
+ void close()
+ {
+ if (fh_ != nullptr)
{
- if (fh_ != nullptr)
- {
- fclose(fh_);
- }
- fh_ = nullptr;
+ fclose(fh_);
}
-
- /*!
- * \brief RAII destructor.
- *
- * Make sure the filehandle gets closed exactly once.
- */
- ~RAIIFile()
+ fh_ = nullptr;
+ }
+
+ /*!
+ * \brief RAII destructor.
+ *
+ * Make sure the filehandle gets closed exactly once.
+ */
+ ~RAIIFile()
+ {
+ if (fh_ != nullptr)
{
- if (fh_ != nullptr)
- {
- fclose(fh_);
- }
+ fclose(fh_);
}
+ }
- RAIIFile(const RAIIFile&) = delete;
+ RAIIFile(const RAIIFile&) = delete;
- RAIIFile& operator=(const RAIIFile&) = delete;
+ RAIIFile& operator=(const RAIIFile&) = delete;
- RAIIFile(RAIIFile&&) = default;
+ RAIIFile(RAIIFile&&) = default;
- RAIIFile& operator=(RAIIFile&&) = default;
+ RAIIFile& operator=(RAIIFile&&) = default;
- /*!
- * \brief Get the managed filehandle.
- *
- * \return raw pointer to the underlying filehandle.
- */
- FILE* fh() const noexcept
- {
- return fh_;
- }
+ /*!
+ * \brief Get the managed filehandle.
+ *
+ * \return raw pointer to the underlying filehandle.
+ */
+ [[nodiscard]] FILE* fh() const noexcept { return fh_; }
- private:
- /// file handle
- FILE* fh_{nullptr};
+private:
+ /// file handle
+ FILE* fh_{ nullptr };
};
} // end namespace plugin
-#endif //RESTRAINT_SESSIONRESOURCES_H
+#endif // RESTRAINT_SESSIONRESOURCES_H
# The Python module requires the new library we wrote as well as the gmxapi that we found in the top-level
# CMakeLists.txt
-target_link_libraries(gmxapi_extension PRIVATE Gromacs::gmxapi gmxapi_extension_ensemblepotential)
+target_link_libraries(gmxapi_extension PRIVATE
+ Gromacs::gmxapi
+ gmxapi_extension_ensemblepotential
+ gmxapi_extension_test
+ )
-if(GMXAPI_EXTENSION_MASTER_PROJECT)
+if (GMXAPI_EXTENSION_MASTER_PROJECT)
install(TARGETS gmxapi_extension
LIBRARY DESTINATION ${GMXPLUGIN_INSTALL_PATH}
ARCHIVE DESTINATION ${GMXPLUGIN_INSTALL_PATH}
RUNTIME DESTINATION ${GMXPLUGIN_INSTALL_PATH}
)
-endif()
+endif ()
/*! \file
* \brief Provide Python bindings and helper functions for setting up restraint potentials.
*
- * There is currently a lot of boilerplate here that will be generalized and removed in a future version.
- * In the mean time, follow the example for EnsembleRestraint to create the proper helper functions
- * and instantiate the necessary templates.
+ * There is currently a lot of boilerplate here that will be generalized and removed in a future
+ * version. In the mean time, follow the example for EnsembleRestraint to create the proper helper
+ * functions and instantiate the necessary templates.
*
* \author M. Eric Irrgang <ericirrgang@gmail.com>
*/
#include "gmxapi/gmxapi.h"
#include "ensemblepotential.h"
+#include "nullpotential.h"
// Make a convenient alias to save some typing...
namespace py = pybind11;
template<class T>
class PyRestraint : public T, public std::enable_shared_from_this<PyRestraint<T>>
{
- public:
- void bind(py::object object);
-
- using T::name;
-
- /*!
- * \brief
- *
- * T must either derive from gmxapi::MDModule or provide a template specialization for
- * PyRestraint<T>::getModule(). If T derives from gmxapi::MDModule, we can keep a weak pointer
- * to ourself and generate a shared_ptr on request, but std::enable_shared_from_this already
- * does that, so we use it when we can.
- * \return
- */
- std::shared_ptr<gmxapi::MDModule> getModule();
-
- /*!
- * \brief Factory function to get a managed pointer to a new restraint.
- *
- * \tparam ArgsT
- * \param args
- * \return
- */
- template<typename ... ArgsT>
- static std::shared_ptr<PyRestraint<T>> create(ArgsT... args)
- {
- auto newRestraint = std::make_shared<PyRestraint<T>>(args...);
- return newRestraint;
- }
-
- template<typename ... ArgsT>
- explicit PyRestraint(ArgsT... args) :
- T{args...}
- {}
+public:
+ void bind(py::object object);
+
+ using T::name;
+
+ /*!
+ * \brief
+ *
+ * T must either derive from gmxapi::MDModule or provide a template specialization for
+ * PyRestraint<T>::getModule(). If T derives from gmxapi::MDModule, we can keep a weak pointer
+ * to ourself and generate a shared_ptr on request, but std::enable_shared_from_this already
+ * does that, so we use it when we can.
+ * \return
+ */
+ std::shared_ptr<gmxapi::MDModule> getModule();
+
+ /*!
+ * \brief Factory function to get a managed pointer to a new restraint.
+ *
+ * \tparam ArgsT
+ * \param args
+ * \return
+ */
+ template<typename... ArgsT>
+ static std::shared_ptr<PyRestraint<T>> create(ArgsT... args)
+ {
+ auto newRestraint = std::make_shared<PyRestraint<T>>(args...);
+ return newRestraint;
+ }
+ template<typename... ArgsT>
+ explicit PyRestraint(ArgsT... args) : T{ args... }
+ {
+ }
};
/*!
template<class T>
void PyRestraint<T>::bind(py::object object)
{
- PyObject * capsule = object.ptr();
- if (PyCapsule_IsValid(capsule,
- gmxapi::MDHolder::api_name))
+ PyObject* capsule = object.ptr();
+ if (PyCapsule_IsValid(capsule, gmxapi::MDHolder::api_name))
{
- auto holder = static_cast<gmxapi::MDHolder*>(PyCapsule_GetPointer(capsule,
- gmxapi::MDHolder::api_name));
+ auto holder =
+ static_cast<gmxapi::MDHolder*>(PyCapsule_GetPointer(capsule, gmxapi::MDHolder::api_name));
auto workSpec = holder->getSpec();
std::cout << this->name() << " received " << holder->name();
std::cout << " containing spec of size ";
template<class T>
std::shared_ptr<gmxapi::MDModule> PyRestraint<T>::getModule()
{
- auto module = std::make_shared<typename std::enable_if<std::is_base_of<gmxapi::MDModule, T>::value, T>::type>();
+ auto module =
+ std::make_shared<typename std::enable_if<std::is_base_of<gmxapi::MDModule, T>::value, T>::type>();
return module;
}
// New restraints mimicking EnsembleRestraint should specialize getModule() here as above.
//////////////////////////////////////////////////////////////////////////////////////////
+//! Provide the hook used by gmxapi::addSessionRestraint().
+template<>
+std::shared_ptr<gmxapi::MDModule> PyRestraint<plugin::RestraintModule<plugin::NullRestraint>>::getModule()
+{
+ return shared_from_this();
+}
-////////////////////
-// Begin MyRestraint
-/*!
- * \brief No-op restraint class for testing and demonstration.
- */
-class MyRestraint
+class EnsembleRestraintBuilder
{
- public:
- static const char* docstring;
+public:
+ explicit EnsembleRestraintBuilder(const py::object& element)
+ {
+ name_ = py::cast<std::string>(element.attr("name"));
+ assert(!name_.empty());
- static std::string name()
- { return "MyRestraint"; };
-};
+ // It looks like we need some boilerplate exceptions for plugins so we have something to
+ // raise if the element is invalid.
+ assert(py::hasattr(element, "params"));
-template<>
-std::shared_ptr<gmxapi::MDModule> PyRestraint<MyRestraint>::getModule()
-{
- auto module = std::make_shared<gmxapi::MDModule>();
- return module;
-}
+ // Params attribute should be a Python list
+ py::dict parameter_dict = element.attr("params");
+ // \todo Check for the presence of these dictionary keys to avoid hard-to-diagnose error.
+ // Get positional parameters.
+ py::list sites = parameter_dict["sites"];
+ for (auto&& site : sites)
+ {
+ siteIndices_.emplace_back(py::cast<int>(site));
+ }
-// Raw string will have line breaks and indentation as written between the delimiters.
-const char* MyRestraint::docstring =
- R"rawdelimiter(Some sort of custom potential.
-)rawdelimiter";
-// end MyRestraint
-//////////////////
+ auto nbins = py::cast<size_t>(parameter_dict["nbins"]);
+ auto binWidth = py::cast<double>(parameter_dict["binWidth"]);
+ auto minDist = py::cast<double>(parameter_dict["min_dist"]);
+ auto maxDist = pybind11::cast<double>(parameter_dict["max_dist"]);
+ auto experimental = pybind11::cast<std::vector<double>>(parameter_dict["experimental"]);
+ auto nSamples = pybind11::cast<unsigned int>(parameter_dict["nsamples"]);
+ auto samplePeriod = pybind11::cast<double>(parameter_dict["sample_period"]);
+ auto nWindows = pybind11::cast<unsigned int>(parameter_dict["nwindows"]);
+ auto k = pybind11::cast<double>(parameter_dict["k"]);
+ auto sigma = pybind11::cast<double>(parameter_dict["sigma"]);
+
+ auto params = plugin::makeEnsembleParams(
+ nbins, binWidth, minDist, maxDist, experimental, nSamples, samplePeriod, nWindows, k, sigma);
+ params_ = std::move(*params);
+
+ // Note that if we want to grab a reference to the Context or its communicator, we can get
+ // it here through element.workspec._context. We need a more general API solution, but this
+ // code is in the Python bindings code, so we know we are in a Python Context.
+ assert(py::hasattr(element, "workspec"));
+ auto workspec = element.attr("workspec");
+ assert(py::hasattr(workspec, "_context"));
+ context_ = workspec.attr("_context");
+ }
+ /*!
+ * \brief Add node(s) to graph for the work element.
+ *
+ * \param graph networkx.DiGraph object still evolving in gmx.context.
+ *
+ * \todo This may not follow the latest graph building protocol as described.
+ */
+ void build(const py::object& graph)
+ {
+ if (!subscriber_)
+ {
+ return;
+ }
+ else
+ {
+ if (!py::hasattr(subscriber_, "potential"))
+ throw gmxapi::ProtocolError("Invalid subscriber");
+ }
-class EnsembleRestraintBuilder
-{
- public:
- explicit EnsembleRestraintBuilder(py::object element)
+ // Restraints do not currently add any new nodes to the graph, so we
+ // mark this standard 'graph' argument unused.
+ (void)graph;
+
+ // Temporarily subvert things to get quick-and-dirty solution for testing.
+ // Need to capture Python communicator and pybind syntax in closure so EnsembleResources
+ // can just call with matrix arguments.
+
+ // This can be replaced with a subscription and delayed until launch, if necessary.
+ if (!py::hasattr(context_, "ensemble_update"))
{
- name_ = py::cast<std::string>(element.attr("name"));
- assert(!name_.empty());
+ throw gmxapi::ProtocolError("context does not have 'ensemble_update'.");
+ }
+ // make a local copy of the Python object so we can capture it in the lambda
+ auto update = context_.attr("ensemble_update");
+ // Make a callable with standardizeable signature.
+ const std::string name{ name_ };
+ auto functor = [update, name](const plugin::Matrix<double>& send, plugin::Matrix<double>* receive)
+ { update(send, receive, py::str(name)); };
+
+ // To use a reduce function on the Python side, we need to provide it with a Python buffer-like object,
+ // so we will create one here. Note: it looks like the SharedData element will be useful after all.
+ auto resources = std::make_shared<plugin::Resources>(std::move(functor));
+
+ auto potential = PyRestraint<plugin::RestraintModule<plugin::EnsembleRestraint>>::create(
+ name_, siteIndices_, params_, resources);
+
+ auto subscriber = subscriber_;
+ py::list potentialList = subscriber.attr("potential");
+ potentialList.append(potential);
+ };
+
+ /*!
+ * \brief Accept subscription of an MD task.
+ *
+ * \param subscriber Python object with a 'potential' attribute that is a Python list.
+ *
+ * During build, an object is added to the subscriber's self.potential, which is then bound with
+ * system.add_potential(potential) during the subscriber's launch()
+ */
+ void addSubscriber(const py::object& subscriber)
+ {
+ assert(py::hasattr(subscriber, "potential"));
+ subscriber_ = subscriber;
+ };
- // It looks like we need some boilerplate exceptions for plugins so we have something to
- // raise if the element is invalid.
- assert(py::hasattr(element,
- "params"));
+ py::object subscriber_;
+ py::object context_;
+ std::vector<int> siteIndices_;
- // Params attribute should be a Python list
- py::dict parameter_dict = element.attr("params");
- // \todo Check for the presence of these dictionary keys to avoid hard-to-diagnose error.
+ plugin::ensemble_input_param_type params_;
- // Get positional parameters.
- py::list sites = parameter_dict["sites"];
- for (auto&& site : sites)
- {
- siteIndices_.emplace_back(py::cast<int>(site));
- }
-
- auto nbins = py::cast<size_t>(parameter_dict["nbins"]);
- auto binWidth = py::cast<double>(parameter_dict["binWidth"]);
- auto minDist = py::cast<double>(parameter_dict["min_dist"]);
- auto maxDist = pybind11::cast<double>(parameter_dict["max_dist"]);
- auto experimental = pybind11::cast<std::vector<double>>(parameter_dict["experimental"]);
- auto nSamples = pybind11::cast<unsigned int>(parameter_dict["nsamples"]);
- auto samplePeriod = pybind11::cast<double>(parameter_dict["sample_period"]);
- auto nWindows = pybind11::cast<unsigned int>(parameter_dict["nwindows"]);
- auto k = pybind11::cast<double>(parameter_dict["k"]);
- auto sigma = pybind11::cast<double>(parameter_dict["sigma"]);
-
- auto params = plugin::makeEnsembleParams(nbins,
- binWidth,
- minDist,
- maxDist,
- experimental,
- nSamples,
- samplePeriod,
- nWindows,
- k,
- sigma);
- params_ = std::move(*params);
-
- // Note that if we want to grab a reference to the Context or its communicator, we can get it
- // here through element.workspec._context. We need a more general API solution, but this code is
- // in the Python bindings code, so we know we are in a Python Context.
- assert(py::hasattr(element,
- "workspec"));
- auto workspec = element.attr("workspec");
- assert(py::hasattr(workspec,
- "_context"));
- context_ = workspec.attr("_context");
+ std::string name_;
+};
+
+class NullRestraintBuilder
+{
+public:
+ explicit NullRestraintBuilder(const py::object& element)
+ {
+ name_ = py::cast<std::string>(element.attr("name"));
+
+ if (name_.empty())
+ {
+ throw gmxapi::ProtocolError("Restraint must provide a *name*.");
}
+ if (!py::hasattr(element, "params"))
+ {
+ throw gmxapi::ProtocolError("Invalid WorkflowElement. (missing *params*)");
+ }
+
+ // Params attribute should be a Python list
+ py::dict parameter_dict = element.attr("params");
- /*!
- * \brief Add node(s) to graph for the work element.
- *
- * \param graph networkx.DiGraph object still evolving in gmx.context.
- *
- * \todo This may not follow the latest graph building protocol as described.
- */
- void build(py::object graph)
+ // Get positional parameters.
+ py::list sites = parameter_dict["sites"];
+ for (auto&& site : sites)
{
- if (!subscriber_)
- {
- return;
- }
- else
- {
- if (!py::hasattr(subscriber_, "potential")) throw gmxapi::ProtocolError("Invalid subscriber");
- }
+ siteIndices_.emplace_back(py::cast<int>(site));
+ }
- // Restraints do not currently add any new nodes to the graph, so we
- // mark this standard 'graph' argument unused.
- (void) graph;
+ params_ = plugin::makeNullParams(std::vector<int>(siteIndices_));
- // Temporarily subvert things to get quick-and-dirty solution for testing.
- // Need to capture Python communicator and pybind syntax in closure so EnsembleResources
- // can just call with matrix arguments.
+ // Note that if we want to grab a reference to the Context or its communicator, we can get
+ // it here through element.workspec._context. We need a more general API solution, but this
+ // code is in the Python bindings code, so we know we are in a Python Context.
+ assert(py::hasattr(element, "workspec"));
+ auto workspec = element.attr("workspec");
+ assert(py::hasattr(workspec, "_context"));
+ context_ = workspec.attr("_context");
+ }
- // This can be replaced with a subscription and delayed until launch, if necessary.
- if (!py::hasattr(context_, "ensemble_update"))
- {
- throw gmxapi::ProtocolError("context does not have 'ensemble_update'.");
- }
- // make a local copy of the Python object so we can capture it in the lambda
- auto update = context_.attr("ensemble_update");
- // Make a callable with standardizeable signature.
- const std::string name{name_};
- auto functor = [update, name](const plugin::Matrix<double>& send,
- plugin::Matrix<double>* receive) {
- update(send,
- receive,
- py::str(name));
- };
-
- // To use a reduce function on the Python side, we need to provide it with a Python buffer-like object,
- // so we will create one here. Note: it looks like the SharedData element will be useful after all.
- auto resources = std::make_shared<plugin::Resources>(std::move(functor));
-
- auto potential = PyRestraint<plugin::RestraintModule<plugin::EnsembleRestraint>>::create(name_,
- siteIndices_,
- params_,
- resources);
-
- auto subscriber = subscriber_;
- py::list potentialList = subscriber.attr("potential");
- potentialList.append(potential);
-
- };
-
- /*!
- * \brief Accept subscription of an MD task.
- *
- * \param subscriber Python object with a 'potential' attribute that is a Python list.
- *
- * During build, an object is added to the subscriber's self.potential, which is then bound with
- * system.add_potential(potential) during the subscriber's launch()
- */
- void addSubscriber(py::object subscriber)
+ /*!
+ * \brief Add node(s) to graph for the work element.
+ *
+ * \param graph networkx.DiGraph object still evolving in gmx.context.
+ *
+ * \todo This may not follow the latest graph building protocol as described.
+ */
+ void build([[maybe_unused]] const py::object& graph)
+ {
+ if (!subscriber_)
{
- assert(py::hasattr(subscriber,
- "potential"));
- subscriber_ = subscriber;
- };
+ return;
+ }
+ else
+ {
+ if (!py::hasattr(subscriber_, "potential"))
+ throw gmxapi::ProtocolError("Invalid subscriber");
+ }
- py::object subscriber_;
- py::object context_;
- std::vector<int> siteIndices_;
+ // This restraint does not need any session resources.
+ auto null_resources = std::make_shared<plugin::Resources>(
+ [](const plugin::Matrix<double>&, plugin::Matrix<double>*) {});
+
+ auto potential = PyRestraint<plugin::RestraintModule<plugin::NullRestraint>>::create(
+ name_, siteIndices_, params_, std::move(null_resources));
+
+ auto subscriber = subscriber_;
+ py::list potentialList = subscriber.attr("potential");
+ potentialList.append(potential);
+ };
+
+ /*!
+ * \brief Accept subscription of an MD task builder.
+ *
+ * \param subscriber Python object with a 'potential' attribute that is a Python list.
+ *
+ * During build, an object is added to the subscriber's self.potential, which is then bound with
+ * system.add_potential(potential) during the subscriber's launch()
+ */
+ void addSubscriber(const py::object& subscriber)
+ {
+ assert(py::hasattr(subscriber, "potential"));
+ subscriber_ = subscriber;
+ };
+
+ py::object subscriber_;
+ py::object context_;
+ std::vector<int> siteIndices_;
- plugin::ensemble_input_param_type params_;
+ plugin::NullRestraint::input_param_type params_;
- std::string name_;
+ std::string name_;
};
-namespace {
+
+namespace
+{
/*!
* \brief Factory function to create a new builder for use during Session launch.
return builder;
}
+/*!
+ * \brief Factory function to create a new builder for use during Session launch.
+ *
+ * \param element WorkElement provided through Context
+ * \return ownership of new builder object
+ */
+std::unique_ptr<NullRestraintBuilder> createNullRestraintBuilder(const py::object& element)
+{
+ using std::make_unique;
+ auto builder = make_unique<NullRestraintBuilder>(element);
+ return builder;
}
+} // namespace
+
+
////////////////////////////////////////////////////////////////////////////////////////////
// New potentials modeled after EnsembleRestraint should define a Builder class and define a
// factory function here, following the previous two examples. The factory function should be
// The first argument is the name of the module when importing to Python. This should be the same as the name specified
// as the OUTPUT_NAME for the shared object library in the CMakeLists.txt file. The second argument, 'm', can be anything
// but it might as well be short since we use it to refer to aspects of the module we are defining.
-PYBIND11_MODULE(myplugin, m) {
+PYBIND11_MODULE(myplugin, m)
+{
m.doc() = "sample plugin"; // This will be the text of the module's docstring.
// Matrix utility class (temporary). Borrowed from http://pybind11.readthedocs.io/en/master/advanced/pycpp/numpy.html#arrays
- py::class_<plugin::Matrix<double>, std::shared_ptr<plugin::Matrix<double>>>(m,
- "Matrix",
- py::buffer_protocol())
- .def_buffer([](plugin::Matrix<double>& matrix) -> py::buffer_info {
- return py::buffer_info(
- matrix.data(), /* Pointer to buffer */
- sizeof(double), /* Size of one scalar */
- py::format_descriptor<double>::format(), /* Python struct-style format descriptor */
- 2, /* Number of dimensions */
- {matrix.rows(), matrix.cols()}, /* Buffer dimensions */
- {sizeof(double) * matrix.cols(), /* Strides (in bytes) for each index */
- sizeof(double)}
- );
- });
+ py::class_<plugin::Matrix<double>, std::shared_ptr<plugin::Matrix<double>>>(
+ m, "Matrix", py::buffer_protocol())
+ .def_buffer(
+ [](plugin::Matrix<double>& matrix) -> py::buffer_info
+ {
+ return py::buffer_info(
+ matrix.data(), /* Pointer to buffer */
+ sizeof(double), /* Size of one scalar */
+ py::format_descriptor<double>::format(), /* Python struct-style format descriptor */
+ 2, /* Number of dimensions */
+ { matrix.rows(), matrix.cols() }, /* Buffer dimensions */
+ { sizeof(double) * matrix.cols(), /* Strides (in bytes) for each index */
+ sizeof(double) });
+ });
//////////////////////////////////////////////////////////////////////////
// Begin EnsembleRestraint
//
// Define Builder to be returned from ensemble_restraint Python function defined further down.
- pybind11::class_<EnsembleRestraintBuilder> ensembleBuilder(m,
- "EnsembleBuilder");
- ensembleBuilder.def("add_subscriber",
- &EnsembleRestraintBuilder::addSubscriber);
- ensembleBuilder.def("build",
- &EnsembleRestraintBuilder::build);
+ pybind11::class_<EnsembleRestraintBuilder> ensembleBuilder(m, "EnsembleBuilder");
+ ensembleBuilder.def("add_subscriber", &EnsembleRestraintBuilder::addSubscriber);
+ ensembleBuilder.def("build", &EnsembleRestraintBuilder::build);
// Get more concise name for the template instantiation...
using PyEnsemble = PyRestraint<plugin::RestraintModule<plugin::EnsembleRestraint>>;
// Export a Python class for our parameters struct
- py::class_<plugin::EnsembleRestraint::input_param_type> ensembleParams(m, "EnsembleRestraintParams");
- m.def("make_ensemble_params",
- &plugin::makeEnsembleParams);
+ py::class_<plugin::EnsembleRestraint::input_param_type> ensembleParams(
+ m, "EnsembleRestraintParams");
+ m.def("make_ensemble_params", &plugin::makeEnsembleParams);
// API object to build.
py::class_<PyEnsemble, std::shared_ptr<PyEnsemble>> ensemble(m, "EnsembleRestraint");
// EnsembleRestraint can only be created via builder for now.
- ensemble.def("bind",
- &PyEnsemble::bind,
- "Implement binding protocol");
+ ensemble.def("bind", &PyEnsemble::bind, "Implement binding protocol");
/*
* To implement gmxapi_workspec_1_0, the module needs a function that a Context can import that
- * produces a builder that translates workspec elements for session launching. The object returned
- * by our function needs to have an add_subscriber(other_builder) method and a build(graph) method.
- * The build() method returns None or a launcher. A launcher has a signature like launch(rank) and
- * returns None or a runner.
+ * produces a builder that translates workspec elements for session launching. The object
+ * returned by our function needs to have an add_subscriber(other_builder) method and a
+ * build(graph) method. The build() method returns None or a launcher. A launcher has a
+ * signature like launch(rank) and returns None or a runner.
*/
// Generate the name operation that will be used to specify elements of Work in gmxapi workflows.
// WorkElements will then have namespace: "myplugin" and operation: "ensemble_restraint"
m.def("ensemble_restraint",
- [](const py::object element) { return createEnsembleBuilder(element); });
+ [](const py::object& element) { return createEnsembleBuilder(element); });
//
// End EnsembleRestraint
///////////////////////////////////////////////////////////////////////////
-
-
-
+ ///////////////////////////////////////////////////////////////////////////
+ // Begin NullRestraint
+ //
+ // Define the Builder Python interface.
+ py::class_<NullRestraintBuilder> nullRestraintBuilder(m, "NullBuilder");
+ nullRestraintBuilder.def("add_subscriber", &NullRestraintBuilder::addSubscriber);
+ nullRestraintBuilder.def("build", &NullRestraintBuilder::build);
+
+ // Get concise name for the template instantiation.
+ using PyNullR = PyRestraint<plugin::RestraintModule<plugin::NullRestraint>>;
+
+ // Export a Python class for our data structure.
+ py::class_<plugin::NullRestraint::input_param_type> nullParams(m, "NullRestraintParams");
+ m.def("make_null_params", &plugin::makeNullParams);
+
+ // Describe the Python API object that is built.
+ py::class_<PyNullR, std::shared_ptr<PyNullR>> nullRestraint(m, "NullRestraint");
+ nullRestraint.def("bind", &PyNullR::bind, "Implement binding protocol.");
+ // We need a protocol for interacting with pluggable extension code and its data.
+ // See #3038, #3133, #3145, #4079.
+ nullRestraint.def(
+ "count",
+ [](PyNullR* restraint)
+ {
+ return plugin::count(
+ dynamic_cast<plugin::NullRestraint*>(restraint->getRestraint().get())->data_);
+ });
+
+ // Export the factory method that will resolve for
+ // {namespace: "myplugin", operation: "null_restraint"} in a gmxapi WorkElement.
+ m.def("null_restraint",
+ [](const py::object& element) { return createNullRestraintBuilder(element); });
+ //
+ // End NullRestraint
+ ///////////////////////////////////////////////////////////////////////////
}
#include "pybind11/pybind11.h"
#include "pybind11/stl.h"
-#endif //GMXAPI_SAMPLE_RESTRAINT_EXPORT_PLUGIN_H
+#endif // GMXAPI_SAMPLE_RESTRAINT_EXPORT_PLUGIN_H
import logging
import os
+try:
+ import mpi4py.MPI as _MPI
+except (ImportError, ModuleNotFoundError):
+ _MPI = None
+
import gmxapi as gmx
from gmxapi.simulation.context import Context
from gmxapi.simulation.workflow import WorkElement, from_tpr
assert myplugin
+@pytest.mark.usefixtures("cleandir")
+def test_binding_protocol(spc_water_box, mdrun_kwargs):
+ """Test that gmxapi successfully attaches MD plugins."""
+ import myplugin
+
+ if _MPI is not None:
+ _size = _MPI.COMM_WORLD.Get_size()
+ _rank = _MPI.COMM_WORLD.Get_rank()
+ else:
+ _size = 1
+ _rank = 0
+
+ tpr_filename = spc_water_box
+ logger.info("Testing plugin potential with input file {}".format(os.path.abspath(tpr_filename)))
+
+ assert gmx.version.api_is_at_least(0, 2, 1)
+ md = from_tpr([tpr_filename] * _size, append_output=False, **mdrun_kwargs)
+
+ potential = WorkElement(namespace="myplugin",
+ operation="null_restraint",
+ params={'sites': [1, 4]})
+ potential.name = "null restraint"
+ md.add_dependency(potential)
+
+ context = Context(md)
+
+ with context as session:
+ session.run()
+
+ # See also #3038, #3145, #4079
+ assert isinstance(context.potentials, list)
+ assert len(context.potentials) > 0
+ for restraint in context.potentials:
+ if isinstance(restraint, myplugin.NullRestraint):
+ assert restraint.count() > 1
+
+
@pytest.mark.usefixtures("cleandir")
def test_ensemble_potential_nompi(spc_water_box, mdrun_kwargs):
"""Test ensemble potential without an ensemble.
"""
tpr_filename = spc_water_box
- print("Testing plugin potential with input file {}".format(os.path.abspath(tpr_filename)))
+ logger.info("Testing plugin potential with input file {}".format(os.path.abspath(tpr_filename)))
assert gmx.version.api_is_at_least(0, 0, 5)
md = from_tpr([tpr_filename], append_output=False, **mdrun_kwargs)
endif()
if(GMXAPI_MASTER_PROJECT)
- # TODO: Retain compatibility with libgmxapi 0.1 and back down the requirement.
- find_package(gmxapi 0.2.0 REQUIRED
+ find_package(gmxapi 0.2.1 REQUIRED
HINTS "$ENV{GROMACS_DIR}"
)
endif()
# to the `check` target. Normal usage is to first install the Python package,
# then run `pytest` on the `tests` directory. Refer to gmxapi package documentation.
if(NOT GMXAPI_MASTER_PROJECT)
- add_subdirectory(test)
+ if (BUILD_TESTING)
+ add_subdirectory(test)
+ endif()
endif()
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
// Export system container class
py::class_<System, std::shared_ptr<System>> system(m, "MDSystem");
- system.def("launch",
- [](System* system, std::shared_ptr<PyContext> context) {
- auto newSession = system->launch(context->get());
- return newSession;
- },
- "Launch the configured workflow in the provided context.");
+ system.def(
+ "launch",
+ [](System* system, std::shared_ptr<PyContext> context) {
+ auto work = gmxapi::getWork(*system->get());
+ auto newSession = context->launch(*work);
+ return newSession;
+ },
+ "Launch the configured workflow in the provided context.");
// Module-level function
m.def("from_tpr",
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*/
#include "pycontext.h"
+#include "gmxapi/exceptions.h"
#include "gmxapi/gmxapi.h"
#include "gmxapi/md.h"
+#include "gmxapi/session.h"
+#include "gmxapi/status.h"
namespace py = pybind11;
std::shared_ptr<gmxapi::Session> PyContext::launch(const gmxapi::Workflow& work)
{
assert(context_);
- return context_->launch(work);
+ std::shared_ptr<gmxapi::Session> session = nullptr;
+
+ // TODO: gmxapi::Workflow, gmxapi::MDWorkSpec, and gmxapi::MDModule need sensible consolidation.
+ session = gmxapi::launchSession(context_.get(), work);
+ if (!session)
+ {
+ throw gmxapi::ProtocolError("Context::launch() expected to produce non-null session.");
+ }
+
+ for (auto&& module : workNodes_->getModules())
+ {
+ // TODO: This should be the job of the launching code that produces the Session.
+ // Configure the restraints in a restraint manager made available to the session launcher.
+ auto status = gmxapi::addSessionRestraint(session.get(), module);
+ }
+
+ return session;
}
std::shared_ptr<gmxapi::MDWorkSpec> PyContext::getSpec() const
assert(workNodes_);
}
-void PyContext::addMDModule(pybind11::object force_object)
+void PyContext::addMDModule(const pybind11::object& force_object) const
{
// If force_object has a bind method, give it a PyCapsule with a pointer
// to our C++ object.
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
std::shared_ptr<gmxapi::Session> launch(const gmxapi::Workflow& work);
std::shared_ptr<gmxapi::Context> get() const;
- void addMDModule(pybind11::object forceProvider);
+ void addMDModule(const pybind11::object& forceProvider) const;
/*!
* \brief Borrow shared ownership of the System's container of associated modules.
/*
- * define GMX_USE_RDTSCP to use the serializing rdtscp instruction instead of rdtsc.
+ * define GMX_USE_RDTSCP=1 to use the serializing rdtscp instruction instead of rdtsc.
* This is only supported on newer Intel/AMD hardware, but provides better accuracy.
*/
typedef __int64 tMPI_Cycles_t;
static __inline tMPI_Cycles_t tMPI_Cycles_read(void)
{
-#ifdef GMX_USE_RDTSCP
+#if GMX_USE_RDTSCP
unsigned int ui;
return __rdtscp(&ui);
#else
/* Here AWH can be extended to work on other coordinates than pull. */
};
-/*! \brief Checks whether any dimension uses pulling as a coordinate provider.
+/*! \brief Checks whether any dimension uses the given coordinate provider type.
*
* \param[in] awhBiasParams The bias params to check.
- * \returns true if any dimension of the bias is linked to pulling.
+ * \param[in] awhCoordProvider The type of coordinate provider
+ * \returns true if any dimension of the bias is linked to the given provider
*/
-static bool anyDimUsesPull(const AwhBiasParams& awhBiasParams)
+static bool anyDimUsesProvider(const AwhBiasParams& awhBiasParams,
+ const AwhCoordinateProviderType awhCoordProvider)
{
- return std::any_of(
- awhBiasParams.dimParams().begin(), awhBiasParams.dimParams().end(), [](const auto& awhDimParam) {
- return awhDimParam.coordinateProvider() == AwhCoordinateProviderType::Pull;
- });
+ return std::any_of(awhBiasParams.dimParams().begin(),
+ awhBiasParams.dimParams().end(),
+ [&awhCoordProvider](const auto& awhDimParam) {
+ return awhDimParam.coordinateProvider() == awhCoordProvider;
+ });
}
-/*! \brief Checks whether any dimension uses pulling as a coordinate provider.
+/*! \brief Checks whether any dimension uses the given coordinate provider type.
*
* \param[in] awhParams The AWH params to check.
- * \returns true if any dimension of awh is linked to pulling.
+ * \param[in] awhCoordProvider The type of coordinate provider
+ * \returns true if any dimension of awh is linked to the given provider type.
*/
-static bool anyDimUsesPull(const AwhParams& awhParams)
+static bool anyDimUsesProvider(const AwhParams& awhParams, const AwhCoordinateProviderType awhCoordProvider)
{
return std::any_of(awhParams.awhBiasParams().begin(),
awhParams.awhBiasParams().end(),
- [](const auto& awhBiasParam) { return anyDimUsesPull(awhBiasParam); });
+ [&awhCoordProvider](const auto& awhBiasParam) {
+ return anyDimUsesProvider(awhBiasParam, awhCoordProvider);
+ });
}
/*! \brief Checks whether any dimension uses pulling as a coordinate provider.
numFepLambdaStates_(numFepLambdaStates),
fepLambdaState_(fepLambdaState)
{
- if (anyDimUsesPull(awhParams))
+ if (anyDimUsesProvider(awhParams, AwhCoordinateProviderType::Pull))
{
GMX_RELEASE_ASSERT(inputRecord.pull != nullptr, "With AWH we should have pull parameters");
GMX_RELEASE_ASSERT(pull_work != nullptr,
if (fplog != nullptr)
{
please_cite(fplog, "Lindahl2014");
+
+ if (anyDimUsesProvider(awhParams, AwhCoordinateProviderType::FreeEnergyLambda))
+ {
+ please_cite(fplog, "Lundborg2021");
+ }
}
if (haveBiasSharingWithinSimulation(awhParams))
void Awh::registerAwhWithPull(const AwhParams& awhParams, pull_t* pull_work)
{
- GMX_RELEASE_ASSERT(!anyDimUsesPull(awhParams) || pull_work, "Need a valid pull object");
+ GMX_RELEASE_ASSERT(!anyDimUsesProvider(awhParams, AwhCoordinateProviderType::Pull) || pull_work,
+ "Need a valid pull object");
for (const auto& biasParam : awhParams.awhBiasParams())
{
if (dd->constraint_comm)
{
- dc->ga2la->clear();
+ dc->ga2la->clearAndResizeHashTable();
}
}
{
if (dd->vsite_comm)
{
- dd->ga2la_vsite->clear();
+ dd->ga2la_vsite->clearAndResizeHashTable();
}
}
}
}
- //! Clear all the entries in the list.
- void clear()
+ /*! \brief Clear all the entries in the list.
+ *
+ * \param[in] resizeHashTable When true the hash table is optimized based on the current number of entries stored
+ */
+ void clear(const bool resizeHashTable)
{
if (usingDirect_)
{
entry.cell = -1;
}
}
+ else if (resizeHashTable)
+ {
+ data_.hashed.clearAndResizeHashTable();
+ }
else
{
data_.hashed.clear();
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
return nullptr;
}
- /*! \brief Clear all the entries in the list
- *
- * Also optimizes the size of the table based on the current
- * number of elements stored.
- */
+ //! Clear all the entries in the list
void clear()
{
- const int oldNumElements = numElements_;
-
for (hashEntry& entry : table_)
{
entry.key = -1;
}
startIndexForSpaceForListEntry_ = bucket_count();
numElements_ = 0;
+ }
+
+ /*! \brief Clear all the entries in the list and resizes the hash table
+ *
+ * Optimizes the size of the hash table based on the current
+ * number of elements stored.
+ */
+ void clearAndResizeHashTable()
+ {
+ const int oldNumElements = numElements_;
+
+ clear();
/* Resize the hash table when the occupation is far from optimal.
* Do not resize with 0 elements to avoid minimal size when clear()
if (!keepLocalAtomIndices)
{
/* Clear the whole list without the overhead of searching */
- ga2la.clear();
+ ga2la.clear(true);
}
else
{
state_change_natoms(state_local, comm->atomRanges.numHomeAtoms());
/* Rebuild all the indices */
- dd->ga2la->clear();
+ dd->ga2la->clear(false);
ncgindex_set = 0;
wallcycle_sub_stop(wcycle, WallCycleSubCounter::DDGrid);
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
EXPECT_LT(map.bucket_count(), 128);
// Check that the table size is double #elements after clear()
- map.clear();
+ map.clearAndResizeHashTable();
EXPECT_EQ(map.bucket_count(), 128);
// Check that calling clear() a second time does not resize
- map.clear();
+ map.clearAndResizeHashTable();
EXPECT_EQ(map.bucket_count(), 128);
map.insert(2, 'b');
EXPECT_EQ(map.bucket_count(), 128);
// Check that calling clear with 1 elements sizes down
- map.clear();
+ map.clearAndResizeHashTable();
EXPECT_LT(map.bucket_count(), 128);
}
cogIndices_.clear();
cogs_.clear();
numAtomsPerCog_.clear();
- globalToLocalMap_.clear();
+ globalToLocalMap_.clearAndResizeHashTable();
}
} // namespace gmx
153,
2020,
"114107" },
+ { "Lundborg2021",
+ "M. Lundborg, J. Lidmar, B. Hess",
+ "The accelerated weight histogram method for alchemical free energy calculations",
+ "J. Chem. Phys.",
+ 154,
+ 2021,
+ "204103" },
};
#define NSTR static_cast<int>(asize(citedb))
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff