DPD Exercises: Additional theory
================================
.. _DPD_Ex1Theory:
Theory to DPD Exercise 1: DPD, hydrophobicity and parameterisation
------------------------------------------------------------------
Groot and Warren's [Groot1997]_ form of conservative interaction:
.. math :: \mathbf{F}_{ij}^{C} = \left\{ \begin{matrix} A_{ij}\left( 1 - \frac{r_{ij}}{r_{c}} \right){\widehat{\mathbf{r}}}_{ij} & (r_{ij} < r_{c}) \\ 0 & (r_{ij} \geq r_{c}) \\ \end{matrix} \right.\ ,
provides a quadratic potential and the following equation of state:
.. math :: p = \rho k_B T + \alpha A_{ij} \rho^2 ~~~~~ \left(\rho>2\right)
:label: eostheory
where :math:`\alpha \approx 0.101r_{c}^{4}` and :math:`\rho` is the numerical particle density in our DPD simulation. This quadratic equation of state is not particularly realistic - a cubic equation of state would be better - but it is still possible to use it by considering the *compressibility* of fluids. By differentiating equation (1) with respect to the particle density but keeping the temperature constant, we can obtain the rescaled reciprocal of the isothermal compressibility :math:`\kappa_T`:
.. math :: \kappa^{-1} = \frac{V_{mol} N_m}{k_B T \kappa_T } = \frac{M_r N_m}{\rho_m R T \kappa_T} = \frac{1}{k_B T} \left( \frac{\partial p}{\partial \rho} \right)_T = 1 + \frac{2 \alpha A_{ij} \rho}{k_B T},
:label: compress
where :math:`N_m` is the number of molecules per bead we wish to use, and :math:`V_{mol} = \frac{M_r}{N_A \rho_m}` is the (physical) volume of a molecule of our fluid, which can be found using the molecular mass :math:`M_r` (mass per mole), Avogadro's constant :math:`N_A` (number per mole) and the density of the material :math:`\rho_m` (mass per unit volume). Putting this relationship into equation :eq:`compress` allows us to replace the quantity :math:`k_B N_A` with the more manageable universal gas constant :math:`R`.
Equation :eq:`compress` allows us to find the appropriate value of :math:`A_{ij}` for a given species. For instance, if we want to use water at room temperature (298 K) and one molecule of water per bead (:math:`N_m = 1`), assuming the required physical properties are:
- :math:`\kappa_T = 4.533 \times 10^{-10}` Pa :sup:`-1`
- :math:`M_r = 18.01528 \times 10^{-3}` kg mol :sup:`-1`
- :math:`\rho_m = 996.95` kg m :sup:`-3`
we can therefore find that the molecular volume is :math:`V_{mol} \approx 3.001 \times 10^{-29}` m :sup:`3` and the reciprocal of dimensionless isothermal compressibility is :math:`\kappa^{-1} \approx 16.09`. Rearranging equation (2), we can find that the required conservative force parameter for water is:
.. math :: A_{ij} = \frac{\left(\kappa^{-1} - 1 \right) k_B T}{2 \alpha \rho} \approx \frac{75 k_B T}{\rho}.
Given that the quantity :math:`k_B T` (the energy scale, often referred to as 'temperature') is usually set to 1 in internal units for DPD simulations and the numerical bead density is often set to :math:`\rho=3`, this gives the *very* frequently used value of :math:`A_{ij} = 25`. [#]_
.. [#] Note that we have not yet specified the mass, length or time scales of the simulation: we can do so in this case based on the desired numbers of molecules per bead and the numerical bead density. For instance, if we specify one water molecule per bead and a bead density of 3, this gives a mass per bead of:
:math:`m = \frac{N_m M_r}{N_A} \approx 2.99 \times 10^{-26}` kg,
while a unit volume equals
:math:`r_c^3 = \rho N_m V_{mol} = \frac{\rho N_m M_r}{N_A \rho_m} \approx 9.002 \times 10^{-29}` m :math:`^{3}`,
so the length scale (interaction cutoff distance) :math:`r_c \approx 4.48 \times 10^{-10}` m. We can then obtain the time scale from
:math:`\tau = r_c \sqrt{\frac{m}{k_B T}}`
which in this case is approximately :math:`1.21 \times 10^{-12}` s. A DPD timestep can typically be up to :math:`\Delta t = 0.05\tau \approx 6.04 \times 10^{-14}` s (0.0604 ps), approximately two orders of magnitude larger than typical timesteps for MD simulations (around 1 fs, :math:`10^{-15}` s).
Alternatively, the time scale can be obtained by comparing the self-diffusivity of water beads (dependent on the dissipative force parameter :math:`\gamma`) with the experimental value [Groot2001]_: for :math:`\gamma = 4.5` and one water molecule per bead, this leads to
:math:`\tau \approx 1.41 \times 10^{-11}` s.
The above parameterisation works for single fluid species, but we often want to use DPD to model how different species interact with each other. One approach we can use is to map the conservative force parameter between bead species to the Flory-Huggins solution theory of polymers. From the assumption that beads of species A and B have the same size and molecules are made up of those beads, we can express a Gibbs free energy of mixing based on the volume fractions of each species (:math:`\phi_A`, :math:`\phi_B`) and the number of beads per molecule (:math:`n_A`, :math:`n_B`):
.. math :: \frac{\Delta G^{mix}}{k_B T} = \frac{\phi_A}{n_A} \ln \phi_A + \frac{\phi_B}{n_B} \ln \phi_B + \chi^{\text{AB}} \phi_A \phi_B, ~~~~~ (3)
where we introduce a parameter :math:`\chi^{\text{AB}}` that represents the *non-ideal* part of the mixing energy. Groot and Warren [Groot1997]_ found that by making the assumption that two species (A and B) interact among themselves in the same way (i.e. :math:`A_{ij}^{\text{AA}} = A_{ij}^{\text{BB}}`), the :math:`\chi^{\text{AB}}`-parameter can be given in terms of the interaction parameter between the two species by the following predicted relationship:
.. math :: \chi^{\text{AB}} = \left( A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}} \right) \frac{2\alpha\rho}{k_B T}.
:label: chiAtheory
The value of :math:`A_{ij}^{\text{AB}}` determines how well the two components will mix together: the higher the value, the more the components will repel each other and thus the less miscible they become.
In reality, while :math:`\chi^{\text{AB}}` is found to be linearly proportional to :math:`\left( A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}} \right)` in DPD simulations, the proportionality constant is generally not equal to :math:`2\alpha\rho/k_B T`. You can work out what the proportionality constant actually is for a given particle density by carrying out a series of simulations with different values of :math:`A_{ij}^{\text{AB}}`, determining the volume fraction of one of the species (e.g. :math:`\phi_{A}` for species A) in
bead-separated regions and using the relationship for monomers (assuming :math:`n_A = n_B = 1`):
.. math :: \chi^{\text{AB}} = \frac{\ln\lbrack(1 - \phi_{A})/\phi_{A}\rbrack}{1 - 2\phi_{A}}.
:label: chiphitheory
Return to :ref:`DPD_Exercise1` to complete the exercise.
.. _DPD_Ex3Theory:
Theory to DPD Exercise 3: Transport properties of DPD fluids
------------------------------------------------------------
One important aspect of the DPD thermostat is its ability to control temperature while ensuring the hydrodynamics of the fluid represented in the simulation are correct. This makes DPD suitable for studying systems with applied flow fields such as pressure-driven flow or linear shear.
The dissipative force parameter :math:`\gamma` in the DPD thermostat is one way we can control the dynamic viscosity :math:`\mu`, an important macroscopic property of a fluid. There are limitations, however - assuming the usual linear function of distance in the random force [#]_ and **no conservative forces**, we can obtain the following relationship for the viscosity [Marsh1997]_:
.. math :: \mu \approx \frac{45 k_B T}{4 \pi \gamma r_c^3} + \frac{2 \pi \gamma \rho^2 r_c^5}{1575}
:label: viscDPDtheory
which, on the surface, is pretty complicated and probably not really suited for modelling flows of dense liquids. Adding the conservative forces will change this relationship, but it might still be difficult to obtain a realistic viscosity for a liquid with the DPD thermostat. (**Note**: This does *not* mean the DPD thermostat cannot be used for liquid systems. For instance, if you are more concerned about equilibrium structures rather than detailed hydrodynamics, the DPD thermostat is perfectly adequate.)
.. [#] Referring to equations :eq:`dissipative`-:eq:`screen` in :ref:`DPD_Intro`, :math:`w^R (r_{ij}) = 1 - \frac{r_{ij}}{r_c}` for :math:`r_{ij} < r_c` and :math:`w^R (r_{ij}) = 0` for :math:`r_{ij} \ge r_c`.
One way around this is to use a different pairwise thermostat, and at least a few have been proposed since DPD came into being. Two of these that are particularly interesting to us for flow-based systems are the Lowe-Andersen [Lowe1999]_ and Stoyanov-Groot [Stoyanov2005]_ thermostats. These are both based on a pairwise version of the Andersen thermostat: a random selection of particle pairs have their relative velocities replaced with values chosen from a Maxwell-Boltzmann distribution for the required temperature. The Stoyanov thermostat also applies additional forces based on instantaneous temperatures to other particle pairs. Both of these thermostats use a collision frequency :math:`\Gamma` as their main parameter - the fluid viscosity is thought to be linearly proportional to this parameter - and they are applied after conservative forces have been integrated.
To actually determine the viscosity of a fluid, we need to determine the relationship between shear stress and shear rate (velocity gradient), which we can work out by carrying out some simulations. We can make use of a special boundary condition available in DL_MESO_DPD. The Lees-Edwards periodic boundary condition [Lees1972]_ can apply a linear shear flow field to a particle simulation by changing the velocities and tangential position of a particle when it move through the boundary. By applying a constant velocity gradient throughout the box, we can obtain a constant shear rate for the system.
We can calculate a stress tensor for the system at a given timestep by summing products of components of pairwise forces and vectors between particle pairs, i.e.
.. math :: \sigma_{\alpha \beta} = \sum_i \left(m_i v_{i,\alpha} v_{i,\beta} + \sum_{j>i} F_{ij,\alpha} r_{ij,\beta} \right), ~~~~~ (\alpha, \beta = x, y, z)
which can be averaged over a reasonably large number of timesteps to obtain representatives values. For a system with shear flow, only one of the symmetric off-diagonal terms (:math:`\sigma_{\alpha \beta} = \sigma_{\beta \alpha}`, :math:`\alpha \ne \beta`) is required. For instance, if the :math:`x`-component of velocity is varying along :math:`y` (so the velocity gradient is :math:`\frac{\partial v_x}{\partial y}`), we will need to look at either :math:`\sigma_{yx}` or :math:`\sigma_{xy}`.
The viscosity of a fluid is equal to the derivative of shear stress :math:`\tau = |\sigma_{yx}|` with respect to shear rate :math:`\dot{\gamma} = \left|\frac{\partial v_x}{\partial y}\right|`:
.. math :: \mu = \frac{d \tau}{d \dot{\gamma}},
:label: viscositytheory
which is generally a function of the shear rate. For simple fluids, often described as *Newtonian*, this value will be a constant value regardless of shear rate.
Other more complex fluids (e.g. polymer suspensions) might have viscosities that increase or decrease as the shear rate increases, or even require a minimum shear stress to be applied before it starts to flow: these describe shear thickening (dilatant), shear thinning (pseudo plastic) and Bingham plastic fluids. Some of these effects may be observed with DPD simulations, but we are not covering this today.
Return to :ref:`DPD_Exercise3` to complete the exercise.
.. [Groot2001] RD Groot and KL Rabone, Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants, *Biophysical Journal*, **81**, p. 725-736, 2001, doi: 10.1016/S0006-3495(01)75737-2.
.. [Marsh1997] CA Marsh, G Backx and MH Ernst, Static and dynamic properties of dissipative particle dynamics, *Physical Review E*, **56**, p. 1676-1691, 1997, doi: 10.1103/PhysRevE.56.1676.