## Abstract

Fluctuations of liquids at the scales where the hydrodynamic and atomistic descriptions overlap are considered. The importance of these fluctuations for atomistic motions is discussed and examples of their accurate modelling with a multi-space–time-scale fluctuating hydrodynamics scheme are provided. To resolve microscopic details of liquid systems, including biomolecular solutions, together with macroscopic fluctuations in space–time, a novel hybrid atomistic–fluctuating hydrodynamics approach is introduced. For a smooth transition between the atomistic and continuum representations, an analogy with two-phase hydrodynamics is used that leads to a strict preservation of macroscopic mass and momentum conservation laws. Examples of numerical implementation of the new hybrid approach for the multiscale simulation of liquid argon in equilibrium conditions are provided.

## 1. Introduction

Whether liquids should be described as a continuum or as a collection of atoms depends on the space and time scales being considered. One description replaces the other when the scales gradually decrease from macroscopic (metres, seconds) to microscopic (ångströms, picoseconds). There is a range of scales (nanometres, nanoseconds) at which both descriptions are valid. At these scales, the fluid dynamics acquires non-negligible fluctuations while in the atomistic dynamics collective phenomena involving hundreds of atoms become important. Although this fact is not new, only recently has it become possible to investigate the processes at these overlapping scales, either computationally (large-scale high-performance molecular dynamics (MD), fluctuating hydrodynamics) or experimentally (the whole field of ‘nano-’ sciences).

However, there is still a conceptual gap in the theoretical description of these phenomena, as the continuous and atomistic theories do not interconvert into each other smoothly. Statistical mechanics converges to macroscopic laws only in the ‘macroscopic’ limit and a rigorous mathematical transition from Hamiltonian particle dynamics to hydrodynamics that is governed by the Navier–Stokes (NS) equations [1] still does not exist [2]. The same is true for going down the scales: the fluctuations of fluid dynamics do not become atoms in the limit.

Resolving microscopic details within a macroscopic flow is important for many practical applications. One such application is the rapid change of chemical environment around the protein molecule that forces it to fold. In one practical example, a micromixer was used to diffuse the protein to and from a hydrodynamically focused free-stream flow in a microchamber with conditions that very quickly change the concentration of a denaturant [3]. A typical mixer would have triple flow (three inlet channels, one outlet). A similar idea was applied [4] with a microchannel composed of two inlets and one outlet, where membrane fabrication was achieved by controlled shear stress and an interfacial cross-linking reaction between two streams of fluid. In another two-branch inlet and single outlet system, Ivorra *et al.* [5] performed a geometry optimization study to control the mixing time of the denaturing agent to control peptide folding. Most of the simulations reported in the literature for this type of problem are limited to rather simple modelling based on the Stokes equation with account taken of a passive concentration field. This does not include microscopic fluctuations despite the fact that for these processes (for example, protein folding) the molecular details are critical [6]. MD simulations can be used here. They are very effective for modelling nano-fluids and complex flows [7] but remain far too expensive even with the rapid development of supercomputers. For example, the characteristic time step in MD is about a femtosecond, which is nine orders of magnitude smaller in comparison with the characteristic mixing times of the best microfluidic mixer designs in the literature [5,8]. It is not only impossible to model a realistic device with pure MD simulation, but it is also useless because atomic resolution is required in very limited regions of space and it is infeasible to handle such a large amount of information.

Examples of compressing this information via a scale transformation from atomistic to hydrodynamic are abundant and go under a common theme of ‘coarse-graining’ in the MD literature. There are several approaches that aim at gradually lumping individual atoms together into larger quasi-particles that faithfully reproduce the behaviour of the original system [9]. Another coarse-graining method is based on representing the molecule as a string of beads and is mostly used in modelling polymers [10].

A different motivation lies behind the direction of research where the effect of continuum flow dynamics is investigated. For example, this is the case when the collective motion of molecules in pure liquid or biomolecular solutions is of interest, such as the dynamics of proteins in water. For the former, the glassy dynamics is the subject of intensive research, where large domains including hundreds of molecules are found to move as a whole with very different dynamics of molecules between the domains [11]. Similarly, the coherent behaviour of large conglomerates of water atoms, which form microscopic ‘water vortices’ on the surface of proteins, was reported by Umezawa *et al.* [12]. This is consistent with the hypothesis [13] that the entire dynamics of proteins is ‘encoded’ in the collective dynamics of water. Recently, the latter hypothesis has been indirectly confirmed by Nerukh & Karabasov [14], who demonstrated that the fluctuations of the water continuum around a peptide molecule are crucial for the dynamics of the peptide itself. One possible way to address the problem of collective phenomena in complex atomistic systems is to use a hybrid MD–continuum hydrodynamics (CH) modelling framework.

Attempts to use hybrid MD-CH methods in the literature have been known since the 1990s, where one of the most popular approaches is based on geometrical domain decomposition. The latter approach is based on using a fast continuum fluid dynamics solver in a large part of the solution domain and limiting the expensive MD model to a small region. The continuum–atomistic domains are connected through some boundary condition, and depending on this connection interface state variables schemes and flux coupling schemes are introduced.

In the state variables schemes [15–17], the coupling between the MD and CH regions is established through a finite-size overlapping zone while ensuring the conservation of bulk mass and momentum fluxes. For coupling of steady states, efficient algorithms based on variants of the alternating Schwarz approach [18–20] have been developed. Having a finite overlap region is important in these methods to allow the MD and CH solutions to relax before they are coupled together, avoiding sharp density oscillations. In the flux coupling schemes, the idea is to bridge the atomistic and hydrodynamic parts of the solution through a flux interaction exchanged through a control interface. As a continuum model for unsteady MD-CH coupling based on the flux-based approach, several hybrid schemes in the literature use the fluctuating hydrodynamics equations of Landau and Lifshitz [21–24]. Landau–Lifshitz fluctuating hydrodynamics (LL-FH) is a statistical model of thermal fluid fluctuations at the microscale. Prior to coupling with fully atomistic calculations, the LL-FH model has already been very popular for modelling colloids and suspensions in the framework of the lattice-Boltzmann method [25] and dissipative particle dynamics [26].

In the following, we present our results on (i) developing the methodology that allows natural change of space and time scales for fluid dynamics, correctly taking into account microscopic fluctuations, (ii) demonstrating the decisive role of water density fluctuations in the dynamics of a biomolecule at the atomistic level and (iii) conceptually combining (i) and (ii) by developing a novel framework for hybrid fluid dynamics–MD modelling of liquid systems, which has the properties of pure fluid dynamics and MD in the limits and at the same time correctly reproduces the fluctuations in the overlapping region of the scales.

## 2. Microscopic fluctuations in liquids

### (a) Fluctuating hydrodynamics model for macroscopic flow modelling

The LL-FH equations are a generalization of the deterministic NS equations to microscopic flow modelling by accounting for fluctuating stochastic sources in the momentum and energy equations, as described below:
2.1
Here is the specific energy density, *C*_{V}=*R*/(*γ*−1) is the heat capacity at constant volume, *R*=*k*_{B}/*m* is the gas constant, and terms on the right-hand side consist of an averaged part and a stochastic fluctuation part denoted with tilde.

The average term on the right-hand side of equations (2.1) is the standard stress tensor
2.2
and heat flow
2.3
where *η* and *ζ* are the shear and bulk viscosity, *D* is the spatial dimension, *κ* is the coefficient of thermal conductivity and *T* is the temperature.

As for the deterministic NS equations, the equations for time-averaged quantities are closed using an equation of state, which can be obtained from MD for instance [23].

The fluctuating part of the LL-FH equations corresponds to a continuum description of the thermal velocity fluctuations, the statistics of which is governed by the fluctuation–dissipation theorem [21]. The stochastic stress tensor is described as a random Gaussian matrix with zero mean and covariance, given by the formula
2.4
Using this correlation, the stochastic stress tensor can be expressed explicitly as [27]
2.5
where **G**_{αβ} is a random Gaussian matrix with zero mean and covariance 〈**G**_{αβ}**G**_{γδ}〉=*δ*_{αγ}*δ*_{βδ}, is a random symmetric matrix with zero trace, **E**_{αβ} is the identity matrix and tr[**G**] is the trace of a matrix **G**_{αβ}. The stochastic heat flow also has a zero average, while the covariance component is determined by the expression
2.6
or explicitly,
2.7
where **G**_{α} is a random Gaussian vector with zero mean and covariance 〈**G**_{α}**G**_{β}〉=*δ*_{αβ}. Furthermore, the stochastic stress tensor and stochastic heat flows are not correlated.

The LL-FH model allows accurate modelling of statistical properties of the atomistic system, such as the standard deviation of mass, momentum and kinetic energy, for control volumes as small as 5–10 atoms per cell size [24]. In the limit of large control volumes, the LL-FH model naturally transitions to the conventional NS equations.

The LL-FH equations can be efficiently solved with Eulerian methods, which have reached a mature state in computational fluid dynamics (CFD), in the framework of the time asynchronous relative dimension in space (TARDIS) scheme [28] for example. The idea of TARDIS is the extension of modern multi-time-step multi-grid-size algorithms developed in the CFD community [29,30] to multiscale modelling in fluid dynamics when the size of the system is itself a parameter that influences the governing equations as for the case of the LL-FH equations. TARDIS is based on concurrent transformation of the governing continuum fluid dynamics equations in space–time and then solving the transformed equations on a uniform Cartesian grid with the corresponding causality conditions at the grid interfaces. For example, in one dimension, let us consider the coordinate *x*=*x*_{0} that corresponds to the interface that separates large physical scales in space and time *l*_{s}(*x*)=*L*_{s} and *t*_{s}(*x*)=*T*_{s} for *x*≥*x*_{0}+Δ*x*/2 from the small ones *l*_{s}(*x*)=*l*_{s} and *t*_{s}(*x*)=*t*_{s} for *x*<*x*_{0}−Δ*x*/2. Then the following space–time transformation rule is considered:
2.8
and the TARDIS algorithm solves the governing equations in the transformed space–time variables and . By using the scale-adjustable zoom-in function (*l*_{s}(*x*)/*l*_{s})^{−1}=(*t*_{s}(*x*)/*t*_{s})^{−1}≡*f*(*x*), TARDIS can be explicitly tailored to the desired space–time scale mapping of the underlying theoretical model.

As an illustration of TARDIS's application to resolving nanoscale thermal velocity fluctuations in the region of interest of microscale-size channel flow on a grid of 200×200 control volumes, figure 1 shows the velocity profile through two cross sections along the height of the channel. These plots show that the fluctuations are most evident in the centre of the domain, corresponding to the small scales (CS2). In the large-scale part of the domain (CS1), the fluctuations become almost indistinguishable from the background Poiseuille flow profile. Figure 2 shows the results of TARDIS application for modelling of thermal density fluctuations in water at equilibrium conditions when the spatial length scale of the control volume ranges from eight atoms in the centre of the domain to 2190 atoms near the boundaries. For this calculation, the TARDIS grid is 256×10 control volumes. The problem considered corresponds to the grand canonical ensemble of volume *V* _{c} (size of the control volume) and temperature *T*, for which the following linear relationship between the standard deviation of density fluctuations and the inverse of isothermal sound speed *c*_{T} (1510 m s^{−1} for water at normal conditions) is valid [23]:
2.9
For the entire range of scales considered, the standard deviation of density that is computed from the LL-FH solution satisfies the theoretical relation within a few per cent error, which is within the statistical noise level (figure 3).

### (b) The role of density fluctuations in the dynamics of biomolecules: molecular dynamics results

We have simulated a system consisting of a dialanine zwitterion in a bath of explicit water. Classical MD has been used with the SPC (simple point charge) water model and the united-atom force field GROMOS for the peptide. The GROMACS software package was used for integrating the molecular equations of motion. The shape of the peptide is defined by the two dihedral angles (figure 4*a*). By calculating the ensemble-averaged density of water around the peptide (figure 4*b*), we aimed at answering the question: to what extent is the dynamics of water density connected with the dynamics of the peptide (quantified by its dihedral angles *ϕ* and *ψ*)?

To answer the question, we first calculated the fluctuations of the water density field *ρ*(**x**,*t*) around its average value . We then analysed the dynamical correlations between the dihedral angles *ϕ*, *ψ* and the water density using the linear stochastic estimation (LSE) model, which was originally developed for visualizing coherent structures in turbulent flows [31] and the identification of noise sources in turbulent jets [32]. In this approach, an approximation to *ρ*′ is found that minimizes the error [33], where are the angles averaged over all the frames at time *t*, *α*,*β* are stationary fields, and quantifies the density fluctuations correlated with the angles. The system of linear equations is solved to calculate the constants *α*,*β*:
2.10
and
2.11
Figure 4*c* shows the original density fluctuations field *ρ*′(**x**,*t*) on the right together with its part correlated with the peptide conformations on the left. The correlations are very strong in the first hydration shell of the biomolecule and remain non-negligible at a significant distance from it. Moreover, as we have shown [14], these correlations are qualitatively different in various regimes of the peptide's dynamics: they are strong at the moments of conformational transitions (when the shape of the peptide changes significantly) and weak at the times of metastable fluctuations of the peptide around the same average shape. Thus, the microscopic fluctuations not only are non-negligible, but also constitute the essence of the molecular phenomenon.

Summarizing, these above subsections clearly demonstrate that it is possible to accurately calculate the fluctuations of the liquid at the scales where the hydrodynamic and atomistic descriptions overlap and these fluctuations are important for atomistic motions. This leads to the question of how both continuum and atomistic descriptions can be implemented within the same multiscale framework, as discussed in §4.

## 3. Coupling of the Landau–Lifshitz fluctuating hydrodynamics and molecular dynamics equations based on a hydrodynamic analogy

### (a) Theory

Consistent coupling of the LL-FH model with MD equations remains a problem since the existing LL-FH–MD approaches [22] share the drawbacks of many other hybrid approaches, where a combination of artificial repulsive wall boundary condition based on a constrained particle dynamics concept and/or a particle reservoir with periodic boundaries is used. The latter prevents particles from drifting away, but may introduce spurious correlations and limits the geometry of the MD region. One of the questions that emerges for coupling LL-FH and MD is the choice of a hybrid modelling strategy that allows a smooth reduction of the number of degrees of freedom from fine scales (MD model) to large scales (LL-FH model). Essentially, this is the problem of developing a suitable low-order model for reducing the complexity of the MD model to the LL-FH framework based on Gaussian statistics in the macroscopic limit. On the other hand, mathematical modelling approaches based on low-order models have been used in continuum fluid dynamics for several decades. A classic example is the Lighthill acoustic analogy, which was introduced to bridge the scale difference that spans three to four orders of magnitude between the sound waves in the range of audible frequencies and the turbulent flow structures that generate sound. Since the pioneering work of Lighthill [34] various hybrid methods of this kind have been introduced in the literature. They all share a general idea to exactly rearrange the governing NS equations into the form of non-homogeneous linear equations and a nonlinear source at the right-hand side. For most advanced approaches of this type [35–37], the nonlinear source is directly related to the properties of fine-scale solution (space–time scales of the turbulence). In this work, the model based on a hydrodynamic analogy with two-phase flows is proposed for coupling the LL-FH and MD equations.

Let us consider a nominally two-phase liquid. The two phases are a Lagrangian and an Eulerian representation of the same chemical substance, which correspond to the atomistic (MD) and the continuum (LL-FH) model, respectively (figure 5). The phases are immersed into each other as ‘fine grains’, the surface tension effects are irrelevant, and both phases simultaneously occupy the same control volume. The partial concentrations of the MD ‘phase’ and the LL-FH ‘phase’ are equal to *s* and 1−*s*, respectively, where *s* is a parameter of the model, 0≤*s*≤1. In general, *s* is a user-defined scale function of space and time, which controls how much atomistic information is required in a particular region of the simulation domain.

For simplicity, let us consider the case of liquids for which the adiabatic heat ratio approaches unity, and the first two equations for mass and momentum in the LL-FH equation (2.1) can be solved separately from the energy equation [23]. Following the standard approach in two-phase modelling [38], the conservation laws for the mass of each phase are given by
3.1
and
3.2
Here *ρ*_{p}=*m*_{p}/*V* is the effective density of MD particle *p* that occupies elementary control volume *V* , corresponding to the grand canonical ensemble for the LL-FH model, *ρ* is LL-FH phase density, *m*_{p} is particle mass, **u**_{p} is the MD velocity, is the average velocity of the mixture , where *u*_{i} is the velocity of the continuum LL-FH phase, *N* is the number of particles per control volume *V* and *J*^{(ρ)} is the mass source/sink term which describes the phase transition. The latter describes the exchange of mass between the MD and LL-FH models, which needs to be defined. It is easy to see that (3.1) and (3.2) automatically satisfy the conservation law of the total mass of the liquid, .

In a similar manner, the momentum equations of the two phases are introduced as
3.3
and
3.4
where *F*_{ip} is the MD force exerted on particle *p* due to the pair potential interactions, for which an explicit expression is not required in the model as will be discussed later, and *J*^{(u)} is the LL-FH/MD momentum exchange term. Similar to *J*^{(ρ)}, *J*^{(u)} is a hybrid model quantity to be defined. By combining (3.3) and (3.4), it can be seen that Newton's second law, which equates the change of the total momentum , where to the force applied, is satisfied if for each control volume.

To close the model, the source/sink terms in the mass and momentum equations *J*^{(ρ)} and *J*^{(u)} must be defined. These terms are specified so that the conserved mass and momentum per each control volume *V* , and are driven towards the macroscopic averages obtained from the MD simulation in the hybrid MD/LL-FH zone, 0<*s*<1:
3.5
Here , and *L*^{(ρ)}(⋅) and *L*^{(u)}(⋅) are some operators, which play the role of LL-FH/MD boundary conditions to drive the corresponding deviations to the prescribed values within the hybrid MD/LL-FH zone 0<*s*<1 and return zero value for *s*=0 and *s*=1. For *s*=1 and no MD particles, it is easy to see that the solution of (3.3) converges to the continuum LL-FH model.

For equilibrium, when there are no external macroscopic forces, the continuum force term is equal to zero in accordance with the fluctuation–dissipation theorem. For non-equilibrium, this continuum force term is equal to the external force.

For the current implementation, several models for *L*^{(ρ)}(⋅) and *L*^{(u)}(⋅) have been tested, which correspond to driving the difference between the continuum average state and the MD average to zero using either (soft) diffusive or (hard) exponential forcing, (3.6) and (3.6*a*), respectively:
3.6
3.6a
where *α*,*β*>0 are numerical adjustable parameters.

The soft diffusive forcing (3.6) leads to the best results in terms of the hybrid algorithm's robustness and accuracy for a good of range of parameters 50<*α*,*β*<2000 (in reduced units specified in table 1) and 0<*s*<1.

From (3.5) and (3.6), the MD/LL-FH phase source terms *J*^{(ρ)} and *J*^{(u)} in the continuum FF-LH equations (3.1) and (3.3) are obtained as
3.7
To complete the hybrid model description (3.1)–(3.4), one needs to modify the classical MD equations d**x**_{p}/d*t*=**u**_{p},d**u**_{p}/d*t*=**F**_{p}, where **x**_{p}, **u**_{p} and **F**_{p} are the coordinate, velocity and pair potential force applied to the particle, so that there is an exact correspondence between (3.7), which satisfies the governing equations for the LL-FH phase (3.1), equation (3.3) and the sink sources in the macroscopic conservation laws for MD particles, (3.2) and (3.4).

This is achieved by the following modification of the MD equations:
3.8
The above relations for the modified MD particle velocity and acceleration, d**x**_{p}/d*t* and d*u*_{ip}/d*t*, have been derived starting from the conservation of mass and the corresponding equation for velocity and acceleration of the MD particles per control volume
3.9
and then exactly rearranging these equations to the form of macroscopic conservation laws (3.2) and (3.4) where the sink terms are given by (3.8).

Notably, in comparison with constrained particle dynamics models used in the literature [17], the current model (i) preserves a strict conservation of mass and momentum and (ii) in the continuum limit *s*→1, it correctly captures the mean and the variance of thermal fluctuations due to its convergence to the LL-FH equations.

For numerical implementation, the system of equations (3.1), (3.2) and (3.4) is solved with an explicit time-marching scheme where the MD and LL-FH solution parts are updated concurrently.

### (b) Numerical implementation results

For the current proof of concept study, liquid argon has been simulated for equilibrium conditions at temperature *T*=2.4 and mean density *ρ*=0.607. The reduced units are used and the main numerical parameters are summarized in table 1.

To initialize the solution, 50 000 MD iterations are first performed with a Nosé–Hoover thermostat [39,40] to reach the prescribed macroscopic parameters, after which the thermostat is turned off and the hybrid model simulations are conducted. There are about 2 000 000 MD time steps performed in total and 1 000 000 MD time steps are typically required to reach statistical convergence within 5%. The exchange between the continuum phase (FH) and the MD phase solutions is performed at every 10th MD step (MD step is equal to 1/10th of FH step) and statistics are collected every 50 time steps as in [24] for faster convergence. To investigate the hybrid model accuracy for multi-time stepping, the same model has also been tested by reducing the FH step so that one MD steps correspond to one FH step. The accuracy of the hybrid model in predicting the statistics of fluctuations remains fairly insensitive to the change of the FH step.

The simulations of the current hybrid LL-FH/MD model have been performed for several values of the model parameters, 50<*α*,*β*<2000 and 0<*s*<1. In all cases, a two-dimensional periodic box solution domain is considered. Figure 6 shows a typical behaviour of the standard deviation (std) of conservative velocity variable for *i*=1 that corresponds to the *x*-velocity component as a function of the number of MD iterations. Two representative values of the hybrid model parameter are considered. One corresponds to the case when the hybrid model is mostly driven by macroscopic dynamics (*s*=0.8). The other corresponds to the case when the hybrid model is mostly driven by MD dynamics (*s*=0.1). For reference, the pure atomistic MD solution and the pure continuum LL-FH solution of the same problem are shown on the same graph.

It can be seen that the hybrid model solution converges to the reference value of the pure MD solution and that of the pure LL-FH solution within the statistical uncertainty for both cases. For the solutions shown, *α*,*β*=50, but it has also been checked that a similar convergence to the correct limit can be obtained for *α*,*β*=2000.

To investigate the hybrid model's capability to deal with the situation when the hybrid system is not in equilibrium initially, the following test case has been considered.

The previous system of liquid argon is simulated when the continuum part of the solution is initially given a certain mean flow velocity drift while the MD part of the solution is initialized with no drift. Figure 7 shows the behaviour of the two parts of the solution for *s*=0.5. Having started from two different initial average velocities, the two parts of the hybrid model converge in the middle. This behaviour is completely consistent with the preservation of the mean flow rate that is expected in the model.

## 4. Conclusion

A new self-consistent hybrid MD-CH modelling framework based on an analogy with two-phase hydrodynamics is developed. For this framework, the MD equations are constrained so that their macroscopic behaviour satisfies the boundary conditions in the hybrid LL-FH/MD zone with a strict conservation of macroscopic mass and momentum throughout the solution domain.

The course-graining parameter, *s*, is introduced into the governing equations by an analogy with the two-phase hydrodynamics to permit a very natural treatment of the scale differences in the MD and the LL-FH representations. The LL-FH model for macroscopic flow description allows the statistics of thermal velocity fluctuations to be captured for a wide range of control volumes from 8 to 2190 atoms per cell length. Its implementation can be as efficient as that of the conventional deterministic CFD methods by using non-uniform Eulerian grids and asynchronous time stepping. Microscopic fluctuations that can be captured by the LL-FH model not only are non-negligible, but also constitute the essence of the molecular phenomenon, as demonstrated for modelling of a dialanine zwitterion in a bath of explicit water.

The new hybrid LL-FH/MD method is implemented for the simulation of liquid argon at equilibrium conditions. The method is shown to produce a converged solution that correctly captures macroscopic parameters such as standard deviation of thermal fluctuations and the mean flow velocity. This is achieved without introducing any ad-hoc treatment of MD particles such as ‘rescaling’ of the MD solution.

At present, the MD part of the current hybrid LL-FH/MD method is limited to a box domain where the average number of MD particles is maintained the same in all control volumes including the part of the computational domain that is dominated by continuum dynamics. Future work will be devoted to extending the current implementation to multi-time-scale modelling of the MD particles. One possibility for multiscale treatment of MD particles can include their gradual coarse-graining in the continuum part of the domain (*s*→1) so that the individual atoms there gradually become ‘hard spheres’, their response time to the pair potential force increases accordingly, and their dynamics is governed by the continuum LL-FH field.

Once fully implemented, it is expected that the new hybrid LL-FH/MD approach will become an efficient tool for multiscale simulations of biomolecules in liquids such as the simulation of hydrodynamic control of protein conformations in microfluidic devices.

## Funding statement

This work has been supported in the framework of the G8 Research Councils Initiative on Multilateral Research Funding (Engineering and Physical Sciences Research Council grant no. EP/J004308/1) and Russian Foundation for Basic Research (11–07–93938-G8_a).

## Acknowledgements

One of the authors is grateful to the Royal Society of London for their continuing support.

## Footnotes

One contribution of 13 to a Theme Issue ‘Multi-scale systems in fluids and soft matter: approaches, numerics and applications’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.