Brownian Dynamics within the existing VV loop

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},

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)},
 

doc/sphinx/system_setup.rst

@@ -267,25 +267,25 @@ 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.
 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

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.},

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

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

src/config/features.def

@@ -19,6 +19,7 @@ MASS
 EXCLUSIONS
 BOND_CONSTRAINT
 LANGEVIN_PER_PARTICLE
+BROWNIAN_PER_PARTICLE
 COLLISION_DETECTION
 METADYNAMICS
 NPT

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

src/core/global.cpp

@@ -168,7 +168,26 @@ const std::unordered_map<int, Datafield> 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;

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