diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index 9f7a3e9e9e5..8637d82f132 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/doc/sphinx/inter_non-bonded.rst @@ -559,60 +559,80 @@ DPD interaction Feature ``DPD`` required. This is a special interaction that is to be used in conjunction with the -`Dissipative Particle Dynamics (DPD)` thermostat, for a general description -of the algorithm see there. The parameters can be set via:: +:ref:`Dissipative Particle Dynamics (DPD)` thermostat. +The parameters can be set via:: system.non_bonded_inter[type1, type2].dpd.set_params(**kwargs) -This command defines a velocity-dependent interaction -between particles of the types ``type1`` and ``type2``. For a description of the input arguments -see :class:`espressomd.interactions.DPDInteraction`. The interaction -only has an effect if the DPD thermostat activated, and acts according to the -temperature of the thermostat. +This command defines an interaction between particles of the types ``type1`` and ``type2`` +that contains velocity-dependent friction and noise according +to a temperature set by :py:meth:`espressomd.thermostat.Thermostat.set_dpd()`. +The parameters for the interaction are -The interaction consists of a dissipative force :math:`\vec{F}_{ij}^{D}` and -a random force :math:`\vec{F}_{ij}^R`, and is decomposed into a component -parallel and perpendicular to the distance vector of the particle pair :math:`\vec{F}_{ij}`. -The parameters for the two parts can be chosen independently. -The force contributions of the parallel part are +* ``gamma`` +* ``weight_function`` +* ``r_cut`` +* ``trans_gamma`` +* ``trans_weight_function`` +* ``trans_r_cut`` -.. math:: \vec{F}_{ij}^{D} = -\zeta w^D (r_{ij}) (\hat{r}_{ij} \cdot \vec{v}_{ij}) \hat{r}_{ij} +which will be explained below. The interaction +only has an effect if the DPD thermostat is activated. -for the dissipative force and +The interaction consists of a friction force :math:`\vec{F}_{ij}^{D}` and +a random force :math:`\vec{F}_{ij}^R` added to the other interactions +between particles :math:`i` and :math:`j`. It is decomposed into a component +parallel and perpendicular to the distance vector :math:`\vec{r}_{ij}` of the particle pair . +The friction force contributions of the parallel part are -.. math:: \vec{F}_{ij}^R = \sigma w^R (r_{ij}) \Theta_{ij} \hat{r}_{ij} +.. math:: \vec{F}_{ij}^{D} = -\gamma_\parallel w_\parallel (r_{ij}) (\hat{r}_{ij} \cdot \vec{v}_{ij}) \hat{r}_{ij} -for the random force. Here :math:`w^D` and :math:`w^R` are weight functions that -can be specified via the weight_function parameter of the interaction. The dissipative -and random weight function are related by the dissipation-fluctuation theorem: +for the dissipative force and -.. math:: (\sigma w^R (r_{ij}))^2=\zeta w^D (r_{ij}) \text{k}_\text{B} T +.. math:: \vec{F}_{ij}^R = \sqrt{2 k_B T \gamma_\parallel w_\parallel (r_{ij}) } \eta_{ij}(t) \hat{r}_{ij} -The possible values for weight_function are 0 and 1, corresponding to the -order of :math:`w^D`: +for the random force. This introduces the friction coefficient :math:`\gamma_\parallel` (parameter ``gamma``) and the weight function +:math:`w_\parallel`. The thermal energy :math:`k_B T` is not set by the interaction, +but by the DPD thermostat (:py:meth:`espressomd.thermostat.Thermostat.set_dpd()`) +to be equal for all particles. The weight function can be specified via the ``weight_function`` switch. +The possible values for ``weight_function`` are 0 and 1, corresponding to the +order of :math:`w_\parallel`: .. math:: - w^D (r_{ij}) = ( w^R (r_{ij})) ^2 = - \left\{ + w_\parallel (r_{ij}) = \left\{ \begin{array}{clcr} 1 & , \; \text{weight_function} = 0 \\ - {( 1 - \frac{r_{ij}}{r_c}} )^2 & , \; \text{weight_function} = 1 + {( 1 - \frac{r_{ij}}{r^\text{cut}_\parallel}} )^2 & , \; \text{weight_function} = 1 \end{array} \right. + +Both weight functions are set to zero for :math:`r_{ij}>r^\text{cut}_\parallel` (parameter ``r_cut``). -For the perpendicular part, the dissipative force is calculated by +The random force has the properties -.. math:: \vec{F}_{ij}^{D} = -\zeta w^D (r_{ij}) (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{v}_{ij} +.. math:: <\eta_{ij}(t)> = 0 , <\eta_{ij}^\alpha(t)\eta_{kl}^\beta(t')> = \delta_{\alpha\beta} \delta_{ik}\delta_{jl}\delta(t-t') -The random force by +and is numerically discretized to a random number :math:`\overline{\eta}` for each spatial +component for each particle pair drawn from a uniform distribution +with properties -.. math:: \vec{F}_{ij}^R = \sigma w^R (r_{ij}) (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{\Theta}_{ij} +.. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt + +For the perpendicular part, the dissipative and random force are calculated analogously -The parameters define the strength of the friction and the cutoff in the -same way as above. Note: This interaction does *not* conserve angular +.. math:: \vec{F}_{ij}^{D} = -\gamma_\bot w^D (r_{ij}) (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{v}_{ij} +.. math:: \vec{F}_{ij}^R = \sqrt{2 k_B T \gamma_\bot w_\bot (r_{ij})} (I-\hat{r}_{ij}\otimes\hat{r}_{ij}) \cdot \vec{\eta}_{ij} + +and introduce the second set of parameters prefixed with ``trans_``. +For :math:`w_\bot (r_{ij})` (parameter ``trans_weight_function``) +the same options are available as for :math:`w_\parallel (r_{ij})`. + +Note: This interaction does *not* conserve angular momentum. +A more detailed description of the interaction can be found in :cite:`soddeman03a`. + .. _Thole correction: Thole correction diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst index 7a4476aa738..0de5f1aecdb 100644 --- a/doc/sphinx/introduction.rst +++ b/doc/sphinx/introduction.rst @@ -548,7 +548,7 @@ report so to the developers. +--------------------------------+------------------------+------------------+------------+ | Langevin Thermostat | Core | Core | Yes | +--------------------------------+------------------------+------------------+------------+ -| Isotropic NPT | None | Single | Yes | +| Isotropic NPT | Experimental | None | Yes | +--------------------------------+------------------------+------------------+------------+ | Quaternion Integrator | Core | Good | Yes | +--------------------------------+------------------------+------------------+------------+ diff --git a/doc/sphinx/lb.rst b/doc/sphinx/lb.rst index 9e858a8a977..0a56763c113 100644 --- a/doc/sphinx/lb.rst +++ b/doc/sphinx/lb.rst @@ -121,53 +121,60 @@ load in the particles with the correct forces, and use:: sys.integrator.run(steps=number_of_steps, reuse_forces=True) -upon the first call to :ref:`run `. This causes the +upon the first call ``integrator.run``. This causes the old forces to be reused and thus conserves momentum. -.. _LB as a thermostat: +.. _Interpolating velocities: -LB as a thermostat ------------------- +Interpolating velocities +------------------------ -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. Furthermore a seed has to be given for the -thermalization of the particle coupling. The magnitude of the fricitional coupling can be adjusted by -the parameter ``gamma``. -To enable the LB thermostat, use:: +To get interpolated velocity values between lattice nodes, the function:: - sys.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5) + lb.get_interpolated_velocity(pos = [1.1,1.2,1.3]) + +with a single position ``pos`` as an argument can be used. +For the GPU fluid :class:`espressomd.lb.LBFluidGPU` +also :py:meth:`espressomd.lb.LBFluidGPU.get_interpolated_fluid_velocity_at_positions()` +is available, which expects a numpy array of positions as an argument. + +By default, the interpolation is done linearly between the nearest 8 LB nodes, +but for the GPU implementation also a quadratic scheme involving 27 nodes is implemented +(see eqs. 297 and 301 in :cite:`duenweg08a`). +You can choose by calling +one of:: + + lb.set_interpolation_order('linear') + lb.set_interpolation_order('quadratic') + +A note on boundaries: -The LBM implementation in |es| uses Ahlrichs and Dünweg's point coupling -method to couple MD particles the LB fluid. This coupling consists of a -frictional and a random force, similar to the :ref:`Langevin thermostat`: +Both interpolation schemes don't take into account the physical location of the boundaries +(e.g. in the middle between two nodes for a planar wall) but will use the boundary node slip velocity +at the node position. This means that every interpolation involving at least one +boundary node will introduce an error. -.. math:: \vec{F} = -\gamma \left(\vec{v}-\vec{u}\right) + \vec{F}_R. +.. _Coupling LB to a MD simulation: -The momentum acquired by the particles is then transferred back to the -fluid using a linear interpolation scheme, to preserve total momentum. -In the GPU implementation the force can alternatively be interpolated -using a three point scheme which couples the particles to the nearest 27 -LB nodes. This can be called using "lbfluid 3pt" and is described in -Dünweg and Ladd by equation 301 :cite:`duenweg08a`. +Coupling LB to a MD simulation +------------------------------ -The frictional force tends to decrease the relative -velocity between the fluid and the particle whereas the random forces -are chosen so large that the average kinetic energy per particle -corresponds to the given temperature, according to a fluctuation -dissipation theorem. No other thermostatting mechanism is necessary -then. Please switch off any other thermostat before starting the LB -thermostatting mechanism. +MD particles can be coupled to a LB fluid through frictional coupling. The friction force -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. + .. math:: F_{i,\text{frict}} = - \gamma (v_i(t)-u(x_i(t),t)) + +depends on the particle velocity :math:`v` and the fluid velocity :math:`u`. It acts both +on the particle and the fluid (in opposite direction). Because the fluid is also affected, +multiple particles can interact via hydrodynamic interactions. As friction in molecular systems is +accompanied by fluctuations, the particle-fluid coupling has to be activated through +the :ref:`LB thermostat` (See more detailed description there). A short example is:: + + sys.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5) -.. note:: Coupling between LB and MD only happens if the LB thermostat is set with a :math:`\gamma \ge 0.0`. +where ``lbf`` is an instance of either :class:`espressomd.lb.LBFluid` or :class:`espressomd.lb.LBFluidGPU`, +``gamma`` the friction coefficient and ``seed`` the seed for the random number generator involved +in the thermalization. -Regarding the unit of the temperature, please refer to -Section :ref:`On units`. .. _Reading and setting properties of single lattice nodes: diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index 89232041ac1..f83bbb57ace 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -3,25 +3,53 @@ Running the simulation ====================== -.. _Integrator: - -Integrator ----------- - To run the integrator call the method :meth:`espressomd.integrate.Integrator.run`:: - system.integrator.run(steps=number_of_steps, recalc_forces=False, reuse_forces=False) + 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. +should perform. The two main integration schemes of |es| are the Velocity Verlet algorithm +and an adaption of the Velocity Verlet algorithms to simulate an NPT ensemble. +A steepest descent implementation is also available. + +.. _Velocity Verlet Algorithm: + +Velocity Verlet algorithm +------------------------- + +:func:`espressomd.integrate.Integrator.set_vv` + +The equations of motion for the trajectory of pointlike particles read + +.. math:: \dot v_i(t) = F_i(\{x_j\},v_i,t)/m_i \\ \dot x_i(t) = v_i(t), + +where :math:`x_i`, :math:`v_i`, :math:`m_i` are position, velocity and mass of +particle :math:`i` and :math:`F_i(\{x_j\},v_i,t)` the forces acting on it. +These forces comprise all interactions with other particles and external fields +as well as non-deterministic contributions described in :ref:`Thermostats`. + +For numerical integration, this equation is discretized to the following steps (:cite:`rapaport04` eqs. 3.5.8 - 3.5.10): + +1. Calculate the velocity at the half step + +.. math:: v(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2 -|es| uses the Velocity Verlet algorithm for the integration of the equations -of motion. +2. Calculate the new position + +.. math:: x(t+dt) = x(t) + v(t+dt/2) dt + +3. Calculate the force based on the new position + +.. math:: F = F(x(t+dt), v(t+dt/2), t+dt) + +4. Calculate the new velocity + +.. 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 -forces. That is, they are computed once in the middle of the time step, -but used twice, at the beginning and end. In the first time +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 to compute them before the first time step. That has two consequences: first, random forces are redrawn, resulting in a narrower distribution @@ -49,51 +77,115 @@ 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. -.. _Run steepest descent minimization: +.. _Isotropic NPT thermostat integrator: -Run steepest descent minimization ---------------------------------- +Isotropic NPT thermostat +------------------------ -:func:`espressomd.integrate.Integrator.set_steepest_descent` +:py:func:`~espressomd.integrate.Integrator.set_isotropic_npt` +As the NPT thermostat alters the way the equations of motion are integrated, it is +discussed here and only a brief summary is given in :ref:`Thermostats`. +To activate the NPT integrator, use :py:func:`~espressomd.integrate.Integrator.set_isotropic_npt` +with parameters: -This feature is used to propagate each particle by a small distance parallel to the force acting on it. -When only conservative forces for which a potential exists are in use, this is equivalent to a steepest descent energy minimization. -A common application is removing overlap between randomly placed particles. + * ``ext_pressure``: (float) The external pressure as float variable. + * ``piston``: (float) The mass of the applied piston as float variable. + * ``direction``: [int, int, int] Flags to enable/disable box dimensions to be subject to fluctuations. By default, all directions are enabled. -Please note that the behavior is undefined if a thermostat is activated. -It runs a simple steepest descent algorithm: +Additionally, a NPT thermostat has to be set by :py:func:`~espressomd.thermostat.Thermostat.set_npt()` +with parameters: -Iterate + * ``kT``: (float) Thermal energy of the heat bath + * ``gamma0``: (float) Friction coefficient of the bath + * ``gammav``: (float) Artificial friction coefficient for the volume fluctuations. -.. math:: p_i = p_i + \min(\texttt{gamma} \times F_i, \texttt{max_displacement}), +A code snippet would look like:: -while the maximal force is bigger than ``f_max`` or for at most ``max_steps`` times. The energy -is relaxed by ``gamma``, while the change per coordinate per step is limited to ``max_displacement``. -The combination of ``gamma`` and ``max_displacement`` can be used to get a poor man's adaptive update. -Rotational degrees of freedom are treated similarly: each particle is -rotated around an axis parallel to the torque acting on the particle. -Please be aware of the fact that this needs not to converge to a local -minimum in periodic boundary conditions. Translational and rotational -coordinates that are fixed using the ``fix`` and ``rotation`` attribute of particles are not altered. + import espressomd -Usage example:: + system = espressomd.System() + system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0) + system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) - system.integrator.set_steepest_descent( - f_max=0, gamma=0.1, max_displacement=0.1) - system.integrator.run(20) - system.integrator.set_vv() # to switch back to velocity verlet +The physical meaning of these parameters is described below: + +The relaxation towards a desired pressure :math:`P` (parameter ``ext_pressure``) +is enabled by treating the box +volume :math:`V` as a degree of freedom with corresponding momentum :math:`\Pi = Q\dot{V}`, +where :math:`Q` (parameter ``piston``) is an artificial piston mass. +Which box dimensions are affected to change the volume can be controlled by a list of +boolean flags for parameter ``direction``. +An additional energy :math:`H_V = 1/(2Q)\Pi + PV` +associated with the volume is postulated. This results in a "force" on the box such that + +.. math:: \dot{\Pi} = \mathcal{P} - P + +where + +.. math:: \mathcal{P} = \frac{1}{Vd} \sum_{i,j} f_{ij}x_{ij} + \frac{1}{Vd} \sum_i m_i v_i^2 + +Here :math:`\mathcal{P}` is the instantaneous pressure, :math:`d` the dimension of the system (number of flags set by ``direction``), :math:`f_{ij}` the short range interaction force between particles :math:`i` and :math:`j` and :math:`x_{ij}= x_j - x_i`. + +In addition to this deterministic force, a friction :math:`-\frac{\gamma^V}{Q}\Pi(t)` and noise :math:`\sqrt{k_B T \gamma^V} \eta(t)` are added for the box +volume dynamics and the particle dynamics. . This introduces three new parameters: +The friction coefficient for the box :math:`\gamma^V` (parameter ``gammav``), the friction coefficient of the particles :math:`\gamma^0` (parameter ``gamma0``) and the thermal energy :math:`k_BT` (parameter ``kT``). +For a discussion of these terms and their discretisation, see :ref:`Langevin thermostat`, which uses the same approach, but only for particles. +As a result of box geometry changes, the particle positions and velocities have to be rescaled +during integration. +The discretisation consists of the following steps (see :cite:`kolb99a` for a full derivation of the algorithm): +1. Calculate the particle velocities at the half step +.. math:: v'(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2 -.. _Integrating rotational degrees of freedom: +2. Calculate the instantaneous pressure and "volume momentum" -Integrating rotational degrees of freedom ------------------------------------------ -When the feature ``ROTATION`` is compiled in, Particles not only have a position, but also an orientation. -Just as a force on a particle leads to an increase in linear velocity, a torque on a particle leads to an increase in angular velocity. The rotational degrees of freedom are also integrated using a velocity Verlet scheme. +.. math:: \mathcal{P} = \mathcal{P}(x(t),V(t),f(x(t)), v'(t+dt/2)) +.. math:: \Pi(t+dt/2) = \Pi(t) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t) dt/2 + \sqrt{k_B T \gamma^V dt} \overline{\eta} + +3. Calculate box volume and scaling parameter :math:`L` at half step and full step, scale the simulation box accordingly + +.. math:: V(t+dt/2) = V(t) + \frac{\Pi(t+dt/2)}{Q} dt/2 +.. math:: L(t+dt/2) = V(t+dt/2)^{1/d} +.. math:: V(t+dt) = V(t+dt/2) + \frac{\Pi(t+dt/2)}{Q} dt/2 +.. math:: L(t+dt) = V(t+dt)^{1/d} + +4. Update particle positions and scale velocities + +.. math:: x(t+dt) = \frac{L(t+dt)}{L(t)} \left[ x(t) + \frac{L^2(t)}{L^2(t+dt/2)} v(t+dt/2) dt \right] +.. math:: v(t+dt/2) = \frac{L(t)}{L(t+dt)} v'(t+dt/2) + +5. Calculate forces, instantaneous pressure and "volume momentum" + +.. math:: F = F(x(t+dt),v(t+dt/2),t) +.. math:: \mathcal{P} = \mathcal{P}(x(t+dt),V(t+dt),f(x(t+dt)), v(t+dt/2)) +.. math:: \Pi(t+dt) = \Pi(t+dt/2) + (\mathcal{P}-P) dt/2 -\frac{\gamma^V}{Q}\Pi(t+dt/2) dt/2 + \sqrt{k_B T \gamma^V dt} \overline{\eta} + +6. Update the velocities + +.. math:: v(t+dt) = v(t+dt/2) + \frac{F(t+dt)}{m} dt/2 + +Notes: + +* The NPT algorithm is only tested for all 3 directions enabled for scaling. Usage of ``direction`` is considered an experimental feature. +* In step 4, only those coordinates are scaled for which ``direction`` is set. +* The for the instantaneous pressure the same limitations of applicability hold as described in :ref:`Pressure`. +* The particle forces :math:`F` include interactions as well as a friction and noise term analogous to the terms in the :ref:`Langevin thermostat`. +* The particle forces 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`. + 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. @@ -129,3 +221,36 @@ The following particle properties are related to rotation: * :attr:`espressomd.particle_data.ParticleHandle.rotation` * :attr:`espressomd.particle_data.ParticleHandle.torque_lab` +.. _Steepest descent: + +Steepest descent +---------------- + +:func:`espressomd.integrate.Integrator.set_steepest_descent` + +This feature is used to propagate each particle by a small distance parallel to the force acting on it. +When only conservative forces for which a potential exists are in use, this is equivalent to a steepest descent energy minimization. +A common application is removing overlap between randomly placed particles. + +Please note that the behavior is undefined if a thermostat is activated. +It runs a simple steepest descent algorithm: + +Iterate + +.. math:: p_i = p_i + \min(\texttt{gamma} \times F_i, \texttt{max_displacement}), + +while the maximal force is bigger than ``f_max`` or for at most ``max_steps`` times. The energy +is relaxed by ``gamma``, while the change per coordinate per step is limited to ``max_displacement``. +The combination of ``gamma`` and ``max_displacement`` can be used to get a poor man's adaptive update. +Rotational degrees of freedom are treated similarly: each particle is +rotated around an axis parallel to the torque acting on the particle. +Please be aware of the fact that this needs not to converge to a local +minimum in periodic boundary conditions. Translational and rotational +coordinates that are fixed using the ``fix`` and ``rotation`` attribute of particles are not altered. + +Usage example:: + + system.integrator.set_steepest_descent( + f_max=0, gamma=0.1, max_displacement=0.1) + system.integrator.run(20) + system.integrator.set_vv() # to switch back to velocity verlet diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 52c4fbcedf5..9419d8e4a97 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -140,7 +140,7 @@ Details about the cell system can be obtained by :meth:`espressomd.System().cell Domain decomposition ~~~~~~~~~~~~~~~~~~~~ -Invoking :py:attr:`~espressomd.cellsystem.CellSystem.set_domain_decomposition` +Invoking :py:meth:`~espressomd.cellsystem.CellSystem.set_domain_decomposition` selects the domain decomposition cell scheme, using Verlet lists for the calculation of the interactions. If you specify ``use_verlet_lists=False``, only the domain decomposition is used, but not the Verlet lists. :: @@ -166,7 +166,7 @@ calculate all pair interactions. N-squared ~~~~~~~~~ -Invoking :py:attr:`~espressomd.cellsystem.CellSystem.set_n_square` +Invoking :py:meth:`~espressomd.cellsystem.CellSystem.set_n_square` selects the very primitive nsquared cellsystem, which calculates the interactions for all particle pairs. Therefore it loops over all particles, giving an unfavorable computation time scaling of @@ -205,7 +205,7 @@ odd number of nodes, if with multiple processors at all. Layered cell system ~~~~~~~~~~~~~~~~~~~ -Invoking :py:attr:`~espressomd.cellsystem.CellSystem.set_layered` +Invoking :py:meth:`~espressomd.cellsystem.CellSystem.set_layered` selects the layered cell system, which is specifically designed for the needs of the MMM2D algorithm. Basically it consists of a nsquared algorithm in x and y, but a domain decomposition along z, i.e. the @@ -251,7 +251,7 @@ Langevin thermostat ~~~~~~~~~~~~~~~~~~~ In order to activate the Langevin thermostat the member function -:py:attr:`~espressomd.thermostat.Thermostat.set_langevin` of the thermostat +:py:meth:`~espressomd.thermostat.Thermostat.set_langevin` of the thermostat class :class:`espressomd.thermostat.Thermostat` has to be invoked. Best explained in an example:: @@ -262,18 +262,31 @@ Best explained in an example:: therm.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 consists of a friction and noise term coupled -via the fluctuation-dissipation theorem. The friction term is a function -of the particle velocities. By specifying the diffusion coefficient for -the particle becomes -.. math:: D = \frac{\text{temperature}}{\text{gamma}}. +The Langevin thermostat is based on an extension of Newton's equation of motion to -The relaxation time is given by :math:`\text{gamma}/\text{MASS}`, with -``MASS`` the particle's mass. For a more detailed explanation, refer to -:cite:`grest86a`. An anisotropic diffusion coefficient tensor is available to -simulate anisotropic colloids (rods, etc.) properly. It can be enabled by the -feature ``PARTICLE_ANISOTROPY``. +.. 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 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 @@ -283,43 +296,94 @@ using a fixed seed) or not (by using a randomized seed). 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 diffusion coefficient for the rotation is set to the +specified then the friction coefficient for the rotation is set to the same value as that for the translation. - -A separate rotational diffusion coefficient can be set by inputting -``gamma_rotate``. This also allows one to properly match the translational and -rotational diffusion coefficients of a sphere. Feature ``ROTATIONAL_INERTIA`` -enables an anisotropic rotational diffusion coefficient tensor through -corresponding friction coefficients. - -Finally, the two options allow one to switch the translational and rotational +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 the rotational degrees of freedom and +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 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. +Detailes 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. + +Furthermore a ``seed`` has to be given for the +thermalization of the particle coupling. The magnitude of the fricitional coupling can be adjusted by +the parameter ``gamma``. +To enable the LB thermostat, use:: + + sys.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) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -To realize a complete DPD fluid model, three parts are needed: +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. -The DPD thermostat can be invoked by the function: -:py:attr:`espressomd.thermostat.Thermostat.set_dpd` -which takes :math:`k_\mathrm{B} T` as the only argument. +The temperature is set via +:py:meth:`espressomd.thermostat.Thermostat.set_dpd` +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 see -there. +: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 @@ -328,15 +392,6 @@ 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 samples/dpd.py. -DPD adds a velocity dependent dissipative force and a random force -to the conservative pair forces. - -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. - 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 @@ -369,10 +424,7 @@ and the parameters: * ``ext_pressure``: (float) The external pressure as float variable. * ``piston``: (float) The mass of the applied piston as float variable. -This thermostat is based on the Andersen thermostat (see -:cite:`andersen80a,mann05d`) and will thermalize the box -geometry. It will only do isotropic changes of the box. -See this code snippet for the two commands:: +For example:: import espressomd @@ -380,6 +432,8 @@ See this code snippet for the two commands:: system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0) system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) +For an explanation of the algorithm involved, see :ref:`Isotropic NPT thermostat` + 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. @@ -389,7 +443,7 @@ behavior, please contribute to solving the problem. CUDA ---- -:py:attr:`~espressomd.cuda_init.CudaInitHandle()` command can be used to choose the GPU for all subsequent +:py:meth:`~espressomd.cuda_init.CudaInitHandle()` command can be used to choose the GPU for all subsequent GPU-computations. Note that due to driver limitations, the GPU cannot be changed anymore after the first GPU-using command has been issued, for example ``lbfluid``. If you do not choose the GPU manually before that, diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index 9a1a201176d..00321264741 100644 --- a/doc/sphinx/zrefs.bib +++ b/doc/sphinx/zrefs.bib @@ -1,5 +1,15 @@ % This file was created with JabRef 2.9. % Encoding: US-ASCII +@article{ahlrichs99, +author = {Ahlrichs,Patrick and Dünweg,Burkhard }, +title = {Simulation of a single polymer chain in solution by combining lattice Boltzmann and molecular dynamics}, +journal = {The Journal of Chemical Physics}, +volume = {111}, +number = {17}, +pages = {8225-8239}, +year = {1999}, +doi = {10.1063/1.480156} +} @ARTICLE{allen06, author = {Allen, Michael P.}, @@ -364,20 +374,6 @@ @article{Gompper1996 pages = {1305--1320} } -@ARTICLE{grest86a, - author = {Gary S. Grest and Kurt Kremer}, - title = {Molecular dynamics simulation for polymers in the presence of a heat - bath}, - journal = {Phys. Rev. A}, - year = {1986}, - volume = {33}, - pages = {3628--31}, - number = {5}, - doi = {10.1103/PhysRevA.33.3628}, - owner = {olenz}, - timestamp = {2007.06.13} -} - @ARTICLE{hickey10a, author = {Hickey, Owen A. and Holm, Christian and Harden, James L. and Slater, Gary W.}, @@ -680,6 +676,18 @@ @ARTICLE{Nikunen03 pages = {407} } +@article{omelyan98, +author = {Omelyan,Igor P. }, +title = {On the numerical integration of motion for rigid polyatomics: The modified quaternion approach}, +journal = {Computers in Physics}, +volume = {12}, +number = {1}, +pages = {97-103}, +year = {1998}, +doi = {10.1063/1.168642}, +URL = {https://aip.scitation.org/doi/abs/10.1063/1.168642} +} + @ARTICLE{pasichnyk04a, author = {Igor Pasichnyk and Burkhard D{\"u}nweg}, title = {Coulomb interactions via local dynamics: {A} molecular-dynamics algorithm}, @@ -765,6 +773,16 @@ @ARTICLE{ramirez10a unique-id = {{ISI:000283359300008}} } +@book{rapaport04, + author = {Rapaport, D. C.}, + title = {The Art of Molecular Dynamics Simulation}, + year = {2004}, + isbn = {0521825687, 9780521825689}, + edition = {2nd}, + publisher = {Cambridge University Press}, + address = {New York, NY, USA}, +} + @article{reed92a, title={{M}onte {C}arlo study of titration of linear polyelectrolytes}, author={Reed, Christopher E and Reed, Wayne F},