From 6e3b0862e6a5b5244ee197161fcd31f67f415435 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Thu, 13 Jun 2019 17:50:23 +0200 Subject: [PATCH 1/8] VV and LD explained, no papers yet --- doc/sphinx/running.rst | 130 ++++++++++++++++++++++-------------- doc/sphinx/system_setup.rst | 59 +++++++++------- 2 files changed, 114 insertions(+), 75 deletions(-) diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index 89232041ac1..cfe9f2fb095 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -3,25 +3,52 @@ 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 basic integration scheme of |es| is the Velocity Verlet algorithm, +but 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 TODO:ref:thermostats. + +For numerical integration, this equation is discretized to the following steps: + +1. Calculate the velocity at the half step -|es| uses the Velocity Verlet algorithm for the integration of the equations -of motion. +.. math:: v(t+dt/2) = v(t) + F(x(t),v(t-dt/2),t)/m dt/2 + +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) + 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 +76,19 @@ 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: - -Run steepest descent minimization ---------------------------------- - -: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 +.. _Rotational degrees of freedom and particle anisotropy: +Rotational degrees of freedom and particle anisotropy +----------------------------------------------------- -.. _Integrating rotational degrees of freedom: +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. -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. +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 TODO:ref:"On the numerical integration of motion for rigid polyatomics: +The modified quaternion approach", Omeylan, Igor (1998), Eq. 12." + 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 +124,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 2393cb97692..60060319298 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -263,18 +263,23 @@ 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 the particles 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(t) + \eta(t). + +Here, :math:`f_i` are all deterministic forces from interactions, :math:`\gamma` the friction coefficient and :math:`\eta` a random, "thermal" force. +The random force mimics collisions of the particle with surrounding solvent molecules at temperature :math:`T` and satisfies + +.. math:: <\eta(t)> = 0 , <\eta_i(t)\eta_j(t')> = 2 k_B T \gamma \delta_{ij}\delta(t-t') + +(:math:`<\cdot>` denotes the ensemble average). + +In the Velocity Verlet discretisation scheme, this thermostat only enters in the force calculation. The random process :math:`\eta(t)` is realized by drawing an uncorrelated random number :math:`\eta*` from a uniform distribution for each component of all the particle forces. It is then scaled and shifted such that the mean and variance match the required properties + +.. math:: <\eta*> = 0 , <\eta*\eta*> = 2 k_B T \gamma \delta_{ij}/dt + +TODO:find paper/book that uses this approach 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 @@ -282,27 +287,33 @@ can be omitted in subsequent calls of ``set_langevin()``. 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). +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. + +The diffusion coeffictient :math:`D` of the particle can be obtained by the Einstein-Smoluchowski relation + +.. math:: D = \frac{k_B T}{\gamma}. + +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`. + 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. -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. +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.). + + .. _Dissipative Particle Dynamics (DPD): From 350ad4c1959e2db01a684ee76e753e4dcaabc096 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Fri, 14 Jun 2019 11:00:34 +0200 Subject: [PATCH 2/8] found some refs --- doc/sphinx/running.rst | 2 ++ doc/sphinx/system_setup.rst | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index cfe9f2fb095..8dd7999b775 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -46,6 +46,8 @@ For numerical integration, this equation is discretized to the following steps: .. math:: v(t+dt) = v(t+dt/2) + F(x(t+dt),t+dt)/m dt/2 +TODO:ref:Rapaport 3.5.8-3.5.10 or https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet + 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 diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 60060319298..155d9507897 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -275,7 +275,7 @@ The random force mimics collisions of the particle with surrounding solvent mol (:math:`<\cdot>` denotes the ensemble average). -In the Velocity Verlet discretisation scheme, this thermostat only enters in the force calculation. The random process :math:`\eta(t)` is realized by drawing an uncorrelated random number :math:`\eta*` from a uniform distribution for each component of all the particle forces. It is then scaled and shifted such that the mean and variance match the required properties +In the Velocity Verlet discretisation scheme, this thermostat only enters in the force calculation. The random process :math:`\eta(t)` is dicretized by drawing an uncorrelated random number :math:`\eta*` from a uniform distribution for each component of all the particle forces. It is then scaled and shifted such that the mean and variance match the required properties TODO:ref:eqs12-14 in The Nature of Folded States of Globular Proteins J. D. HONEYCUTT’ and D. THlRUMALAl .. math:: <\eta*> = 0 , <\eta*\eta*> = 2 k_B T \gamma \delta_{ij}/dt From 9f745c509d35e0addf77282d3baa80fbf0020051 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Wed, 19 Jun 2019 17:13:35 +0200 Subject: [PATCH 3/8] DPD, NPT, LB now with references --- doc/sphinx/inter_non-bonded.rst | 78 +++++++++++++-------- doc/sphinx/lb.rst | 70 +++++++++---------- doc/sphinx/running.rst | 113 +++++++++++++++++++++++++++--- doc/sphinx/system_setup.rst | 118 +++++++++++++++++++++++--------- doc/sphinx/zrefs.bib | 32 +++++++++ 5 files changed, 302 insertions(+), 109 deletions(-) diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index 9f7a3e9e9e5..28eb20a68f0 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/doc/sphinx/inter_non-bonded.rst @@ -559,60 +559,78 @@ 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:attr:`espressomd.thermostat.Thermostat.set_dpd()` +. The parameters 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 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 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/lb.rst b/doc/sphinx/lb.rst index 42fbb517561..b5091d95103 100644 --- a/doc/sphinx/lb.rst +++ b/doc/sphinx/lb.rst @@ -121,53 +121,53 @@ 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:attr:`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:: -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`: + lb.set_interpolation_order('linear') + lb.set_interpolation_order('quadratic') -.. 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 8dd7999b775..8e79241f5d8 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -9,8 +9,9 @@ To run the integrator call the method 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 basic integration scheme of |es| is the Velocity Verlet algorithm, -but a steepest descent implementation is also available. +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: @@ -26,13 +27,13 @@ The equations of motion for the trajectory of pointlike particles read 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 TODO:ref:thermostats. +as well as non-deterministic contributions described in :ref:`Thermostats`. -For numerical integration, this equation is discretized to the following steps: +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) + F(x(t),v(t-dt/2),t)/m dt/2 +.. math:: v(t+dt/2) = v(t) + \frac{F(x(t),v(t-dt/2),t)}{m} dt/2 2. Calculate the new position @@ -44,9 +45,7 @@ For numerical integration, this equation is discretized to the following steps: 4. Calculate the new velocity -.. math:: v(t+dt) = v(t+dt/2) + F(x(t+dt),t+dt)/m dt/2 - -TODO:ref:Rapaport 3.5.8-3.5.10 or https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet +.. 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 in step 1. That is, they are computed once in step 3, @@ -78,7 +77,102 @@ 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 thermostat integrator: + +Isotropic NPT thermostat +------------------------ + +: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: + + * ``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. + +Additionally, a NPT thermostat has to be set by :py:func:`~espressomd.thermostat.Thermostat.set_npt()` +with parameters: + + * ``kT``: (float) Thermal energy of the heat bath + * ``gamma0``: (float) Friction coefficient of the bath + * ``gammav``: (float) Artificial friction coefficient for the volume fluctuations. + +A code snippet would look like:: + + import espressomd + + system = espressomd.System() + system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0) + system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0) + +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 + +2. Calculate the instantaneous pressure and "volume momentum" + +.. 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 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: @@ -88,8 +182,7 @@ 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 TODO:ref:"On the numerical integration of motion for rigid polyatomics: -The modified quaternion approach", Omeylan, Igor (1998), Eq. 12." +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. diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 155d9507897..4c75f8d289a 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -264,22 +264,29 @@ Best explained in an example:: As explained before the temperature is set as thermal energy :math:`k_\mathrm{B} T`. -The Langevin thermostat is based on an extension of the particles equation of motion to +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(t) + \eta(t). +.. 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 friction coefficient and :math:`\eta` a random, "thermal" force. -The random force mimics collisions of the particle with surrounding solvent molecules at temperature :math:`T` and satisfies +Here, :math:`f_i` are all deterministic forces from interactions, +:math:`\gamma` the friction coefficient and :math:`\eta` a random, "thermal" force. +The random force mimics collisions of the particle with surrounding solvent molecules +at temperature :math:`T` and satisfies -.. math:: <\eta(t)> = 0 , <\eta_i(t)\eta_j(t')> = 2 k_B T \gamma \delta_{ij}\delta(t-t') +.. 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). +(:math:`<\cdot>` denotes the ensemble average and :math:`\alpha,\beta` are spatial coordinates). -In the Velocity Verlet discretisation scheme, this thermostat only enters in the force calculation. The random process :math:`\eta(t)` is dicretized by drawing an uncorrelated random number :math:`\eta*` from a uniform distribution for each component of all the particle forces. It is then scaled and shifted such that the mean and variance match the required properties TODO:ref:eqs12-14 in The Nature of Folded States of Globular Proteins J. D. HONEYCUTT’ and D. THlRUMALAl +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. -.. math:: <\eta*> = 0 , <\eta*\eta*> = 2 k_B T \gamma \delta_{ij}/dt - -TODO:find paper/book that uses this approach +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 @@ -287,11 +294,6 @@ can be omitted in subsequent calls of ``set_langevin()``. 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). -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. - The diffusion coeffictient :math:`D` of the particle can be obtained by the Einstein-Smoluchowski relation .. math:: D = \frac{k_B T}{\gamma}. @@ -311,8 +313,55 @@ 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 anisotropic diffusion of anisotropic colloids (rods, etc.). +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): @@ -320,18 +369,29 @@ The keywords ``gamma`` and ``gamma_rotate`` can be specified as a scalar, or, wi 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: +The temperature is set via :py:attr:`espressomd.thermostat.Thermostat.set_dpd` -which takes :math:`k_\mathrm{B} T` as the only argument. +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 @@ -340,15 +400,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 @@ -381,10 +432,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 @@ -392,6 +440,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. diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index 9a1a201176d..30336b1dc58 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.}, @@ -680,6 +690,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 +787,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}, From cf45b9f3dff1eaa1df37f4aa43e95653b30d0de8 Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Tue, 25 Jun 2019 08:55:53 +0200 Subject: [PATCH 4/8] NPT direction, typo --- doc/sphinx/inter_non-bonded.rst | 2 +- doc/sphinx/running.rst | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index 28eb20a68f0..b9309806370 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/doc/sphinx/inter_non-bonded.rst @@ -577,7 +577,7 @@ to a temperature set by :py:attr:`espressomd.thermostat.Thermostat.set_dpd()` * ``trans_r_cut`` which will be explained below. The interaction -only has an effect if the DPD thermostat activated. +only has an effect if the DPD thermostat is activated. 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 diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index 8e79241f5d8..ad8f0898d61 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -170,7 +170,8 @@ The discretisation consists of the following steps (see :cite:`kolb99a` for a fu Notes: -* The for the instantaneous pressure the same limitations of applicability hold as described in :ref:`Pressure` +* 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. From f233b4aef022bd1302dbaff95e40baa23b67cceb Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Tue, 25 Jun 2019 10:26:09 +0200 Subject: [PATCH 5/8] DPD temp, lb interp bounds, sphinx :meth: --- doc/sphinx/inter_non-bonded.rst | 8 +++++--- doc/sphinx/lb.rst | 9 ++++++++- doc/sphinx/system_setup.rst | 12 ++++++------ 3 files changed, 19 insertions(+), 10 deletions(-) diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst index b9309806370..8637d82f132 100644 --- a/doc/sphinx/inter_non-bonded.rst +++ b/doc/sphinx/inter_non-bonded.rst @@ -566,8 +566,8 @@ The parameters can be set via:: 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:attr:`espressomd.thermostat.Thermostat.set_dpd()` -. The parameters are +to a temperature set by :py:meth:`espressomd.thermostat.Thermostat.set_dpd()`. +The parameters for the interaction are * ``gamma`` * ``weight_function`` @@ -592,7 +592,9 @@ for the dissipative force and .. math:: \vec{F}_{ij}^R = \sqrt{2 k_B T \gamma_\parallel w_\parallel (r_{ij}) } \eta_{ij}(t) \hat{r}_{ij} for the random force. This introduces the friction coefficient :math:`\gamma_\parallel` (parameter ``gamma``) and the weight function -:math:`w_\parallel`. The weight function can be specified via the ``weight_function`` switch. +: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`: diff --git a/doc/sphinx/lb.rst b/doc/sphinx/lb.rst index b5091d95103..10bf7018911 100644 --- a/doc/sphinx/lb.rst +++ b/doc/sphinx/lb.rst @@ -135,7 +135,7 @@ To get interpolated velocity values between lattice nodes, the function:: with a single position ``pos`` as an argument can be used. For the GPU fluid :class:`espressomd.lb.LBFluidGPU` -also :py:attr:`espressomd.lb.LBFluidGPU.get_interpolated_fluid_velocity_at_positions()` +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, @@ -146,6 +146,13 @@ one of:: lb.set_interpolation_order('linear') lb.set_interpolation_order('quadratic') + +A note on boundaries: + +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. .. _Coupling LB to a MD simulation: diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 4c75f8d289a..8e1ad53b53f 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -141,7 +141,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. :: @@ -167,7 +167,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 @@ -206,7 +206,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 @@ -252,7 +252,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:: @@ -380,7 +380,7 @@ interaction between the particles that make up the fluid, see :ref:`DPD interaction`, and a repulsive conservative force. The temperature is set via -:py:attr:`espressomd.thermostat.Thermostat.set_dpd` +:py:meth:`espressomd.thermostat.Thermostat.set_dpd` which takes ``kT`` as the only argument. The friction coefficients and cutoff are controlled via the @@ -451,7 +451,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, From 1fe1af9ffe981532dcb4bfcd4b23937efeef9e2d Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Tue, 25 Jun 2019 10:57:09 +0200 Subject: [PATCH 6/8] langevin: rem diffusion coeff --- doc/sphinx/system_setup.rst | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 8e1ad53b53f..626b232f022 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -269,7 +269,9 @@ 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 friction coefficient and :math:`\eta` a random, "thermal" force. +:math:`\gamma` the bare friction coefficient and :math:`\eta` a random, "thermal" force. +The friction term accounts for dissipation in a surrounding fluid (the actual force +response of a particle will also depend on the interaction with other particles). The random force mimics collisions of the particle with surrounding solvent molecules at temperature :math:`T` and satisfies @@ -294,14 +296,6 @@ can be omitted in subsequent calls of ``set_langevin()``. 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). -The diffusion coeffictient :math:`D` of the particle can be obtained by the Einstein-Smoluchowski relation - -.. math:: D = \frac{k_B T}{\gamma}. - -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`. - 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 From 115f7ab516f0d6b244ac86e9db5fdb2fa67526bf Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Tue, 25 Jun 2019 14:48:18 +0200 Subject: [PATCH 7/8] reintroduced diffusion coeff, removed ref --- doc/sphinx/system_setup.rst | 12 +++++++++--- doc/sphinx/zrefs.bib | 14 -------------- 2 files changed, 9 insertions(+), 17 deletions(-) diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 626b232f022..f04519b6748 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -270,15 +270,21 @@ The Langevin thermostat is based on an extension of Newton's equation of motion 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 (the actual force -response of a particle will also depend on the interaction with other particles). -The random force mimics collisions of the particle with surrounding solvent molecules +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 a system where particles only interact with the implicitly +modelled surrounding fluid, the diffusion constant :math:`D` +is determined only by the bare friction and the temperature through +the Einstein-Smoluchowski relation + +.. math:: D = \frac{k_B T}{\gamma}. + 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 diff --git a/doc/sphinx/zrefs.bib b/doc/sphinx/zrefs.bib index 30336b1dc58..00321264741 100644 --- a/doc/sphinx/zrefs.bib +++ b/doc/sphinx/zrefs.bib @@ -374,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.}, From 47494fee516889d9c7ec45545978983a3aa2844a Mon Sep 17 00:00:00 2001 From: Christoph Lohrmann Date: Thu, 25 Jul 2019 13:01:56 +0200 Subject: [PATCH 8/8] doc: marked NPT as experimental, rem diff coeff (again) --- doc/sphinx/introduction.rst | 2 +- doc/sphinx/running.rst | 1 + doc/sphinx/system_setup.rst | 7 ------- 3 files changed, 2 insertions(+), 8 deletions(-) diff --git a/doc/sphinx/introduction.rst b/doc/sphinx/introduction.rst index 85c0ac28efa..c708f84d506 100644 --- a/doc/sphinx/introduction.rst +++ b/doc/sphinx/introduction.rst @@ -541,7 +541,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/running.rst b/doc/sphinx/running.rst index ad8f0898d61..f83bbb57ace 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -170,6 +170,7 @@ The discretisation consists of the following steps (see :cite:`kolb99a` for a fu 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`. diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index f04519b6748..7641ccabd31 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -278,13 +278,6 @@ at temperature :math:`T` and satisfies (:math:`<\cdot>` denotes the ensemble average and :math:`\alpha,\beta` are spatial coordinates). -In a system where particles only interact with the implicitly -modelled surrounding fluid, the diffusion constant :math:`D` -is determined only by the bare friction and the temperature through -the Einstein-Smoluchowski relation - -.. math:: D = \frac{k_B T}{\gamma}. - 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