diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst index 7c84c22f939..a6f699890d9 100644 --- a/doc/sphinx/advanced_methods.rst +++ b/doc/sphinx/advanced_methods.rst @@ -1894,7 +1894,7 @@ The call of ``add_reaction`` define the insertion :math:`\mathrm{\emptyset \to t Multiple reactions for the insertions of different types can be added to the same ``WidomInsertion`` instance. Measuring the excess chemical potential using the insertion method is done via calling ``widom.measure_excess_chemical_potential(0)``. If another particle insertion is defined, then the excess chemical potential for this insertion can be measured by calling ``widom.measure_excess_chemical_potential(1)``. -Be aware that the implemented method only works for the canonical ensemble. If the numbers of particles fluctuate (i.e. in a semi grand canonical simulation) one has to adapt the formulas from which the excess chemical potential is calculated! This is not implemented. Also in a isobaric-isothermal simulation (NPT) the corresponding formulas for the excess chemical potentials need to be adapted. This is not implemented. +Be aware that the implemented method only works for the canonical ensemble. If the numbers of particles fluctuate (i.e. in a semi grand canonical simulation) one has to adapt the formulas from which the excess chemical potential is calculated! This is not implemented. Also in a isobaric-isothermal simulation (NpT) the corresponding formulas for the excess chemical potentials need to be adapted. This is not implemented. The implementation can also deal with the simultaneous insertion of multiple particles and can therefore measure the change of excess free energy of multiple particles like e.g.: diff --git a/doc/sphinx/electrostatics.rst b/doc/sphinx/electrostatics.rst index b4417a44b03..d39e7db170e 100644 --- a/doc/sphinx/electrostatics.rst +++ b/doc/sphinx/electrostatics.rst @@ -318,7 +318,7 @@ MMM1D is used with:: where the prefactor :math:`C` is defined in Eqn. :eq:`coulomb_prefactor`. MMM1D Coulomb method for systems with periodicity (0 0 1). Needs the -nsquared cell system (see section :ref:`Cellsystems`). The first form sets parameters +N-squared cell system (see section :ref:`Cellsystems`). The first form sets parameters manually. The switch radius determines at which xy-distance the force calculation switches from the near to the far formula. The Bessel cutoff does not need to be specified as it is automatically determined from the @@ -336,7 +336,7 @@ MMM1D on GPU :class:`espressomd.electrostatics.MMM1DGPU` MMM1D is also available in a GPU implementation. Unlike its CPU -counterpart, it does not need the nsquared cell system. +counterpart, it does not need the N-squared cell system. :: diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index 8efbf28b043..8ca57ab431b 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -355,9 +355,9 @@ General features In addition, there are switches that enable additional features in the integrator or thermostat: -- ``NPT`` Enables an on-the-fly NPT integration scheme. +- ``NPT`` Enables an on-the-fly NpT integration scheme. - .. seealso:: :ref:`Isotropic NPT thermostat` + .. seealso:: :ref:`Isotropic NpT thermostat` - ``ENGINE`` diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst index a633689d304..d955488f03b 100644 --- a/doc/sphinx/introduction.rst +++ b/doc/sphinx/introduction.rst @@ -448,7 +448,7 @@ report so to the developers using the instructions in :ref:`Contributing`. +--------------------------------+------------------------+------------------+------------+ | Langevin Thermostat | Core | Core | Yes | +--------------------------------+------------------------+------------------+------------+ -| Isotropic NPT | Experimental | None | Yes | +| Isotropic NpT | Experimental | None | Yes | +--------------------------------+------------------------+------------------+------------+ | Quaternion Integrator | Core | Good | Yes | +--------------------------------+------------------------+------------------+------------+ diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index 83911aa0f97..773325d982b 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -3,20 +3,80 @@ Running the simulation ====================== +.. _Particle integration and propagation: + +Particle integration and propagation +------------------------------------ + +The main integration scheme of |es| is the velocity Verlet algorithm. +A steepest descent algorithm is used to minimize the system. + +Additional integration schemes are available, which can be coupled to +thermostats to enable Langevin dynamics, Brownian dynamics, Stokesian dynamics, +dissipative particle dynamics, and simulations in the NpT ensemble. + +.. _Rotational degrees of freedom and particle anisotropy: + +Rotational degrees of freedom and particle anisotropy +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +When the feature ``ROTATION`` is compiled in, particles not only have a position, but also an orientation that changes with an angular velocity. A torque on a particle leads to a change in angular velocity depending on the particles rotational inertia. The property :attr:`espressomd.particle_data.ParticleHandle.rinertia` has to be specified as the three eigenvalues of the particles rotational inertia tensor. + +The rotational degrees of freedom are also integrated using a velocity Verlet scheme. +The implementation is based on a quaternion representation of the particle orientation and described in :cite:`omelyan98` with quaternion components indexing made according to the formalism :math:`q = a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k}` :cite:`allen2017`. + +When the Langevin thermostat is enabled, the rotational degrees of freedom are also thermalized. + +Whether or not rotational degrees of freedom are propagated, is controlled on a per-particle and per-axis level, where the axes are the Cartesian axes of the particle in its body-fixed frame. +It is important to note that starting from version 4.0 and unlike in earlier versions of |es|, the particles' rotation is disabled by default. +In this way, just compiling in the ``ROTATION`` feature no longer changes the physics of the system. + +The rotation of a particle is controlled via the :attr:`espressomd.particle_data.ParticleHandle.rotation` property. E.g., the following code adds a particle with rotation enabled on the x axis:: + + import espressomd + system = espressomd.System(box_l=[1, 1, 1]) + system.part.add(pos=(0, 0, 0), rotation=(1, 0, 0)) + +Notes: + +* The orientation of a particle is stored as a quaternion in the :attr:`espressomd.particle_data.ParticleHandle.quat` property. For a value of (1,0,0,0), the body and space frames coincide. +* The space-frame direction of the particle's z-axis in its body frame is accessible through the :attr:`espressomd.particle_data.ParticleHandle.director` property. +* Any other vector can be converted from body to space fixed frame using the :meth:`espressomd.particle_data.ParticleHandle.convert_vector_body_to_space` method. +* When ``DIPOLES`` are compiled in, the particles dipole moment is always co-aligned with the z-axis in the body-fixed frame. +* Changing the particles dipole moment and director will re-orient the particle such that its z-axis in space frame is aligned parallel to the given vector. No guarantees are made for the other two axes after setting the director or the dipole moment. + + +The following particle properties are related to rotation: + +* :attr:`espressomd.particle_data.ParticleHandle.dip` +* :attr:`espressomd.particle_data.ParticleHandle.director` +* :attr:`espressomd.particle_data.ParticleHandle.ext_torque` +* :attr:`espressomd.particle_data.ParticleHandle.gamma_rot` +* :attr:`espressomd.particle_data.ParticleHandle.gamma_rot` +* :attr:`espressomd.particle_data.ParticleHandle.omega_body` +* :attr:`espressomd.particle_data.ParticleHandle.omega_lab` +* :attr:`espressomd.particle_data.ParticleHandle.quat` +* :attr:`espressomd.particle_data.ParticleHandle.rinertia` +* :attr:`espressomd.particle_data.ParticleHandle.rotation` +* :attr:`espressomd.particle_data.ParticleHandle.torque_lab` + + +.. _Integrators: + +Integrators +----------- + To run the integrator call the method -:meth:`espressomd.integrate.Integrator.run`:: +:meth:`system.integrate.run() `:: system.integrator.run(number_of_steps, recalc_forces=False, reuse_forces=False) -where ``number_of_steps`` is the number of time steps the integrator -should perform. The two main integration schemes of |es| are the Velocity Verlet algorithm -and an adaption of the Velocity Verlet algorithm to simulate an NpT ensemble. -A steepest descent implementation is also available. +where ``number_of_steps`` is the number of time steps the integrator should perform. .. _Velocity Verlet Algorithm: Velocity Verlet algorithm -------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^ :meth:`espressomd.integrate.IntegratorHandle.set_vv` @@ -47,7 +107,7 @@ For numerical integration, this equation is discretized to the following steps ( .. math:: v(t+dt) = v(t+dt/2) + \frac{F(x(t+dt),t+dt)}{m} dt/2 -Note that this implementation of the Velocity Verlet algorithm reuses +Note that this implementation of the velocity Verlet algorithm reuses forces in step 1. That is, they are computed once in step 3, but used twice, in step 4 and in step 1 of the next iteration. In the first time step after setting up, there are no forces present yet. Therefore, |es| has @@ -77,10 +137,10 @@ would like to recompute the forces, despite the fact that they are already correctly calculated. To this aim, the option ``recalc_forces`` can be used to enforce force recalculation. -.. _Isotropic NPT integrator: +.. _Isotropic NpT integrator: -Isotropic NPT integrator ------------------------- +Isotropic NpT integrator +^^^^^^^^^^^^^^^^^^^^^^^^ :meth:`espressomd.integrate.IntegratorHandle.set_isotropic_npt` @@ -185,55 +245,10 @@ Notes: * The particle forces :math:`F` include interactions as well as a friction (:math:`\gamma^0`) and noise term (:math:`\sqrt{k_B T \gamma^0 dt} \overline{\eta}`) analogous to the terms in the :ref:`Langevin thermostat`. * The particle forces are only calculated in step 5 and then reused in step 1 of the next iteration. See :ref:`Velocity Verlet Algorithm` for the implications of that. -.. _Rotational degrees of freedom and particle anisotropy: - -Rotational degrees of freedom and particle anisotropy ------------------------------------------------------ - -When the feature ``ROTATION`` is compiled in, particles not only have a position, but also an orientation that changes with an angular velocity. A torque on a particle leads to a change in angular velocity depending on the particles rotational inertia. The property :attr:`espressomd.particle_data.ParticleHandle.rinertia` has to be specified as the three eigenvalues of the particles rotational inertia tensor. - -The rotational degrees of freedom are also integrated using a velocity Verlet scheme. -The implementation is based on a quaternion representation of the particle orientation and described in :cite:`omelyan98` with quaternion components indexing made according to the formalism :math:`q = a + b\mathbf{i} + c\mathbf{j} + d\mathbf{k}` :cite:`allen2017`. - -When the Langevin thermostat is enabled, the rotational degrees of freedom are also thermalized. - -Whether or not rotational degrees of freedom are propagated, is controlled on a per-particle and per-axis level, where the axes are the Cartesian axes of the particle in its body-fixed frame. -It is important to note that starting from version 4.0 and unlike in earlier versions of |es|, the particles' rotation is disabled by default. -In this way, just compiling in the ``ROTATION`` feature no longer changes the physics of the system. - -The rotation of a particle is controlled via the :attr:`espressomd.particle_data.ParticleHandle.rotation` property. E.g., the following code adds a particle with rotation enabled on the x axis:: - - import espressomd - system = espressomd.System(box_l=[1, 1, 1]) - system.part.add(pos=(0, 0, 0), rotation=(1, 0, 0)) - -Notes: - -* The orientation of a particle is stored as a quaternion in the :attr:`espressomd.particle_data.ParticleHandle.quat` property. For a value of (1,0,0,0), the body and space frames coincide. -* The space-frame direction of the particle's z-axis in its body frame is accessible through the :attr:`espressomd.particle_data.ParticleHandle.director` property. -* Any other vector can be converted from body to space fixed frame using the :meth:`espressomd.particle_data.ParticleHandle.convert_vector_body_to_space` method. -* When ``DIPOLES`` are compiled in, the particles dipole moment is always co-aligned with the z-axis in the body-fixed frame. -* Changing the particles dipole moment and director will re-orient the particle such that its z-axis in space frame is aligned parallel to the given vector. No guarantees are made for the other two axes after setting the director or the dipole moment. - - -The following particle properties are related to rotation: - -* :attr:`espressomd.particle_data.ParticleHandle.dip` -* :attr:`espressomd.particle_data.ParticleHandle.director` -* :attr:`espressomd.particle_data.ParticleHandle.ext_torque` -* :attr:`espressomd.particle_data.ParticleHandle.gamma_rot` -* :attr:`espressomd.particle_data.ParticleHandle.gamma_rot` -* :attr:`espressomd.particle_data.ParticleHandle.omega_body` -* :attr:`espressomd.particle_data.ParticleHandle.omega_lab` -* :attr:`espressomd.particle_data.ParticleHandle.quat` -* :attr:`espressomd.particle_data.ParticleHandle.rinertia` -* :attr:`espressomd.particle_data.ParticleHandle.rotation` -* :attr:`espressomd.particle_data.ParticleHandle.torque_lab` - .. _Steepest descent: Steepest descent ----------------- +^^^^^^^^^^^^^^^^ :meth:`espressomd.integrate.IntegratorHandle.set_steepest_descent` @@ -265,7 +280,7 @@ Usage example:: system.integrator.set_vv() # to switch back to velocity Verlet Using a custom convergence criterion -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +"""""""""""""""""""""""""""""""""""" The ``f_max`` parameter can be set to zero to prevent the integrator from halting when a specific force/torque is reached. The integration can then @@ -309,11 +324,18 @@ The correct forces need to be re-calculated after running the integration:: system.integrator.run(0, recalc_forces=True) # re-calculate forces from virtual sites system.integrator.set_vv() +.. _Brownian Dynamics: + +Brownian Dynamics +^^^^^^^^^^^^^^^^^ + +Brownian Dynamics integrator :cite:`schlick2010`. +See details in :ref:`Brownian thermostat`. .. _Stokesian Dynamics: Stokesian Dynamics ------------------- +^^^^^^^^^^^^^^^^^^ .. note:: @@ -365,7 +387,7 @@ thermalize the system, the SD thermostat needs to be activated (see .. _Important_SD: Important -~~~~~~~~~ +""""""""" The particles must be prevented from overlapping. It is mathematically allowed for the particles to overlap to a certain degree. However, once the distance @@ -381,3 +403,314 @@ The near field (so-called lubrication) correction is planned. For now, Stokesian Dynamics provides a good approximation of the hydrodynamics in dilute systems where the average distance between particles is several sphere diameters. + + +.. _Thermostats: + +Thermostats +----------- + +To add a thermostat, call the appropriate setter:: + + system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=41) + +The different thermostats available in |es| will be described in the following +subsections. + +You may combine different thermostats at your own risk by turning them on +one by one. The list of active thermostats can be cleared at any time with +:py:meth:`system.thermostat.turn_off() `. +Not all combinations of thermostats are allowed, though (see +:py:func:`espressomd.thermostat.AssertThermostatType` for details). +Some integrators only work with a specific thermostat and throw an +error otherwise. Note that there is only one temperature for all +thermostats, although for some thermostats like the Langevin thermostat, +particles can be assigned individual temperatures. + +Since |es| does not enforce a particular unit system, it cannot know about +the current value of the Boltzmann constant. Therefore, when specifying +the temperature of a thermostat, you actually do not define the +temperature, but the value of the thermal energy :math:`k_B T` in the +current unit system (see the discussion on units, Section :ref:`On units`). + +All thermostats have a ``seed`` argument that controls the state of the random +number generator (Philox Counter-based RNG). This seed is required on first +activation of a thermostat, unless stated otherwise. It can be omitted in +subsequent calls of the method that activates the same thermostat. The random +sequence also depends on the thermostats counters that are +incremented after each integration step. + +.. _Langevin thermostat: + +Langevin thermostat +^^^^^^^^^^^^^^^^^^^ + +In order to activate the Langevin thermostat the member function +:py:meth:`~espressomd.thermostat.Thermostat.set_langevin` of the thermostat +class :class:`espressomd.thermostat.Thermostat` has to be invoked. +Best explained in an example:: + + import espressomd + system = espressomd.System(box_l=[1, 1, 1]) + system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=41) + +As explained before the temperature is set as thermal energy :math:`k_\mathrm{B} T`. + +The Langevin thermostat is based on an extension of Newton's equation of motion to + +.. 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. +The friction term accounts for dissipation in a surrounding fluid whereas +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). + +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. +The distribution of :math:`\overline{\eta}` is uniform and satisfies + +.. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt + +If the feature ``ROTATION`` is compiled in, the rotational degrees of freedom are +also coupled to the thermostat. If only the first two arguments are +specified then the friction coefficient for the rotation is set to the +same value as that for the translation. +A separate rotational friction coefficient can be set by inputting +``gamma_rotate``. The two options allow one to switch the translational and rotational +thermalization on or off separately, maintaining the frictional behavior. This +can be useful, for instance, in high Péclet number active matter systems, where +one wants to thermalize only the rotational degrees of freedom while +translational degrees of freedom are affected 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 anisotropic diffusion of anisotropic colloids (rods, etc.). + +Using the Langevin thermostat, it is possible to set a temperature and a +friction coefficient for every particle individually via the feature +``LANGEVIN_PER_PARTICLE``. Consult the reference of the ``part`` command +(chapter :ref:`Setting up particles`) for information on how to achieve this. + +.. _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(box_l=[1, 1, 1]) + 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: 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. + +.. _Isotropic NpT thermostat: + +Isotropic NpT thermostat +^^^^^^^^^^^^^^^^^^^^^^^^ + +This feature allows to simulate an (on average) homogeneous and isotropic system in the NpT ensemble. +In order to use this feature, ``NPT`` has to be defined in the :file:`myconfig.hpp`. +Activate the NpT thermostat with the command :py:meth:`~espressomd.thermostat.Thermostat.set_npt` +and setup the integrator for the NpT ensemble with :py:meth:`~espressomd.integrate.IntegratorHandle.set_isotropic_npt`. + +For example:: + + import espressomd + + system = espressomd.System(box_l=[1, 1, 1]) + system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=41) + system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) + +For an explanation of the algorithm involved, see :ref:`Isotropic NpT integrator`. + +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. + +.. _Dissipative Particle Dynamics (DPD): + +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 +encoded in a dissipative interaction between particles :cite:`soddeman03a`. + +To realize a complete DPD fluid model in |es|, three parts are needed: +the DPD thermostat, which controls the temperate, a dissipative interaction +between the particles that make up the fluid, see :ref:`DPD interaction`, +and a repulsive conservative force, see :ref:`Hat interaction`. + +The temperature is set via +:py:meth:`espressomd.thermostat.Thermostat.set_dpd` +which takes ``kT`` and ``seed`` as arguments. + +The friction coefficients and cutoff are controlled via the +:ref:`DPD interaction` on a per type-pair basis. + +The friction (dissipative) and noise (random) term are coupled via the +fluctuation-dissipation theorem. The friction term is a function of the +relative velocity of particle pairs. The DPD thermostat is better for +dynamics than the Langevin thermostat, since it mimics hydrodynamics in +the system. + +As a conservative force any interaction potential can be used, +see :ref:`Isotropic non-bonded interactions`. A common choice is +a force ramp which is implemented as :ref:`Hat interaction`. + +A complete example of setting up a DPD fluid and running it +to sample the equation of state can be found in :file:`/samples/dpd.py`. + +When using a Lennard-Jones interaction, :math:`{r_\mathrm{cut}} = +2^{\frac{1}{6}} \sigma` is a good value to choose, so that the +thermostat acts on the relative velocities between nearest neighbor +particles. Larger cutoffs including next nearest neighbors or even more +are unphysical. + +Boundary conditions for DPD can be introduced by adding the boundary +as a particle constraint, and setting a velocity and a type on it, see +:class:`espressomd.constraints.Constraint`. Then a +:ref:`DPD interaction` with the type can be defined, which acts as a +boundary condition. + +.. _LB thermostat: + +Lattice-Boltzmann thermostat +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The :ref:`Lattice-Boltzmann` thermostat acts similar to the :ref:`Langevin thermostat` in that the governing equation for particles is + +.. 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. + +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. +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. + +The magnitude of the frictional coupling can be adjusted by the +parameter ``gamma``. To enable the LB thermostat, use:: + + import espressomd + import espressomd.lb + system = espressomd.System(box_l=[1, 1, 1]) + lbf = espressomd.lb.LBFluid(agrid=1, dens=1, visc=1, tau=0.01) + system.actors.add(lbf) + system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5) + +No other thermostatting mechanism is necessary +then. Please switch off any other thermostat before starting the LB +thermostatting mechanism. + +The LBM implementation provides a fully thermalized LB fluid, all +nonconserved modes, including the pressure tensor, fluctuate correctly +according to the given temperature and the relaxation parameters. All +fluctuations can be switched off by setting the temperature to 0. + +.. note:: Coupling between LB and MD only happens if the LB thermostat is set with a :math:`\gamma \ge 0.0`. + +.. _Stokesian thermostat: + +Stokesian thermostat +^^^^^^^^^^^^^^^^^^^^ + +.. note:: + + Requires ``STOKESIAN_DYNAMICS`` external feature, enabled with + ``-DWITH_STOKESIAN_DYNAMICS=ON``. + +In order to thermalize a Stokesian Dynamics simulation, the SD thermostat +needs to be activated via:: + + import espressomd + system = espressomd.System(box_l=[1.0, 1.0, 1.0]) + system.periodicity = [False, False, False] + system.time_step = 0.01 + system.cell_system.skin = 0.4 + system.part.add(pos=[0, 0, 0], rotation=[1, 0, 0], ext_force=[0, 0, -1]) + system.thermostat.set_stokesian(kT=1.0, seed=43) + system.integrator.set_stokesian_dynamics(viscosity=1.0, radii={0: 1.0}) + system.integrator.run(100) + +where ``kT`` denotes the desired temperature of the system, and ``seed`` the +seed for the random number generator. diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index b6bf671bad6..59b070115c5 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -145,7 +145,7 @@ range should be much smaller than the total system size, leaving out all interactions between non-adjacent cells can mean a tremendous speed-up. Moreover, since for constant interaction range, the number of particles in a cell depends only on the density. The number of interactions is -therefore of the order N instead of order :math:`N^2` if one has to +therefore of the order :math:`N` instead of order :math:`N^2` if one has to calculate all pair interactions. .. _N-squared: @@ -154,7 +154,7 @@ N-squared ~~~~~~~~~ Invoking :py:meth:`~espressomd.cellsystem.CellSystem.set_n_square` -selects the very primitive nsquared cellsystem, which calculates +selects the very primitive N-squared cellsystem, which calculates the interactions for all particle pairs. Therefore it loops over all particles, giving an unfavorable computation time scaling of :math:`N^2`. However, algorithms like MMM1D or the plain Coulomb @@ -165,7 +165,7 @@ interactions. :: system.cell_system.set_n_square() -In a multiple processor environment, the nsquared cellsystem uses a +In a multiple processor environment, the N-squared cellsystem uses a simple particle balancing scheme to have a nearly equal number of particles per CPU, :math:`n` nodes have :math:`m` particles, and :math:`p-n` nodes have :math:`m+1` particles, such that @@ -184,318 +184,9 @@ this node is twice as high. For 3 processors, the interactions are 0-0, 1-1, 2-2, 0-1, 1-2, 0-2. Of these interactions, node 0 treats 0-0 and 0-2, node 1 treats 1-1 and 0-1, and node 2 treats 2-2 and 1-2. -Therefore it is highly recommended that you use nsquared only with an +Therefore it is highly recommended that you use N-squared only with an odd number of nodes, if with multiple processors at all. -.. _Thermostats: - -Thermostats ------------ - -The thermostat can be controlled by the class :class:`espressomd.thermostat.Thermostat`. -The different thermostats available in |es| will be described in the following -subsections. - -You may combine different thermostats at your own risk by turning them on -one by one. The list of active thermostats can be cleared at any time with -:py:meth:`system.thermostat.turn_off() `. -Not all combinations of thermostats are allowed, though (see -:py:func:`espressomd.thermostat.AssertThermostatType` for details). -Some integrators only work with a specific thermostat and throw an -error otherwise. Note that there is only one temperature for all -thermostats, although for some thermostats like the Langevin thermostat, -particles can be assigned individual temperatures. - -Since |es| does not enforce a particular unit system, it cannot know about -the current value of the Boltzmann constant. Therefore, when specifying -the temperature of a thermostat, you actually do not define the -temperature, but the value of the thermal energy :math:`k_B T` in the -current unit system (see the discussion on units, Section :ref:`On units`). - -All thermostats have a ``seed`` argument that controls the state of the random -number generator (Philox Counter-based RNG). This seed is required on first -activation of a thermostat, unless stated otherwise. It can be omitted in -subsequent calls of the method that activates the same thermostat. The random -sequence also depends on the thermostats counters that are -incremented after each integration step. - -.. _Langevin thermostat: - -Langevin thermostat -~~~~~~~~~~~~~~~~~~~ - -In order to activate the Langevin thermostat the member function -:py:meth:`~espressomd.thermostat.Thermostat.set_langevin` of the thermostat -class :class:`espressomd.thermostat.Thermostat` has to be invoked. -Best explained in an example:: - - import espressomd - system = espressomd.System(box_l=[1, 1, 1]) - system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=41) - -As explained before the temperature is set as thermal energy :math:`k_\mathrm{B} T`. - -The Langevin thermostat is based on an extension of Newton's equation of motion to - -.. 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. -The friction term accounts for dissipation in a surrounding fluid whereas -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). - -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. -The distribution of :math:`\overline{\eta}` is uniform and satisfies - -.. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt - -If the feature ``ROTATION`` is compiled in, the rotational degrees of freedom are -also coupled to the thermostat. If only the first two arguments are -specified then the friction coefficient for the rotation is set to the -same value as that for the translation. -A separate rotational friction coefficient can be set by inputting -``gamma_rotate``. The two options allow one to switch the translational and rotational -thermalization on or off separately, maintaining the frictional behavior. This -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 affected 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 anisotropic diffusion of anisotropic colloids (rods, etc.). - -Using the Langevin thermostat, it is possible to set a temperature and a -friction coefficient for every particle individually via the feature -``LANGEVIN_PER_PARTICLE``. Consult the reference of the ``part`` command -(chapter :ref:`Setting up particles`) for information on how to achieve this. - -.. _LB thermostat: - -Lattice-Boltzmann thermostat -~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -The :ref:`Lattice-Boltzmann` thermostat acts similar to the :ref:`Langevin thermostat` in that the governing equation for particles is - -.. 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. - -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. -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. - -The magnitude of the frictional coupling can be adjusted by the -parameter ``gamma``. To enable the LB thermostat, use:: - - import espressomd - import espressomd.lb - system = espressomd.System(box_l=[1, 1, 1]) - lbf = espressomd.lb.LBFluid(agrid=1, dens=1, visc=1, tau=0.01) - system.actors.add(lbf) - system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5) - -No other thermostatting mechanism is necessary -then. Please switch off any other thermostat before starting the LB -thermostatting mechanism. - -The LBM implementation provides a fully thermalized LB fluid, all -nonconserved modes, including the pressure tensor, fluctuate correctly -according to the given temperature and the relaxation parameters. All -fluctuations can be switched off by setting the temperature to 0. - -.. note:: Coupling between LB and MD only happens if the LB thermostat is set with a :math:`\gamma \ge 0.0`. - - -.. _Dissipative Particle Dynamics (DPD): - -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 -encoded in a dissipative interaction between particles :cite:`soddeman03a`. - -To realize a complete DPD fluid model in |es|, three parts are needed: -the DPD thermostat, which controls the temperate, a dissipative interaction -between the particles that make up the fluid, see :ref:`DPD interaction`, -and a repulsive conservative force, see :ref:`Hat interaction`. - -The temperature is set via -:py:meth:`espressomd.thermostat.Thermostat.set_dpd` -which takes ``kT`` and ``seed`` as arguments. - -The friction coefficients and cutoff are controlled via the -:ref:`DPD interaction` on a per type-pair basis. For details -see there. - -The friction (dissipative) and noise (random) term are coupled via the -fluctuation-dissipation theorem. The friction term is a function of the -relative velocity of particle pairs. The DPD thermostat is better for -dynamics than the Langevin thermostat, since it mimics hydrodynamics in -the system. - -As a conservative force any interaction potential can be used, -see :ref:`Isotropic non-bonded interactions`. A common choice is -a force ramp which is implemented as :ref:`Hat interaction`. - -A complete example of setting up a DPD fluid and running it -to sample the equation of state can be found in :file:`/samples/dpd.py`. - -When using a Lennard-Jones interaction, :math:`{r_\mathrm{cut}} = -2^{\frac{1}{6}} \sigma` is a good value to choose, so that the -thermostat acts on the relative velocities between nearest neighbor -particles. Larger cutoffs including next nearest neighbors or even more -are unphysical. - -Boundary conditions for DPD can be introduced by adding the boundary -as a particle constraint, and setting a velocity and a type on it, see -:class:`espressomd.constraints.Constraint`. Then a -:ref:`DPD interaction` with the type can be defined, which acts as a -boundary condition. - -.. _Isotropic NPT thermostat: - -Isotropic NPT thermostat -~~~~~~~~~~~~~~~~~~~~~~~~ - -This feature allows to simulate an (on average) homogeneous and isotropic system in the NPT ensemble. -In order to use this feature, ``NPT`` has to be defined in the :file:`myconfig.hpp`. -Activate the NPT thermostat with the command :py:meth:`~espressomd.thermostat.Thermostat.set_npt` -and setup the integrator for the NPT ensemble with :py:meth:`~espressomd.integrate.IntegratorHandle.set_isotropic_npt`. - -For example:: - - import espressomd - - system = espressomd.System(box_l=[1, 1, 1]) - system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=41) - system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) - -For an explanation of the algorithm involved, see :ref:`Isotropic NPT integrator`. - -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(box_l=[1, 1, 1]) - 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. - -.. _Stokesian thermostat: - -Stokesian thermostat -~~~~~~~~~~~~~~~~~~~~ - -.. note:: - - Requires ``STOKESIAN_DYNAMICS`` external feature, enabled with - ``-DWITH_STOKESIAN_DYNAMICS=ON``. - -In order to thermalize a Stokesian Dynamics simulation, the SD thermostat -needs to be activated via:: - - import espressomd - system = espressomd.System(box_l=[1.0, 1.0, 1.0]) - system.periodicity = [False, False, False] - system.time_step = 0.01 - system.cell_system.skin = 0.4 - system.part.add(pos=[0, 0, 0], rotation=[1, 0, 0], ext_force=[0, 0, -1]) - system.thermostat.set_stokesian(kT=1.0, seed=43) - system.integrator.set_stokesian_dynamics(viscosity=1.0, radii={0: 1.0}) - system.integrator.run(100) - -where ``kT`` denotes the desired temperature of the system, and ``seed`` the -seed for the random number generator. - .. _CUDA: diff --git a/doc/tutorials/charged_system/charged_system-1.ipynb b/doc/tutorials/charged_system/charged_system-1.ipynb index 485150305b3..d743b8731db 100644 --- a/doc/tutorials/charged_system/charged_system-1.ipynb +++ b/doc/tutorials/charged_system/charged_system-1.ipynb @@ -506,7 +506,7 @@ " * resets the system clock\n", "\n", "**Hints:**\n", - "* The relevant parts of the documentation can be found here: [Thermostats](http://espressomd.org/html/doc/system_setup.html#thermostats), [Particle List](http://espressomd.org/html/doc/espressomd.html#espressomd.particle_data.ParticleList),\n", + "* The relevant parts of the documentation can be found here: [Thermostats](http://espressomd.org/html/doc/running.html#thermostats), [ParticleList](http://espressomd.org/html/doc/espressomd.html#espressomd.particle_data.ParticleList),\n", "[AutoUpdateAccumulators](http://espressomd.org/html/doc/espressomd.html#espressomd.accumulators.AutoUpdateAccumulators),\n", "[System properties](http://espressomd.org/html/doc/espressomd.html#module-espressomd.system)" ]