diff --git a/doc/doxygen/bibliography.bib b/doc/doxygen/bibliography.bib index e30b3f9fd6b..df73a90672e 100644 --- a/doc/doxygen/bibliography.bib +++ b/doc/doxygen/bibliography.bib @@ -351,6 +351,17 @@ @Article{neumann85b doi = {10.1063/1.448553}, } +@Book{Pottier2010, + title={{Nonequilibrium Statistical Physics}}, + subtitle={{Linear Irreversible Processes}}, + author={Pottier, No\"{e}lle}, + year={2010}, + series={Oxford Graduate Texts}, + publisher={OUP Oxford}, + isbn={9780199556885}, + url={https://global.oup.com/academic/product/nonequilibrium-statistical-physics-9780199556885}, +} + @Article{reed92a, title = {{Monte Carlo study of titration of linear polyelectrolytes}}, author = {Reed, C. E. and Reed, W. F.}, @@ -361,6 +372,18 @@ @Article{reed92a doi = {10.1063/1.462145}, } +@Book{schlick2010, +title = {{Molecular Modeling and Simulation: An Interdisciplinary Guide}}, +author = {Schlick, Tamar}, +series = {Interdisciplinary Applied Mathematics}, +volume = {21}, +year = {2010}, +publisher = {Springer New York}, +address = {New York, NY}, +isbn = {978-1-4419-6350-5}, +doi = {10.1007/978-1-4419-6351-2}, +} + @Article{smith94c, title={{The reaction ensemble method for the computer simulation of chemical and phase equilibria. I. Theory and basic examples}}, author={Smith, W. R. and Triska, B.}, @@ -398,7 +421,7 @@ @book{allen2017 title={Computer simulation of liquids}, author={Allen, Michael P and Tildesley, Dominic J}, year={2017}, - publisher={Oxford university press}, + publisher={Oxford University Press}, url={https://global.oup.com/academic/product/computer-simulation-of-liquids-9780198803201}, doi={10.1093/oso/9780198803195.001.0001}, isbn={9780198803195}, diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index 7fabba2fdb5..a4214956f67 100644 --- a/doc/sphinx/analysis.rst +++ b/doc/sphinx/analysis.rst @@ -411,7 +411,7 @@ are calculated is beyond the scope of this manual. Three body potentials are implemented following the procedure in Ref. :cite:`thompson09a`. A different formula is used to calculate contribution from electrostatic interactions. For -electrostatic interactions in P3M, the :math:`k`-space contribution is implemented according to :cite:`essmann1995smooth`. +electrostatic interactions in P3M, the :math:`k`-space contribution is implemented according to :cite:`essmann95a`. The implementation of the Coulomb P3M pressure is tested against LAMMPS. Four-body dihedral potentials are not included. Except of @@ -439,7 +439,7 @@ The instantaneous virial stress tensor is calculated by where the notation is the same as for the pressure. The superscripts :math:`k` and :math:`l` correspond to the components in the tensors and vectors. -If electrostatic interactions are present then also the coulombic parts of the stress tensor need to be calculated. If P3M is present, then the instantaneous stress tensor is added to the above equation in accordance with :cite:`essmann1995smooth` : +If electrostatic interactions are present then also the coulombic parts of the stress tensor need to be calculated. If P3M is present, then the instantaneous stress tensor is added to the above equation in accordance with :cite:`essmann95a` : .. math :: p^\text{Coulomb, P3M}_{(k,l)} =p^\text{Coulomb, P3M, dir}_{(k,l)} + p^\text{Coulomb, P3M, rec}_{(k,l)}, diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index aa84982c890..85e71acce01 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -267,15 +267,15 @@ The Langevin thermostat is based on an extension of Newton's equation of motion .. math:: m_i \dot{v}_i(t) = f_i(\{x_j\},v_i,t) - \gamma v_i(t) + \sqrt{2\gamma k_B T} \eta_i(t). -Here, :math:`f_i` are all deterministic forces from interactions, -:math:`\gamma` the bare friction coefficient and :math:`\eta` a random, "thermal" force. +Here, :math:`f_i` are all deterministic forces from interactions, +:math:`\gamma` the bare friction coefficient and :math:`\eta` a random, "thermal" force. The friction term accounts for dissipation in a surrounding fluid whereas -the random force mimics collisions of the particle with solvent molecules +the random force mimics collisions of the particle with solvent molecules at temperature :math:`T` and satisfies .. math:: <\eta(t)> = 0 , <\eta^\alpha_i(t)\eta^\beta_j(t')> = \delta_{\alpha\beta} \delta_{ij}\delta(t-t') -(:math:`<\cdot>` denotes the ensemble average and :math:`\alpha,\beta` are spatial coordinates). +(:math:`<\cdot>` denotes the ensemble average and :math:`\alpha,\beta` are spatial coordinates). In the |es| implementation of the Langevin thermostat, the additional terms only enter in the force calculation. @@ -283,9 +283,9 @@ This reduces the accuracy of the Velocity Verlet integrator by one order in :math:`dt` because forces are now velocity-dependent. The random process :math:`\eta(t)` is discretized by drawing an uncorrelated random number -:math:`\overline{\eta}` for each component of all the particle forces. +:math:`\overline{\eta}` for each component of all the particle forces. The distribution of :math:`\overline{\eta}` is uniform and satisfies - + .. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt The keyword ``seed`` controls the state of the random number generator (Philox @@ -305,9 +305,9 @@ can be useful, for instance, in high Péclet number active matter systems, where one only wants to thermalize only the rotational degrees of freedom and translational motion is effected by the self-propulsion. -The keywords ``gamma`` and ``gamma_rotate`` can be specified as a scalar, -or, with feature ``PARTICLE_ANISOTROPY`` compiled in, as the three eigenvalues -of the respective friction coefficient tensor. This is enables the simulation of +The keywords ``gamma`` and ``gamma_rotate`` can be specified as a scalar, +or, with feature ``PARTICLE_ANISOTROPY`` compiled in, as the three eigenvalues +of the respective friction coefficient tensor. This is enables the simulation of the anisotropic diffusion of anisotropic colloids (rods, etc.). Using the Langevin thermostat, it is possible to set a temperature and a @@ -324,17 +324,17 @@ The :ref:`Lattice-Boltzmann` thermostat acts similar to the :ref:`Langevin therm .. math:: m_i \dot{v}_i(t) = f_i(\{x_j\},v_i,t) - \gamma (v_i(t)-u(x_i(t),t)) + \sqrt{2\gamma k_B T} \eta_i(t). -where :math:`u(x,t)` is the fluid velocity at position :math:`x` and time :math:`t`. -To preserve momentum, an equal and opposite friction force and random force act on the fluid. +where :math:`u(x,t)` is the fluid velocity at position :math:`x` and time :math:`t`. +To preserve momentum, an equal and opposite friction force and random force act on the fluid. -Numerically the fluid velocity is determined from the lattice-Boltzmann node velocities +Numerically the fluid velocity is determined from the lattice-Boltzmann node velocities by interpolating as described in :ref:`Interpolating velocities`. -The backcoupling of friction forces and noise to the fluid is also done by distributing those forces amongst the nearest LB nodes. +The backcoupling of friction forces and noise to the fluid is also done by distributing those forces amongst the nearest LB nodes. Details for both the interpolation and the force distribution can be found in :cite:`ahlrichs99` and :cite:`duenweg08a`. The LB fluid can be used to thermalize particles, while also including their hydrodynamic interactions. The LB thermostat expects an instance of either :class:`espressomd.lb.LBFluid` or :class:`espressomd.lb.LBFluidGPU`. -Temperature is set via the ``kT`` argument of the LB fluid. +Temperature is set via the ``kT`` argument of the LB fluid. Furthermore a ``seed`` has to be given for the thermalization of the particle coupling. The magnitude of the frictional coupling can be adjusted by @@ -362,8 +362,8 @@ Dissipative Particle Dynamics (DPD) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The DPD thermostat adds friction and noise to the particle -dynamics like the :ref:`Langevin thermostat`, but these -are not applied to every particle individually but instead +dynamics like the :ref:`Langevin thermostat`, but these +are not applied to every particle individually but instead encoded in a dissipative interaction between particles :cite:`soddeman03a`. To realize a complete DPD fluid model in |es|, three parts are needed: @@ -376,7 +376,7 @@ The temperature is set via which takes ``kT`` as the only argument. The friction coefficients and cutoff are controlled via the -:ref:`DPD interaction` on a per type-pair basis. For details +:ref:`DPD interaction` on a per type-pair basis. For details see there. The friction (dissipative) and noise (random) term are coupled via the @@ -428,6 +428,86 @@ Be aware that this feature is neither properly examined for all systems nor is it maintained regularly. If you use it and notice strange behavior, please contribute to solving the problem. +.. _Brownian thermostat: + +Brownian thermostat +~~~~~~~~~~~~~~~~~~~ + +Brownian thermostat is a formal name of a thermostat enabling the +Brownian Dynamics feature (see :cite:`schlick2010`) which implies +a propagation scheme involving systematic and thermal parts of the +classical Ermak-McCammom's (see :cite:`ermak78a`) +Brownian Dynamics. Currently it is implemented without +hydrodynamic interactions, i.e. +with a diagonal diffusion tensor. +The hydrodynamic interactions feature will be available later +as a part of the present Brownian Dynamics or +implemented separately within the Stokesian Dynamics. + +In order to activate the Brownian thermostat, the member function +:py:attr:`~espressomd.thermostat.Thermostat.set_brownian` of the thermostat +class :class:`espressomd.thermostat.Thermostat` has to be invoked. +The system integrator should be also changed. +Best explained in an example:: + + import espressomd + system = espressomd.System() + system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=41) + system.integrator.set_brownian_dynamics() + +where ``gamma`` (hereinafter :math:`\gamma`) is a viscous friction coefficient. +In terms of the Python interface and setup, the Brownian thermostat is very +similar to the :ref:`Langevin thermostat`. The feature +``BROWNIAN_PER_PARTICLE`` is used to control the per-particle +temperature and the friction coefficient setup. The major differences are +its internal integrator implementation and other temporal constraints. +The integrator is still a symplectic Velocity Verlet-like one. +It is implemented via a viscous drag part and a random walk of both the position and +velocity. Due to a nature of the Brownian Dynamics method, its time step :math:`\Delta t` +should be large enough compared to the relaxation time +:math:`m/\gamma` where :math:`m` is the particle mass. +This requirement is just a conceptual one +without specific implementation technical restrictions. +Note that with all similarities of +Langevin and Brownian Dynamics, the Langevin thermostat temporal constraint +is opposite. A velocity is restarting from zero at every step. +Formally, the previous step velocity at the beginning of the the :math:`\Delta t` interval +is dissipated further +and does not contribute to the end one as well as to the positional random walk. +Another temporal constraint +which is valid for both Langevin and Brownian Dynamics: conservative forces +should not change significantly over the :math:`\Delta t` interval. + +The viscous terminal velocity :math:`\Delta v` and corresponding positional +step :math:`\Delta r` are fully driven by conservative forces :math:`F`: + +.. math:: \Delta r = \frac{F \cdot \Delta t}{\gamma} + +.. math:: \Delta v = \frac{F}{\gamma} + +A positional random walk variance of each coordinate :math:`\sigma_p^2` +corresponds to a diffusion within the Wiener process: + +.. math:: \sigma_p^2 = 2 \frac{kT}{\gamma} \cdot \Delta t + +Each velocity component random walk variance :math:`\sigma_v^2` is defined by the heat +component: + +.. math:: \sigma_v^2 = \frac{kT}{m} + +Note that the velocity random walk is propagated from zero at each step. + +A rotational motion is implemented similarly. +Note: the rotational Brownian dynamics implementation is compatible with particles which have +the isotropic moment of inertia tensor only. Otherwise, the viscous terminal angular velocity +is not defined, i.e. it has no constant direction over the time. + +The keyword ``seed`` controls the state of the random number generator (Philox +Counter-based RNG) and is required on first activation of the thermostat. It +can be omitted in subsequent calls of ``set_brownian()``. It is the user's +responsibility to decide whether the thermostat should be deterministic (by +using a fixed seed) or not (by using a randomized seed). + .. _CUDA: CUDA diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index c5beac32854..82d4d9b7ce2 100644 --- a/doc/sphinx/zrefs.bib +++ b/doc/sphinx/zrefs.bib @@ -993,6 +993,30 @@ @ARTICLE{wolff04a timestamp = {2007.06.14} } +@book{schlick2010, +address = {New York, NY}, +author = {Schlick, Tamar}, +doi = {10.1007/978-1-4419-6351-2}, +isbn = {978-1-4419-6350-5}, +publisher = {Springer New York}, +series = {Interdisciplinary Applied Mathematics}, +title = {{Molecular Modeling and Simulation: An Interdisciplinary Guide}}, +volume = {21}, +year = {2010}, +} + +@article{ermak78a, + title={Brownian dynamics with hydrodynamic interactions}, + author={Ermak, Donald L and McCammon, J Andrew}, + journal={J. Chem. Phys.}, + volume={69}, + number={4}, + pages={1352--1360}, + year={1978}, + publisher={AIP}, + doi={10.1063/1.436761}, +} + @Article{cerda08d, title = {{P3M} algorithm for dipolar interactions.}, author = {Juan J. Cerd\`{a} and V. Ballenegger and O. Lenz and Ch. Holm}, @@ -1006,7 +1030,7 @@ @Article{cerda08d timestamp = {2008.01.14} } -@article{essmann1995smooth, +@article{essmann95a, title={A smooth particle mesh Ewald method}, author={Essmann, Ulrich and Perera, Lalith and Berkowitz, Max L and Darden, Tom and Lee, Hsing and Pedersen, Lee G}, journal={J. Chem. Phys.}, diff --git a/maintainer/configs/maxset.hpp b/maintainer/configs/maxset.hpp index aac46dbfe72..7ea4862c139 100644 --- a/maintainer/configs/maxset.hpp +++ b/maintainer/configs/maxset.hpp @@ -33,6 +33,7 @@ #define BOND_CONSTRAINT #define COLLISION_DETECTION #define LANGEVIN_PER_PARTICLE +#define BROWNIAN_PER_PARTICLE #define SWIMMER_REACTIONS #define NPT diff --git a/maintainer/configs/no_rotation.hpp b/maintainer/configs/no_rotation.hpp index 5412eadf2fc..8e2fd495f81 100644 --- a/maintainer/configs/no_rotation.hpp +++ b/maintainer/configs/no_rotation.hpp @@ -30,6 +30,7 @@ #define MASS #define EXTERNAL_FORCES #define LANGEVIN_PER_PARTICLE +#define BROWNIAN_PER_PARTICLE #define BOND_CONSTRAINT #define NPT #define DPD diff --git a/src/config/features.def b/src/config/features.def index 8bfe30874a3..a16d9319e4d 100644 --- a/src/config/features.def +++ b/src/config/features.def @@ -19,6 +19,7 @@ MASS EXCLUSIONS BOND_CONSTRAINT LANGEVIN_PER_PARTICLE +BROWNIAN_PER_PARTICLE COLLISION_DETECTION METADYNAMICS NPT diff --git a/src/config/myconfig-default.hpp b/src/config/myconfig-default.hpp index da600142b02..54465b344f7 100644 --- a/src/config/myconfig-default.hpp +++ b/src/config/myconfig-default.hpp @@ -33,6 +33,7 @@ #define PARTICLE_ANISOTROPY #define EXTERNAL_FORCES #define LANGEVIN_PER_PARTICLE +#define BROWNIAN_PER_PARTICLE #define BOND_CONSTRAINT #define NPT #define DPD diff --git a/src/core/global.cpp b/src/core/global.cpp index cde3fe5e04f..c6dc5802fc3 100644 --- a/src/core/global.cpp +++ b/src/core/global.cpp @@ -168,7 +168,26 @@ const std::unordered_map fields{ "n_thermalized_bonds"}}, /* 56 from thermalized_bond.cpp */ {FIELD_FORCE_CAP, {&force_cap, Datafield::Type::DOUBLE, 1, "force_cap"}}, {FIELD_THERMO_VIRTUAL, - {&thermo_virtual, Datafield::Type::BOOL, 1, "thermo_virtual"}}}; + {&thermo_virtual, Datafield::Type::BOOL, 1, "thermo_virtual"}}, +#ifndef PARTICLE_ANISOTROPY + {FIELD_BROWNIAN_GAMMA_ROTATION, + {&brownian.gamma_rotation, Datafield::Type::DOUBLE, 1, + "brownian.gamma_rotation"}}, /* 57 from thermostat.cpp */ +#else + {FIELD_BROWNIAN_GAMMA_ROTATION, + {brownian.gamma_rotation.data(), Datafield::Type::DOUBLE, 3, + "brownian.gamma_rotation"}}, /* 57 from thermostat.cpp */ +#endif +#ifndef PARTICLE_ANISOTROPY + {FIELD_BROWNIAN_GAMMA, + {&brownian.gamma, Datafield::Type::DOUBLE, 1, + "brownian.gamma"}}, /* 58 from thermostat.cpp */ +#else + {FIELD_BROWNIAN_GAMMA, + {brownian.gamma.data(), Datafield::Type::DOUBLE, 3, + "brownian.gamma"}}, /* 58 from thermostat.cpp */ +#endif // PARTICLE_ANISOTROPY +}; std::size_t hash_value(Datafield const &field) { using boost::hash_range; diff --git a/src/core/global.hpp b/src/core/global.hpp index d2c28952e86..30e5417b11d 100644 --- a/src/core/global.hpp +++ b/src/core/global.hpp @@ -98,7 +98,11 @@ enum Fields { FIELD_THERMALIZEDBONDS, FIELD_FORCE_CAP, FIELD_THERMO_VIRTUAL, - FIELD_SWIMMING_PARTICLES_EXIST + FIELD_SWIMMING_PARTICLES_EXIST, + /** index of \ref BrownianThermostat::gamma */ + FIELD_BROWNIAN_GAMMA, + /** index of \ref BrownianThermostat::gamma_rotation */ + FIELD_BROWNIAN_GAMMA_ROTATION, }; #endif diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 62bfdcfa5ad..7ea4e93a683 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -57,11 +57,13 @@ #include "thermostat.hpp" #include "virtual_sites.hpp" +#include "integrators/brownian_inline.hpp" #include "integrators/steepest_descent.hpp" #include "integrators/velocity_verlet_inline.hpp" #include "integrators/velocity_verlet_npt.hpp" #include +#include #include #include @@ -126,6 +128,10 @@ bool integrator_step_1(ParticleRange &particles) { velocity_verlet_npt_step_1(particles); break; #endif + case INTEG_METHOD_BD: + // the Ermak-McCammon's Brownian Dynamics requires a single step + // so, just skip here + break; default: throw std::runtime_error("Unknown value for integ_switch"); } @@ -146,6 +152,10 @@ void integrator_step_2(ParticleRange &particles) { velocity_verlet_npt_step_2(particles); break; #endif + case INTEG_METHOD_BD: + // the Ermak-McCammon's Brownian Dynamics requires a single step + brownian_dynamics_propagator(particles); + break; default: throw std::runtime_error("Unknown value for INTEG_SWITCH"); } @@ -305,6 +315,9 @@ void philox_counter_increment() { if (thermo_switch & THERMO_LANGEVIN) { langevin_rng_counter_increment(); } + if (thermo_switch & THERMO_BROWNIAN) { + brownian_rng_counter_increment(); + } if (thermo_switch & THERMO_DPD) { #ifdef DPD dpd_rng_counter_increment(); @@ -402,6 +415,11 @@ void integrate_set_nvt() { mpi_bcast_parameter(FIELD_INTEG_SWITCH); } +void integrate_set_bd() { + integ_switch = INTEG_METHOD_BD; + mpi_bcast_parameter(FIELD_INTEG_SWITCH); +} + #ifdef NPT int integrate_set_npt_isotropic(double ext_pressure, double piston, bool xdir_rescale, bool ydir_rescale, diff --git a/src/core/integrate.hpp b/src/core/integrate.hpp index ee25c85ecf3..1723c1b0eca 100644 --- a/src/core/integrate.hpp +++ b/src/core/integrate.hpp @@ -30,6 +30,7 @@ #define INTEG_METHOD_NPT_ISO 0 #define INTEG_METHOD_NVT 1 #define INTEG_METHOD_STEEPEST_DESCENT 2 +#define INTEG_METHOD_BD 3 /** Switch determining which integrator to use. */ extern int integ_switch; @@ -124,6 +125,9 @@ int integrate_set_steepest_descent(double f_max, double gamma, int max_steps, /** @brief Set the velocity Verlet integrator for the NVT ensemble. */ void integrate_set_nvt(); +/** @brief Set the Brownian Dynamics integrator. */ +void integrate_set_bd(); + /** @brief Set the velocity Verlet integrator modified for the NpT ensemble * with isotropic rescaling. * diff --git a/src/core/integrators/brownian_inline.hpp b/src/core/integrators/brownian_inline.hpp new file mode 100644 index 00000000000..27d055bceb7 --- /dev/null +++ b/src/core/integrators/brownian_inline.hpp @@ -0,0 +1,571 @@ +/* + * Copyright (C) 2010-2019 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 . + */ +/** \file */ + +#ifndef BROWNIAN_INLINE_HPP +#define BROWNIAN_INLINE_HPP + +#include "random.hpp" +#include "thermostat.hpp" + +extern BrownianThermostat brownian; + +/** Propagate position: viscous drag driven by conservative forces. + * From eq. (14.39) in @cite Schlick2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +inline void bd_drag(Particle &p, double dt) { + // The friction tensor Z from the Eq. (14.31) of Schlick2010: + Thermostat::GammaType local_gamma; + +#ifdef BROWNIAN_PER_PARTICLE + if (p.p.gamma >= Thermostat::GammaType{}) { + local_gamma = p.p.gamma; + } else +#endif + { + local_gamma = brownian.gamma; + } + + bool aniso_flag; // particle anisotropy flag + +#ifdef PARTICLE_ANISOTROPY + // Particle frictional isotropy check. + aniso_flag = + (local_gamma[0] != local_gamma[1]) || (local_gamma[1] != local_gamma[2]); +#endif + + Utils::Vector3d force_body; + Utils::Vector3d delta_pos_body, delta_pos_lab; + +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + force_body = convert_vector_space_to_body(p, p.f.f); + } +#endif + + for (int j = 0; j < 3; j++) { + // Second (deterministic) term of the Eq. (14.39) of Schlick2010. + // Only a conservative part of the force is used here +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + delta_pos_body[j] = force_body[j] * dt / (local_gamma[j]); + } else { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.r.p[j] += p.f.f[j] * dt / (local_gamma[j]); + } + } +#else +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.r.p[j] += p.f.f[j] * dt / (local_gamma); + } +#endif // PARTICLE_ANISOTROPY + } +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + delta_pos_lab = convert_vector_body_to_space(p, delta_pos_body); + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.r.p[j] += delta_pos_lab[j]; + } + } + } +#endif // PARTICLE_ANISOTROPY +} + +/** Set the terminal velocity driven by the conservative forces drag. + * From eq. (14.34) in @cite Schlick2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +inline void bd_drag_vel(Particle &p, double dt) { + // The friction tensor Z from the eq. (14.31) of Schlick2010: + Thermostat::GammaType local_gamma; + +#ifdef BROWNIAN_PER_PARTICLE + if (p.p.gamma >= Thermostat::GammaType{}) { + local_gamma = p.p.gamma; + } else +#endif + { + local_gamma = brownian.gamma; + } + + bool aniso_flag; // particle anisotropy flag + +#ifdef PARTICLE_ANISOTROPY + // Particle frictional isotropy check. + aniso_flag = + (local_gamma[0] != local_gamma[1]) || (local_gamma[1] != local_gamma[2]); +#endif + + Utils::Vector3d force_body; + Utils::Vector3d vel_body, vel_lab; + +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + force_body = convert_vector_space_to_body(p, p.f.f); + } +#endif + + for (int j = 0; j < 3; j++) { + // First (deterministic) term of the eq. (14.34) of Schlick2010 taking + // into account eq. (14.35). Only conservative part of the force is used + // here. NOTE: velocity is assigned here and propagated by thermal part + // further on top of it +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + vel_body[j] = force_body[j] / (local_gamma[j]); + } else { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.m.v[j] = p.f.f[j] / (local_gamma[j]); + } +#ifdef EXTERNAL_FORCES + else { + p.m.v[j] = 0.0; + } +#endif + } +#else // PARTICLE_ANISOTROPY +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.m.v[j] = p.f.f[j] / (local_gamma); + } +#endif // PARTICLE_ANISOTROPY + } + +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + vel_lab = convert_vector_body_to_space(p, vel_body); + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.m.v[j] = vel_lab[j]; + } +#ifdef EXTERNAL_FORCES + else { + p.m.v[j] = 0.0; + } +#endif + } + } +#endif // PARTICLE_ANISOTROPY +} + +/** Propagate the positions: random walk part. + * From eq. (14.37) in @cite Schlick2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +void bd_random_walk(Particle &p, double dt) { + // skip the translation thermalizing for virtual sites unless enabled + extern bool thermo_virtual; + if (p.p.is_virtual && !thermo_virtual) + return; + // first, set defaults + Thermostat::GammaType brown_sigma_pos_inv = brownian.sigma_pos_inv; + + // Override defaults if per-particle values for T and gamma are given +#ifdef BROWNIAN_PER_PARTICLE + auto const constexpr brown_temp_coeff = 2.0; + + if (p.p.gamma >= Thermostat::GammaType{}) { + // Is a particle-specific temperature also specified? + if (p.p.T >= 0.) { + if (p.p.T > 0.0) { + brown_sigma_pos_inv = sqrt(p.p.gamma / (brown_temp_coeff * p.p.T)); + } else { + // just an indication of the infinity + brown_sigma_pos_inv = brownian.gammatype_nan; + } + } else + // default temperature but particle-specific gamma + if (temperature > 0.0) { + brown_sigma_pos_inv = sqrt(p.p.gamma / (brown_temp_coeff * temperature)); + } else { + brown_sigma_pos_inv = brownian.gammatype_nan; + } + } // particle-specific gamma + else { + // No particle-specific gamma, but is there particle-specific temperature + if (p.p.T >= 0.) { + if (p.p.T > 0.0) { + brown_sigma_pos_inv = sqrt(brownian.gamma / (brown_temp_coeff * p.p.T)); + } else { + // just an indication of the infinity + brown_sigma_pos_inv = brownian.gammatype_nan; + } + } else { + // default values for both + brown_sigma_pos_inv = brownian.sigma_pos_inv; + } + } +#endif /* BROWNIAN_PER_PARTICLE */ + + bool aniso_flag; // particle anisotropy flag + +#ifdef PARTICLE_ANISOTROPY + // Particle frictional isotropy check. + aniso_flag = (brown_sigma_pos_inv[0] != brown_sigma_pos_inv[1]) || + (brown_sigma_pos_inv[1] != brown_sigma_pos_inv[2]); +#else + aniso_flag = false; +#endif // PARTICLE_ANISOTROPY + + Utils::Vector3d delta_pos_body, delta_pos_lab; + + // Eq. (14.37) is factored by the Gaussian noise (12.22) with its squared + // magnitude defined in the second eq. (14.38), Schlick2010. + Utils::Vector3d noise = Random::v_noise_g( + brownian.rng_counter->value(), p.p.identity); + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { +#ifndef PARTICLE_ANISOTROPY + if (brown_sigma_pos_inv > 0.0) { + delta_pos_body[j] = (1.0 / brown_sigma_pos_inv) * sqrt(dt) * noise[j]; + } else { + delta_pos_body[j] = 0.0; + } +#else + if (brown_sigma_pos_inv[j] > 0.0) { + delta_pos_body[j] = + (1.0 / brown_sigma_pos_inv[j]) * sqrt(dt) * noise[j]; + } else { + delta_pos_body[j] = 0.0; + } +#endif // PARTICLE_ANISOTROPY + } + } + +#ifdef PARTICLE_ANISOTROPY + if (aniso_flag) { + delta_pos_lab = convert_vector_body_to_space(p, delta_pos_body); + } +#endif + + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + p.r.p[j] += aniso_flag ? delta_pos_lab[j] : delta_pos_body[j]; + } + } +} + +/** Determine the velocities: random walk part. + * From eq. (10.2.16) in @cite Pottier2010. + * @param[in, out] p Particle + * @param[in] dt Time interval + */ +inline void bd_random_walk_vel(Particle &p, double dt) { + // skip the translation thermalizing for virtual sites unless enabled + extern bool thermo_virtual; + if (p.p.is_virtual && !thermo_virtual) + return; + // Just a square root of kT, see eq. (10.2.17) and comments in 2 paragraphs + // afterwards, Pottier2010 + double brown_sigma_vel; + + // Override defaults if per-particle values for T and gamma are given +#ifdef BROWNIAN_PER_PARTICLE + auto const constexpr brown_temp_coeff = 1.0; + // Is a particle-specific temperature specified? + if (p.p.T >= 0.) { + brown_sigma_vel = sqrt(brown_temp_coeff * p.p.T); + } else { + brown_sigma_vel = brownian.sigma_vel; + } +#else + // defaults + brown_sigma_vel = brownian.sigma_vel; +#endif /* BROWNIAN_PER_PARTICLE */ + + Utils::Vector3d noise = Random::v_noise_g( + brownian.rng_counter->value(), p.identity()); + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + // Random (heat) velocity is added here. It is already initialized in the + // terminal drag part. See eq. (10.2.16) taking into account eq. (10.2.18) + // and (10.2.29), Pottier2010. Note, that the Pottier2010 units system + // (see Eq. (10.1.1) there) has been adapted towards the ESPResSo and the + // referenced above Schlick2010 one, which is defined by the eq. (14.31) + // of Schlick2010. A difference is the mass factor to the friction tensor. + // The noise is Gaussian according to the convention at p. 237 (last + // paragraph), Pottier2010. + p.m.v[j] += brown_sigma_vel * noise[j] / sqrt(p.p.mass); + } + } +} + +#ifdef ROTATION + +/** Propagate quaternions: viscous drag driven by conservative torques. + * An analogy of eq. (14.39) in @cite Schlick2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +void bd_drag_rot(Particle &p, double dt) { + Thermostat::GammaType local_gamma; + +#ifdef BROWNIAN_PER_PARTICLE + if (p.p.gamma_rot >= Thermostat::GammaType{}) { + local_gamma = p.p.gamma_rot; + } else +#endif + { + local_gamma = brownian.gamma_rotation; + } + + Utils::Vector3d dphi = {0.0, 0.0, 0.0}; + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + // only a conservative part of the torque is used here +#ifndef PARTICLE_ANISOTROPY + dphi[j] = p.f.torque[j] * dt / (local_gamma); +#else + dphi[j] = p.f.torque[j] * dt / (local_gamma[j]); +#endif // ROTATIONAL_INERTIA + } + } // j + dphi = mask(p.p.rotation, dphi); + double dphi_m = dphi.norm(); + if (dphi_m) { + Utils::Vector3d dphi_u; + dphi_u = dphi / dphi_m; + local_rotate_particle_body(p, dphi_u, dphi_m); + } +} + +/** Set the terminal angular velocity driven by the conservative torques drag. + * An analogy of the 1st term of eq. (14.34) in @cite Schlick2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +void bd_drag_vel_rot(Particle &p, double dt) { + Thermostat::GammaType local_gamma; + +#ifdef BROWNIAN_PER_PARTICLE + if (p.p.gamma_rot >= Thermostat::GammaType{}) { + local_gamma = p.p.gamma_rot; + } else +#endif + { + local_gamma = brownian.gamma_rotation; + } + + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (p.p.ext_flag & COORD_FIXED(j)) { + p.m.omega[j] = 0.0; + } else +#endif + { + // only conservative part of the force is used here + // NOTE: velocity is assigned here and propagated by thermal part further + // on top of it +#ifndef PARTICLE_ANISOTROPY + p.m.omega[j] = p.f.torque[j] / (local_gamma); +#else + p.m.omega[j] = p.f.torque[j] / (local_gamma[j]); +#endif // ROTATIONAL_INERTIA + } + } + p.m.omega = mask(p.p.rotation, p.m.omega); +} + +/** Propagate the quaternions: random walk part. + * An analogy of eq. (14.37) in @cite Schlick2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +void bd_random_walk_rot(Particle &p, double dt) { + // first, set defaults + Thermostat::GammaType brown_sigma_pos_inv = brownian.sigma_pos_rotation_inv; + + // Override defaults if per-particle values for T and gamma are given +#ifdef BROWNIAN_PER_PARTICLE + auto const constexpr brown_temp_coeff = 2.0; + + if (p.p.gamma_rot >= Thermostat::GammaType{}) { + // Is a particle-specific temperature also specified? + if (p.p.T >= 0.) { + if (p.p.T > 0.0) { + brown_sigma_pos_inv = sqrt(p.p.gamma_rot / (brown_temp_coeff * p.p.T)); + } else { + // just an indication of the infinity + brown_sigma_pos_inv = brownian.gammatype_nan; + } + } else if (temperature > 0.) { + // Default temperature but particle-specific gamma + brown_sigma_pos_inv = + sqrt(p.p.gamma_rot / (brown_temp_coeff * temperature)); + } else { + // just an indication of the infinity + brown_sigma_pos_inv = brownian.gammatype_nan; + } + } // particle-specific gamma + else { + // No particle-specific gamma, but is there particle-specific temperature + if (p.p.T >= 0.) { + if (p.p.T > 0.0) { + brown_sigma_pos_inv = + sqrt(brownian.gamma_rotation / (brown_temp_coeff * p.p.T)); + } else { + // just an indication of the infinity + brown_sigma_pos_inv = brownian.gammatype_nan; + } + } else { + // Defaut values for both + brown_sigma_pos_inv = brownian.sigma_pos_rotation_inv; + } + } +#endif /* BROWNIAN_PER_PARTICLE */ + + Utils::Vector3d dphi = {0.0, 0.0, 0.0}; + Utils::Vector3d noise = Random::v_noise_g( + brownian.rng_counter->value(), p.p.identity); + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { +#ifndef PARTICLE_ANISOTROPY + if (brown_sigma_pos_inv > 0.0) { + dphi[j] = noise[j] * (1.0 / brown_sigma_pos_inv) * sqrt(dt); + } else { + dphi[j] = 0.0; + } +#else + if (brown_sigma_pos_inv[j] > 0.0) { + dphi[j] = noise[j] * (1.0 / brown_sigma_pos_inv[j]) * sqrt(dt); + } else { + dphi[j] = 0.0; + } +#endif // ROTATIONAL_INERTIA + } + } + dphi = mask(p.p.rotation, dphi); + // making the algorithm to be independ on an order of the rotations + double dphi_m = dphi.norm(); + if (dphi_m) { + Utils::Vector3d dphi_u; + dphi_u = dphi / dphi_m; + local_rotate_particle_body(p, dphi_u, dphi_m); + } +} + +/** Determine the angular velocities: random walk part. + * An analogy of eq. (10.2.16) in @cite Pottier2010. + * @param[in,out] p Particle + * @param[in] dt Time interval + */ +void bd_random_walk_vel_rot(Particle &p, double dt) { + double brown_sigma_vel; + + // Override defaults if per-particle values for T and gamma are given +#ifdef BROWNIAN_PER_PARTICLE + auto const constexpr brown_temp_coeff = 1.0; + // Is a particle-specific temperature specified? + if (p.p.T >= 0.) { + brown_sigma_vel = sqrt(brown_temp_coeff * p.p.T); + } else { + brown_sigma_vel = brownian.sigma_vel_rotation; + } +#else + // set defaults + brown_sigma_vel = brownian.sigma_vel_rotation; +#endif /* BROWNIAN_PER_PARTICLE */ + + Utils::Vector3d domega; + Utils::Vector3d noise = Random::v_noise_g( + brownian.rng_counter->value(), p.p.identity); + for (int j = 0; j < 3; j++) { +#ifdef EXTERNAL_FORCES + if (!(p.p.ext_flag & COORD_FIXED(j))) +#endif + { + // velocity is added here. It is already initialized in the terminal drag + // part. + domega[j] = brown_sigma_vel * noise[j] / sqrt(p.p.rinertia[j]); + } + } + domega = mask(p.p.rotation, domega); + p.m.omega += domega; +} +#endif // ROTATION + +inline void brownian_dynamics_propagator(const ParticleRange &particles) { + for (auto &p : particles) { + // Don't propagate translational degrees of freedom of vs +#ifdef VIRTUAL_SITES + extern bool thermo_virtual; + if (!(p.p.is_virtual) or thermo_virtual) +#endif + { + bd_drag(p, time_step); + bd_drag_vel(p, time_step); + bd_random_walk(p, time_step); + bd_random_walk_vel(p, time_step); + /* Verlet criterion check */ + if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) + set_resort_particles(Cells::RESORT_LOCAL); + } +#ifdef ROTATION + if (!p.p.rotation) + continue; + convert_torque_to_body_frame_apply_fix(p); + bd_drag_rot(p, time_step); + bd_drag_vel_rot(p, time_step); + bd_random_walk_rot(p, time_step); + bd_random_walk_vel_rot(p, time_step); +#endif // ROTATION + } + sim_time += time_step; +} + +#endif // BROWNIAN_INLINE_HPP diff --git a/src/core/random.hpp b/src/core/random.hpp index fcf91b7701b..08f02aad059 100644 --- a/src/core/random.hpp +++ b/src/core/random.hpp @@ -30,6 +30,7 @@ #include #include +#include #include #include @@ -50,13 +51,17 @@ enum class RNGSalt : uint64_t { PARTICLES, LANGEVIN, LANGEVIN_ROT, + BROWNIAN_WALK, + BROWNIAN_INC, + BROWNIAN_ROT_INC, + BROWNIAN_ROT_WALK, SALT_DPD, THERMALIZED_BOND }; namespace Random { /** - * @brief 3d uniform vector noise. + * @brief get 4 random uint 64 from the Philox RNG * * This uses the Philox PRNG, the state is controlled * by the counter, the salt and two keys. @@ -66,7 +71,8 @@ namespace Random { * */ template -Utils::Vector3d v_noise(uint64_t counter, int key1, int key2 = 0) { +Utils::Vector philox_4_uint64s(uint64_t counter, int key1, + int key2 = 0) { using rng_type = r123::Philox4x64; using ctr_type = rng_type::ctr_type; @@ -78,7 +84,24 @@ Utils::Vector3d v_noise(uint64_t counter, int key1, int key2 = 0) { auto const id2 = static_cast(key2); const key_type k{id1, id2}; - auto const noise = rng_type{}(c, k); + auto const res = rng_type{}(c, k); + return {res[0], res[1], res[2], res[3]}; +} + +/** + * @brief 3d uniform vector noise. + * + * 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 +Utils::Vector3d v_noise(uint64_t counter, int key1, int key2 = 0) { + + auto const noise = philox_4_uint64s(counter, key1, key2); using Utils::uniform; return Utils::Vector3d{uniform(noise[0]), uniform(noise[1]), @@ -86,6 +109,48 @@ Utils::Vector3d v_noise(uint64_t counter, int key1, int key2 = 0) { Utils::Vector3d::broadcast(0.5); } +/** @brief Generator for Gaussian random 3d vector. + * + * 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 + * random numbers are not zero (approx one in 2^64). To avoid this case, + * such numbers are replaced by std::numeric_limits::min() + * This breaks statistics in rare cases but allows for consistent RNG + * counters across MPI ranks. + * + * @param counter counter for the random number generation + * @param key1 key for random number generation + * @param key2 key for random number generation + * @tparam salt (decorrelates different thermostat types) + * + * @return 3D vector of Gaussian random numbers. + * + */ +template +inline Utils::Vector3d v_noise_g(uint64_t counter, int key1, int key2 = 0) { + + auto const noise = philox_4_uint64s(counter, key1, key2); + using Utils::uniform; + Utils::Vector4d u{uniform(noise[0]), uniform(noise[1]), uniform(noise[2]), + uniform(noise[3])}; + + // Replace zeros from uniformly distributed numbers (see doc string) + static const double epsilon = std::numeric_limits::min(); + for (int i = 0; i < 4; i++) + if (u[i] < epsilon) + u[i] = epsilon; + + // Box muller transform code adapted from + // https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform + + static const double two_pi = 2.0 * Utils::pi(); + return {sqrt(-2.0 * log(u[0])) * cos(two_pi * u[1]), + sqrt(-2.0 * log(u[0])) * sin(two_pi * u[1]), + sqrt(-2.0 * log(u[2])) * cos(two_pi * u[3])}; +} + extern std::mt19937 generator; extern std::uniform_real_distribution uniform_real_distribution; extern bool user_has_seeded; diff --git a/src/core/rotation.cpp b/src/core/rotation.cpp index f0b044f339f..2b5ec9b045f 100644 --- a/src/core/rotation.cpp +++ b/src/core/rotation.cpp @@ -38,6 +38,7 @@ #ifdef ROTATION #include "integrate.hpp" +#include #include #include @@ -151,11 +152,6 @@ void propagate_omega_quat_particle(Particle &p) { } } -inline void convert_torque_to_body_frame_apply_fix(Particle &p) { - auto const torque = convert_vector_space_to_body(p, p.f.torque); - p.f.torque = mask(p.p.rotation, torque); -} - /** convert the torques to the body-fixed frames and propagate angular * velocities */ void convert_torques_propagate_omega(const ParticleRange &particles) { @@ -201,26 +197,4 @@ void convert_initial_torques(const ParticleRange &particles) { } } -/** Rotate the particle p around the NORMALIZED axis aSpaceFrame by amount phi - */ -void local_rotate_particle(Particle &p, const Utils::Vector3d &axis_space_frame, - const double phi) { - // Rotation turned off entirely? - if (!p.p.rotation) - return; - - // Convert rotation axis to body-fixed frame - auto const axis = - mask(p.p.rotation, convert_vector_space_to_body(p, axis_space_frame)) - .normalize(); - - auto const s = std::sin(phi / 2); - auto const q = - Utils::Vector4d{cos(phi / 2), s * axis[0], s * axis[1], s * axis[2]} - .normalize(); - - // Rotate the particle - p.r.quat = Utils::multiply_quaternions(p.r.quat, q); -} - #endif // ROTATION diff --git a/src/core/rotation.hpp b/src/core/rotation.hpp index b48f0182adf..4018dcfb72d 100644 --- a/src/core/rotation.hpp +++ b/src/core/rotation.hpp @@ -32,6 +32,7 @@ #include "ParticleRange.hpp" #include +#include #include #include @@ -91,8 +92,44 @@ convert_dip_to_quat(const Utils::Vector3d &dip) { #endif -/** Rotate the particle p around the NORMALIZED axis a by amount phi */ -void local_rotate_particle(Particle &p, const Utils::Vector3d &a, double phi); +/** Rotate the particle p around the body-frame defined NORMALIZED axis + * aBodyFrame by amount phi + */ +inline void local_rotate_particle_body(Particle &p, + const Utils::Vector3d &axis_body_frame, + const double phi) { + auto axis = axis_body_frame; + + // Rotation turned off entirely? + if (!p.p.rotation) + return; + + // Convert rotation axis to body-fixed frame + axis = mask(p.p.rotation, axis).normalize(); + + auto const s = std::sin(phi / 2); + auto const q = + Utils::Vector4d{cos(phi / 2), s * axis[0], s * axis[1], s * axis[2]} + .normalize(); + + // Rotate the particle + p.r.quat = Utils::multiply_quaternions(p.r.quat, q); +} + +/** Rotate the particle p around the NORMALIZED axis aSpaceFrame by amount phi + */ +inline void local_rotate_particle(Particle &p, + const Utils::Vector3d &axis_space_frame, + const double phi) { + // Convert rotation axis to body-fixed frame + Utils::Vector3d axis = convert_vector_space_to_body(p, axis_space_frame); + local_rotate_particle_body(p, axis, phi); +} + +inline void convert_torque_to_body_frame_apply_fix(Particle &p) { + auto const torque = convert_vector_space_to_body(p, p.f.torque); + p.f.torque = mask(p.p.rotation, torque); +} #endif // ROTATION #endif diff --git a/src/core/thermostat.cpp b/src/core/thermostat.cpp index 02437d68eda..ce1fc17833a 100644 --- a/src/core/thermostat.cpp +++ b/src/core/thermostat.cpp @@ -41,22 +41,12 @@ bool thermo_virtual = true; using Thermostat::GammaType; -namespace { -/** @name Langevin parameters sentinels. - * These functions return the sentinel value for the Langevin - * parameters, indicating that they have not been set yet. - */ -/*@{*/ -constexpr double sentinel(double) { return -1.0; } -Utils::Vector3d sentinel(Utils::Vector3d) { return {-1.0, -1.0, -1.0}; } -/*@}*/ -} // namespace - GammaType langevin_gamma = sentinel(GammaType{}); GammaType langevin_gamma_rotation = sentinel(GammaType{}); GammaType langevin_pref1; GammaType langevin_pref2; GammaType langevin_pref2_rotation; +BrownianThermostat brownian = {}; /* NPT ISOTROPIC THERMOSTAT */ double nptiso_gamma0 = 0.0; @@ -86,27 +76,54 @@ void mpi_bcast_langevin_rng_counter_slave(const uint64_t 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(); } + void thermo_init_langevin() { langevin_pref1 = -langevin_gamma; langevin_pref2 = sqrt(24.0 * temperature / time_step * langevin_gamma); @@ -134,6 +151,47 @@ void thermo_init_npt_isotropic() { } #endif +/** Brownian thermostat initializer. + * + * Default particle mass is assumed to be unitary in these global parameters. + */ +void thermo_init_brownian() { + /** The heat velocity dispersion corresponds to the Gaussian noise only, + * which is only valid for the BD. Just a square root of kT, see (10.2.17) + * and comments in 2 paragraphs afterwards, @cite Pottier2010. + */ + brownian.sigma_vel = sqrt(temperature); + /** The random walk position dispersion is defined by the second eq. (14.38) + * of @cite Schlick2010. Its time interval factor will be added in the + * Brownian Dynamics functions. Its square root is the standard deviation. + */ + if (temperature > 0.0) { + brownian.sigma_pos_inv = sqrt(brownian.gamma / (2.0 * temperature)); + } else { + brownian.sigma_pos_inv = + brownian.gammatype_nan; // just an indication of the infinity + } +#ifdef ROTATION + /** Note: the BD thermostat assigns the brownian viscous parameters as well. + * They correspond to the friction tensor Z from the eq. (14.31) of + * @cite Schlick2010. + */ + /* If gamma_rotation is not set explicitly, + we use the translational one. */ + if (brownian.gamma_rotation < GammaType{}) { + brownian.gamma_rotation = brownian.gamma; + } + brownian.sigma_vel_rotation = sqrt(temperature); + if (temperature > 0.0) { + brownian.sigma_pos_rotation_inv = + sqrt(brownian.gamma_rotation / (2.0 * temperature)); + } else { + brownian.sigma_pos_rotation_inv = + brownian.gammatype_nan; // just an indication of the infinity + } +#endif // ROTATION +} + void thermo_init() { // Init thermalized bond despite of thermostat if (n_thermalized_bonds) { @@ -152,4 +210,6 @@ void thermo_init() { if (thermo_switch & THERMO_NPT_ISO) thermo_init_npt_isotropic(); #endif + if (thermo_switch & THERMO_BROWNIAN) + thermo_init_brownian(); } diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index 8e85e2182a2..c15c4246896 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -45,6 +45,7 @@ #define THERMO_DPD 2 #define THERMO_NPT_ISO 4 #define THERMO_LB 8 +#define THERMO_BROWNIAN 16 /*@}*/ namespace Thermostat { @@ -58,6 +59,21 @@ using GammaType = double; #endif } // namespace Thermostat +namespace { +/** @name Integrators parameters sentinels. + * These functions return the sentinel value for the Langevin/Brownian + * parameters, indicating that they have not been set yet. + */ +/*@{*/ +constexpr double sentinel(double) { return -1.0; } +constexpr Utils::Vector3d sentinel(Utils::Vector3d) { + return {-1.0, -1.0, -1.0}; +} +constexpr double set_nan(double) { return NAN; } +constexpr Utils::Vector3d set_nan(Utils::Vector3d) { return {NAN, NAN, NAN}; } +/*@}*/ +} // namespace + /************************************************ * exported variables ************************************************/ @@ -93,6 +109,40 @@ extern double nptiso_gammav; /** Langevin RNG counter, used for both translation and rotation. */ extern std::unique_ptr> langevin_rng_counter; +/** %Thermostat for Brownian dynamics. */ +struct BrownianThermostat { +private: + using GammaType = Thermostat::GammaType; + +public: + /** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */ + GammaType gamma = sentinel(GammaType{}); + /** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */ + GammaType gamma_rotation = sentinel(GammaType{}); + /** Inverse of the translational noise standard deviation. + * Stores @f$ \left(\sqrt{2D_{\text{trans}}}\right)^{-1} @f$ with + * @f$ D_{\text{trans}} = k_B T/\gamma_{\text{trans}} @f$ + * the translational diffusion coefficient + */ + GammaType sigma_pos_inv = sentinel(GammaType{}); + /** Inverse of the rotational noise standard deviation. + * Stores @f$ \left(\sqrt{2D_{\text{rot}}}\right)^{-1} @f$ with + * @f$ D_{\text{rot}} = k_B T/\gamma_{\text{rot}} @f$ + * the rotational diffusion coefficient + */ + GammaType sigma_pos_rotation_inv = sentinel(GammaType{}); + /** Sentinel value for divisions by zero. */ + GammaType const gammatype_nan = set_nan(GammaType{}); + /** Translational velocity noise standard deviation. */ + double sigma_vel = 0; + /** Angular velocity noise standard deviation. */ + double sigma_vel_rotation = 0; + /** RNG counter, used for both translation and rotation. */ + std::unique_ptr> rng_counter; +}; + +extern BrownianThermostat brownian; + /************************************************ * functions ************************************************/ @@ -100,11 +150,17 @@ extern std::unique_ptr> langevin_rng_counter; /** Only require seed if rng is not initialized. */ bool langevin_is_seed_required(); +/** 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(); /*@}*/ /** Initialize constants of the thermostat at the start of integration */ diff --git a/src/python/espressomd/globals.pxd b/src/python/espressomd/globals.pxd index d4f0b8b07dd..95ab83edaa7 100644 --- a/src/python/espressomd/globals.pxd +++ b/src/python/espressomd/globals.pxd @@ -36,9 +36,11 @@ cdef extern from "global.hpp": int FIELD_THERMO_VIRTUAL int FIELD_TEMPERATURE int FIELD_LANGEVIN_GAMMA + int FIELD_BROWNIAN_GAMMA int FIELD_SWIMMING_PARTICLES_EXIST IF ROTATION: int FIELD_LANGEVIN_GAMMA_ROTATION + int FIELD_BROWNIAN_GAMMA_ROTATION IF NPT: int FIELD_NPTISO_G0 int FIELD_NPTISO_GV diff --git a/src/python/espressomd/integrate.pxd b/src/python/espressomd/integrate.pxd index a05d8cc6f15..de897544cc1 100644 --- a/src/python/espressomd/integrate.pxd +++ b/src/python/espressomd/integrate.pxd @@ -31,6 +31,7 @@ cdef extern from "integrate.hpp" nogil: cdef int integrate_set_steepest_descent(const double f_max, const double gamma, const int max_steps, const double max_displacement) cdef extern cbool skin_set + cdef void integrate_set_bd() IF NPT: cdef extern from "integrate.hpp" nogil: diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx index 3f5d4bf4afa..8912f468bc5 100644 --- a/src/python/espressomd/integrate.pyx +++ b/src/python/espressomd/integrate.pyx @@ -94,6 +94,13 @@ cdef class IntegratorHandle: """ self._integrator = VelocityVerletIsotropicNPT(*args, **kwargs) + def set_brownian_dynamics(self): + """ + Set the integration method to BD. + + """ + self._integrator = BrownianDynamics() + cdef class Integrator: """ @@ -371,3 +378,31 @@ ELSE: cdef class VelocityVerletIsotropicNPT(Integrator): def __init__(self, *args, **kwargs): raise Exception("NPT not compiled in.") + + +cdef class BrownianDynamics(Integrator): + """ + Brownian Dynamics integrator. + + """ + + def default_params(self): + return {} + + def valid_keys(self): + """All parameters that can be set. + + """ + return {} + + def required_keys(self): + """Parameters that have to be set. + + """ + return {} + + def validate_params(self): + return True + + def _set_params_in_es_core(self): + integrate_set_bd() diff --git a/src/python/espressomd/thermostat.pxd b/src/python/espressomd/thermostat.pxd index c3ece48ea28..860cae7cd1f 100644 --- a/src/python/espressomd/thermostat.pxd +++ b/src/python/espressomd/thermostat.pxd @@ -31,6 +31,7 @@ cdef extern from "thermostat.hpp": int THERMO_LB int THERMO_NPT_ISO int THERMO_DPD + int THERMO_BROWNIAN IF PARTICLE_ANISOTROPY: Vector3d langevin_gamma_rotation @@ -39,10 +40,25 @@ cdef extern from "thermostat.hpp": double langevin_gamma_rotation double langevin_gamma + IF PARTICLE_ANISOTROPY: + ctypedef struct brownian_thermostat_struct "BrownianThermostat": + Vector3d gamma_rotation + Vector3d gamma + ELSE: + ctypedef struct brownian_thermostat_struct "BrownianThermostat": + double gamma_rotation + double gamma + + # links intern C-struct with python object + cdef extern brownian_thermostat_struct brownian + void langevin_set_rng_state(stdint.uint64_t counter) + void brownian_set_rng_state(stdint.uint64_t counter) cbool langevin_is_seed_required() + cbool brownian_is_seed_required() stdint.uint64_t langevin_get_rng_state() + stdint.uint64_t brownian_get_rng_state() cdef extern from "dpd.hpp": IF DPD: diff --git a/src/python/espressomd/thermostat.pyx b/src/python/espressomd/thermostat.pyx index 068808845c0..e2b0798be8e 100644 --- a/src/python/espressomd/thermostat.pyx +++ b/src/python/espressomd/thermostat.pyx @@ -111,6 +111,11 @@ cdef class Thermostat: gammav=thmst["gammav"]) if thmst["type"] == "DPD": self.set_dpd(kT=thmst["kT"], seed=thmst["seed"]) + if thmst["type"] == "BROWNIAN": + self.set_brownian(kT=thmst["kT"], gamma=thmst["gamma"], + gamma_rotation=thmst["gamma_rotation"], + act_on_virtual=thmst["act_on_virtual"], + seed=thmst["seed"]) def get_ts(self): return thermo_switch @@ -145,7 +150,29 @@ cdef class Thermostat: lang_dict["gamma_rotation"] = None thermo_list.append(lang_dict) + if thermo_switch & THERMO_BROWNIAN: + lang_dict = {} + lang_dict["type"] = "BROWNIAN" + lang_dict["kT"] = temperature + lang_dict["act_on_virtual"] = thermo_virtual + lang_dict["seed"] = int(brownian_get_rng_state()) + IF PARTICLE_ANISOTROPY: + lang_dict["gamma"] = [brownian.gamma[0], + brownian.gamma[1], + brownian.gamma[2]] + ELSE: + lang_dict["gamma"] = brownian.gamma + IF ROTATION: + IF PARTICLE_ANISOTROPY: + lang_dict["gamma_rotation"] = [brownian.gamma_rotation[0], + brownian.gamma_rotation[1], + brownian.gamma_rotation[2]] + ELSE: + lang_dict["gamma_rotation"] = brownian.gamma_rotation + ELSE: + lang_dict["gamma_rotation"] = None + thermo_list.append(lang_dict) if thermo_switch & THERMO_LB: lb_dict = {} lb_dict["LB_fluid"] = self._LB_fluid @@ -190,6 +217,13 @@ cdef class Thermostat: ELSE: langevin_gamma = 0. mpi_bcast_parameter(FIELD_LANGEVIN_GAMMA) + global brownian + IF PARTICLE_ANISOTROPY: + for i in range(3): + brownian.gamma[i] = 0. + ELSE: + brownian.gamma = 0. + mpi_bcast_parameter(FIELD_BROWNIAN_GAMMA) global langevin_gamma_rotation IF ROTATION: IF PARTICLE_ANISOTROPY: @@ -198,6 +232,13 @@ cdef class Thermostat: ELSE: langevin_gamma_rotation = 0. mpi_bcast_parameter(FIELD_LANGEVIN_GAMMA_ROTATION) + IF ROTATION: + IF PARTICLE_ANISOTROPY: + for i in range(3): + brownian.gamma_rotation[i] = 0. + ELSE: + brownian.gamma_rotation = 0. + mpi_bcast_parameter(FIELD_BROWNIAN_GAMMA_ROTATION) global thermo_switch thermo_switch = THERMO_OFF @@ -358,6 +399,156 @@ cdef class Thermostat: mpi_bcast_parameter(FIELD_LANGEVIN_GAMMA_ROTATION) return True + @AssertThermostatType(THERMO_BROWNIAN) + def set_brownian(self, kT=None, gamma=None, gamma_rotation=None, + act_on_virtual=False, seed=None): + """Sets the Brownian thermostat. + + Parameters + ----------- + kT : :obj:`float` + Thermal energy of the simulated heat bath. + gamma : :obj:`float` + Contains the friction coefficient of the bath. If the feature + ``PARTICLE_ANISOTROPY`` is compiled in, then ``gamma`` can be a list + of three positive floats, for the friction coefficient in each + cardinal direction. + gamma_rotation : :obj:`float`, optional + The same applies to ``gamma_rotation``, which requires the feature + ``ROTATION`` to work properly. But also accepts three floats + if ``PARTICLE_ANISOTROPY`` is also compiled in. + act_on_virtual : :obj:`bool`, optional + If ``True`` the thermostat will act on virtual sites, default is + ``False``. + seed : :obj:`int` + Initial counter value (or seed) of the philox RNG. + Required on first activation of the Brownian thermostat. + Must be positive. + + """ + + scalar_gamma_def = True + scalar_gamma_rot_def = True + IF PARTICLE_ANISOTROPY: + if hasattr(gamma, "__iter__"): + scalar_gamma_def = False + else: + scalar_gamma_def = True + + IF PARTICLE_ANISOTROPY: + if hasattr(gamma_rotation, "__iter__"): + scalar_gamma_rot_def = False + else: + scalar_gamma_rot_def = True + + if kT is None or gamma is None: + raise ValueError( + "Both, kT and gamma have to be given as keyword args") + utils.check_type_or_throw_except(kT, 1, float, "kT must be a number") + if scalar_gamma_def: + utils.check_type_or_throw_except( + gamma, 1, float, "gamma must be a number") + else: + utils.check_type_or_throw_except( + gamma, 3, float, "diagonal elements of the gamma tensor must be numbers") + if gamma_rotation is not None: + if scalar_gamma_rot_def: + utils.check_type_or_throw_except( + gamma_rotation, 1, float, "gamma_rotation must be a number") + else: + utils.check_type_or_throw_except( + gamma_rotation, 3, float, "diagonal elements of the gamma_rotation tensor must be numbers") + + if scalar_gamma_def: + if float(kT) < 0. or float(gamma) < 0.: + raise ValueError( + "temperature and gamma must be positive numbers") + else: + if float(kT) < 0. or float(gamma[0]) < 0. or float( + gamma[1]) < 0. or float(gamma[2]) < 0.: + raise ValueError( + "temperature and diagonal elements of the gamma tensor must be positive numbers") + if gamma_rotation is not None: + if scalar_gamma_rot_def: + if float(gamma_rotation) < 0.: + raise ValueError( + "gamma_rotation must be positive number") + else: + if float(gamma_rotation[0]) < 0. or float( + gamma_rotation[1]) < 0. or float(gamma_rotation[2]) < 0.: + raise ValueError( + "diagonal elements of the gamma_rotation tensor must be positive numbers") + + # Seed is required if the RNG is not initialized + if seed is None and brownian_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") + brownian_set_rng_state(seed) + + global temperature + temperature = float(kT) + global brownian + IF PARTICLE_ANISOTROPY: + if scalar_gamma_def: + brownian.gamma[0] = gamma + brownian.gamma[1] = gamma + brownian.gamma[2] = gamma + else: + brownian.gamma[0] = gamma[0] + brownian.gamma[1] = gamma[1] + brownian.gamma[2] = gamma[2] + ELSE: + brownian.gamma = float(gamma) + IF ROTATION: + if gamma_rotation is not None: + IF PARTICLE_ANISOTROPY: + if scalar_gamma_rot_def: + brownian.gamma_rotation[0] = gamma_rotation + brownian.gamma_rotation[1] = gamma_rotation + brownian.gamma_rotation[2] = gamma_rotation + else: + brownian.gamma_rotation[0] = gamma_rotation[0] + brownian.gamma_rotation[1] = gamma_rotation[1] + brownian.gamma_rotation[2] = gamma_rotation[2] + ELSE: + if scalar_gamma_rot_def: + brownian.gamma_rotation = gamma_rotation + else: + raise ValueError( + "gamma_rotation must be a scalar since feature PARTICLE_ANISOTROPY is disabled") + else: + IF PARTICLE_ANISOTROPY: + if scalar_gamma_def: + brownian.gamma_rotation[0] = gamma + brownian.gamma_rotation[1] = gamma + brownian.gamma_rotation[2] = gamma + else: + brownian.gamma_rotation[0] = gamma[0] + brownian.gamma_rotation[1] = gamma[1] + brownian.gamma_rotation[2] = gamma[2] + ELSE: + brownian.gamma_rotation = brownian.gamma + + global thermo_switch + thermo_switch = (thermo_switch | THERMO_BROWNIAN) + mpi_bcast_parameter(FIELD_THERMO_SWITCH) + mpi_bcast_parameter(FIELD_TEMPERATURE) + mpi_bcast_parameter(FIELD_BROWNIAN_GAMMA) + + global thermo_virtual + thermo_virtual = act_on_virtual + mpi_bcast_parameter(FIELD_THERMO_VIRTUAL) + + IF ROTATION: + mpi_bcast_parameter(FIELD_BROWNIAN_GAMMA_ROTATION) + return True + @AssertThermostatType(THERMO_LB, THERMO_DPD) def set_lb( self, @@ -461,7 +652,6 @@ cdef class Thermostat: Must be positive. """ - if kT is None: raise ValueError("kT has to be given as keyword args") if not isinstance(kT, float): diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index a42196b2e5a..fef42948fb3 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -48,7 +48,7 @@ endfunction(PYTHON_TEST) # Checkpointing tests: semicolon-separated list of mutually-compatible features. # Separate features with hyphens, use a period to add an optional flag. -foreach(TEST_COMBINATION lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt-minimization;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd) +foreach(TEST_COMBINATION lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt-minimization;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd;lb.off-therm.bd-int.bd) if(${TEST_COMBINATION} MATCHES "\\.gpu") set(TEST_LABELS "gpu") else() @@ -113,6 +113,7 @@ python_test(FILE ek_eof_one_species_y.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE ek_eof_one_species_z.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE exclusions.py MAX_NUM_PROC 2) python_test(FILE langevin_thermostat.py MAX_NUM_PROC 1) +python_test(FILE brownian_dynamics.py MAX_NUM_PROC 1) python_test(FILE nsquare.py MAX_NUM_PROC 4) python_test(FILE virtual_sites_relative.py MAX_NUM_PROC 2) python_test(FILE virtual_sites_tracers.py MAX_NUM_PROC 2) diff --git a/testsuite/python/brownian_dynamics.py b/testsuite/python/brownian_dynamics.py new file mode 100644 index 00000000000..8988f9c5e36 --- /dev/null +++ b/testsuite/python/brownian_dynamics.py @@ -0,0 +1,303 @@ +# +# Copyright (C) 2013-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 . +# +import unittest as ut +import unittest_decorators as utx +import espressomd +import numpy as np +from espressomd.accumulators import Correlator, TimeSeries +from espressomd.observables import ParticleVelocities, ParticleBodyAngularVelocities, ParticlePositions +from tests_common import single_component_maxwell + + +class BrownianDynamics(ut.TestCase): + + """Tests velocity distributions and diffusion for Brownian Dynamics""" + system = espressomd.System(box_l=[1.0, 1.0, 1.0]) + system.cell_system.set_domain_decomposition(use_verlet_lists=True) + system.cell_system.skin = 0 + system.seed = range(system.cell_system.get_state()["n_nodes"]) + system.periodicity = [0, 0, 0] + system.integrator.set_brownian_dynamics() + + @classmethod + def setUpClass(cls): + np.random.seed(42) + + def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT): + """check the recorded particle distributions in velocity against a + histogram with n_bins bins. Drop velocities outside minmax. Check + individual histogram bins up to an accuracy of error_tol against the + analytical result for kT.""" + for i in range(3): + hist = np.histogram( + vel[:, i], range=(-minmax, minmax), bins=n_bins, density=False) + data = hist[0] / float(vel.shape[0]) + bins = hist[1] + for j in range(n_bins): + found = data[j] + expected = single_component_maxwell(bins[j], bins[j + 1], kT) + self.assertLessEqual(abs(found - expected), error_tol) + + def test_00_verify_single_component_maxwell(self): + """Verifies the normalization of the analytical expression.""" + self.assertLessEqual( + abs(single_component_maxwell(-10, 10, 4.) - 1.), 1E-4) + + def check_vel_dist_global_temp(self, recalc_forces, loops): + """Test velocity distribution for global temperature parameters. + + Parameters + ---------- + recalc_forces : :obj:`bool` + True if the forces should be recalculated after every step. + loops : :obj:`int` + Number of sampling loops + """ + N = 200 + system = self.system + system.part.clear() + system.time_step = 1.6 + + # Place particles + system.part.add(pos=np.random.random((N, 3))) + + # Enable rotation if compiled in + if espressomd.has_features("ROTATION"): + system.part[:].rotation = [1, 1, 1] + + kT = 1.1 + gamma = 3.5 + system.thermostat.set_brownian(kT=kT, gamma=gamma, seed=41) + + # Warmup + system.integrator.run(20) + + # Sampling + v_stored = np.zeros((N * loops, 3)) + omega_stored = np.zeros((N * loops, 3)) + for i in range(loops): + system.integrator.run(1, recalc_forces=recalc_forces) + v_stored[i * N:(i + 1) * N, :] = system.part[:].v + if espressomd.has_features("ROTATION"): + omega_stored[i * N:(i + 1) * N, :] = system.part[:].omega_body + + v_minmax = 5 + bins = 4 + error_tol = 0.01 + self.check_velocity_distribution( + v_stored, v_minmax, bins, error_tol, kT) + if espressomd.has_features("ROTATION"): + self.check_velocity_distribution( + omega_stored, v_minmax, bins, error_tol, kT) + + def test_vel_dist_global_temp(self): + """Test velocity distribution for global temperature.""" + self.check_vel_dist_global_temp(False, loops=1500) + + def test_vel_dist_global_temp_initial_forces(self): + """Test velocity distribution for global Brownian parameters, + when using the initial force calculation. + """ + self.check_vel_dist_global_temp(True, loops=170) + + @utx.skipIfMissingFeatures("BROWNIAN_PER_PARTICLE") + def test_05__brownian_per_particle(self): + """Test Brownian dynamics with particle specific kT and gamma. Covers all combinations of + particle specific gamma and temp set or not set. + """ + N = 400 + system = self.system + system.part.clear() + system.time_step = 1.9 + system.part.add(pos=np.random.random((N, 3))) + if espressomd.has_features("ROTATION"): + system.part[:].rotation = [1, 1, 1] + + kT = 0.9 + gamma = 3.2 + gamma2 = 14.3 + kT2 = 1.5 + system.thermostat.set_brownian(kT=kT, gamma=gamma, seed=41) + # Set different kT on 2nd half of particles + system.part[int(N / 2):].temp = kT2 + # Set different gamma on half of the particles (overlap over both kTs) + if espressomd.has_features("PARTICLE_ANISOTROPY"): + system.part[int(N / 4):int(3 * N / 4)].gamma = 3 * [gamma2] + else: + system.part[int(N / 4):int(3 * N / 4)].gamma = gamma2 + + system.integrator.run(50) + loops = 300 + + v_kT = np.zeros((int(N / 2) * loops, 3)) + v_kT2 = np.zeros((int(N / 2 * loops), 3)) + + if espressomd.has_features("ROTATION"): + omega_kT = np.zeros((int(N / 2) * loops, 3)) + omega_kT2 = np.zeros((int(N / 2 * loops), 3)) + + for i in range(loops): + system.integrator.run(1) + v_kT[int(i * N / 2):int((i + 1) * N / 2), + :] = system.part[:int(N / 2)].v + v_kT2[int(i * N / 2):int((i + 1) * N / 2), + :] = system.part[int(N / 2):].v + + if espressomd.has_features("ROTATION"): + omega_kT[int(i * N / 2):int((i + 1) * N / 2), :] = \ + system.part[:int(N / 2)].omega_body + omega_kT2[int(i * N / 2):int((i + 1) * N / 2), :] = \ + system.part[int(N / 2):].omega_body + v_minmax = 5 + bins = 4 + error_tol = 0.012 + self.check_velocity_distribution(v_kT, v_minmax, bins, error_tol, kT) + self.check_velocity_distribution(v_kT2, v_minmax, bins, error_tol, kT2) + + if espressomd.has_features("ROTATION"): + self.check_velocity_distribution( + omega_kT, v_minmax, bins, error_tol, kT) + self.check_velocity_distribution( + omega_kT2, v_minmax, bins, error_tol, kT2) + + def setup_diff_mass_rinertia(self, p): + if espressomd.has_features("MASS"): + p.mass = 0.5 + if espressomd.has_features("ROTATION"): + p.rotation = [1, 1, 1] + # Make sure rinertia does not change diff coeff + if espressomd.has_features("ROTATIONAL_INERTIA"): + p.rinertia = [0.4, 0.4, 0.4] + + def test_msd_global_temp(self): + """Tests diffusion via MSD for global gamma and temeprature""" + + gamma = 9.4 + kT = 0.37 + dt = 0.5 + + system = self.system + system.part.clear() + p = system.part.add(pos=(0, 0, 0), id=0) + system.time_step = dt + system.thermostat.set_brownian(kT=kT, gamma=gamma, seed=42) + system.cell_system.skin = 0.4 + + pos_obs = ParticlePositions(ids=(p.id,)) + + c_pos = Correlator(obs1=pos_obs, tau_lin=16, tau_max=100., delta_N=10, + corr_operation="square_distance_componentwise", + compress1="discard1") + system.auto_update_accumulators.add(c_pos) + + system.integrator.run(500000) + + c_pos.finalize() + + # Check MSD + msd = c_pos.result() + system.auto_update_accumulators.clear() + + def expected_msd(x): + return 2. * kT / gamma * x + + for i in range(2, 6): + np.testing.assert_allclose( + msd[i, 2:5], expected_msd(msd[i, 0]), rtol=0.02) + + @utx.skipIfMissingFeatures("VIRTUAL_SITES") + def test_07__virtual(self): + system = self.system + system.time_step = 0.01 + system.part.clear() + + virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) + physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + + system.thermostat.set_brownian( + kT=0, gamma=1, gamma_rotation=1., act_on_virtual=False, seed=41) + + system.integrator.run(1) + + np.testing.assert_almost_equal(np.copy(virtual.v), [1, 0, 0]) + np.testing.assert_almost_equal(np.copy(physical.v), [0, 0, 0]) + + system.part.clear() + virtual = system.part.add(pos=[0, 0, 0], virtual=True, v=[1, 0, 0]) + physical = system.part.add(pos=[0, 0, 0], virtual=False, v=[1, 0, 0]) + + system.thermostat.set_brownian( + kT=0, gamma=1, gamma_rotation=1., act_on_virtual=True, seed=41) + system.integrator.run(1) + + np.testing.assert_almost_equal(np.copy(virtual.v), [0, 0, 0]) + np.testing.assert_almost_equal(np.copy(physical.v), [0, 0, 0]) + + def test_08__noise_correlation(self): + """Checks that the Brownian noise is uncorrelated""" + + system = self.system + system.part.clear() + system.time_step = 0.01 + system.cell_system.skin = 0.1 + system.part.add(id=(0, 1), pos=np.zeros((2, 3))) + vel_obs = ParticleVelocities(ids=system.part[:].id) + vel_series = TimeSeries(obs=vel_obs) + system.auto_update_accumulators.add(vel_series) + if espressomd.has_features("ROTATION"): + system.part[:].rotation = (1, 1, 1) + omega_obs = ParticleBodyAngularVelocities(ids=system.part[:].id) + omega_series = TimeSeries(obs=omega_obs) + system.auto_update_accumulators.add(omega_series) + + kT = 3.2 + system.thermostat.set_brownian(kT=kT, gamma=2.1, seed=17) + steps = int(1e6) + system.integrator.run(steps) + system.auto_update_accumulators.clear() + + # test translational noise correlation + vel = np.array(vel_series.time_series()) + for i in range(6): + for j in range(i, 6): + corrcoef = np.dot(vel[:, i], vel[:, j]) / steps / kT + if i == j: + self.assertAlmostEqual(corrcoef, 1.0, delta=0.04) + else: + self.assertLessEqual(np.abs(corrcoef), 0.04) + + # test rotational noise correlation + if espressomd.has_features("ROTATION"): + omega = np.array(omega_series.time_series()) + for i in range(6): + for j in range(6): + corrcoef = np.dot(omega[:, i], omega[:, j]) / steps / kT + if i == j: + self.assertAlmostEqual(corrcoef, 1.0, delta=0.04) + else: + self.assertLessEqual(np.abs(corrcoef), 0.04) + # translational and angular velocities should be independent + for i in range(6): + for j in range(6): + corrcoef = np.dot(vel[:, i], omega[:, j]) / steps / kT + self.assertLessEqual(np.abs(corrcoef), 0.04) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/rotational-diffusion-aniso.py b/testsuite/python/rotational-diffusion-aniso.py index f5dafbf28f5..9492e87248b 100644 --- a/testsuite/python/rotational-diffusion-aniso.py +++ b/testsuite/python/rotational-diffusion-aniso.py @@ -22,10 +22,11 @@ import tests_common -@utx.skipIfMissingFeatures(["PARTICLE_ANISOTROPY", +@utx.skipIfMissingFeatures(["ROTATION", "PARTICLE_ANISOTROPY", "ROTATIONAL_INERTIA", "DIPOLES"]) class RotDiffAniso(ut.TestCase): longMessage = True + round_error_prec = 1E-14 # Handle for espresso system system = espressomd.System(box_l=[1.0, 1.0, 1.0]) system.cell_system.skin = 5.0 @@ -43,11 +44,15 @@ class RotDiffAniso(ut.TestCase): def setUp(self): self.system.time = 0.0 self.system.part.clear() + if "BROWNIAN_DYNAMICS" in espressomd.features(): + self.system.thermostat.turn_off() + # the default integrator is supposed implicitly + self.system.integrator.set_nvt() - def rot_diffusion_param_setup(self, n): + def add_particles_setup(self, n): """ - Setup the parameters for the rotational diffusion - test check_rot_diffusion(). + Adding particles according to the + previously set parameters. Parameters ---------- @@ -56,18 +61,21 @@ def rot_diffusion_param_setup(self, n): """ - # Time - # The time step should be less than t0 ~ mass / gamma - self.system.time_step = 3E-3 + for ind in range(n): + part_pos = np.random.random(3) * self.box + self.system.part.add(rotation=(1, 1, 1), id=ind, + pos=part_pos) + self.system.part[ind].rinertia = self.J + if espressomd.has_features("ROTATION"): + self.system.part[ind].omega_body = [0.0, 0.0, 0.0] - # Space - box = 10.0 - self.system.box_l = box, box, box - self.system.periodicity = [0, 0, 0] + def set_anisotropic_param(self): + """ + Select parameters for anisotropic particles. + + """ # NVT thermostat - # Just some temperature range to cover by the test: - self.kT = np.random.uniform(1.0, 1.5) # Note: here & hereinafter specific variations in the random parameter # ranges are related to the test execution duration to achieve the # required statistical averages faster. The friction gamma_global should @@ -84,7 +92,7 @@ def rot_diffusion_param_setup(self, n): # eq. (10.2.26) [N. Pottier, doi:10.1007/s10955-010-0114-6 (2010)]. self.gamma_global = 1E2 * np.random.uniform(0.35, 1.05, (3)) - # Particles + # Particles' properties # As far as the problem characteristic time is t0 ~ J / gamma # and the Langevin equation finite-difference approximation is stable # only for time_step << t0, it is needed to set the moment of inertia @@ -94,11 +102,46 @@ def rot_diffusion_param_setup(self, n): # too much of the CPU time: the in silico time should clock over the # t0. self.J = np.random.uniform(1.5, 16.5, (3)) - for ind in range(n): - part_pos = np.random.random(3) * box - self.system.part.add(rotation=(1, 1, 1), id=ind, rinertia=self.J, - pos=part_pos) - self.system.part[ind].omega_body = [0.0, 0.0, 0.0] + + def set_isotropic_param(self): + """ + Select parameters for isotropic particles. + + Parameters + ---------- + + """ + + # NVT thermostat + # see the comments in set_anisotropic_param() + self.gamma_global[0] = 1E2 * np.random.uniform(0.35, 1.05) + self.gamma_global[1] = self.gamma_global[0] + self.gamma_global[2] = self.gamma_global[0] + # Particles' properties + # see the comments in set_anisotropic_param() + self.J[0] = np.random.uniform(1.5, 16.5) + self.J[1] = self.J[0] + self.J[2] = self.J[0] + + def rot_diffusion_param_setup(self): + """ + Setup the parameters for the rotational diffusion + test check_rot_diffusion(). + + """ + + # Time + # The time step should be less than t0 ~ mass / gamma + self.system.time_step = 3E-3 + + # Space + self.box = 10.0 + self.system.box_l = 3 * [self.box] + self.system.periodicity = [0, 0, 0] + + # NVT thermostat + # Just some temperature range to cover by the test: + self.kT = np.random.uniform(0.5, 1.5) def check_rot_diffusion(self, n): """ @@ -239,6 +282,11 @@ def check_rot_diffusion(self, n): if i != j: D1D1 += D[i] * D[j] D1D1 /= 6.0 + # Technical workaround of a digital arithmetic issue for isotropic + # particle + if np.absolute((D0**2 - D1D1) / (D0**2 + D1D1) + ) < self.round_error_prec: + D1D1 *= (1.0 - 2.0 * self.round_error_prec) # Eq. (32) [Perrin1936]. dcosjj2_validate = 1. / 3. + (1. / 3.) * (1. + (D - D0) / (2. * np.sqrt(D0**2 - D1D1))) \ * np.exp(-6. * (D0 - np.sqrt(D0**2 - D1D1)) * self.system.time) \ @@ -302,14 +350,42 @@ def check_rot_diffusion(self, n): 'rotational diffusion is too large: {2}' .format(i, j, dcosij2_dev[i, j])) + # Langevin Dynamics / Anisotropic def test_case_00(self): - n = 300 - self.rot_diffusion_param_setup(n) + n = 800 + self.rot_diffusion_param_setup() + self.set_anisotropic_param() + self.add_particles_setup(n) self.system.thermostat.set_langevin( kT=self.kT, gamma=self.gamma_global, seed=42) # Actual integration and validation run self.check_rot_diffusion(n) + # Langevin Dynamics / Isotropic + def test_case_01(self): + n = 800 + self.rot_diffusion_param_setup() + self.set_isotropic_param() + self.add_particles_setup(n) + self.system.thermostat.set_langevin( + kT=self.kT, gamma=self.gamma_global, seed=42) + # Actual integration and validation run + self.check_rot_diffusion(n) + + if "BROWNIAN_DYNAMICS" in espressomd.features(): + # Brownian Dynamics / Isotropic + def test_case_10(self): + n = 800 + self.system.thermostat.turn_off() + self.rot_diffusion_param_setup() + self.set_isotropic_param() + self.add_particles_setup(n) + self.system.thermostat.set_brownian( + kT=self.kT, gamma=self.gamma_global, seed=42) + self.system.integrator.set_brownian_dynamics() + # Actual integration and validation run + self.check_rot_diffusion(n) + if __name__ == '__main__': ut.main() diff --git a/testsuite/python/rotational_inertia.py b/testsuite/python/rotational_inertia.py index f3f08d8bd45..31779801d42 100644 --- a/testsuite/python/rotational_inertia.py +++ b/testsuite/python/rotational_inertia.py @@ -32,10 +32,6 @@ class RotationalInertia(ut.TestCase): L_0_lab = np.zeros((3)) L_lab = np.zeros((3)) - def convert_vec_body_to_space(self, part, vec): - A = tests_common.rotation_matrix_quat(self.system, part) - return np.dot(A.transpose(), vec) - # Angular momentum def L_body(self, part): return self.system.part[part].omega_body[:] * \ @@ -44,11 +40,13 @@ def L_body(self, part): # Set the angular momentum def set_L_0(self, part): L_0_body = self.L_body(part) - self.L_0_lab = self.convert_vec_body_to_space(part, L_0_body) + self.L_0_lab = tests_common.convert_vec_body_to_space( + self.system, part, L_0_body) def set_L(self, part): L_body = self.L_body(part) - self.L_lab = self.convert_vec_body_to_space(part, L_body) + self.L_lab = tests_common.convert_vec_body_to_space( + self.system, part, L_body) def test_stability(self): self.system.part.clear() diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py index deb226d42f7..377e3018436 100644 --- a/testsuite/python/save_checkpoint.py +++ b/testsuite/python/save_checkpoint.py @@ -138,6 +138,8 @@ # set thermostat if 'THERM.LANGEVIN' in modes: system.thermostat.set_langevin(kT=1.0, gamma=2.0, seed=42) + 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) elif 'THERM.DPD' in modes and has_features('DPD'): @@ -151,6 +153,8 @@ max_displacement=0.01) elif 'INT.NVT' in modes: system.integrator.set_nvt() + elif 'INT.BD' in modes: + system.integrator.set_brownian_dynamics() # set minimization if 'MINIMIZATION' in modes: steepest_descent(system, f_max=1, gamma=10, max_steps=0, diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index 7c4e623862d..d07d0660d88 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -157,7 +157,22 @@ def test_thermostat_Langevin(self): self.assertEqual(thmst['type'], 'LANGEVIN') self.assertEqual(thmst['kT'], 1.0) self.assertEqual(thmst['seed'], 42) - np.testing.assert_array_equal(thmst['gamma'], np.array(3 * [2.0])) + self.assertFalse(thmst['act_on_virtual']) + np.testing.assert_array_equal(thmst['gamma'], 3 * [2.0]) + if espressomd.has_features('ROTATION'): + np.testing.assert_array_equal(thmst['gamma_rotation'], 3 * [2.0]) + + @ut.skipIf('THERM.BD' not in modes, + 'Brownian thermostat not in modes') + def test_thermostat_Brownian(self): + thmst = system.thermostat.get_state()[0] + self.assertEqual(thmst['type'], 'BROWNIAN') + self.assertEqual(thmst['kT'], 1.0) + self.assertEqual(thmst['seed'], 42) + self.assertFalse(thmst['act_on_virtual']) + np.testing.assert_array_equal(thmst['gamma'], 3 * [2.0]) + if espressomd.has_features('ROTATION'): + np.testing.assert_array_equal(thmst['gamma_rotation'], 3 * [2.0]) @utx.skipIfMissingFeatures('DPD') @ut.skipIf('THERM.DPD' not in modes, 'DPD thermostat not in modes') @@ -210,6 +225,13 @@ def test_integrator_VV(self): params = integ.get_params() self.assertEqual(params, {}) + @ut.skipIf('INT.BD' not in modes, 'BD integrator not in modes') + def test_integrator_BD(self): + integ = system.integrator.get_state() + self.assertIsInstance(integ, espressomd.integrate.BrownianDynamics) + params = integ.get_params() + self.assertEqual(params, {}) + @utx.skipIfMissingFeatures('LENNARD_JONES') @ut.skipIf('LJ' not in modes, "Skipping test due to missing mode.") def test_non_bonded_inter(self): diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index 72782efcf22..ce015a733f7 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -167,6 +167,11 @@ def transform_vel_from_cartesian_to_polar_coordinates(pos, vel): (pos[0] * vel[1] - pos[1] * vel[0]) / (pos[0]**2 + pos[1]**2), vel[2]]) +def convert_vec_body_to_space(system, part, vec): + A = rotation_matrix_quat(system, part) + return np.dot(A.transpose(), vec) + + def rotation_matrix(axis, theta): """ Return the rotation matrix associated with counterclockwise rotation about