From 3228ee4ca9f4023d8c50864bbab1aa7805247e72 Mon Sep 17 00:00:00 2001 From: RubenChM Date: Thu, 23 Apr 2026 09:16:10 +0200 Subject: [PATCH 1/6] Apply fixes from #49 --- include/aggregation/multiscalar_time_series.hpp | 1 + src/aggregation/scalar_time_series.cpp | 1 + src/geometry/spline_curve_1D.cpp | 1 + src/path-finding/molecular_path.cpp | 2 +- src/statistics/kernel_function.cpp | 2 +- 5 files changed, 5 insertions(+), 2 deletions(-) diff --git a/include/aggregation/multiscalar_time_series.hpp b/include/aggregation/multiscalar_time_series.hpp index 454aa7c3..34f10826 100644 --- a/include/aggregation/multiscalar_time_series.hpp +++ b/include/aggregation/multiscalar_time_series.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include "aggregation/scalar_time_series.hpp" diff --git a/src/aggregation/scalar_time_series.cpp b/src/aggregation/scalar_time_series.cpp index 8600ead3..71d2e03a 100644 --- a/src/aggregation/scalar_time_series.cpp +++ b/src/aggregation/scalar_time_series.cpp @@ -22,6 +22,7 @@ // THE SOFTWARE. +#include #include "aggregation/scalar_time_series.hpp" diff --git a/src/geometry/spline_curve_1D.cpp b/src/geometry/spline_curve_1D.cpp index ecfa870c..bf488cbc 100644 --- a/src/geometry/spline_curve_1D.cpp +++ b/src/geometry/spline_curve_1D.cpp @@ -24,6 +24,7 @@ #include #include +#include #include diff --git a/src/path-finding/molecular_path.cpp b/src/path-finding/molecular_path.cpp index 70c1e983..cd9ecfb5 100644 --- a/src/path-finding/molecular_path.cpp +++ b/src/path-finding/molecular_path.cpp @@ -27,7 +27,7 @@ #include #include #include - +#include #include #include diff --git a/src/statistics/kernel_function.cpp b/src/statistics/kernel_function.cpp index be58d524..dc46a30e 100644 --- a/src/statistics/kernel_function.cpp +++ b/src/statistics/kernel_function.cpp @@ -23,7 +23,7 @@ #include - +#include #include "statistics/kernel_function.hpp" From 642342d9ab4fde9e474af9e6e83e2d2ae6b2b520 Mon Sep 17 00:00:00 2001 From: RubenChM Date: Thu, 23 Apr 2026 09:40:27 +0200 Subject: [PATCH 2/6] Update Gromacs library installation notes #46 --- README.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index a4029143..334e2d19 100644 --- a/README.md +++ b/README.md @@ -16,10 +16,11 @@ Prior to installing CHAP, make sure that you have the following libraries and to 3. The [Boost][Boost] C++ libraries, which on Ubuntu can be installed using `sudo apt-get install libboost-all-dev`. Boost algorithms are used in CHAP to solve some root finding and optimisation problems. 4. The CBLAS and LAPACKE linear algebra libraries. On Ubuntu, the easiest way to obtain these is by typing `sudo apt-get install libblas-dev libatlas-base-dev libopenblas-dev liblapacke-dev`. The linear algebra libraries are used in CHAP's spline interpolation. 5. The `libgromacs` library of the [Gromacs][Gromacs] molecular dynamics engine in version 2016 or higher. Comprehensive installation instructions for Gromacs can be found [here][Gromacs-install]. -Please note that for using Gromacs as a library, the underlying FFTW libray +Please note that for using Gromacs as a library, the underlying FFTW library may **not** be installed automatically, i.e. you need to set `-DGMX_BUILD_OWN_FFTW=OFF` when running CMake during the Gromacs installation. +For Gromacs>2019, you also need `-DGMX_INSTALL_LEGACY_API=ON` to get the required headers for CHAP. CHAP also depends on [RapidJSON](http://rapidjson.org/), but this is included as a header-only library, and on [GTest][GTest], but this is downloaded and installed automatically by CMake, so you don't need to do anything about either of these (you will however need Internet access when installing CHAP). From 735f59b278711c388f20cff2d77916ca21fb24dd Mon Sep 17 00:00:00 2001 From: RubenChM Date: Thu, 23 Apr 2026 10:07:28 +0200 Subject: [PATCH 3/6] Update CMake minimum version to 3.5 and Google Test URL #15 --- CMakeLists.txt | 9 +++++---- README.md | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d680ff05..344934f3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,7 +23,7 @@ # minimum version of cmake required for building: -cmake_minimum_required(VERSION 3.2 FATAL_ERROR) +cmake_minimum_required(VERSION 3.5 FATAL_ERROR) # set project parameters: project(CHAP LANGUAGES CXX C VERSION 0.9.1) @@ -98,7 +98,8 @@ include(ExternalProject) # Google test as external project: ExternalProject_Add( googletest - URL https://github.com/google/googletest/archive/release-1.7.0.zip + URL https://github.com/google/googletest/archive/refs/tags/release-1.12.1.zip + DOWNLOAD_EXTRACT_TIMESTAMP TRUE # Disable install step INSTALL_COMMAND "" ) @@ -107,8 +108,8 @@ ExternalProject_Add( ExternalProject_Get_Property(googletest source_dir binary_dir) # set include and library path variables: -set(GTEST_INCLUDE_DIR ${source_dir}/include) -set(GTEST_LIBRARY_PATH ${binary_dir}/${CMAKE_FIND_LIBRARY_PREFIXES}gtest.a) +set(GTEST_INCLUDE_DIR ${source_dir}/googletest/include) +set(GTEST_LIBRARY_PATH ${binary_dir}/lib/${CMAKE_FIND_LIBRARY_PREFIXES}gtest.a) set(GTEST_LIBRARY gtest) include_directories(${GTEST_INCLUDE_DIR}) diff --git a/README.md b/README.md index 334e2d19..fc4b0dd5 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ CHAP is a tool for the functional annotation of ion channel structures written i Prior to installing CHAP, make sure that you have the following libraries and tools installed: -1. The [CMake][CMake] tool in version 3.2 or higher. This will typically be available through your system's package manager. For example, on Ubuntu you can install CMake by typing `sudo apt-get install cmake`. CMake is used to check the availability of libraries and compilers on your system and will ensure that CHAP is installed properly. +1. The [CMake][CMake] tool in version 3.5 or higher. This will typically be available through your system's package manager. For example, on Ubuntu you can install CMake by typing `sudo apt-get install cmake`. CMake is used to check the availability of libraries and compilers on your system and will ensure that CHAP is installed properly. 2. A C++ compiler that supports the `C++11` standard. A popular choice is the [GNU Compiler Collection][GCC], which on Ubuntu can be obtained by typing `sudo apt-get install gcc`. 3. The [Boost][Boost] C++ libraries, which on Ubuntu can be installed using `sudo apt-get install libboost-all-dev`. Boost algorithms are used in CHAP to solve some root finding and optimisation problems. 4. The CBLAS and LAPACKE linear algebra libraries. On Ubuntu, the easiest way to obtain these is by typing `sudo apt-get install libblas-dev libatlas-base-dev libopenblas-dev liblapacke-dev`. The linear algebra libraries are used in CHAP's spline interpolation. From f886e5fee8437076f6fca3191b308ed4b987b27a Mon Sep 17 00:00:00 2001 From: RubenChM Date: Thu, 23 Apr 2026 11:48:23 +0200 Subject: [PATCH 4/6] Update for GROMACS API > 2020 #11 #34 --- README.md | 4 +- .../residue_information_provider.hpp | 1 + .../geometry/abstract_cubic_spline_interp.hpp | 2 +- include/geometry/abstract_spline_curve.hpp | 2 +- include/geometry/spline_curve_1D.hpp | 2 +- include/geometry/spline_curve_3D.hpp | 2 +- include/gromacs/math/3dtransforms.h | 80 ++ include/gromacs/random/seed.h | 108 ++ include/gromacs/random/threefry.h | 954 ++++++++++++++++++ .../gromacs/random/uniformrealdistribution.h | 293 ++++++ include/io/colour.hpp | 2 +- include/io/pdb_io.hpp | 5 +- include/io/wavefront_mtl_io.hpp | 2 +- include/io/wavefront_obj_io.hpp | 2 +- include/path-finding/abstract_path_finder.hpp | 2 - .../abstract_probe_path_finder.hpp | 1 - .../inplane_optimised_probe_path_finder.hpp | 2 - include/path-finding/molecular_path.hpp | 2 +- .../naive_cylindrical_path_finder.hpp | 2 +- include/path-finding/vdw_radius_provider.hpp | 1 + .../chap_trajectory_analysis.hpp | 5 +- .../residue_information_provider.cpp | 4 +- src/geometry/spline_curve_1D.cpp | 2 +- src/io/molecular_path_obj_exporter.cpp | 2 +- src/io/pdb_io.cpp | 10 +- src/main.cpp | 3 +- .../inplane_optimised_probe_path_finder.cpp | 2 +- src/path-finding/molecular_path.cpp | 2 +- src/path-finding/vdw_radius_provider.cpp | 4 +- .../chap_trajectory_analysis.cpp | 43 +- test/path-finding/ut_molecular_path.cpp | 2 +- 31 files changed, 1489 insertions(+), 59 deletions(-) create mode 100644 include/gromacs/math/3dtransforms.h create mode 100644 include/gromacs/random/seed.h create mode 100644 include/gromacs/random/threefry.h create mode 100644 include/gromacs/random/uniformrealdistribution.h diff --git a/README.md b/README.md index fc4b0dd5..7a696474 100644 --- a/README.md +++ b/README.md @@ -15,12 +15,12 @@ Prior to installing CHAP, make sure that you have the following libraries and to 2. A C++ compiler that supports the `C++11` standard. A popular choice is the [GNU Compiler Collection][GCC], which on Ubuntu can be obtained by typing `sudo apt-get install gcc`. 3. The [Boost][Boost] C++ libraries, which on Ubuntu can be installed using `sudo apt-get install libboost-all-dev`. Boost algorithms are used in CHAP to solve some root finding and optimisation problems. 4. The CBLAS and LAPACKE linear algebra libraries. On Ubuntu, the easiest way to obtain these is by typing `sudo apt-get install libblas-dev libatlas-base-dev libopenblas-dev liblapacke-dev`. The linear algebra libraries are used in CHAP's spline interpolation. -5. The `libgromacs` library of the [Gromacs][Gromacs] molecular dynamics engine in version 2016 or higher. Comprehensive installation instructions for Gromacs can be found [here][Gromacs-install]. +5. The `libgromacs` library of the [Gromacs][Gromacs] molecular dynamics engine in version 2020 or higher. Comprehensive installation instructions for Gromacs can be found [here][Gromacs-install]. Please note that for using Gromacs as a library, the underlying FFTW library may **not** be installed automatically, i.e. you need to set `-DGMX_BUILD_OWN_FFTW=OFF` when running CMake during the Gromacs installation. -For Gromacs>2019, you also need `-DGMX_INSTALL_LEGACY_API=ON` to get the required headers for CHAP. +You also need `-DGMX_INSTALL_LEGACY_API=ON` to get the required headers for CHAP. CHAP also depends on [RapidJSON](http://rapidjson.org/), but this is included as a header-only library, and on [GTest][GTest], but this is downloaded and installed automatically by CMake, so you don't need to do anything about either of these (you will however need Internet access when installing CHAP). diff --git a/include/analysis-setup/residue_information_provider.hpp b/include/analysis-setup/residue_information_provider.hpp index 92cac954..be1080f6 100644 --- a/include/analysis-setup/residue_information_provider.hpp +++ b/include/analysis-setup/residue_information_provider.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include "external/rapidjson/document.h" diff --git a/include/geometry/abstract_cubic_spline_interp.hpp b/include/geometry/abstract_cubic_spline_interp.hpp index 3c200849..ba28ecb1 100644 --- a/include/geometry/abstract_cubic_spline_interp.hpp +++ b/include/geometry/abstract_cubic_spline_interp.hpp @@ -27,7 +27,7 @@ #include -#include +#include #include "geometry/basis_spline.hpp" diff --git a/include/geometry/abstract_spline_curve.hpp b/include/geometry/abstract_spline_curve.hpp index 33163b5a..f7b0858a 100644 --- a/include/geometry/abstract_spline_curve.hpp +++ b/include/geometry/abstract_spline_curve.hpp @@ -27,7 +27,7 @@ #include -#include +#include #include "geometry/bspline_basis_set.hpp" diff --git a/include/geometry/spline_curve_1D.hpp b/include/geometry/spline_curve_1D.hpp index 9ec9e793..7f63033b 100644 --- a/include/geometry/spline_curve_1D.hpp +++ b/include/geometry/spline_curve_1D.hpp @@ -30,7 +30,7 @@ #include -#include +#include #include "geometry/abstract_spline_curve.hpp" diff --git a/include/geometry/spline_curve_3D.hpp b/include/geometry/spline_curve_3D.hpp index 1c13cc7d..e01e33f4 100644 --- a/include/geometry/spline_curve_3D.hpp +++ b/include/geometry/spline_curve_3D.hpp @@ -30,7 +30,7 @@ #include #include -#include +#include #include "geometry/abstract_spline_curve.hpp" diff --git a/include/gromacs/math/3dtransforms.h b/include/gromacs/math/3dtransforms.h new file mode 100644 index 00000000..4446fbd8 --- /dev/null +++ b/include/gromacs/math/3dtransforms.h @@ -0,0 +1,80 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 1991- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ +#ifndef GMX_MATH_3DTRANSFORMS_H +#define GMX_MATH_3DTRANSFORMS_H + +#include + +#include "gromacs/utility/real.h" +#include "gromacs/utility/vectypes.h" + +/** Index for the fourth dimension for `vec4`. */ +#define WW 3 + +/*! \brief + * 4D vector type used in 3D transformations. + * + * In \Gromacs, only a limited set of 3D transformations are used, and all of + * them operate on coordinates, so the fourth element is assumed to be one and + * ignored in all contexts. + */ +typedef real vec4[4]; + +/*! \brief + * 4D matrix type used in 3D transformations. + */ +typedef real mat4[4][4]; + +void gmx_mat4_copy(mat4 a, mat4 b); + +void gmx_mat4_transform_point(mat4 m, const rvec x, vec4 v); + +/*! \brief + * Computes the product of two `mat4` matrices as A = B * C. + * + * Note that the order of operands is different from mmul() in vec.h! + */ +void gmx_mat4_mmul(mat4 A, mat4 B, mat4 C); + +void gmx_mat4_init_unity(mat4 m); + +void gmx_mat4_init_rotation(int axis, real angle, mat4 A); + +void gmx_mat4_init_translation(real tx, real ty, real tz, mat4 A); + +void gmx_mat4_print(FILE* fp, const char* s, mat4 A); + +void gmx_vec4_print(FILE* fp, const char* s, vec4 a); + +#endif diff --git a/include/gromacs/random/seed.h b/include/gromacs/random/seed.h new file mode 100644 index 00000000..b85e3e23 --- /dev/null +++ b/include/gromacs/random/seed.h @@ -0,0 +1,108 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2015- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ + +/*! \file + * \brief Random seed and domain utilities + * + * This file contains utilities to create true random seeds from the system, + * and logic to keep track of different random domains for random engines such + * as ThreeFry that can take a second seed value. + * + * \author Erik Lindahl + * \inpublicapi + * \ingroup module_random + */ + +#ifndef GMX_RANDOM_SEED_H +#define GMX_RANDOM_SEED_H + +#include + +#include + +#include "gromacs/utility/basedefinitions.h" + +namespace gmx +{ + +/*! \brief Return 64 random bits from the random device, suitable as seed. + * + * If the internal random device output is smaller than 64 bits, this routine + * will use multiple calls internally until we have 64 bits of random data. + * + * \return 64-bit unsigned integer with random bits. + */ +uint64_t makeRandomSeed(); + +/*! \brief Random device + * + * For now this is identical to the standard library, but since we use + * the GROMACS random module for all other random engines and distributions + * it is convenient to have this too in the same module. + */ +typedef std::random_device RandomDevice; + +/*! \brief Enumerated values for fixed part of random seed (domain) + * + * Random numbers are used in many places in GROMACS, and to avoid identical + * streams the random seeds should be different. Instead of keeping track of + * several different user-provided seeds, it is better to use the fact that + * generators like ThreeFry take two 64-bit keys, and combine a general + * user-provided 64-bit random seed with a second constant value from this list + * to make each stream guaranteed unique. + * + * \note There is no reason to go overboard with adding options; we only + * need to guarantee different streams for cases that might be present + * simultaneously in a single simulation. As an example, two different + * integrators (or thermostats) can reuse the same domain. + * \note When you do add options, leave some space between the values so + * you can group new options with old ones without changing old values. + */ +enum class RandomDomain +{ + Other = 0x00000000, //!< Generic - stream uniqueness is not important + MaxwellVelocities = 0x00001000, //!< Veolcity assignment from Maxwell distribution + TestParticleInsertion = 0x00002000, //!< Test particle insertion + UpdateCoordinates = 0x00003000, //!< Particle integrators + UpdateConstraints = 0x00004000, //!< Second integrator step for constraints + Thermostat = 0x00005000, //!< Stochastic temperature coupling + Barostat = 0x00006000, //!< Stochastic pressure coupling + ReplicaExchange = 0x00007000, //!< Replica exchange metropolis moves + ExpandedEnsemble = 0x00008000, //!< Expanded ensemble lambda moves + AwhBiasing = 0x00009000 //!< AWH biasing reference value moves +}; + +} // namespace gmx + +#endif // GMX_RANDOM_SEED_H diff --git a/include/gromacs/random/threefry.h b/include/gromacs/random/threefry.h new file mode 100644 index 00000000..d7768645 --- /dev/null +++ b/include/gromacs/random/threefry.h @@ -0,0 +1,954 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2015- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ +/*! \defgroup module_random Random engines and distributions (random) + * \ingroup group_utilitymodules + * \brief + * Provides efficient and portable random generators and distributions + * + *

Basic Use

+ * + * \Gromacs relies on random numbers in several different modules, and in + * particular for methods that influence the integration we both require the + * generation to be very fast and the resulting numbers of high quality. + * In addition, it is highly desirable that we generate the same trajectories + * in parallel as for a single-core run. + * + * To realize this, we have implemented the ThreeFry2x64 counter-based random + * engine. In contrast to a normal random engine that is seeded and then keeps + * an internal state, ThreeFry2x64 is derived from cryptographic applications + * where we use a key to turn a highly regular counter int a stream of random + * numbers. This makes it possible to quickly set the counter in the random + * engine based e.g. on the timestep and atom index, and get the same random + * numbers regardless of parallelization. + * + * The TreeFry2x64 engine has been implemented to be fully compatible with + * standard C++11 random engines. There is a gmx::ThreeFry2x64General class that + * allows full control over the accuracy (more iterations means higher quality), + * and gmx::ThreeFry2x64 and gmx::ThreeFry2x64Fast that are specialized to 20 + * and 13 iterations, respectively. With 20 iterations this engine passes all + * tests in the standard BigCrush test, and with 13 iterations only a single + * test fails (in comparision, Mersenne Twister fails two). + * + * All these engines take a template parameter that specifies the number of + * bits to reserve for an internal counter. This is based on an idea of + * John Salmon, and it makes it possible to set your external counter based + * on two simple values (usually timestep and particle index), but then it is + * still possible to draw more than one value for this external counter since + * the internal counter increments. If you run out of internal counter space + * the class will throw an exception to make sure you don't silently end up + * with corrupted/overlapping random data. + * + *

But what if I just want a vanilla random number generator?

+ * + * We've thought about that. Just use the gmx::DefaultRandomEngine class and + * forget everything about counters. Initialize the class with a single value + * for the seed (up to 64 bits), and you are good to go. + * + *

Random number distributions

+ * + * The ThreeFry random engine is fully compatible with all distributions from + * the C++11 standard library, but unfortunately implementation differences + * (and bugs) mean you will typically not get the same sequence of numbers from + * two different library implementations. Since this causes problems for our + * unit tests, we prefer to use our own implementations - they should work + * exactly like the corresponding C++11 versions. + * + * The normal distribution is frequently used in integration, and it can be + * a performance bottleneck. To avoid this, we use a special tabulated + * distribution gmx::TabulatedNormalDistribution that provides very high + * performance at the cost of slightly discretized values; the default 14-bit + * table gives us 16,384 unique values, but this has been thoroughly tested to + * be sufficient for all integration usage. + * + * \author Erik Lindahl + */ +/*! \file + * \brief Implementation of the 2x64 ThreeFry random engine + * + * \author Erik Lindahl + * \inpublicapi + * \ingroup module_random + */ + +#ifndef GMX_RANDOM_THREEFRY_H +#define GMX_RANDOM_THREEFRY_H + +#include +#include +#include + +#include "gromacs/math/functions.h" +#include "gromacs/random/seed.h" +#include "gromacs/utility/classhelpers.h" +#include "gromacs/utility/exceptions.h" + +/* + * The GROMACS implementation of the ThreeFry random engine has been + * heavily inspired by the versions proposed to Boost by: + * + * John Salmon, Copyright 2010-2014 by D. E. Shaw Research + * https://github.com/DEShawResearch/Random123-Boost + * + * Thijs van den Berg, Copyright (c) 2014 M.A. (Thijs) van den Berg + * https://github.com/sitmo/threefry + * + * Both of them are covered by the Boost Software License: + * + * Boost Software License - Version 1.0 - August 17th, 2003 + * + * Permission is hereby granted, free of charge, to any person or organization + * obtaining a copy of the software and accompanying documentation covered by + * this license (the "Software") to use, reproduce, display, distribute, + * execute, and transmit the Software, and to prepare derivative works of the + * Software, and to permit third-parties to whom the Software is furnished to + * do so, all subject to the following: + * + * The copyright notices in the Software and this entire statement, including + * the above license grant, this restriction and the following disclaimer, + * must be included in all copies of the Software, in whole or in part, and + * all derivative works of the Software, unless such copies or derivative + * works are solely in the form of machine-executable object code generated by + * a source language processor. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT + * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE + * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, + * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + */ + +namespace gmx +{ + +namespace internal +{ +// Variable-bitfield counters used to increment internal counters as +// part of std::arrays. + +struct highBitCounter +{ + /*! \brief Clear highBits higest bits of ctr, return false if they were non-zero. + * + * This function clears the space required for the internal counters, + * and returns true if they were correctly zero when calling, false otherwise. + * + * \tparam UIntType Integer type to use for each word in counter + * \tparam words Number of UIntType words in counter + * \tparam highBits Number of bits to check. The template parameter makes it + * possible to optimize this extensively at compile time. + * \param ctr Reference to counter to check and clear. + */ + template + static bool checkAndClear(std::array* ctr) + { + const std::size_t bitsPerWord = std::numeric_limits::digits; + const std::size_t bitsTotal = bitsPerWord * words; + + static_assert(highBits <= bitsTotal, "High bits do not fit in counter."); + + const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord; + const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord; + const UIntType lastWordOne = static_cast(1) << lastWordLowBitIdx; + const UIntType mask = lastWordOne - 1; + + bool isClear = true; + + for (unsigned int i = words - 1; i > lastWordIdx; --i) + { + if ((*ctr)[i]) + { + isClear = false; + (*ctr)[i] = 0; + } + } + if (highBits > 0 && (*ctr)[lastWordIdx] >= lastWordOne) + { + isClear = false; + (*ctr)[lastWordIdx] &= mask; + } + return isClear; + } + + /*! \brief Increment the internal counter in highBits by one + * + * \tparam UIntType Integer type to use for each word in counter + * \tparam words Number of UIntType words in counter + * \tparam highBits Number of bits reserved for the internal counter. + * \param ctr Reference to the counter value to increment. + * + * \throws InternalError if internal counter space is exhausted. + * + * This routine will work across the word boundaries for any number + * of internal counter bits that fits in the total counter. + */ + template + static void increment(std::array* ctr) + { + const std::size_t bitsPerWord = std::numeric_limits::digits; + const std::size_t bitsTotal = bitsPerWord * words; + + static_assert(highBits <= bitsTotal, "High bits do not fit in counter."); + + const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord; + const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord; + const UIntType lastWordOne = static_cast(1) << lastWordLowBitIdx; + + // For algorithm & efficiency reasons we need to store the internal counter in + // the same array as the user-provided counter, so we use the higest bits, possibly + // crossing several words. + // + // To have the computer help us with the dirty carry arithmetics we store the bits + // in the internal counter part in normal fashion, but the internal counter words in + // reverse order; the highest word of the total counter array (words-1) is thus + // the least significant part of the internal counter (if it spans several words). + // + // The incrementation works as follows: + // + // 0) If the index of the least significant internal counter word is larger + // than words-1, there was never any space. + // 1) If the internal counter spans more than one word, we must have one or + // more internal counter words that correspond entirely to the this counter. + // Start with the least significant one (words-1) and increment it. + // If the new value is not zero we did not loop around (no carry), so everything + // is good, and we are done - return! + // If the new value is zero, we need to move the carry result to the next word, + // so we just continue the loop until we have gone through all words that + // are internal-counter-only. + // 2) After the loop, there is stuff remaining to add, and by definition there + // is some internal counter space in the next word, but the question + // is if we have exhausted it. We already created a constant that corresponds + // to the bit that represents '1' for the internal counter part of this word. + // When we add this constant it will not affect the user-counter-part at all, + // and if we exhaust the internal counter space the high bits will cause the entire + // word to wrap around, and the result will be smaller than the bit we added. + // If this happens we throw, otherwise we're done. + // + // Since all constants will be evaluated at compile time, this entire routine + // will usually be reduced to simply incrementing a word by a constant, and throwing + // if the result is smaller than the constant. + + if (lastWordIdx >= words) + { + GMX_THROW(InternalError( + "Cannot increment random engine defined with 0 internal counter bits.")); + } + + for (unsigned int i = words - 1; i > lastWordIdx; --i) + { + (*ctr)[i]++; + if ((*ctr)[i]) + { + return; // No carry means we are done + } + } + (*ctr)[lastWordIdx] += lastWordOne; + if ((*ctr)[lastWordIdx] < lastWordOne) + { + GMX_THROW(InternalError("Random engine stream ran out of internal counter space.")); + } + } + + /*! \brief Increment the internal counter in highBits by a value. + * + * \tparam UIntType Integer type to use for each word in counter + * \tparam words Number of UIntType words in counter + * \tparam highBits Number of bits reserved for the internal counter. + * \param ctr Reference to the counter to increment. + * \param addend Value to add to internal. + * + * \throws InternalError if internal counter space is exhausted. + * + * This routine will work across the word boundaries for any number + * of internal counter bits that fits in the total counter. + */ + template + static void increment(std::array* ctr, UIntType addend) + { + const std::size_t bitsPerWord = std::numeric_limits::digits; + const std::size_t bitsTotal = bitsPerWord * words; + + static_assert(highBits <= bitsTotal, "High bits do not fit in counter."); + + const std::size_t lastWordIdx = (bitsTotal - highBits) / bitsPerWord; + const std::size_t lastWordLowBitIdx = (bitsTotal - highBits) % bitsPerWord; + const UIntType lastWordOne = static_cast(1) << lastWordLowBitIdx; + const UIntType lastWordMaxVal = (~static_cast(0)) >> lastWordLowBitIdx; + + if (lastWordIdx >= words) + { + GMX_THROW(InternalError( + "Cannot increment random engine defined with 0 internal counter bits.")); + } + + for (unsigned int i = words - 1; i > lastWordIdx; --i) + { + (*ctr)[i] += addend; + addend = ((*ctr)[i] < addend); // 1 is the carry! + if (addend == 0) + { + return; + } + } + + if (addend > lastWordMaxVal) + { + GMX_THROW(InternalError("Random engine stream ran out of internal counter space.")); + } + addend *= lastWordOne; + + (*ctr)[lastWordIdx] += addend; + + if ((*ctr)[lastWordIdx] < addend) + { + GMX_THROW(InternalError("Random engine stream ran out of internal counter space.")); + } + } +}; +} // namespace internal + +/*! \brief General implementation class for ThreeFry counter-based random engines. + * + * This class is used to implement several different ThreeFry2x64 random engines + * differing in the number of rounds executed in and the number of bits reserved + * for the internal counter. It is compatible with C++11 random engines, and + * can be used e.g. with all random distributions from the standard library. + * + * ThreeFry is a counter-based rather than state-based random engine. This + * means that we seed it with a "key", after which we can get the + * N:th random number in a sequence (specified by a counter) directly. This + * means we are guaranteed the same sequence of numbers even when running in + * parallel if using e.g. step and atom index as counters. + * + * However, it is also useful to be able to use it as a normal random engine, + * for instance if you need more than 2 64-bit random values for a specific + * counter value, not to mention where you just need good normal random numbers. + * To achieve this, this implementation uses John Salmon's idea of reserving + * a couple of the highest bits in the user-provided counter for an internal + * counter. For instance, if reserving 3 bits, this means you get a stream of + * 8 iterations (each with 2 random values) after every restart. If you call + * the engine after these bits have been exhausted, it will throw an + * exception to make sure you don't get overlapping streams by mistake. + * Reserving 3 bits also means you can only use 64-3=61 bits of the highest + * word when restarting (i.e., setting) the counters. + * + * This version also supports using internalCounterBits=0. In this case the + * random engine will be able to return a single counter round, i.e. 2 64-bit + * values for ThreeFry2x64, after which an exception is thrown. In this case no + * high bits are reserved, which means the class implements the raw ThreeFry2x64 + * random function. + * + * \tparam rounds The number of encryption iterations used when generating. + * This can in principle be any value, but 20 rounds has been + * shown to pass all BigCrush random tests, and with 13 rounds + * only one fails. This is a very stringent test, and the + * standard Mersenne Twister engine fails two, so 13 rounds + * should be a perfectly fine balance in most cases. + * \tparam internalCounterBits + * Number of high bits in the user-provided counter reserved + * for the internal counter. The number of values the engine + * can return after each restart will be + * words*2^internalCounterBits. + */ +template +class ThreeFry2x64General +{ + // While this class will formally work with any value for rounds, there is + // no reason to go lower than 13, and this might help catch some typos. + // If we find a reason to use lower values in the future, or if you simply + // want to test, this assert can safely be removed. + static_assert(rounds >= 13, + "You should not use less than 13 encryption rounds for ThreeFry2x64."); + +public: + // result_type must be lower case to be compatible with C++11 standard library + + /*! \brief Integer type for output. */ + typedef uint64_t result_type; + /*! \brief Use array for counter & key states so it is allocated on the stack */ + typedef std::array counter_type; + +private: + /*! \brief Rotate value left by specified number of bits + * + * \param i Value to rotate (result_type, which should be 64-bit). + * \param bits Number of bits to rotate i. + * + * \return Input value rotated 'bits' left. + */ + result_type rotLeft(result_type i, unsigned int bits) + { + return (i << bits) | (i >> (std::numeric_limits::digits - bits)); + } + + /*! \brief Perform encryption step for ThreeFry2x64 algorithm + * + * It performs the encryption step of the standard ThreeFish symmetric-key + * tweakable block cipher, which is the core of the ThreeFry random + * engine. The number of encryption rounds is specified by the class + * template parameter 'rounds'. + * + * \param key Reference to key value + * \param ctr Counter value to use + * + * \return Newly encrypted 2x64 block, according to the class template parameters. + */ + counter_type generateBlock(const counter_type& key, const counter_type& ctr) + { + const unsigned int rotations[] = { 16, 42, 12, 31, 16, 32, 24, 21 }; + counter_type x = ctr; + + result_type ks[3] = { 0x0, 0x0, 0x1bd11bdaa9fc1a22 }; + + // This is actually a pretty simple routine that merely executes the + // for-block specified further down 'rounds' times. However, both + // clang and gcc have problems unrolling and replacing rotations[r%8] + // with constants, so we unroll the first 20 iterations manually. + + if (rounds > 0) + { + ks[0] = key[0]; + ks[2] ^= key[0]; + x[0] = x[0] + key[0]; + ks[1] = key[1]; + ks[2] ^= key[1]; + x[1] = x[1] + key[1]; + x[0] += x[1]; + x[1] = rotLeft(x[1], 16); + x[1] ^= x[0]; + } + if (rounds > 1) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 42); + x[1] ^= x[0]; + } + if (rounds > 2) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 12); + x[1] ^= x[0]; + } + if (rounds > 3) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 31); + x[1] ^= x[0]; + x[0] += ks[1]; + x[1] += ks[2] + 1; + } + if (rounds > 4) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 16); + x[1] ^= x[0]; + } + if (rounds > 5) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 32); + x[1] ^= x[0]; + } + if (rounds > 6) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 24); + x[1] ^= x[0]; + } + if (rounds > 7) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 21); + x[1] ^= x[0]; + x[0] += ks[2]; + x[1] += ks[0] + 2; + } + if (rounds > 8) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 16); + x[1] ^= x[0]; + } + if (rounds > 9) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 42); + x[1] ^= x[0]; + } + if (rounds > 10) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 12); + x[1] ^= x[0]; + } + if (rounds > 11) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 31); + x[1] ^= x[0]; + x[0] += ks[0]; + x[1] += ks[1] + 3; + } + if (rounds > 12) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 16); + x[1] ^= x[0]; + } + if (rounds > 13) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 32); + x[1] ^= x[0]; + } + if (rounds > 14) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 24); + x[1] ^= x[0]; + } + if (rounds > 15) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 21); + x[1] ^= x[0]; + x[0] += ks[1]; + x[1] += ks[2] + 4; + } + if (rounds > 16) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 16); + x[1] ^= x[0]; + } + if (rounds > 17) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 42); + x[1] ^= x[0]; + } + if (rounds > 18) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 12); + x[1] ^= x[0]; + } + if (rounds > 19) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], 31); + x[1] ^= x[0]; + x[0] += ks[2]; + x[1] += ks[0] + 5; + } + + for (unsigned int r = 20; r < rounds; r++) + { + x[0] += x[1]; + x[1] = rotLeft(x[1], rotations[r % 8]); + x[1] ^= x[0]; + if (((r + 1) & 3) == 0) + { + unsigned int r4 = (r + 1) >> 2; + x[0] += ks[r4 % 3]; + x[1] += ks[(r4 + 1) % 3] + r4; + } + } + return x; + } + +public: + //! \brief Smallest value that can be returned from random engine. +#if !defined(_MSC_VER) + static constexpr +#else + // Avoid constexpr bug in MSVC 2015, note that max() below does work + static +#endif + result_type + min() + { + return std::numeric_limits::min(); + } + + //! \brief Largest value that can be returned from random engine. + static constexpr result_type max() { return std::numeric_limits::max(); } + + /*! \brief Construct random engine with 2x64 key values + * + * This constructor takes two values, and should only be used with + * the 2x64 implementations. + * + * \param key0 Random seed in the form of a 64-bit unsigned value. + * \param domain Random domain. This is used to guarantee that different + * applications of a random engine inside the code get different + * streams of random numbers, without requiring the user + * to provide lots of random seeds. Pick a value from the + * RandomDomain class, or RandomDomain::Other if it is + * not important. In the latter case you might want to use + * \ref gmx::DefaultRandomEngine instead. + * + * \note The random domain is really another 64-bit seed value. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. + */ + //NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + ThreeFry2x64General(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) + { + seed(key0, domain); + } + + /*! \brief Construct random engine from 2x64-bit unsigned integers + * + * This constructor assigns the raw 128 bit key data from unsigned integers. + * It is meant for the case when you want full control over the key, + * for instance to compare with reference values of the ThreeFry + * function during testing. + * + * \param key0 First word of key/random seed. + * \param key1 Second word of key/random seed. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. To test arbitrary values, use 0 internal counter bits. + */ + //NOLINTNEXTLINE(cppcoreguidelines-pro-type-member-init) + ThreeFry2x64General(uint64_t key0, uint64_t key1) { seed(key0, key1); } + + /*! \brief Seed 2x64 random engine with two 64-bit key values + * + * \param key0 First word of random seed, in the form of 64-bit unsigned values. + * \param domain Random domain. This is used to guarantee that different + * applications of a random engine inside the code get different + * streams of random numbers, without requiring the user + * to provide lots of random seeds. Pick a value from the + * RandomDomain class, or RandomDomain::Other if it is + * not important. In the latter case you might want to use + * \ref gmx::DefaultRandomEngine instead. + * + * \note The random domain is really another 64-bit seed value. + * + * Re-initialized the seed similar to the counter constructor. + * Same rules apply: The highest few bits of the last word are + * reserved to encode the number of internal counter bits, but + * to save the user the trouble of making sure these are zero + * when using e.g. a random device, we just ignore them. + */ + void seed(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) + { + seed(key0, static_cast(domain)); + } + + /*! \brief Seed random engine from 2x64-bit unsigned integers + * + * This assigns the raw 128 bit key data from unsigned integers. + * It is meant for the case when you want full control over the key, + * for instance to compare with reference values of the ThreeFry + * function during testing. + * + * \param key0 First word of key/random seed. + * \param key1 Second word of key/random seed. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. To test arbitrary values, use 0 internal counter bits. + */ + void seed(uint64_t key0, uint64_t key1) + { + const unsigned int internalCounterBitsBits = + (internalCounterBits > 0) ? (StaticLog2::value + 1) : 0; + + key_ = { { key0, key1 } }; + + if (internalCounterBits > 0) + { + internal::highBitCounter::checkAndClear(&key_); + internal::highBitCounter::increment( + &key_, internalCounterBits - 1); + } + restart(0, 0); + } + + /*! \brief Restart 2x64 random engine counter from 2 64-bit values + * + * \param ctr0 First word of new counter, in the form of 64-bit unsigned values. + * \param ctr1 Second word of new counter + * + * Restarting the engine with a new counter is extremely fast with ThreeFry64, + * and basically just consists of storing the counter value, so you should + * use this liberally in your innermost loops to restart the engine with + * e.g. the current step and atom index as counter values. + * + * \throws InternalError if any of the highest bits that are reserved + * for the internal part of the counter are set. The number of + * reserved bits is to the last template parameter to the class. + */ + void restart(uint64_t ctr0 = 0, uint64_t ctr1 = 0) + { + + counter_ = { { ctr0, ctr1 } }; + if (!internal::highBitCounter::checkAndClear(&counter_)) + { + GMX_THROW(InternalError( + "High bits of counter are reserved for the internal stream counter.")); + } + block_ = generateBlock(key_, counter_); + index_ = 0; + } + + /*! \brief Generate the next random number + * + * This will return the next stored 64-bit value if one is available, + * and otherwise generate a new block, update the internal counters, and + * return the first value while storing the others. + * + * \throws InternalError if the internal counter space is exhausted. + */ + result_type operator()() + { + if (index_ >= c_resultsPerCounter_) + { + internal::highBitCounter::increment(&counter_); + block_ = generateBlock(key_, counter_); + index_ = 0; + } + return block_[index_++]; + } + + /*! \brief Skip next n random numbers + * + * Moves the internal random stream for the give key/counter value + * n positions forward. The count is based on the number of random values + * returned, such that skipping 5 values gives exactly the same result as + * drawing 5 values that are ignored. + * + * \param n Number of values to jump forward. + * + * \throws InternalError if the internal counter space is exhausted. + */ + void discard(uint64_t n) + { + index_ += n % c_resultsPerCounter_; + n /= c_resultsPerCounter_; + + if (index_ > c_resultsPerCounter_) + { + index_ -= c_resultsPerCounter_; + n++; + } + + // Make sure the state is the same as if we came to this counter and + // index by natural generation. + if (index_ == 0 && n > 0) + { + index_ = c_resultsPerCounter_; + n--; + } + internal::highBitCounter::increment(&counter_, n); + block_ = generateBlock(key_, counter_); + } + + /*! \brief Return true if two ThreeFry2x64 engines are identical + * + * \param x Instance to compare with. + * + * This routine should return true if the two engines will generate + * identical random streams when drawing. + */ + bool operator==(const ThreeFry2x64General& x) const + { + // block_ is uniquely specified by key_ and counter_. + return (key_ == x.key_ && counter_ == x.counter_ && index_ == x.index_); + } + + /*! \brief Return true of two ThreeFry2x64 engines are not identical + * + * \param x Instance to compare with. + * + * This routine should return true if the two engines will generate + * different random streams when drawing. + */ + bool operator!=(const ThreeFry2x64General& x) const + { + return !operator==(x); + } + +private: + /*! \brief Number of results returned for each invocation of the block generation */ + static const unsigned int c_resultsPerCounter_ = + static_cast(sizeof(counter_type) / sizeof(result_type)); + + /*! \brief ThreeFry2x64 key, i.e. the random seed for this stream. + * + * The highest few bits of the key are replaced to encode the value of + * internalCounterBits, in order to make all streams unique. + */ + counter_type key_; + + /*! \brief ThreeFry2x64 total counter. + * + * The highest internalCounterBits are reserved for an internal counter + * so that the combination of a key and counter provides a stream that + * returns 2*2^internalCounterBits (ThreeFry2x64) random 64-bit values before + * the internal counter space is exhausted and an exception is thrown. + */ + counter_type counter_; + /*! \brief The present block encrypted from values of key and counter. */ + counter_type block_; + /*! \brief Index of the next value in block_ to return from random engine */ + unsigned int index_; + + GMX_DISALLOW_COPY_AND_ASSIGN(ThreeFry2x64General); +}; + + +/*! \brief ThreeFry2x64 random engine with 20 iteractions. + * + * \tparam internalCounterBits, default 64. + * + * This class provides very high quality random numbers that pass all + * BigCrush tests, it works with two 64-bit values each for keys and + * counters, and is most efficient when we only need a few random values + * before restarting the counters with new values. + */ +template +class ThreeFry2x64 : public ThreeFry2x64General<20, internalCounterBits> +{ +public: + /*! \brief Construct ThreeFry random engine with 2x64 key values, 20 rounds. + * + * \param key0 Random seed in the form of a 64-bit unsigned value. + * \param domain Random domain. This is used to guarantee that different + * applications of a random engine inside the code get different + * streams of random numbers, without requiring the user + * to provide lots of random seeds. Pick a value from the + * RandomDomain class, or RandomDomain::Other if it is + * not important. In the latter case you might want to use + * \ref gmx::DefaultRandomEngine instead. + * + * \note The random domain is really another 64-bit seed value. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. + */ + ThreeFry2x64(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) : + ThreeFry2x64General<20, internalCounterBits>(key0, domain) + { + } + + /*! \brief Construct random engine from 2x64-bit unsigned integers, 20 rounds + * + * This constructor assigns the raw 128 bit key data from unsigned integers. + * It is meant for the case when you want full control over the key, + * for instance to compare with reference values of the ThreeFry + * function during testing. + * + * \param key0 First word of key/random seed. + * \param key1 Second word of key/random seed. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. To test arbitrary values, use 0 internal counter bits. + */ + ThreeFry2x64(uint64_t key0, uint64_t key1) : + ThreeFry2x64General<20, internalCounterBits>(key0, key1) + { + } +}; + +/*! \brief ThreeFry2x64 random engine with 13 iteractions. + * + * \tparam internalCounterBits, default 64. + * + * This class provides relatively high quality random numbers that only + * fail one BigCrush test, and it is a bit faster than the 20-round version. + * It works with two 64-bit values each for keys and counters, and is most + * efficient when we only need a few random values before restarting + * the counters with new values. + */ +template +class ThreeFry2x64Fast : public ThreeFry2x64General<13, internalCounterBits> +{ +public: + /*! \brief Construct ThreeFry random engine with 2x64 key values, 13 rounds. + * + * \param key0 Random seed in the form of a 64-bit unsigned value. + * \param domain Random domain. This is used to guarantee that different + * applications of a random engine inside the code get different + * streams of random numbers, without requiring the user + * to provide lots of random seeds. Pick a value from the + * RandomDomain class, or RandomDomain::Other if it is + * not important. In the latter case you might want to use + * \ref gmx::DefaultRandomEngine instead. + * + * \note The random domain is really another 64-bit seed value. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. + */ + ThreeFry2x64Fast(uint64_t key0 = 0, RandomDomain domain = RandomDomain::Other) : + ThreeFry2x64General<13, internalCounterBits>(key0, domain) + { + } + + /*! \brief Construct ThreeFry random engine from 2x64-bit unsigned integers, 13 rounds. + * + * This constructor assigns the raw 128 bit key data from unsigned integers. + * It is meant for the case when you want full control over the key, + * for instance to compare with reference values of the ThreeFry + * function during testing. + * + * \param key0 First word of key/random seed. + * \param key1 Second word of key/random seed. + * + * \throws InternalError if the high bits needed to encode the number of counter + * bits are nonzero. To test arbitrary values, use 0 internal counter bits. + */ + ThreeFry2x64Fast(uint64_t key0, uint64_t key1) : + ThreeFry2x64General<13, internalCounterBits>(key0, key1) + { + } +}; + + +/*! \brief Default fast and accurate random engine in Gromacs + * + * This engine will return 2*2^64 random results using the default + * gmx::RandomDomain::Other stream, and can be initialized with a single + * seed argument without having to remember empty template angle brackets. + */ +typedef ThreeFry2x64Fast<> DefaultRandomEngine; + +} // namespace gmx + +#endif // GMX_RANDOM_THREEFRY_H diff --git a/include/gromacs/random/uniformrealdistribution.h b/include/gromacs/random/uniformrealdistribution.h new file mode 100644 index 00000000..a51e6b2e --- /dev/null +++ b/include/gromacs/random/uniformrealdistribution.h @@ -0,0 +1,293 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright 2015- The GROMACS Authors + * and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. + * Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * https://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at https://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out https://www.gromacs.org. + */ + +/*! \file + * \brief The uniform real distribution + * + * Portable version of the uniform real that generates the same sequence + * on all platforms. Since stdlibc++ and libc++ provide different sequences + * we prefer this one so unit tests produce the same values on all platforms. + * + * \author Erik Lindahl + * \inpublicapi + * \ingroup module_random + */ + +#ifndef GMX_RANDOM_UNIFORMREALDISTRIBUTION_H +#define GMX_RANDOM_UNIFORMREALDISTRIBUTION_H + +#include + +#include +#include +#include +#include + +#include "gromacs/math/functions.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/classhelpers.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/real.h" + +/* + * The portable version of the uniform real distribution (to make sure we get + * the same values on all platforms) has been modified from the LLVM libcxx + * headers, distributed under the MIT license: + * + * Copyright (c) The LLVM compiler infrastructure + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN + * THE SOFTWARE. + */ + +namespace gmx +{ + +/*! \brief Generate a floating-point value with specified number of random bits + * + * \tparam RealType Floating-point type to generate + * \tparam Bits Number of random bits to generate + * \tparam Rng Random number generator class + * + * \param g Random number generator to use + * + * This implementation avoids the bug in libc++ and stdlibc++ (which is due + * to the C++ standard being unclear) where 1.0 can be returned occasionally. + * + */ +template +RealType generateCanonical(Rng& g) +{ + // No point in using more bits than fit in RealType + const uint64_t digits = std::numeric_limits::digits; + const uint64_t realBits = std::min(digits, static_cast(Bits)); + const uint64_t log2R = std::numeric_limits::digits; + uint64_t k = realBits / log2R + (realBits % log2R != 0) + (realBits == 0); + // Note that Rng::max and Rng::min are typically an integer type. + // Only unsigned integer types can express the range using the + // same type. Converting to RealType before computing the range + // would work but we have no need for that. + static_assert(std::is_unsigned_v && std::is_unsigned_v, + "Rng::max and Rng::min must be unsigned"); + RealType r = RealType(Rng::max() - Rng::min()) + RealType(1); + RealType s = g() - Rng::min(); + RealType base = r; + RealType result; + + for (uint64_t i = 1; i < k; ++i) + { + s += RealType(g() - Rng::min()) * base; + base *= r; + } + result = s / base; + + // This implementation is specified by the C++ standard, but unfortunately it + // has a bug where 1.0 can be generated occasionally due to the limited + // precision of floating point, while 0.0 is only generated half as often as + // it should. We "solve" both these issues by swapping 1.0 for 0.0 when it happens. + // + // See: + // https://llvm.org/bugs/show_bug.cgi?id=18767 + // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176 + // + // Note that we prefer not to use the gcc 'fix' of looping until the result + // is smaller than 1.0, since that breaks the strict specification of the + // number of times the rng will be called. + // + // This can only happen when we ask for the same number of bits that fit + // in RealType, so by checking for that we avoid the extra code in all other + // cases. If you are worried about it: Use RealType=double with 32 bits. + // + if (realBits == digits && result == 1.0) + { + result = 0.0; + } + return result; +} + + +/*! \brief Uniform real distribution + * + * The C++ standard library does provide this distribution, but even + * though they all sample from the correct distribution different standard + * library implementations appear to return different sequences of numbers + * for the same random number generator. To make it easier to use GROMACS + * unit tests that depend on random numbers we have our own implementation. + * + * \tparam RealType Floating-point type, real by default in GROMACS. + */ +template +class UniformRealDistribution +{ +public: + /*! \brief Type of values returned */ + typedef RealType result_type; + + /*! \brief Uniform real distribution parameters */ + class param_type + { + /*! \brief Lower end of range (inclusive) */ + result_type a_; + /*! \brief Upper end of range (exclusive) */ + result_type b_; + + public: + /*! \brief Reference back to the distribution class */ + typedef UniformRealDistribution distribution_type; + + /*! \brief Construct parameter block + * + * \param a Lower end of range (inclusive) + * \param b Upper end of range (exclusive) + */ + explicit param_type(result_type a = 0.0, result_type b = 1.0) : a_(a), b_(b) + { + GMX_RELEASE_ASSERT(a < b, "The uniform real distribution requires a + result_type operator()(Rng& g) + { + return (*this)(g, param_); + } + + /*! \brief Return value from uniform real distribution with given parameters + * + * \tparam Rng Random engine class + * + * \param g Random engine + * \param param Parameters to use + */ + template + result_type operator()(Rng& g, const param_type& param) + { + result_type r = generateCanonical::digits>(g); + return (param.b() - param.a()) * r + param.a(); + } + + /*! \brief Return the lower range uniform real distribution */ + result_type a() const { return param_.a(); } + + /*! \brief Return the upper range of the uniform real distribution */ + result_type b() const { return param_.b(); } + + /*! \brief Return the full parameter class of the uniform real distribution */ + param_type param() const { return param_; } + + /*! \brief Smallest value that can be returned from uniform real distribution */ + result_type min() const { return a(); } + + /*! \brief Largest value that can be returned from uniform real distribution */ + result_type max() const { return b(); } + + /*! \brief True if two uniform real distributions will produce the same values. + * + * \param x Instance to compare with. + */ + bool operator==(const UniformRealDistribution& x) const { return param_ == x.param_; } + + /*! \brief True if two uniform real distributions will produce different values. + * + * \param x Instance to compare with. + */ + bool operator!=(const UniformRealDistribution& x) const { return !operator==(x); } + +private: + /*! \brief Internal value for parameters, can be overridden at generation time. */ + param_type param_; + + GMX_DISALLOW_COPY_AND_ASSIGN(UniformRealDistribution); +}; + +} // namespace gmx + +#endif // GMX_RANDOM_UNIFORMREALDISTRIBUTION_H diff --git a/include/io/colour.hpp b/include/io/colour.hpp index 118f380c..8cd413ec 100644 --- a/include/io/colour.hpp +++ b/include/io/colour.hpp @@ -29,7 +29,7 @@ #include #include -#include +#include #include "external/rapidjson/document.h" diff --git a/include/io/pdb_io.hpp b/include/io/pdb_io.hpp index 71f4c994..bc7ac320 100644 --- a/include/io/pdb_io.hpp +++ b/include/io/pdb_io.hpp @@ -30,6 +30,7 @@ #include #include +#include #include #include @@ -58,8 +59,8 @@ class PdbStructure // data required for writing PDB file: t_atoms atoms_; - rvec *coords_; - int ePBC_; + const rvec *coords_; + PbcType ePBC_; matrix box_; }; diff --git a/include/io/wavefront_mtl_io.hpp b/include/io/wavefront_mtl_io.hpp index 49b1919e..5eb42eba 100644 --- a/include/io/wavefront_mtl_io.hpp +++ b/include/io/wavefront_mtl_io.hpp @@ -29,7 +29,7 @@ #include #include -#include +#include #include diff --git a/include/io/wavefront_obj_io.hpp b/include/io/wavefront_obj_io.hpp index e294ca3f..16273c64 100644 --- a/include/io/wavefront_obj_io.hpp +++ b/include/io/wavefront_obj_io.hpp @@ -32,7 +32,7 @@ #include #include -#include +#include #include diff --git a/include/path-finding/abstract_path_finder.hpp b/include/path-finding/abstract_path_finder.hpp index 3588eccf..8d053dbc 100644 --- a/include/path-finding/abstract_path_finder.hpp +++ b/include/path-finding/abstract_path_finder.hpp @@ -27,8 +27,6 @@ #include -#include - #include "path-finding/molecular_path.hpp" diff --git a/include/path-finding/abstract_probe_path_finder.hpp b/include/path-finding/abstract_probe_path_finder.hpp index 0a922fbe..4182e992 100644 --- a/include/path-finding/abstract_probe_path_finder.hpp +++ b/include/path-finding/abstract_probe_path_finder.hpp @@ -27,7 +27,6 @@ #include -#include #include #include "path-finding/abstract_path_finder.hpp" diff --git a/include/path-finding/inplane_optimised_probe_path_finder.hpp b/include/path-finding/inplane_optimised_probe_path_finder.hpp index ac73037f..7c73be00 100644 --- a/include/path-finding/inplane_optimised_probe_path_finder.hpp +++ b/include/path-finding/inplane_optimised_probe_path_finder.hpp @@ -29,8 +29,6 @@ #include #include -#include - #include "path-finding/abstract_probe_path_finder.hpp" diff --git a/include/path-finding/molecular_path.hpp b/include/path-finding/molecular_path.hpp index 15ae7e40..06287736 100644 --- a/include/path-finding/molecular_path.hpp +++ b/include/path-finding/molecular_path.hpp @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include #include diff --git a/include/path-finding/naive_cylindrical_path_finder.hpp b/include/path-finding/naive_cylindrical_path_finder.hpp index daf4d7b9..3aeba518 100644 --- a/include/path-finding/naive_cylindrical_path_finder.hpp +++ b/include/path-finding/naive_cylindrical_path_finder.hpp @@ -26,7 +26,7 @@ #define NAIVE_CYLINDRICAL_PATH_FINDER_HPP #include -#include +#include #include "path-finding/abstract_path_finder.hpp" diff --git a/include/path-finding/vdw_radius_provider.hpp b/include/path-finding/vdw_radius_provider.hpp index 4c74ba04..7237385d 100644 --- a/include/path-finding/vdw_radius_provider.hpp +++ b/include/path-finding/vdw_radius_provider.hpp @@ -30,6 +30,7 @@ #include #include +#include #include #include "external/rapidjson/document.h" diff --git a/include/trajectory-analysis/chap_trajectory_analysis.hpp b/include/trajectory-analysis/chap_trajectory_analysis.hpp index ffd8c0a4..317fd0d2 100644 --- a/include/trajectory-analysis/chap_trajectory_analysis.hpp +++ b/include/trajectory-analysis/chap_trajectory_analysis.hpp @@ -30,7 +30,10 @@ #include #include -#include +#include +#include +#include +#include #include "analysis-setup/residue_information_provider.hpp" diff --git a/src/analysis-setup/residue_information_provider.cpp b/src/analysis-setup/residue_information_provider.cpp index 88de2e33..4c5973ee 100644 --- a/src/analysis-setup/residue_information_provider.cpp +++ b/src/analysis-setup/residue_information_provider.cpp @@ -47,7 +47,7 @@ ResidueInformationProvider::nameFromTopology( const gmx::TopologyInformation &top) { // get list of all atoms (including residue information): - t_atoms atoms = top.topology() -> atoms; + const t_atoms atoms = *(top.atoms()); // loop over all residues: for(int i = 0; i < atoms.nres; i++) @@ -67,7 +67,7 @@ ResidueInformationProvider::chainFromTopology( const gmx::TopologyInformation &top) { // get list of all atoms (includingb residue information): - t_atoms atoms = top.topology() -> atoms; + const t_atoms atoms = *(top.atoms()); // loop over all residues: for(int i = 0; i < atoms.nres; i++) diff --git a/src/geometry/spline_curve_1D.cpp b/src/geometry/spline_curve_1D.cpp index bf488cbc..1cfd659b 100644 --- a/src/geometry/spline_curve_1D.cpp +++ b/src/geometry/spline_curve_1D.cpp @@ -230,7 +230,7 @@ SplineCurve1D::minimum(const std::pair &lim) { // internal parameters: real maxSampleDist = 0.1; - boost::uintmax_t maxIter = 100; + std::uintmax_t maxIter = 100; // draw sample of values and find minimum: real length = lim.second - lim.first; diff --git a/src/io/molecular_path_obj_exporter.cpp b/src/io/molecular_path_obj_exporter.cpp index 1cfb311b..0dc18084 100644 --- a/src/io/molecular_path_obj_exporter.cpp +++ b/src/io/molecular_path_obj_exporter.cpp @@ -28,7 +28,7 @@ #include #include -#include +#include #include #include "io/molecular_path_obj_exporter.hpp" diff --git a/src/io/pdb_io.cpp b/src/io/pdb_io.cpp index 6e73e121..f9bb8801 100644 --- a/src/io/pdb_io.cpp +++ b/src/io/pdb_io.cpp @@ -40,14 +40,14 @@ PdbStructure::fromTopology( const gmx::TopologyInformation &top) { // retrieve coordinates and box matrix from topology: - top.getTopologyConf(&coords_, box_); + coords_ = gmx::as_rvec_array(top.x().data()); + top.getBox(box_); // retrieve list of atoms in topology: - t_topology *topol = top.topology(); - atoms_ = topol -> atoms; + atoms_ = *(top.atoms()); // retrieve periodic BC: - ePBC_ = top.ePBC(); + ePBC_ = top.pbcType(); } @@ -69,7 +69,7 @@ PdbStructure::setPoreFacing( // fill with proper values (occup and bfac assigned zero by default): for(size_t i = 0; i < atoms_.nr; i++) { - atoms_.pdbinfo[i].type = epdbATOM; + atoms_.pdbinfo[i].type = PdbRecordType::Atom; atoms_.pdbinfo[i].atomnr = i + 1; atoms_.pdbinfo[i].altloc = ' '; atoms_.pdbinfo[i].occup = 0.0; diff --git a/src/main.cpp b/src/main.cpp index 4b387338..449ad06f 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -26,6 +26,7 @@ #include "config/back_matter.hpp" #include "config/front_matter.hpp" +#include #include "trajectory-analysis/chap_trajectory_analysis.hpp" using namespace gmx; @@ -45,7 +46,7 @@ int main(int argc, char **argv) argc++; // run trajectory analysis: - int status = TrajectoryAnalysisCommandLineRunner::runAsMain(argc, argv); + int status = gmx::TrajectoryAnalysisCommandLineRunner::runAsMain(argc, argv); // print back matter: BackMatter::print(); diff --git a/src/path-finding/inplane_optimised_probe_path_finder.cpp b/src/path-finding/inplane_optimised_probe_path_finder.cpp index 2f59bb4a..d87b2657 100644 --- a/src/path-finding/inplane_optimised_probe_path_finder.cpp +++ b/src/path-finding/inplane_optimised_probe_path_finder.cpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include "optim/simulated_annealing_module.hpp" #include "optim/nelder_mead_module.hpp" diff --git a/src/path-finding/molecular_path.cpp b/src/path-finding/molecular_path.cpp index cd9ecfb5..10c0f1bf 100644 --- a/src/path-finding/molecular_path.cpp +++ b/src/path-finding/molecular_path.cpp @@ -587,7 +587,7 @@ MolecularPath::minRadius() { // internal parameters: real maxSampleDist = 0.1; - boost::uintmax_t maxIter = 100; + std::uintmax_t maxIter = 100; // draw radius samples along path: int nSamples = std::ceil(length_/maxSampleDist); diff --git a/src/path-finding/vdw_radius_provider.cpp b/src/path-finding/vdw_radius_provider.cpp index a3a7389b..97c41a4a 100644 --- a/src/path-finding/vdw_radius_provider.cpp +++ b/src/path-finding/vdw_radius_provider.cpp @@ -28,7 +28,7 @@ //#include #include #include - +#include #include "path-finding/vdw_radius_provider.hpp" @@ -182,7 +182,7 @@ VdwRadiusProvider::vdwRadiiForTopology(const gmx::TopologyInformation &top, std::vector mappedIds) { // get list of all atoms: - t_atoms atoms = top.topology() -> atoms; + t_atoms atoms = *(top.atoms()); // sanity check: int maxId = (*std::max_element(mappedIds.begin(), mappedIds.end())); diff --git a/src/trajectory-analysis/chap_trajectory_analysis.cpp b/src/trajectory-analysis/chap_trajectory_analysis.cpp index 94bf89f0..067ee3d2 100644 --- a/src/trajectory-analysis/chap_trajectory_analysis.cpp +++ b/src/trajectory-analysis/chap_trajectory_analysis.cpp @@ -25,8 +25,11 @@ #include #include +#include +#include +#include #include -#include +#include #include "trajectory-analysis/chap_trajectory_analysis.hpp" @@ -201,7 +204,7 @@ ChapTrajectoryAnalysis::initOptions( const char * const allowedPathFindingMethod[] = {"cylindrical", "inplane_optim"}; pfMethod_ = ePathFindingMethodInplaneOptimised; - options -> addOption(EnumOption("pf-method") + options -> addOption(LegacyEnumOption("pf-method") .enumValue(allowedPathFindingMethod) .store(&pfMethod_) .description("Path finding method. The default " @@ -221,7 +224,7 @@ ChapTrajectoryAnalysis::initOptions( "hole_xplor", "user"}; pfVdwRadiusDatabase_ = eVdwRadiusDatabaseHoleSimple; - options -> addOption(EnumOption("pf-vdwr-database") + options -> addOption(LegacyEnumOption("pf-vdwr-database") .enumValue(allowedVdwRadiusDatabase) .store(&pfVdwRadiusDatabase_) .description("Database of van-der-Waals radii to be " @@ -250,7 +253,7 @@ ChapTrajectoryAnalysis::initOptions( const char * const allowedPathAlignmentMethod[] = {"none", "ipp"}; pfPathAlignmentMethod_ = ePathAlignmentMethodIpp; - options -> addOption(EnumOption("pf-align-method") + options -> addOption(LegacyEnumOption("pf-align-method") .enumValue(allowedPathAlignmentMethod) .store(&pfPathAlignmentMethod_) .description("Method for aligning pathway " @@ -387,7 +390,7 @@ ChapTrajectoryAnalysis::initOptions( const char * const allowedDensityEstimationMethod[] = {"histogram", "kernel"}; deMethod_ = eDensityEstimatorKernel; - options -> addOption(EnumOption("de-method") + options -> addOption(LegacyEnumOption("de-method") .enumValue(allowedDensityEstimationMethod) .store(&deMethod_) .description("Method used for estimating the " @@ -444,7 +447,7 @@ ChapTrajectoryAnalysis::initOptions( "memprotmd", "user"}; hydrophobicityDatabase_ = eHydrophobicityDatabaseWimleyWhite1996; - options -> addOption(EnumOption("hydrophob-database") + options -> addOption(LegacyEnumOption("hydrophob-database") .enumValue(allowedHydrophobicityDatabase) .store(&hydrophobicityDatabase_) .description("Database of hydrophobicity scale for " @@ -482,19 +485,9 @@ ChapTrajectoryAnalysis::initAnalysis( const TrajectoryAnalysisSettings& /*settings*/, const TopologyInformation &top) { - // the following code ensures compatibility across Gromacs versions: - // (there was an API break between 2016 and 2018) - #if GROMACS_VERSION_MAJOR>=2018 - - std::unique_ptr topologyPointer(new gmx_mtop_t); - *topologyPointer = *top.mtop(); - - #elif GROMACS_VERSION_MAJOR>=2016 - - std::unique_ptr topologyPointer(new t_topology); - *topologyPointer = *top.topology(); - - #endif + // GROMACS 202x exposes gmx_mtop_t as non-copyable; use the topology owned by + // TopologyInformation and cast away const for legacy APIs that still require it. + gmx_mtop_t* topologyPointer = const_cast(top.mtop()); // set path name of NDX file: @@ -628,20 +621,20 @@ ChapTrajectoryAnalysis::initAnalysis( if( customNdxFileName_.size() != 0 ) { gmx_ana_indexgrps_init(&poreIdxGroups, - topologyPointer.get(), + topologyPointer, customNdxFileName_.c_str()); } else { gmx_ana_indexgrps_init(&poreIdxGroups, - topologyPointer.get(), + topologyPointer, NULL); } // create selections as defined above: poreMappingSelCal_ = poreMappingSelCol_.parseFromString(poreMappingSelCalString)[0]; poreMappingSelCog_ = poreMappingSelCol_.parseFromString(poreMappingSelCogString)[0]; - poreMappingSelCol_.setTopology(topologyPointer.get(), 0); + poreMappingSelCol_.setTopology(topologyPointer, 0); poreMappingSelCol_.setIndexGroups(poreIdxGroups); poreMappingSelCol_.compile(); @@ -677,13 +670,13 @@ ChapTrajectoryAnalysis::initAnalysis( if( customNdxFileName_.size() != 0 ) { gmx_ana_indexgrps_init(&solvIdxGroups, - topologyPointer.get(), + topologyPointer, customNdxFileName_.c_str()); } else { gmx_ana_indexgrps_init(&solvIdxGroups, - topologyPointer.get(), + topologyPointer, NULL); } @@ -694,7 +687,7 @@ ChapTrajectoryAnalysis::initAnalysis( solvMappingSelCog_ = solvMappingSelCol_.parseFromString(solvMappingSelCogString)[0]; // compile the selections: - solvMappingSelCol_.setTopology(topologyPointer.get(), 0); + solvMappingSelCol_.setTopology(topologyPointer, 0); solvMappingSelCol_.setIndexGroups(solvIdxGroups); solvMappingSelCol_.compile(); diff --git a/test/path-finding/ut_molecular_path.cpp b/test/path-finding/ut_molecular_path.cpp index ddea2c05..183371b1 100644 --- a/test/path-finding/ut_molecular_path.cpp +++ b/test/path-finding/ut_molecular_path.cpp @@ -24,7 +24,7 @@ #include -#include +#include #include "path-finding/molecular_path.hpp" From 5f99668335bb39f90403cf9b062caa5d58ae8805 Mon Sep 17 00:00:00 2001 From: RubenChM Date: Thu, 23 Apr 2026 12:03:07 +0200 Subject: [PATCH 5/6] Add scripts to build with conda enviroment in preparation for #18 --- .gitignore | 1 + CMakeLists.txt | 35 +++++++++++++++++++++++++++++++++++ README.md | 17 +++++++++++++---- conda/build.sh | 14 ++++++++++++++ conda/env.yml | 17 +++++++++++++++++ 5 files changed, 80 insertions(+), 4 deletions(-) create mode 100644 conda/build.sh create mode 100644 conda/env.yml diff --git a/.gitignore b/.gitignore index f0d1997a..81e92a7a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ # ignored subdirecotires: +.vscode/* build/* bin/* devdoc/* diff --git a/CMakeLists.txt b/CMakeLists.txt index 344934f3..a577388a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -85,6 +85,41 @@ endif() find_package(GROMACS 2016 REQUIRED) gromacs_check_double(GMX_DOUBLE) gromacs_check_compiler(CXX) + +# GROMACS Conda package exports carry a non-relocatable absolute libm path +# from the build environment. Drop that stale path and keep portable math linking. +if(TARGET Gromacs::libgromacs) + get_target_property(_gmx_iface_libs Gromacs::libgromacs INTERFACE_LINK_LIBRARIES) + if(_gmx_iface_libs) + set(_gmx_iface_libs_sanitized "") + set(_gmx_removed_abs_libm FALSE) + set(_gmx_has_portable_m FALSE) + foreach(_lib IN LISTS _gmx_iface_libs) + if(_lib MATCHES ".*/libm\\.so(\\.[0-9]+)*$") + set(_gmx_removed_abs_libm TRUE) + else() + if(_lib STREQUAL "m") + set(_gmx_has_portable_m TRUE) + endif() + list(APPEND _gmx_iface_libs_sanitized "${_lib}") + endif() + endforeach() + if(_gmx_removed_abs_libm) + if(NOT _gmx_has_portable_m) + list(APPEND _gmx_iface_libs_sanitized m) + endif() + set_target_properties(Gromacs::libgromacs PROPERTIES + INTERFACE_LINK_LIBRARIES "${_gmx_iface_libs_sanitized}" + ) + message(STATUS "Patched GROMACS imported target: replaced absolute libm.so dependency with portable math link") + endif() + unset(_gmx_iface_libs_sanitized) + unset(_gmx_removed_abs_libm) + unset(_gmx_has_portable_m) + endif() + unset(_gmx_iface_libs) +endif() + include_directories(${GROMACS_INCLUDE_DIRS}) add_definitions(${GROMACS_DEFINITIONS}) diff --git a/README.md b/README.md index 7a696474..183c31b3 100644 --- a/README.md +++ b/README.md @@ -6,8 +6,8 @@ CHAP is a tool for the functional annotation of ion channel structures written in C++. See the website under [www.channotation.org][CHANNOTATION] for a full documentation including installation instructions and usage examples. If you have any questions please use the GitHub [Issue Tracker](https://github.com/channotation/chap/issues). - -## Prerequisites ## +## Building from source ## +### Prerequisites ### Prior to installing CHAP, make sure that you have the following libraries and tools installed: @@ -25,7 +25,7 @@ You also need `-DGMX_INSTALL_LEGACY_API=ON` to get the required headers for CHAP CHAP also depends on [RapidJSON](http://rapidjson.org/), but this is included as a header-only library, and on [GTest][GTest], but this is downloaded and installed automatically by CMake, so you don't need to do anything about either of these (you will however need Internet access when installing CHAP). -## Installation ## +### Installation ### For a minimal install of CHAP create a `build` directory parallel to the source tree and from there run `cmake`, `make`, `make check`, and `make install`. @@ -58,7 +58,16 @@ which should bring up an online help for using CHAP. [GTest]: https://github.com/google/googletest [CHANNOTATION]: http://www.channotation.org +## Building with Conda ## + +If you have conda installed, you can create a conda environment with all dependencies and install chap using the following commands: + +```bash +conda env create -f conda/env.yml +conda activate chap-env +bash conda/build.sh +``` -## New ## +## Running with Singularity ## If you have ubuntu and you want to run a precompiled chap version with a singularity image, follow these instructions [here](https://github.com/channotation/chap/blob/singularity_branch/docs/_docs/getting-started/index.md) diff --git a/conda/build.sh b/conda/build.sh new file mode 100644 index 00000000..828c1c00 --- /dev/null +++ b/conda/build.sh @@ -0,0 +1,14 @@ +#!/bin/bash + +mkdir build +cd build +# Ensure installed binaries can find shared libs in conda-style prefixes without requiring LD_LIBRARY_PATH. +cmake -DCMAKE_INSTALL_PREFIX=${CONDA_PREFIX} \ + -DCMAKE_INSTALL_RPATH="\$ORIGIN/../lib;\$ORIGIN/../lib.AVX2_256;\$ORIGIN/../../lib;\$ORIGIN/../../lib.AVX2_256" \ + -DCMAKE_BUILD_WITH_INSTALL_RPATH=OFF \ + -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=OFF \ + .. +make +make check +make install +cp ${CONDA_PREFIX}/chap/bin/* ${CONDA_PREFIX}/bin diff --git a/conda/env.yml b/conda/env.yml new file mode 100644 index 00000000..dfd2980e --- /dev/null +++ b/conda/env.yml @@ -0,0 +1,17 @@ +name: chap-env +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - make + - gcc + - cmake >=3.5 + - boost-cpp + - libcblas * *_netlib # Only for build + - liblapacke * *_netlib # Only for build + - libtmglib # Only for build + - gromacs + - rapidjson + - intel-compute-runtime + From 353c300e1b0f1d1bc84ed2792bad2dd81c34469e Mon Sep 17 00:00:00 2001 From: RubenChM Date: Thu, 23 Apr 2026 13:29:44 +0200 Subject: [PATCH 6/6] Add legal documentation for vendored GROMACS headers --- include/gromacs/legal.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 include/gromacs/legal.md diff --git a/include/gromacs/legal.md b/include/gromacs/legal.md new file mode 100644 index 00000000..33534333 --- /dev/null +++ b/include/gromacs/legal.md @@ -0,0 +1,15 @@ +## Vendored GROMACS Headers ## + +CHAP currently vendors selected GROMACS-derived headers to allow its compilation with newer GROMACS versions (>=2018) that have changed their API. This includes are mostly reimplementations of random number generators of other libraries with additional upstream notices: + +1. `random/threefry.h` (includes notice for the Boost Software License) +2. `random/uniformrealdistribution.h` (includes notice for The LLVM compiler infrastructure MIT-licensed implementation) + +To stay compliant when vendoring these files, use the following checklist: + +1. Place vendored files in a dedicated folder, for example: `/third_party/gromacs/`. +2. Do not modify these files unless you are prepared to distribute those modifications under the GNU Lesser General Public License (LGPL) v2.1 or later, consistent with GROMACS licensing terms. +3. In your About/Legal material, include this exact statement: "This software uses the ThreeFry implementation from GROMACS, licensed under LGPL v2.1." +4. Provide a link to the GROMACS source code: https://github.com/gromacs/gromacs + +When redistributing, preserve all original license and attribution notices present in the vendored files, including the additional Boost and LLVM-related notices where applicable. \ No newline at end of file