Additional details on DPD¶
Scaling DPD calculations¶
Dissipative Particle Dynamics (DPD) simulations are often described in terms of three unit quantities:
- Mass (of a bead), \(m\)
- Length (interaction cutoff distance), \(r_c\)
- Energy (system temperature), \(k_B T\)
each of which are normally set to 1 in DPD calculations and can be defined based on what we want to model. Starting with the representation of a particular chemical species as DPD particles or ‘beads’, we can obtain our mass unit from the contents of a particular particle. The length unit is also obtained from the particle contents, along with the numerical density of the DPD calculation \(\rho\), while the energy unit is normally set to the product of the operating temperature and the Boltzmann constant \(k_B\).
For the example of modelling water at room temperature (298 K) with \(N_m\) molecules in each DPD bead, given that the molecular mass of water is \(M_r = 18.01528 \times 10^{-3}\) kg mol\(^{-1}\), our mass unit is equal to:
while the energy unit is equal to:
A ‘real-life’ density for water of \(\rho_m = 996.95\) kg m\(^{-3}\) gives the physical volume of this bead as:
If we use a numerical density \(\rho = 3\) in our simulation [1], the unit volume (related to the interaction cutoff \(r_c\)) becomes:
and thus the length unit
There are two possible approaches to obtaining a time unit \(\tau_0\) for the simulation. The first is to apply dimensional analysis and make use of the three units defined above:
The second is to compare the self-diffusivity of the DPD beads - measured by calculating mean-squared displacements from simulations - with experimental values for the fluid we wish to model [Groot2001]. In our case of modelling water beads each with \(N_m\) molecules and setting the dissipative force parameter \(\gamma = 4.5\), the time unit is found to be:
Parameterising DPD interactions¶
Groot and Warren’s [Groot1997] form of conservative interaction:
is able to specify both an equation of state and different degrees of repulsion between particle (bead) species. The values of \(A_{ij}\) for different bead pairs can be chosen or fitted to give specific behaviours or simulation outcomes. An example of choosing \(A_{ij}\) values can include specifying degrees of hydrophobicity for different bead species, e.g. those found in amphiphilic molecules with hydrophilic and hydrophobic sections.
Experimental properties or outcomes that can be used to fit \(A_{ij}\) values can include expected phases or structures (e.g. micelles, vesicles, bilayers) at given concentrations, critical micelle concentrations (CMCs) and/or partition coefficients (ratios of concentrations in two immiscible phases) [Anderson2017].
Alternative and additional approaches for interactions between DPD particles include many-body DPD, bond interactions (stretching, fixed-length holonomic constraints, angles, dihedrals etc.) and electrostatics.
Like-like interactions¶
The equation of state resulting from the Groot-Warren conservative interaction has been found from viral theorem to be:
where \(\alpha \approx 0.101 r_c^{4}\) when the bead density \(r_c^3 \rho > 2\). While it is not especially realistic in describing a real fluid’s behaviour, we can still use it to parameterise interactions between particles of the same species by deriving a relationship with isothermal compressibility \(\kappa_T\), the rate of change of density with pressure at constant temperature. The derivative of (1) with respect to density at constant temperature gives the following relationship:
Applying this to water at room temperature - which has an isothermal compressibility of \(\kappa_T = 4.533 \times 10^{-10}\) Pa\(^{-1}\) - with one molecule per bead (\(N_m = 1\)), we can determine a reciprocal of dimensionless isothermal compressibility of \(\kappa^{-1} \approx 16.09\). Rearrangement of the above equation with this value then gives the conservative force parameter as:
Setting the bead density as \(\rho = 3 r_c^{-3}\) leads to \(A_{ij} = 25 k_B T r_c^{-1}\), which is very commonly used in DPD simulations. Alternatively, if \(N_m = 3\) (three water molecules per bead), the conservative force parameter would end up being \(A_{ij} = 78 k_B T r_c^{-1}\) [Groot2001].
Flory-Huggins solution theory¶
To parameterise conservative force parameters \(A_{ij}^{AB}\) between different particle species (A and B), we can use an approach based on the Flory-Huggins solution theory of polymers [Groot1997]. Assuming beads of species A and B are of equal size and molecules are made from these beads, we can express a Gibbs free energy of mixing based on the volume fractions of each species (\(\phi_A\), \(\phi_B\)) and the number of beads per molecule (\(n_A\), \(n_B\)):
where we introduce a parameter \(\chi^{\text{AB}}\) to represent the non-ideal part of the mixing energy: the higher the value, the more the components will repel each other. Groot and Warren found that the \(\chi^{\text{AB}}\)-parameter can be given in terms of the interaction parameter between the two species by the following predicted relationship:
provided that the two species (A and B) interact among themselves in the same way (i.e. \(A_{ij}^{\text{AA}} = A_{ij}^{\text{BB}}\)) [2]. The value of \(A_{ij}^{\text{AB}}\) thus determines how well the two components will mix together.
In real DPD simulations, while the linear relationship between \(\chi^{\text{AB}}\) and \(\left( A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}} \right)\) holds, the proportionality constant is generally not equal to \(2\alpha\rho/k_B T\). We can find the actual proportionality constant for a given particle density from series of simulations with separating beads or molecules, each with different values of \(A_{ij}^{\text{AB}}\). The volume fraction of one species (e.g. \(\phi_{A}\) for species A) in separated regions can be used to determine \(\chi^{\text{AB}}\):
Values of \(\chi^{\text{AB}}\), either experimentally determined or estimated from atomistic molecular dynamics simulations, can then be used to determine the required values of \(A_{ij}^{\text{AB}}\).
Infinite dilution activity coefficients¶
Another approach to parameterise DPD interactions between different particle species is to match up \(A_{ij}^{\text{AB}}\) with infinite dilution activity coefficients \(\gamma^{\infty}\) [Vishnyakov2013]. The activity coefficient describes non-ideal mixing behaviour and the effect this has on chemical potential compared to the value for a pure component \(i\):
with the infinite dilution activity coefficient \(\gamma^{\infty}_i\) as the limiting value when the mole fraction \(x_i\) approaches zero. This value can be obtained for real systems either from experimental values or by using predictive models based on functional groups (e.g. UNIFAC).
To find the relationship between \(A_{ij}^{\text{AB}}\) and \(\gamma^{\infty}\), we can carry out Widom trial insertions of each bead type into simulation configurations of pure components. The excess chemical potentials \(\mu_i - \mu_i^{\circ}\) of inserting molecules of A and B into pure A are replaced with logarithms of block-averaged exponentials of inserted energies, thus:
where \(n_i\) is the number of beads in a molecule of component \(i\).
For Groot-Warren interactions, the block-averages in the above expression end up being directly relatable to conservative force parameters and thus give linear relationships between \(\ln \gamma^{\infty}_B\) and \(A_{ij}^{AB} - A_{ij}^{AA}\). Unlike the Flory-Huggins approach above, no assumption is made about the like-like parameters \(A_{ij}^{AA}\) and \(A_{ij}^{BB}\), which can differ for the various particle types.
Many-body DPD¶
While Groot-Warren interactions can be used for many mesoscopic simulation systems, they are not always sufficient to accurately describe the thermodynamic behaviour of real systems. In particular, their resulting quadratic equation of state (1) cannot describe fluids with multiple coexisting phases, e.g. liquid and vapour.
One approach that can provide better thermodynamics at the mesoscale is ‘multibody’ or many-body DPD [Pagonabarraga2001] [Trofimov2002], which directly specifies the free energy of fluids. The total free energy of a system can be expressed as the sum of free energies per particle \(\psi\), which are dependent on localised densities \(\tilde{\rho}\) that can be calculated for each particle, i.e.
The free energy per particle can be split into ideal, bulk and excess terms. Of these contributions, only the excess part of the free energy is involved in effective (conservative) interactions between particles, which can be expressed in a pairwise form:
where \(w^{C}\) is a function of distance. The form of excess free energy per particle \(\psi^{ex}\) has an effect on the resulting equation of state for a single fluid:
The localised density can be calculated by summing up contributions for another function of distance:
which can be restricted to another cutoff distance, i.e. \(r_{ij} < r_d\). The two functions are related to each other by:
The ‘standard’ Groot-Warren potential can be recovered by setting the excess free energy per particle to
and the function for localised density calculations to
A simple many-body DPD model to give a cubic equation of state and vapour-liquid coexistence can be established by using an additional density term for the excess free energy [Warren2003] with two different localised densities, i.e.
where \(r_d = r_c\) for \(w^{\rho} \left(r_{ij}\right)\) when calculating the first localised density \(\overline{\rho}\). The resulting pairwise conservative force can be simplified to the standard Groot-Warren form with an additional density-based term:
noting that \(A_{ij} \equiv A\) can be specified for different pairs of species and set to negative values to provide attractions between particles, while the repulsive parameter \(B\) must be kept constant for all species to ensure the forces are conservative [Warren2013a]. These two parameters can be chosen or fitted for e.g. both isothermal compressibility (as per standard Groot-Warren interactions) and interfacial tension.
Inclusion of electrostatics¶
The long-range nature of charge-based interactions can have an effect on systems operating at the mesoscale, and as such there are circumstances when these need to be included in DPD simulations. Techniques used to solve the Poisson equation for molecular dynamics simulations, e.g. Particle-Particle Particle-Mesh (PPPM) and Ewald summation, can be applied to DPD simulations.
One important aspect to consider for DPD electrostatics is the softness of bead-bead interactions. While the hard-core interactions typically used in atomistic and coarse-grained molecular dynamics can prevent oppositely-charged point charges from collapsing on top of each other, this is not automatically guaranteed for e.g. Groot-Warren soft repulsive interactions [3].
To avoid soft ion collapse, point charges are often replaced with charge functions to spread out each charge over a finite volume. These smearing functions are selected to act entirely at shorter ranges, with the potential reverting back to standard Coulombic for point charges at larger distances. This property enables the long-range parts of electrostatic interactions to be calculated using techniques common to molecular dynamics without modifications.
Two approaches to including electostatics in DPD simulations have been applied. The first uses a form of PPPM [Groot2003], which starts by assigning smeared charge distributions to a grid. Each charge \(q_i\) is multiplied by a smearing function \(f(r)\), e.g.
where \(r_e\) is the electrostatic smearing radius, to give a charge concentration \(\rho^{*}\) of cations minus anions, before assigning values to a regular grid representing the system volume. The electrostatic field \(\psi\) is then calculated at all grid points to satisfy the Poisson equation:
where \(p(r)\) is the local polarisability relative to the background fluid and \(\Gamma\) is a dimensionless permittivity coefficient equal to
where \(e\) is the charge on an electron, \(\epsilon_0\) the permittivity of a vacuum and \(\epsilon_r\) is the relative permittivity of the background fluid. For water at room temperature, the permittivity coefficient is related to the number of molecules per water bead, i.e.
The electrostatic field can be solved either by using Fourier transforms or iteratively by real-space successive overdamped relaxation (SOR). The electrostatic force acting on an ion is then found by finding the gradient of the electrostatic field and using the smearing function over grid points within \(r_e\) of the particle centre:
The other approach for DPD electrostatics involves a modified form of the Ewald sum [GonzalezMelchor2006]. The smearing function is now defined as a modification to the Coulombic potential:
where \(q_i\) is the charge (valence) on ion \(i\). Ewald summation splits the potential into real-space (short range) and reciprocal-space (long range) parts, making use of the known Fourier transform for a Gaussian (normal) distribution. Since the charge smearing function is designed to only be applicable at shorter ranges, only the real-space terms for potential (shown below) and corresponding forces need to be modified:
where \(\alpha\) is the Ewald sum real-space convergence parameter. The reciprocal space part of the Ewald sum can be solved using any viable method: standard (analytical) calculations, PPPM, (Smooth) Particle Mesh Ewald etc.
Several choices of smearing function for Ewald sums are possible and have been used in DPD calculations. One choice (originally proposed) is a truncated form of a Slater-type function:
where \(\beta\) is related to a charge decay length. Another choice is to use a Gaussian function [Warren2013b], i.e.
where \(\sigma_G\) is the screening length. If \(\sigma_G = \frac{1}{2 \alpha}\), then all real-space potentials and forces reduce to zero and the calculation can be carried out entirely in reciprocal space.
Transport properties¶
As a method designed to give both temperature control and both local and global momentum conservation (Galilean invariance), DPD as a simulation method is well-suited to examine how chemical species behave when undergoing flow and determine their rheological behaviour.
The rheological behaviour of a species can be determined from how it responds to being deformed, i.e. measuring the shear stress required to apply a particular shear rate or velocity gradient. The relationship between shear stress and shear rate can be determined from DPD simulations by using non-equilibrium MD techniques, such as the application of linear shear (constant velocity gradient).
The random and dissipative forces in the DPD thermostat have an effect on this rheological behaviour. Alternative pairwise thermostats have subsequently been proposed to provide even more realistic hydrodynamic behaviours for denser fluids.
Measuring viscosity using linear shear¶
The dynamic viscosity of a fluid is defined as the ratio of shear stress \(\tau\) (shearing force per unit cross-sectional area) to the rate of shear deformation \(\dot{\gamma}\) (gradient of velocity):
where the shear stress is a result of the applied shear rate.
The relationship between shear stress and shear rate can vary from fluid to fluid. The most common relationship - particularly for simple fluids - assumes the dynamic viscosity \(\mu\) remains constant for all shear rates, describing a Newtonian fluid. 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 start flowing: these describe shear thickening (dilatant), shear thinning (pseudo plastic) and Bingham plastic fluids.
To determine the viscosity of a fluid, we need to find a relationship between applied shear rates (velocity gradients) and resulting shear stresses, which we can work out by running some DPD simulations.
The Lees-Edwards periodic boundary condition [Lees1972] can apply a linear shear flow field to a particle simulation. It works by both changing the velocities and tangential positions of particles that move through the given periodic boundary. For example, to apply wall velocities of \(-\mathbf{v}_{w}\) and \(+\mathbf{v}_{w}\) at the bottom and top boundaries of a box, any particle moving through the top boundary would be shifted by a displacement \(-\mathbf{v}_{w} t\) before re-entering the box - where \(t\) is the time since shearing started - and its velocity modified by \(-2\mathbf{v}_{w}\). Similarly, a particle moving through the bottom boundary would be displaced by \(+\mathbf{v}_{w} t\) and its velocity modified by \(+2\mathbf{v}_{w}\).
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 and dividing by the system volume, i.e.
which can be averaged over a reasonably large number of timesteps to obtain representative values. For a system with shear flow, only one of the symmetric off-diagonal terms (\(\sigma_{\alpha \beta} = \sigma_{\beta \alpha}\), \(\alpha \ne \beta\)) is required, which will be equal to the negative of the shear stress \(\tau\). For instance, if the \(x\)-component of velocity varies along \(y\) (so the velocity gradient is \(\frac{\partial v_x}{\partial y}\)), the shear stress is equal to either \(-\sigma_{yx}\) or \(-\sigma_{xy}\).
We can therefore carry out several simulations with different velocity gradients (by changing the wall velocities \(\pm\mathbf{v}_{w}\)) and find the resulting shear stresses. Plotting the stresses against shear rates will give the relationship between the two, and the gradient of the curve at each point will give the fluid viscosity. If the fluid is Newtonian, the curve will be a straight line crossing through the origin (i.e. zero stress at zero shear rate).
Alternative pairwise thermostats to DPD¶
The dissipative force parameter \(\gamma\) in the DPD thermostat is one way we can control the dynamic viscosity \(\mu\). There are, however, limitations on how much we can control the viscosity using this thermostat.
Assuming the random force has a similar linear dependence to that used by Groot-Warren conservative forces, i.e. for (1)-(4) in Dissipative Particle Dynamics (DPD):
and also not applying conservative forces, we can obtain the following relationship for the viscosity [Marsh1997]:
which is pretty complex and might not really be suited for modelling flows of dense liquids [4]. Adding conservative forces will change this relationship and is likely to increase the viscosity further, 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, particularly if you are more interested in equilibrium structures than the detailed hydrodynamics.
One way around this is to use a different pairwise thermostat. Two particularly interesting choices for flow-based systems are the Lowe-Andersen [Lowe1999] and Stoyanov-Groot [Stoyanov2005] thermostats, both of which are pairwise variants of the Andersen thermostat [Andersen1980] and use a collision frequency \(\Gamma\) as their main parameter.
For each interacting particle pair, we generate a random number \(u_{ij}\) between 0 and 1, and if \(u_{ij} < \Gamma \Delta t\), we select that pair for thermostatting after integration of conservative (interaction) forces. Going through the selected particle pairs in a randomised order and using the most up-to-date values of particle velocities, we apply the following modifications to each pair’s velocities:
where \(\mu_{ij} = \frac{m_i m_j}{m_i + m_j}\) is the reduced mass between the particles and \(v_{ij}^{\circ} = \zeta_{ij} \sqrt{\frac{k_B T}{\mu_{ij}}}\) is a randomised replacement relative velocity. It can be shown that the viscosity for a single species obtained by these thermostats is:
which is clearly simpler than the expression for the DPD thermostat in (3) and can obtain larger viscosities for dense liquids.
The Stoyanov-Groot thermostat adds additional pairwise forces for pairs that are not selected for thermostatting:
where \(\alpha^{T}\) is an additional system-wide coupling parameter for the thermostat, \(w^{T}\) is a weighting function similar to that used for DPD random forces and \(T^{*}\) is an instantaneous temperature calculated using another pairwise weighting function for all interacting particle pairs. These additional pairwise Nosé-Hoover forces can ensure better temperature control, particularly when \(\Gamma\) is low and not many particle pairs are selected for Lowe-Andersen thermostatting.
Footnotes
| [1] | The numerical density shown here has units of \(r_c^{-3}\) (i.e. a number per unit volume). This particular value is frequently used in DPD simulations with Groot-Warren conservative interactions as this is close to the lowest density for which the constant (density-insensitive) equation of state applies. |
| [2] | This follows from the assumption in Flory-Huggins solution theory that all particles (solvent and monomers) are equal in size. |
| [3] | Higher values of \(A_{ij}\) can make the use of point charges in DPD simulations viable, particularly if a larger degree of coarse-graining is in use [Vaiwala2017], but greater care is needed when starting such simulations to ensure oppositely-charged particles do not start too close together. |
| [4] | Dividing the dynamic viscosity by fluid density to give the kinematic viscosity \(\nu\) and dividing this property by the fluid diffusivity \(D\) gives the dimensionless Schmidt number \(\mathrm{Sc}\). For DPD simulations with no conservative forces, this value is typically around 1: fine for gases but generally too low for liquids, whose values are typically around 1000. |
References
| [Groot2001] | (1, 2) 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. |
| [Groot1997] | (1, 2) RD Groot and PB Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, Journal of Chemical Physics, 107, p. 4423–4435, 1997, doi: 10.1063/1.474784. |
| [Vishnyakov2013] | A Vishnyakov, M-T Lee and AV Neimark, Prediction of the critical micelle concentration of nonionic surfactants by dissipative particle dynamics simulations, Journal of Physical Chemistry Letters, 4, p. 797-802, 2013, doi: 10.1021/jz400066k. |
| [Anderson2017] | RL Anderson, DJ Bray, AS Ferrante, MG Noro, IP Stott and PB Warren, Dissipative particle dynamics: systematic parametrization using water-octanol partition coefficients, Journal of Chemical Physics, 147, 094503, 2017. doi: 10.1063/1.4992111. |
| [Pagonabarraga2001] | I Pagonabarraga and D Frenkel, Dissipative particle dynamics for interacting systems, Journal of Chemical Physics, 115, p. 5015-5026, 2001. doi: 10.1063/1.1396848. |
| [Trofimov2002] | SY Trofimov, ELF Ries and MAJ Michels, Thermodynamic consistency in dissipative particle dynamics simulations of strongly nonideal liquids and liquid mixtures, Journal of Chemical Physics, 117, p. 9383-9394, 2002. doi: 10.1063/1.1515774. |
| [Warren2003] | PB Warren, Vapor-liquid coexistence in many-body dissipative particle dynamics, Physical Review E, 68, 066702, 2003. doi: 10.1103/PhysRevE.68.066702. |
| [Warren2013a] | PB Warren, No-go theorem in many-body dissipative particle dynamics, Physical Review E, 87, 045303, 2013. doi: 10.1103/PhysRevE.87.045303. |
| [Groot2003] | RD Groot, Electrostatic interactions in dissipative particle dynamics - simulation of polyelectrolytes and anionic surfactants, Journal of Chemical Physics, 118, p. 11265-11277, 2003, doi: 10.1063/1.1574800. |
| [GonzalezMelchor2006] | M González-Melchor, E Mayoral, ME Velázquez and J Alejandre, Electorstatic interactions in dissipative particle dynamics using the Ewald sums, Journal of Chemical Physics, 125, 224107, 2006, doi: 10.1063/1.2400223. |
| [Warren2013b] | PB Warren, A Vlasov, L Anton and AJ Masters, Screening properties of Gaussian electrolyte models, with application to dissipative particle dynamics, Journal of Chemical Physics, 138, 204907, 2013, doi: 10.1063/1.4807057. |
| [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. |
| [Andersen1980] | HC Andersen, Molecular dynamics simulations at constant pressure and/or temperature, Journal of Chemical Physics, 72, p. 2384-2393, 1980, doi: 10.1063/1.439486. |
| [Vaiwala2017] | R Vaiwala, S Jadhav and R Thaokar, Electrostatic interactions in dissipative particle dynamics - Ewald-like formalism, error analysis, and pressure computation, Journal of Chemical Physics, 146, 124904, 2017, doi: 10.1063/1.4978809. |