Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Documentation of VV, NPT, Thermostats #2939

Merged
merged 9 commits into from
Jul 25, 2019
Merged

Documentation of VV, NPT, Thermostats #2939

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
6e3b086
VV and LD explained, no papers yet
christophlohrmann Jun 13, 2019
350ad4c
found some refs
christophlohrmann Jun 14, 2019
9f745c5
DPD, NPT, LB now with references
christophlohrmann Jun 19, 2019
cf45b9f
NPT direction, typo
christophlohrmann Jun 25, 2019
f233b4a
DPD temp, lb interp bounds, sphinx :meth:
christophlohrmann Jun 25, 2019
1fe1af9
langevin: rem diffusion coeff
christophlohrmann Jun 25, 2019
115f7ab
reintroduced diffusion coeff, removed ref
christophlohrmann Jun 25, 2019
47494fe
doc: marked NPT as experimental, rem diff coeff (again)
christophlohrmann Jul 25, 2019
e64e5f7
Merge branch 'python' into 2680_vv_docu
jngrad Jul 25, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 50 additions & 30 deletions doc/sphinx/inter_non-bonded.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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``
jngrad marked this conversation as resolved.
Show resolved Hide resolved
* ``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
Expand Down
77 changes: 42 additions & 35 deletions doc/sphinx/lb.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <Integrator>`. 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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could you please also mention the issue of interpolating velocities near boundaries?


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:

Expand Down
Loading