diff --git a/doc/doxygen/Doxyfile.in b/doc/doxygen/Doxyfile.in index 05d551aa6ec..237c4221f18 100644 --- a/doc/doxygen/Doxyfile.in +++ b/doc/doxygen/Doxyfile.in @@ -1516,7 +1516,8 @@ INCLUDE_FILE_PATTERNS = *.h *.hpp *.cuh # Use the PREDEFINED tag if you want to use a different macro definition that # overrules the definition found in the source code. -EXPAND_AS_DEFINED = NEW_PARAMLESS_OBSERVABLE +EXPAND_AS_DEFINED = NEW_PARAMLESS_OBSERVABLE \ + NEW_THERMOSTAT # If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then # doxygen's preprocessor will remove all references to function-like macros diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index 01752fa3dc2..e28248d90c4 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -106,7 +106,7 @@ A code snippet would look like:: import espressomd system = espressomd.System() - system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0) + system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=42) system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) The physical meaning of these parameters is described below: @@ -171,6 +171,8 @@ The discretisation consists of the following steps (see :cite:`kolb99a` for a fu .. math:: \mathcal{P} = \mathcal{P}(x(t+dt),V(t+dt),f(x(t+dt)), v(t+dt/2)) .. math:: \Pi(t+dt) = \Pi(t+dt/2) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t+dt/2) dt/2 + \sqrt{k_B T \gamma^V dt} \overline{\eta} + with uncorrelated numbers :math:`\overline{\eta}` drawn from a random uniform process :math:`\eta(t)` + 6. Update the velocities .. math:: v(t+dt) = v(t+dt/2) + \frac{F(t+dt)}{m} dt/2 @@ -180,7 +182,7 @@ Notes: * The NpT algorithm is only tested for all 3 directions enabled for scaling. Usage of ``direction`` is considered an experimental feature. * In step 4, only those coordinates are scaled for which ``direction`` is set. * For the instantaneous pressure, the same limitations of applicability hold as described in :ref:`Pressure`. -* The particle forces :math:`F` include interactions as well as a friction and noise term analogous to the terms in the :ref:`Langevin thermostat`. +* The particle forces :math:`F` include interactions as well as a friction (:math:`\gamma^0`) and noise term (:math:`\sqrt{k_B T \gamma^0 dt} \overline{\eta}`) analogous to the terms in the :ref:`Langevin thermostat`. * The particle forces are only calculated in step 5 and then reused in step 1 of the next iteration. See :ref:`Velocity Verlet Algorithm` for the implications of that. .. _Rotational degrees of freedom and particle anisotropy: diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 85e71acce01..e1a6294bc9f 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -419,7 +419,7 @@ For example:: import espressomd system = espressomd.System() - system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0) + system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=41) system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) For an explanation of the algorithm involved, see :ref:`Isotropic NPT integrator`. diff --git a/samples/visualization_npt.py b/samples/visualization_npt.py index 52998f7ef9f..6edb69fb33c 100644 --- a/samples/visualization_npt.py +++ b/samples/visualization_npt.py @@ -57,7 +57,7 @@ max_displacement=0.1) print("E after minimization:", system.analysis.energy()["total"]) -system.thermostat.set_npt(kT=2.0, gamma0=1.0, gammav=0.01) +system.thermostat.set_npt(kT=2.0, gamma0=1.0, gammav=0.01, seed=42) system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=0.01) diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index c2f978d00e3..32cffb19622 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -320,11 +320,16 @@ void philox_counter_increment() { if (thermo_switch & THERMO_BROWNIAN) { brownian_rng_counter_increment(); } - if (thermo_switch & THERMO_DPD) { +#ifdef NPT + if (thermo_switch & THERMO_NPT_ISO) { + npt_iso_rng_counter_increment(); + } +#endif #ifdef DPD + if (thermo_switch & THERMO_DPD) { dpd_rng_counter_increment(); -#endif } +#endif if (n_thermalized_bonds) thermalized_bond_rng_counter_increment(); } diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index a416e5196b5..199ebbe82ef 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -40,15 +40,15 @@ void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { // Virtual sites are not propagated during integration if (p.p.is_virtual) continue; + auto const noise = friction_therm0_nptiso<2>(npt_iso, p.m.v, p.p.identity); for (int j = 0; j < 3; j++) { if (!(p.p.ext_flag & COORD_FIXED(j))) { if (nptiso.geometry & nptiso.nptgeom_dir[j]) { nptiso.p_vel[j] += Utils::sqr(p.m.v[j] * time_step) * p.p.mass; - p.m.v[j] += 0.5 * time_step / p.p.mass * p.f.f[j] + - friction_therm0_nptiso(npt_iso, p.m.v[j]) / p.p.mass; + p.m.v[j] += (p.f.f[j] * time_step / 2.0 + noise[j]) / p.p.mass; } else // Propagate velocity: v(t+dt) = v(t+0.5*dt) + 0.5*dt * a(t+dt) - p.m.v[j] += 0.5 * time_step * p.f.f[j] / p.p.mass; + p.m.v[j] += p.f.f[j] * time_step / 2.0 / p.p.mass; #ifdef EXTERNAL_FORCES } #endif @@ -171,15 +171,16 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) { for (int j = 0; j < 3; j++) { if (!(p.p.ext_flag & COORD_FIXED(j))) { #ifdef NPT + auto const noise = + friction_therm0_nptiso<1>(npt_iso, p.m.v, p.p.identity); if (integ_switch == INTEG_METHOD_NPT_ISO && (nptiso.geometry & nptiso.nptgeom_dir[j])) { - p.m.v[j] += p.f.f[j] * 0.5 * time_step / p.p.mass + - friction_therm0_nptiso(npt_iso, p.m.v[j]) / p.p.mass; + p.m.v[j] += (p.f.f[j] * time_step / 2.0 + noise[j]) / p.p.mass; nptiso.p_vel[j] += Utils::sqr(p.m.v[j] * time_step) * p.p.mass; } else #endif // Propagate velocities: v(t+0.5*dt) = v(t) + 0.5*dt * a(t) - p.m.v[j] += 0.5 * time_step * p.f.f[j] / p.p.mass; + p.m.v[j] += p.f.f[j] * time_step / 2.0 / p.p.mass; } } } diff --git a/src/core/random.hpp b/src/core/random.hpp index 08f02aad059..955358c0e68 100644 --- a/src/core/random.hpp +++ b/src/core/random.hpp @@ -55,6 +55,9 @@ enum class RNGSalt : uint64_t { BROWNIAN_INC, BROWNIAN_ROT_INC, BROWNIAN_ROT_WALK, + NPTISO0_HALF_STEP1, + NPTISO0_HALF_STEP2, + NPTISOV, SALT_DPD, THERMALIZED_BOND }; @@ -88,9 +91,29 @@ Utils::Vector philox_4_uint64s(uint64_t counter, int key1, return {res[0], res[1], res[2], res[3]}; } +/** + * @brief Uniform noise. + * + * Mean = 0, variance = 1 / 12. + * This uses the Philox PRNG, the state is controlled + * by the counter, the salt and two keys. + * If any of the keys and salt differ, the noise is + * not correlated between two calls along the same counter + * sequence. + * + */ +template double noise(uint64_t counter, int key1, int key2 = 0) { + + auto const noise = philox_4_uint64s(counter, key1, key2); + + using Utils::uniform; + return uniform(noise[0]) - 0.5; +} + /** * @brief 3d uniform vector noise. * + * Mean = 0, variance = 1 / 12. * This uses the Philox PRNG, the state is controlled * by the counter, the salt and two keys. * If any of the keys and salt differ, the noise is @@ -111,7 +134,7 @@ Utils::Vector3d v_noise(uint64_t counter, int key1, int key2 = 0) { /** @brief Generator for Gaussian random 3d vector. * - * Mean = 0, standard deviation = 1.0 + * Mean = 0, standard deviation = 1.0. * Based on the Philox RNG using 4x64 bits. * The Box-Muller transform is used to convert from uniform to normal * distribution. The transform is only valid, if the uniformly distributed @@ -211,8 +234,8 @@ void init_random_seed(int seed); } // namespace Random /** - * @brief Draws a random real number from the uniform distribution in the range - * [0,1) + * @brief Draws a random real number from the uniform distribution in the + * range [0,1) using the Mersenne twister. */ inline double d_random() { using namespace Random; diff --git a/src/core/thermostat.cpp b/src/core/thermostat.cpp index 84e0402b88d..e2e77615f37 100644 --- a/src/core/thermostat.cpp +++ b/src/core/thermostat.cpp @@ -38,63 +38,51 @@ bool thermo_virtual = true; using Thermostat::GammaType; +/** + * @brief Register a thermostat MPI callbacks + * + * @param thermostat The thermostat global variable + * @param thermostat_enum The thermostat enum value + */ +#define REGISTER_THERMOSTAT_CALLBACKS(thermostat, thermostat_enum) \ + void mpi_bcast_##thermostat##_rng_counter_slave(const uint64_t counter) { \ + thermostat.rng_counter = \ + std::make_unique>(counter); \ + } \ + \ + REGISTER_CALLBACK(mpi_bcast_##thermostat##_rng_counter_slave) \ + \ + void mpi_bcast_##thermostat##_rng_counter(const uint64_t counter) { \ + mpi_call(mpi_bcast_##thermostat##_rng_counter_slave, counter); \ + } \ + \ + void thermostat##_rng_counter_increment() { \ + if (thermo_switch & thermostat_enum) \ + thermostat.rng_counter->increment(); \ + } \ + \ + bool thermostat##_is_seed_required() { \ + /* Seed is required if rng is not initialized */ \ + return thermostat.rng_counter == nullptr; \ + } \ + \ + void thermostat##_set_rng_state(const uint64_t counter) { \ + mpi_bcast_##thermostat##_rng_counter(counter); \ + thermostat.rng_counter = \ + std::make_unique>(counter); \ + } \ + \ + uint64_t thermostat##_get_rng_state() { \ + return thermostat.rng_counter->value(); \ + } + LangevinThermostat langevin = {}; BrownianThermostat brownian = {}; IsotropicNptThermostat npt_iso = {}; -void mpi_bcast_langevin_rng_counter_slave(const uint64_t counter) { - langevin.rng_counter = std::make_unique>(counter); -} - -REGISTER_CALLBACK(mpi_bcast_langevin_rng_counter_slave) - -void mpi_bcast_brownian_rng_counter_slave(const uint64_t counter) { - brownian.rng_counter = std::make_unique>(counter); -} - -REGISTER_CALLBACK(mpi_bcast_brownian_rng_counter_slave) - -void mpi_bcast_langevin_rng_counter(const uint64_t counter) { - mpi_call(mpi_bcast_langevin_rng_counter_slave, counter); -} - -void mpi_bcast_brownian_rng_counter(const uint64_t counter) { - mpi_call(mpi_bcast_brownian_rng_counter_slave, counter); -} - -void langevin_rng_counter_increment() { - if (thermo_switch & THERMO_LANGEVIN) - langevin.rng_counter->increment(); -} - -void brownian_rng_counter_increment() { - if (thermo_switch & THERMO_BROWNIAN) - brownian.rng_counter->increment(); -} - -bool langevin_is_seed_required() { - /* Seed is required if rng is not initialized */ - return langevin.rng_counter == nullptr; -} - -bool brownian_is_seed_required() { - /* Seed is required if rng is not initialized */ - return brownian.rng_counter == nullptr; -} - -void langevin_set_rng_state(const uint64_t counter) { - mpi_bcast_langevin_rng_counter(counter); - langevin.rng_counter = std::make_unique>(counter); -} - -void brownian_set_rng_state(const uint64_t counter) { - mpi_bcast_brownian_rng_counter(counter); - brownian.rng_counter = std::make_unique>(counter); -} - -uint64_t langevin_get_rng_state() { return langevin.rng_counter->value(); } - -uint64_t brownian_get_rng_state() { return brownian.rng_counter->value(); } +REGISTER_THERMOSTAT_CALLBACKS(langevin, THERMO_LANGEVIN) +REGISTER_THERMOSTAT_CALLBACKS(brownian, THERMO_BROWNIAN) +REGISTER_THERMOSTAT_CALLBACKS(npt_iso, THERMO_NPT_ISO) void thermo_init() { // Init thermalized bond despite of thermostat diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 1a051f71931..5b73a40184a 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -49,9 +49,6 @@ /*@}*/ namespace Thermostat { - -static auto noise = []() { return (d_random() - 0.5); }; - #ifdef PARTICLE_ANISOTROPY using GammaType = Utils::Vector3d; #else @@ -94,8 +91,13 @@ extern bool thermo_virtual; * parameter structs ************************************************/ +struct BaseThermostat { + /** RNG counter. */ + std::unique_ptr> rng_counter; +}; + /** %Thermostat for Langevin dynamics. */ -struct LangevinThermostat { +struct LangevinThermostat : public BaseThermostat { private: using GammaType = Thermostat::GammaType; @@ -112,9 +114,14 @@ struct LangevinThermostat { } pref_noise_rotation = sigma(temperature, time_step, gamma_rotation); } - /** Calculate the noise standard deviation. */ + /** Calculate the noise prefactor. + * Evaluates the quantity @f$ \sqrt{2 k_B T \gamma / dt} / \sigma_\eta @f$ + * with @f$ \sigma_\eta @f$ the standard deviation of the random uniform + * process @f$ \eta(t) @f$. + */ static GammaType sigma(double kT, double time_step, GammaType const &gamma) { - constexpr auto const temp_coeff = 24.0; + // random uniform noise has variance 1/12 + constexpr auto const temp_coeff = 2.0 * 12.0; return sqrt((temp_coeff * kT / time_step) * gamma); } /** @name Parameters */ @@ -126,21 +133,25 @@ struct LangevinThermostat { /*@}*/ /** @name Prefactors */ /*@{*/ - /** Prefactor for the friction. */ + /** Prefactor for the friction. + * Stores @f$ \gamma_{\text{trans}} @f$. + */ GammaType pref_friction; - /** Prefactor for the translational velocity noise. */ + /** Prefactor for the translational velocity noise. + * Stores @f$ \sqrt{2 k_B T \gamma_{\text{trans}} / dt} / \sigma_\eta @f$. + */ GammaType pref_noise; - /** Prefactor for the angular velocity noise. */ + /** Prefactor for the angular velocity noise. + * Stores @f$ \sqrt{2 k_B T \gamma_{\text{rot}} / dt} / \sigma_\eta @f$. + */ GammaType pref_noise_rotation; /*@}*/ - /** RNG counter, used for both translation and rotation. */ - std::unique_ptr> rng_counter; }; /** %Thermostat for Brownian dynamics. * Default particle mass is assumed to be unitary in these global parameters. */ -struct BrownianThermostat { +struct BrownianThermostat : public BaseThermostat { private: using GammaType = Thermostat::GammaType; @@ -172,12 +183,20 @@ struct BrownianThermostat { sigma_pos_rotation = sigma(temperature, gamma_rotation); #endif // ROTATION } - /** Calculate the noise standard deviation. */ + /** Calculate the noise prefactor. + * Evaluates the quantity @f$ \sqrt{2 k_B T / \gamma} / \sigma_\eta @f$ + * with @f$ \sigma_\eta @f$ the standard deviation of the random gaussian + * process @f$ \eta(t) @f$. + */ static GammaType sigma(double kT, GammaType const &gamma) { constexpr auto const temp_coeff = 2.0; return sqrt(Utils::hadamard_division(temp_coeff * kT, gamma)); } - /** Calculate the noise standard deviation. */ + /** Calculate the noise prefactor. + * Evaluates the quantity @f$ \sqrt{k_B T} / \sigma_\eta @f$ + * with @f$ \sigma_\eta @f$ the standard deviation of the random gaussian + * process @f$ \eta(t) @f$. + */ static double sigma(double kT) { constexpr auto const temp_coeff = 1.0; return sqrt(temp_coeff * kT); @@ -212,42 +231,64 @@ struct BrownianThermostat { */ double sigma_vel_rotation = 0; /*@}*/ - /** RNG counter, used for both translation and rotation. */ - std::unique_ptr> rng_counter; }; /** %Thermostat for isotropic NPT dynamics. */ -struct IsotropicNptThermostat { +struct IsotropicNptThermostat : public BaseThermostat { private: using GammaType = Thermostat::GammaType; public: -#ifdef NPT /** Recalculate prefactors. * Needs to be called every time the parameters are changed. */ void recalc_prefactors(double piston) { +#ifdef NPT assert(piston > 0.0); - pref1 = -gamma0 * 0.5 * time_step; - pref2 = sqrt(12.0 * temperature * gamma0 * time_step); - pref3 = -gammav * (1.0 / piston) * 0.5 * time_step; - pref4 = sqrt(12.0 * temperature * gammav * time_step); - } + auto const half_time_step = time_step / 2.0; + pref_rescale_0 = -gamma0 * half_time_step; + pref_noise_0 = sigma(temperature, gamma0); + pref_rescale_V = -gammav * half_time_step / piston; + pref_noise_V = sigma(temperature, gammav); #endif + } + /** Calculate the noise prefactor. + * Evaluates the quantity @f$ \sqrt{2 k_B T \gamma dt / 2} / \sigma_\eta @f$ + * with @f$ \sigma_\eta @f$ the standard deviation of the random uniform + * process @f$ \eta(t) @f$. + */ + static double sigma(double kT, double gamma) { + // random uniform noise has variance 1/12; the temperature + // coefficient of 2 is canceled out by the half time step + constexpr auto const temp_coeff = 12.0; + return sqrt(temp_coeff * temperature * gamma * time_step); + } /** @name Parameters */ /*@{*/ - /** Friction coefficient @f$ \gamma_0 @f$ */ + /** Friction coefficient of the particles @f$ \gamma^0 @f$ */ double gamma0; - /** Friction coefficient @f$ \gamma_V @f$ */ + /** Friction coefficient for the box @f$ \gamma^V @f$ */ double gammav; /*@}*/ #ifdef NPT /** @name Prefactors */ /*@{*/ - double pref1; - double pref2; - double pref3; - double pref4; + /** %Particle velocity rescaling at half the time step. + * Stores @f$ \gamma^{0}\cdot\frac{dt}{2} @f$. + */ + double pref_rescale_0; + /** %Particle velocity rescaling noise standard deviation. + * Stores @f$ \sqrt{k_B T \gamma^{0} dt} / \sigma_\eta @f$. + */ + double pref_noise_0; + /** Volume rescaling at half the time step. + * Stores @f$ \frac{\gamma^{V}}{Q}\cdot\frac{dt}{2} @f$. + */ + double pref_rescale_V; + /** Volume rescaling noise standard deviation. + * Stores @f$ \sqrt{k_B T \gamma^{V} dt} / \sigma_\eta @f$. + */ + double pref_noise_V; /*@}*/ #endif }; @@ -256,21 +297,20 @@ struct IsotropicNptThermostat { * functions ************************************************/ -/** Only require seed if rng is not initialized. */ -bool langevin_is_seed_required(); +/** + * @brief Register a thermostat public interface + * + * @param thermostat The thermostat name + */ +#define NEW_THERMOSTAT(thermostat) \ + bool thermostat##_is_seed_required(); \ + void thermostat##_rng_counter_increment(); \ + void thermostat##_set_rng_state(uint64_t counter); \ + uint64_t thermostat##_get_rng_state(); -/** Only require seed if rng is not initialized. */ -bool brownian_is_seed_required(); - -/** @name philox functionality: increment, get/set */ -/*@{*/ -void langevin_rng_counter_increment(); -void langevin_set_rng_state(uint64_t counter); -uint64_t langevin_get_rng_state(); -void brownian_rng_counter_increment(); -void brownian_set_rng_state(uint64_t counter); -uint64_t brownian_get_rng_state(); -/*@}*/ +NEW_THERMOSTAT(langevin) +NEW_THERMOSTAT(brownian) +NEW_THERMOSTAT(npt_iso) /** Initialize constants of the thermostat at the start of integration */ void thermo_init(); @@ -278,20 +318,30 @@ void thermo_init(); #ifdef NPT /** Add velocity-dependent noise and friction for NpT-sims to the particle's * velocity - * @param npt_iso Parameters - * @param vj j-component of the velocity - * @return j-component of the noise added to the velocity, also scaled by - * dt (contained in prefactors) + * @tparam step Which half time step to integrate (1 or 2) + * @param npt_iso Parameters + * @param vel particle velocity + * @param p_identity particle identity + * @return noise added to the velocity, already rescaled by + * dt/2 (contained in prefactors) */ -inline double friction_therm0_nptiso(IsotropicNptThermostat const &npt_iso, - double vj) { +template +inline Utils::Vector3d +friction_therm0_nptiso(IsotropicNptThermostat const &npt_iso, + Utils::Vector3d const &vel, int p_identity) { + static_assert(step == 1 or step == 2, "NPT only has 2 integration steps"); + constexpr auto const salt = + (step == 1) ? RNGSalt::NPTISO0_HALF_STEP1 : RNGSalt::NPTISO0_HALF_STEP2; if (thermo_switch & THERMO_NPT_ISO) { - if (npt_iso.pref2 > 0.0) { - return npt_iso.pref1 * vj + npt_iso.pref2 * Thermostat::noise(); + if (npt_iso.pref_noise_0 > 0.0) { + return npt_iso.pref_rescale_0 * vel + + npt_iso.pref_noise_0 * + Random::v_noise(npt_iso.rng_counter->value(), + p_identity); } - return npt_iso.pref1 * vj; + return npt_iso.pref_rescale_0 * vel; } - return 0.0; + return {}; } /** Add p_diff-dependent noise and friction for NpT-sims to \ref @@ -300,13 +350,15 @@ inline double friction_therm0_nptiso(IsotropicNptThermostat const &npt_iso, inline double friction_thermV_nptiso(IsotropicNptThermostat const &npt_iso, double p_diff) { if (thermo_switch & THERMO_NPT_ISO) { - if (npt_iso.pref4 > 0.0) { - return npt_iso.pref3 * p_diff + npt_iso.pref4 * Thermostat::noise(); + if (npt_iso.pref_noise_V > 0.0) { + return npt_iso.pref_rescale_V * p_diff + + npt_iso.pref_noise_V * Random::noise( + npt_iso.rng_counter->value(), 0); } - return npt_iso.pref3 * p_diff; + return npt_iso.pref_rescale_V * p_diff; } return 0.0; } -#endif +#endif // NPT #endif diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index 6b4391b166b..ab3473bad1b 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -47,4 +47,5 @@ unit_test(NAME grid_test SRC grid_test.cpp DEPENDS EspressoCore) unit_test(NAME BoxGeometry_test SRC BoxGeometry_test.cpp DEPENDS EspressoCore) unit_test(NAME LocalBox_test SRC LocalBox_test.cpp DEPENDS EspressoCore) unit_test(NAME thermostats_test SRC thermostats_test.cpp DEPENDS EspressoCore) +unit_test(NAME random_test SRC random_test.cpp DEPENDS EspressoCore) diff --git a/src/core/unit_tests/random_test.cpp b/src/core/unit_tests/random_test.cpp new file mode 100644 index 00000000000..10920f4e26a --- /dev/null +++ b/src/core/unit_tests/random_test.cpp @@ -0,0 +1,133 @@ +/* + * Copyright (C) 2019 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +/* Unit tests for random number generators. */ + +#define BOOST_TEST_MODULE PRNG test +#define BOOST_TEST_DYN_LINK +#include +#include + +#include + +#include "random.hpp" +#include "random_test.hpp" + +BOOST_AUTO_TEST_CASE(test_noise_statistics) { + constexpr size_t const sample_size = 100'000; + constexpr size_t const x = 0, y = 1, z = 2; + constexpr double const tol = 1e-12; + + double value = 1; + std::vector means, variances; + std::vector> covariance; + std::vector> correlation; + std::tie(means, variances, covariance, correlation) = + noise_statistics(std::function()>( + [&value]() -> std::vector { + value *= -1; + return {{Utils::Vector2d{value, -value}}}; + }), + sample_size); + // check pooled mean and variance + BOOST_CHECK_SMALL(std::abs(means[0]), 100 * tol); + BOOST_CHECK_CLOSE(variances[0], 1.0, tol); + // check variance per axis + BOOST_CHECK_CLOSE(covariance[x][x], 1.0, tol); + BOOST_CHECK_CLOSE(covariance[y][y], 1.0, tol); + BOOST_CHECK_CLOSE(covariance[x][y], -1.0, tol); + BOOST_CHECK_EQUAL(covariance[x][y], covariance[y][x]); + // check correlation + BOOST_CHECK_CLOSE(correlation[x][x], 1.0, tol); + BOOST_CHECK_CLOSE(correlation[y][y], 1.0, tol); + BOOST_CHECK_CLOSE(correlation[x][y], -1.0, tol); + BOOST_CHECK_EQUAL(correlation[x][y], correlation[y][x]); +} + +BOOST_AUTO_TEST_CASE(test_noise_uniform_1d) { + constexpr size_t const sample_size = 4'000'000; + + int counter = 0; + std::vector means, variances; + std::vector> covariance; + std::vector> correlation; + std::tie(means, variances, covariance, correlation) = noise_statistics( + std::function()>( + [&counter]() -> std::vector { + return {{Random::noise(counter++, 0)}}; + }), + sample_size); + // check pooled mean and variance + BOOST_CHECK_SMALL(std::abs(means[0]), 2e-4); + BOOST_CHECK_CLOSE(variances[0] * 12.0, 1.0, 0.05); +} + +BOOST_AUTO_TEST_CASE(test_noise_uniform_3d) { + constexpr size_t const sample_size = 4'000'000; + constexpr size_t const x = 0, y = 1, z = 2; + + int counter = 0; + std::vector means, variances; + std::vector> covariance; + std::vector> correlation; + std::tie(means, variances, covariance, correlation) = noise_statistics( + std::function()>( + [&counter]() -> std::vector { + return {{Random::v_noise(counter++, 0)}}; + }), + sample_size); + // check pooled mean and variance + BOOST_CHECK_SMALL(std::abs(means[0]), 2e-4); + BOOST_CHECK_CLOSE(variances[0] * 12.0, 1.0, 0.05); + // check variance per axis + BOOST_CHECK_CLOSE(covariance[x][x] * 12.0, 1.0, 0.2); + BOOST_CHECK_CLOSE(covariance[y][y] * 12.0, 1.0, 0.2); + BOOST_CHECK_CLOSE(covariance[z][z] * 12.0, 1.0, 0.2); + // check correlation + BOOST_CHECK_SMALL(std::abs(correlation[x][y]), 2e-3); + BOOST_CHECK_SMALL(std::abs(correlation[y][z]), 2e-3); + BOOST_CHECK_SMALL(std::abs(correlation[z][x]), 2e-3); +} + +BOOST_AUTO_TEST_CASE(test_noise_gaussian_3d) { + constexpr size_t const sample_size = 4'000'000; + constexpr size_t const x = 0, y = 1, z = 2; + + int counter = 0; + std::vector means, variances; + std::vector> covariance; + std::vector> correlation; + std::tie(means, variances, covariance, correlation) = noise_statistics( + std::function()>( + [&counter]() -> std::vector { + return {{Random::v_noise_g(counter++, 0)}}; + }), + sample_size); + // check pooled mean and variance + BOOST_CHECK_SMALL(std::abs(means[0]), 2e-4); + BOOST_CHECK_CLOSE(variances[0], 1.0, 0.05); + // check variance per axis + BOOST_CHECK_CLOSE(covariance[x][x], 1.0, 0.2); + BOOST_CHECK_CLOSE(covariance[y][y], 1.0, 0.2); + BOOST_CHECK_CLOSE(covariance[z][z], 1.0, 0.2); + // check correlation + BOOST_CHECK_SMALL(std::abs(correlation[x][y]), 2e-3); + BOOST_CHECK_SMALL(std::abs(correlation[y][z]), 2e-3); + BOOST_CHECK_SMALL(std::abs(correlation[z][x]), 2e-3); +} diff --git a/src/core/unit_tests/random_test.hpp b/src/core/unit_tests/random_test.hpp new file mode 100644 index 00000000000..e456fdcdab0 --- /dev/null +++ b/src/core/unit_tests/random_test.hpp @@ -0,0 +1,179 @@ +/* + * Copyright (C) 2019 The ESPResSo project + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef CORE_UNIT_TESTS_RANDOM_TEST_HPP +#define CORE_UNIT_TESTS_RANDOM_TEST_HPP + +/* Helper functions to compute random numbers covariance in a single pass */ + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "random.hpp" + +namespace Utils { +using VariantVectorXd = boost::variant; +} // namespace Utils + +using Utils::VariantVectorXd; + +namespace { + +using Utils::Vector; + +class visitor_size : public boost::static_visitor { +public: + template size_t operator()(Vector const &v) const { + return v.size(); + } + size_t operator()(double v) const { return 1; } +}; + +class visitor_get : public boost::static_visitor { +public: + template + double operator()(Vector const &v, size_t i) const { + return v[i]; + } + double operator()(double v, size_t i) const { + assert(i == 0); + return v; + } +}; + +size_t get_size(VariantVectorXd const &vec) { + return boost::apply_visitor(visitor_size(), vec); +} + +double get_value(VariantVectorXd const &vec, size_t i) { + return boost::apply_visitor( + std::bind(visitor_get(), std::placeholders::_1, i), vec); +} + +template auto square_matrix(size_t N) { + return std::vector>(N, std::vector(N)); +} + +} // namespace + +/** Draw a large sample of 3D vectors from PRNGs and compute statistics. + * Parameter @p noise_function is a generator that returns @f N @f vectors + * of size @f M_i @f. The following statistics are evaluated: @f N @f means + * and @f N @f variances (samples are uncorrelated across axes, so pooling + * them is fine), and a covariance and a correlation matrix of size + * @f \sum M_i @f. + */ +std::tuple, std::vector, + std::vector>, std::vector>> +noise_statistics(std::function()> noise_function, + size_t sample_size) { + + // get size of the arrays and size of the triangular correlation matrix + auto const first_value = noise_function(); + auto const n_vectors = first_value.size(); + std::vector dimensions(n_vectors); + std::transform(first_value.begin(), first_value.end(), dimensions.begin(), + [](auto const &element) { return get_size(element); }); + auto const matrix_dim = std::accumulate(dimensions.begin(), dimensions.end(), + 0, std::plus()); + + // set up boost accumulators + namespace ba = boost::accumulators; + namespace bt = boost::accumulators::tag; + using stat_variance = ba::stats; + using stat_covariance = ba::stats>; + using boost_variance = ba::accumulator_set; + using boost_covariance = ba::accumulator_set; + std::vector acc_variance(n_vectors); + auto acc_covariance = ::square_matrix(matrix_dim); + + // accumulate + for (size_t step = 0; step < sample_size; ++step) { + auto const noise_tuple = noise_function(); + // for each vector, pool the random numbers of all columns + for (size_t vec1 = 0; vec1 < dimensions.size(); ++vec1) { + for (size_t col1 = 0; col1 < dimensions[vec1]; ++col1) { + acc_variance[vec1](::get_value(noise_tuple[vec1], col1)); + } + } + // fill the covariance matrix (upper triangle) + size_t index1 = 0; + for (size_t vec1 = 0; vec1 < dimensions.size(); ++vec1) { + for (size_t col1 = 0; col1 < dimensions[vec1]; ++col1) { + size_t index2 = index1; + for (size_t vec2 = vec1; vec2 < dimensions.size(); ++vec2) { + for (size_t col2 = (vec2 == vec1) ? col1 : 0; col2 < dimensions[vec2]; + ++col2) { + acc_covariance[index1][index2]( + ::get_value(noise_tuple[vec1], col1), + ba::covariate1 = ::get_value(noise_tuple[vec2], col2)); + index2++; + } + } + index1++; + } + } + } + + // compute statistics + std::vector means(n_vectors); + std::vector variances(n_vectors); + for (size_t i = 0; i < n_vectors; ++i) { + means[i] = ba::mean(acc_variance[i]); + variances[i] = ba::variance(acc_variance[i]); + } + auto covariance = ::square_matrix(matrix_dim); + for (size_t i = 0; i < matrix_dim; ++i) { + for (size_t j = i; j < matrix_dim; ++j) { + covariance[i][j] = covariance[j][i] = + ba::covariance(acc_covariance[i][j]); + } + } + auto correlation = ::square_matrix(matrix_dim); + for (size_t i = 0; i < matrix_dim; ++i) { + for (size_t j = i; j < matrix_dim; ++j) { + correlation[i][j] = correlation[j][i] = + covariance[i][j] / sqrt(covariance[i][i] * covariance[j][j]); + } + } + + return std::make_tuple(means, variances, covariance, correlation); +} + +boost::test_tools::predicate_result correlation_almost_equal( + std::vector> const &correlation_matrix, size_t i, + size_t j, double reference, double threshold) { + auto const value = correlation_matrix[i][j]; + auto const diff = std::abs(value - reference); + if (diff > threshold) { + boost::test_tools::predicate_result res(false); + res.message() << "The correlation coefficient M[" << i << "][" << j << "]{" + << value << "} differs from " << reference << " by " << diff + << " (> " << threshold << ")"; + return res; + } + return true; +} +#endif diff --git a/src/core/unit_tests/thermostats_test.cpp b/src/core/unit_tests/thermostats_test.cpp index d66206dc290..55d4e3a3c9c 100644 --- a/src/core/unit_tests/thermostats_test.cpp +++ b/src/core/unit_tests/thermostats_test.cpp @@ -19,12 +19,13 @@ /* Unit tests for thermostats. */ -#define BOOST_TEST_MODULE Particle test +#define BOOST_TEST_MODULE Thermostats test #define BOOST_TEST_DYN_LINK #include #include #include +#include #include #include #include @@ -32,6 +33,7 @@ #include "Particle.hpp" #include "integrators/brownian_inline.hpp" #include "integrators/langevin_inline.hpp" +#include "random_test.hpp" #include "thermostat.hpp" extern double time_step; @@ -209,3 +211,128 @@ BOOST_AUTO_TEST_CASE(test_langevin_dynamics) { } #endif // ROTATION } + +BOOST_AUTO_TEST_CASE(test_noise_statistics) { + time_step = 1.0; + temperature = 2.0; + constexpr size_t const sample_size = 1'000'000; + auto thermostat = thermostat_factory(); + thermostat.rng_counter = std::make_unique>(0); + auto p1 = particle_factory(); + auto p2 = particle_factory(); + p1.p.identity = 0; + p1.p.identity = 1; + + auto const correlation = std::get<3>(noise_statistics( + std::function()>( + [&p1, &p2, &thermostat]() -> std::vector { + thermostat.rng_counter->increment(); + return {{friction_thermo_langevin(thermostat, p1), + -friction_thermo_langevin(thermostat, p1), + friction_thermo_langevin(thermostat, p2)}}; + }), + sample_size)); + for (size_t i = 0; i < correlation.size(); ++i) { + for (size_t j = i; j < correlation.size(); ++j) { + double expected; + if (i == j) { + expected = 1.0; + } else if (i < 3 and j == i + 3) { + expected = -1.0; + } else { + expected = 0.0; + } + BOOST_CHECK(correlation_almost_equal(correlation, i, j, expected, 2e-3)); + } + } +} + +BOOST_AUTO_TEST_CASE(test_brownian_randomness) { + time_step = 1.0; + temperature = 2.0; + constexpr size_t const sample_size = 500'000; + auto const thermostat = thermostat_factory(); + auto p = particle_factory(); +#ifdef ROTATION + p.p.rotation = ROTATION_X | ROTATION_Y | ROTATION_Z; +#endif + + auto const correlation = std::get<3>( + noise_statistics(std::function()>( + [&p, &thermostat]() -> std::vector { + thermostat.rng_counter->increment(); + return {{ + bd_random_walk(thermostat, p, time_step), + bd_random_walk_vel(thermostat, p), +#ifdef ROTATION + bd_random_walk_rot(thermostat, p, time_step), + bd_random_walk_vel_rot(thermostat, p), +#endif + }}; + }), + sample_size)); + for (size_t i = 0; i < correlation.size(); ++i) { + for (size_t j = i + 1; j < correlation.size(); ++j) { + BOOST_CHECK(correlation_almost_equal(correlation, i, j, 0.0, 5e-2)); + } + } +} + +BOOST_AUTO_TEST_CASE(test_langevin_randomness) { + time_step = 1.0; + temperature = 2.0; + constexpr size_t const sample_size = 500'000; + auto const thermostat = thermostat_factory(); + auto p = particle_factory(); + + auto const correlation = std::get<3>(noise_statistics( + std::function()>( + [&p, &thermostat]() -> std::vector { + thermostat.rng_counter->increment(); + return {{ + friction_thermo_langevin(thermostat, p), +#ifdef ROTATION + friction_thermo_langevin_rotation(thermostat, p), +#endif + }}; + }), + sample_size)); + for (size_t i = 0; i < correlation.size(); ++i) { + for (size_t j = i + 1; j < correlation.size(); ++j) { + BOOST_CHECK(correlation_almost_equal(correlation, i, j, 0.0, 5e-2)); + } + } +} + +#ifdef NPT +BOOST_AUTO_TEST_CASE(test_npt_iso_randomness) { + extern int thermo_switch; + thermo_switch |= THERMO_NPT_ISO; + time_step = 1.0; + temperature = 2.0; + constexpr size_t const sample_size = 500'000; + IsotropicNptThermostat thermostat{}; + thermostat.rng_counter = std::make_unique>(0); + thermostat.gamma0 = 2.0; + thermostat.gammav = 0.1; + thermostat.recalc_prefactors(1.0); + auto p = particle_factory(); + + auto const correlation = std::get<3>(noise_statistics( + std::function()>( + [&p, &thermostat]() -> std::vector { + thermostat.rng_counter->increment(); + return {{ + friction_therm0_nptiso<1>(thermostat, p.m.v, 0), + friction_therm0_nptiso<2>(thermostat, p.m.v, 0), + friction_thermV_nptiso(thermostat, 1.5), + }}; + }), + sample_size)); + for (size_t i = 0; i < correlation.size(); ++i) { + for (size_t j = i + 1; j < correlation.size(); ++j) { + BOOST_CHECK(correlation_almost_equal(correlation, i, j, 0.0, 5e-2)); + } + } +} +#endif // NPT diff --git a/src/python/espressomd/thermostat.pxd b/src/python/espressomd/thermostat.pxd index cf9ef7068eb..b7f2f08fccb 100644 --- a/src/python/espressomd/thermostat.pxd +++ b/src/python/espressomd/thermostat.pxd @@ -53,11 +53,14 @@ cdef extern from "thermostat.hpp": void langevin_set_rng_state(stdint.uint64_t counter) void brownian_set_rng_state(stdint.uint64_t counter) + void npt_iso_set_rng_state(stdint.uint64_t counter) cbool langevin_is_seed_required() cbool brownian_is_seed_required() + cbool npt_iso_is_seed_required() stdint.uint64_t langevin_get_rng_state() stdint.uint64_t brownian_get_rng_state() + stdint.uint64_t npt_iso_get_rng_state() cdef extern from "Globals.hpp": # links intern C-struct with python object diff --git a/src/python/espressomd/thermostat.pyx b/src/python/espressomd/thermostat.pyx index ae1de17df7b..0f0568ca4b2 100644 --- a/src/python/espressomd/thermostat.pyx +++ b/src/python/espressomd/thermostat.pyx @@ -108,7 +108,7 @@ cdef class Thermostat: seed=thmst["rng_counter_fluid"]) if thmst["type"] == "NPT_ISO": self.set_npt(kT=thmst["kT"], gamma0=thmst["gamma0"], - gammav=thmst["gammav"]) + gammav=thmst["gammav"], seed=thmst["seed"]) if thmst["type"] == "DPD": self.set_dpd(kT=thmst["kT"], seed=thmst["seed"]) if thmst["type"] == "BROWNIAN": @@ -185,6 +185,7 @@ cdef class Thermostat: npt_dict = {} npt_dict["type"] = "NPT_ISO" npt_dict["kT"] = temperature + npt_dict["seed"] = int(npt_iso_get_rng_state()) npt_dict["gamma0"] = npt_iso.gamma0 npt_dict["gammav"] = npt_iso.gammav npt_dict.update(nptiso) @@ -600,7 +601,7 @@ cdef class Thermostat: IF NPT: @AssertThermostatType(THERMO_NPT_ISO) - def set_npt(self, kT=None, gamma0=None, gammav=None): + def set_npt(self, kT=None, gamma0=None, gammav=None, seed=None): """ Sets the NPT thermostat. @@ -612,6 +613,10 @@ cdef class Thermostat: Friction coefficient of the bath gammav : :obj:`float` Artificial friction coefficient for the volume fluctuations. + seed : :obj:`int` + Initial counter value (or seed) of the philox RNG. + Required on first activation of the Langevin thermostat. + Must be positive. """ @@ -620,6 +625,19 @@ cdef class Thermostat: "kT, gamma0 and gammav have to be given as keyword args") if not isinstance(kT, float): raise ValueError("temperature must be a positive number") + + # Seed is required if the RNG is not initialized + if seed is None and npt_iso_is_seed_required(): + raise ValueError( + "A seed has to be given as keyword argument on first activation of the thermostat") + + if seed is not None: + utils.check_type_or_throw_except( + seed, 1, int, "seed must be a positive integer") + if seed < 0: + raise ValueError("seed must be a positive integer") + npt_iso_set_rng_state(seed) + global temperature temperature = float(kT) global thermo_switch diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index 932e9298a68..bcc011b3c13 100644 --- a/testsuite/python/integrator_npt.py +++ b/testsuite/python/integrator_npt.py @@ -29,7 +29,6 @@ class IntegratorNPT(ut.TestCase): """This compares pressure and compressibility of a LJ system against expected values.""" S = espressomd.System(box_l=[1.0, 1.0, 1.0]) - S.seed = S.cell_system.get_state()['n_nodes'] * [1234] p_ext = 2.0 def setUp(self): @@ -39,8 +38,7 @@ def setUp(self): self.S.time_step = 0.01 self.S.cell_system.skin = 0.25 - data = np.genfromtxt(tests_common.abspath( - "data/npt_lj_system.data")) + data = np.genfromtxt(tests_common.abspath("data/npt_lj_system.data")) # Input format: id pos f for particle in data: @@ -51,7 +49,7 @@ def setUp(self): self.S.non_bonded_inter[0, 0].lennard_jones.set_params( epsilon=1, sigma=1, cutoff=1.12246, shift=0.25) - self.S.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.004) + self.S.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.004, seed=42) self.S.integrator.set_isotropic_npt( ext_pressure=self.p_ext, piston=0.0001) @@ -69,11 +67,10 @@ def test_npt(self): avp /= (n / skip_p) Vs = np.array(ls)**3 - compressibility = pow(np.std(Vs), 2) / np.average(Vs) - print(avp, compressibility) + compressibility = np.var(Vs) / np.average(Vs) - self.assertAlmostEqual(2.0, avp, delta=0.02) - self.assertAlmostEqual(0.2, compressibility, delta=0.02) + self.assertAlmostEqual(avp, 2.0, delta=0.02) + self.assertAlmostEqual(compressibility, 0.32, delta=0.02) if __name__ == "__main__": diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index 377e3018436..c14c8b0a83f 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -141,7 +141,7 @@ elif 'THERM.BD' in modes: system.thermostat.set_brownian(kT=1.0, gamma=2.0, seed=42) elif 'THERM.NPT' in modes and has_features('NPT'): - system.thermostat.set_npt(kT=1.0, gamma0=2.0, gammav=0.1) + system.thermostat.set_npt(kT=1.0, gamma0=2.0, gammav=0.1, seed=42) elif 'THERM.DPD' in modes and has_features('DPD'): system.thermostat.set_dpd(kT=1.0, seed=42) # set integrator diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index d07d0660d88..e9336a61387 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -187,6 +187,7 @@ def test_thermostat_DPD(self): def test_thermostat_NPT(self): thmst = system.thermostat.get_state()[0] self.assertEqual(thmst['type'], 'NPT_ISO') + self.assertEqual(thmst['seed'], 42) self.assertEqual(thmst['gamma0'], 2.0) self.assertEqual(thmst['gammav'], 0.1)