The microscopic quantities a ϵ ( x , y = x ϵ ) (static) or a ϵ ( x , y = x ϵ , t ) (dynamic) may be diffusivities (this example), permeabilities, elasticity moduli, reaction rate constants, etc. They converge to the effective properties when y ≫ 1 . It is worth noticing that homogenization does not aim at averaging the microscopic quantity a ϵ ( x ) , but the underlying conservation laws instead. By considering transport phenomena and due to the intertwining with other microscopic phenomena such as ad/absorption, a ϵ ( x ) should only be envisioned as the property relating the microscale flux j ϵ ( x , t ) = a ϵ ( x ) ∇ u ϵ ( x , t ) with the microscale driving force ∇ u ϵ ( x , t ) .

The corresponding microscopic advection-diffusion problem reads: ∂ u ϵ ( x , y , t ) ∂ t + ∇ ⋅ [ b ( x ϵ ) u ϵ ( x , y , t ) ] = ∇ ⋅ [ a ϵ ( x ϵ ) ∇ u ϵ ( x , y , t ) ] − f ( x ) . (3)

By assuming that a and b vary smoothly over a periodic cell (of volume V ) larger than ϵ c , Eq. 3 homogenizes as: ∂ U ( x ) ∂ t + ∇ ⋅ [ b ¯ ( x ) U ( x ) ] = ∇ ⋅ [ A ( x ) ∇ U ( x ) ] − f ( x ) , (4) where b ¯ is the average of b over V and A is the homogenized operator (stiffness matrix), defined as ∂ ∂ x i [ a i j ∂ ∂ x j ] = ∇ ⋅ [ A ∇ ] , and capturing the effective transport properties along each coordinate.

The simplest case arises from the homogenization of Eq. 3 in one-dimension ( n = 1 ) at steady-state and in the absence of advection, which reads: d d x [ a ϵ ( x ϵ ) d d x u ϵ ] = f   ∈ L 2 ( 0,1 ) with   u ϵ ∈ H 0 1 ( 0,1 ) . (5)

By setting F ( x ) = ∫ 0 x f d x and C ϵ a constant, the double integration of Eq. 5 leads to: ∫ 0 1 d u ϵ d x d x = ∫ 0 1 F ( x ) − C a ϵ ( x ϵ ) d x = 0. (6)

By invoking the separation of scales for ϵ → 0 and that the problem is periodic along y , one gets C = ∫ 0 1 F ( x ) d x and d U d x = lim ϵ → 0 d u ϵ d x . The macroscale solution of the Dirichlet periodic problem reads hence: d d x [ 〈 1 a ϵ 〉 − 1 d d x U ] = f   with   U ∈ H 0 1 ( 0,1 ) . (7)

Thus the homogenized transport coefficient A = 〈 1 a ϵ 〉 − 1 , where 〈 X 〉 represents the volume average of X . The following result can be generalized for arbitrary boundary conditions and dimensions. The homogenized matrix corresponding to A ϵ = [ a 11 ϵ a 12 ϵ a 21 ϵ a 22 ϵ ] (stratified anisotropic 2D media) is [see p 16 of Ref. (Zhikov et al., 1994)]: A = [ 〈 1 a 11 ϵ 〉 − 1 〈 a 12 ϵ a 11 ϵ 〉 〈 1 a 11 ϵ 〉 − 1 〈 a 21 ϵ a 11 ϵ 〉 〈 1 a 11 ϵ 〉 − 1 〈 a 22 ϵ 〉 + 〈 1 a 11 ϵ 〉 − 1 〈 a 21 ϵ a 11 ϵ 〉 〈 a 12 ϵ a 11 ϵ 〉 − 〈 a 12 ϵ a 21 ϵ a 11 ϵ 〉 ] . (8)

The presented results generalize standard approximations for effective media, such as the Maxwell Garnett and the Maxwell-Rayleigh approximations for spherical and non-spherical inclusions. For symmetric positive definite matrix, the following inequalities (known as Voight-Reiss’ inequalities) hold: 〈 A ϵ − 1 〉 − 1 ≤ A ≤ 〈 A ϵ 〉 . (9)

The first equality A = 〈 ( A ϵ ) − 1 〉 − 1 holds only if ∇ × ( A ϵ − 1 ) = 0 , the second A = 〈 A ϵ 〉 only if ∇ ⋅ A ϵ = 0 .

Though the studied transport equation does not include any thermodynamic consideration, particular cases of two-phase systems can be treated with the previous procedure. Soft and hardcore inclusions are coded as A ϵ ( x ) = 0 and A ϵ − 1 ( x ) = 0 , respectively, with nonzero values outside inclusions: ν 1 I ≤ A ϵ ( x ) ≤ ν 2 I with ν 1 , ν 2 > 0 .

Explicit solutions exist only in rare cases such as periodic stratified media, but they cease to operate in large dimensions. Alternatively, the microscale solution u ϵ ( x , y ) can be envisioned as a first-order approximation of the solution of the homogenized problem U ( x ) : u ϵ ( x , y ) ≈ U ( x ) + ϵ δ u ( x , y ) . (10)

Equation 10 accepts a statistical interpretation: u ϵ ( x , y ) ≈ U ( x ) + ( u ϵ ( x , y ) − 〈 u ϵ ( x , y ) 〉 o v e r   y ) , so that ∂ u ϵ ( x , y ) ∂ y = o ( ∂ U ( x ) ∂ x ) . The exact calculation of expansions involves total derivatives. The procedure is accelerated by highlighting the composition properties of the microscale differential operator: Δ ϵ = ∇ ⋅ [ A ϵ ∇ ] = ( ∂ ∂ x i + ϵ − 1 ∂ ∂ y i ) a i j ϵ ( ∂ ∂ x i + ϵ − 1 ∂ ∂ y i ) = ϵ − 2 Δ 1 + ϵ − 1 Δ 2 + Δ 3 , (11) where: Δ 1 = ∂ ∂ y i ( a i j ϵ ∂ ∂ y i ) ; Δ 2 = ∂ ∂ y i ( a i j ϵ ∂ ∂ x j ) + ∂ ∂ x i ( a i j ϵ ∂ ∂ y i ) ; Δ 3 = a i j ϵ ∂ 2 ∂ x j ∂ x j . (12)

The expressions Eq. 11 and Eq. 12 applied to approximation Eq. 10 (by assuming that the macroscale solution is twice differentiable up to the boundary) lead to: Δ ϵ [ U + ϵ δ u ] = ϵ − 1 ( Δ 1 δ u + Δ 2 U ) + ϵ 0 ( Δ 3 U + Δ 2 δ u ) Δ 2 + ϵ Δ 3 δ u . (13)

The separation of scales ( ϵ → 0 ) enforces that Δ 1 δ u = − Δ 2 U , which can be resolved as a periodic boundary value problem to the independent variable y . By denoting its general solution { ϕ k } k = 1... m , one gets δ u ( x , y ) = ϕ k ( y ) ∂ U ∂ x k . This solution introduced in Eq. 13 and averaged over all possible solutions leads to the following homogenized coefficients: a i j = 〈 a i j + a i k ∂ ϕ j ∂ y k 〉 . A second-order approximation can be derived for periodic macroscale solutions U . The results are very general, but it requires that U is of class C 2 . This property fails to be verified if the homogenized medium includes barriers (e.g., cell walls), strong partitioning, or large jumps in transport properties.

In the presence of jumps in a ϵ ( x ) , the homogenizing problem is better solved using a weak formulation of the partial differential equations. The sufficient condition is that u ϵ fulfills Dirichlet or Neumann boundary conditions at jumps. For a sufficiently smooth and periodic a i j ϵ ( x , x ϵ ) , the solution of the microscale problem at steady-state minimizes the standard variational minimization problem: min u ∈ H 0 1 ( V ) [ 1 2 ∫ V ∑ i , j a i j ϵ ( x , x ϵ ) ∂ u ∂ x i ∂ u ∂ x j d V − ∫ V f ( x ) u ( x ) d V ] . (14)

The solution of Eq. 14 is inferred efficiently via a top-down approach sketched in Figure 2C. Standard finite-element techniques (triangular mesh) are used to solve the macroscale problem. The effective macroscale operator herein is approximated by solving the microscale problem only on a small subset of elements taken at quadrature points. This strategy is known as the heterogeneous multi-scale method introduced by W. E. and B. Engquist (Engquist and Huang, 2003) and discussed in Ref. (Weinan et al., 2005). The overall cost does not increase as ϵ decreases and generalizes naturally to time-dependent and nonlinear problems.

Capturing the microscale statistics adequately is essential as the whole approach relies on an approximation of the separation of scales ( ϵ ≪ 1 ) , for which A evolves only smoothly with the position. Due to the independence of microscopic and macroscopic scales, there is no need to apply microscale fluxes with intensities as large as those applied at the macroscale. Indeed, macroscale fluxes describe an evolution of the macroscopic system, whereas microscale ones describe its evolution at some equilibrium states. The concept of equilibrium differs here from strict thermodynamic equilibrium as it assumes an equilibrium of microscale fluxes/forces without enforcing stronger concepts, such as microscopic reversibility or time reversibility. The concept of equilibrium can be approached numerically by introducing a micro time scale and by refining the microscale solutions iteratively until reaching a steady state.

Diffusivities of gases and organic molecules in condensed phases (liquids or solids) can be derived from the second-order stationary displacements of the center-of-mass (CM) of one single diffusant (or averaged over the trajectory of several diffusants) by molecular dynamics simulation. The displacements controlling the translation of CM due to thermal noise (collision with atoms of surrounding molecules) between times τ and τ + t , r CM ( t , τ ) are defined from the positions vectors (3 × 1 vector) of CM, x C M , relatively to the center of mass of the entire system (diffusant + host), x O : r CM ( τ , t ) = x CM ( τ + t ) − x CM ( τ ) − ( x O ( τ + t ) − x O ( τ ) ) . (22)

When the displacements of CM are truly independent, the quadratic displacement increases linearly with the lag operator t with 〈 r C M ( τ , t ) 〉 a l l   t = 0 : d d t 〈 r CM ( τ , t ) r CM ( τ , t ) T 〉 all   τ = 2 D ( t ) , (23) where X T is the transpose of X   . The diffusion matrix D ( t ) is positive-defined, whose eigenvectors code for the main direction of translation on time scale t . It accepts a unique Cholesky decomposition: D ( t ) = B ( t ) Λ ( t ) 1 / 2 ( B ( t ) Λ ( t ) 1 / 2 ) T , (24) where Λ ( t ) is a diagonal matrix containing the eigenvalues and where the columns of B ( t ) are the associated eigenvectors. Displacements fluctuate on the time scale t at a rate d r C M ( τ , t ) d t equal to B ( t ) Λ ( t ) 1 2 X , where X is a random vector, each of whose components follows a reduced centered normal distribution.

Solutes larger than voids exhibit a preferential direction as shown for methoxybenzene (anisole) in Figure 3, so that the random trajectory looks like an extruded tube on large time scales (Figure 3A). The solute seems to jump from one sorption site to the next, separated by distances ranging from 3 to 6 nm (Figure 3B). During the different jump events, the solute progressively loses memory of its previous orientation (Figure 3D). The macroscopic diffusion coefficient eventually emerges as the trace limit of D ( t ) when t → ∞ (note that t = 1 f is the lag operator and not the clock time). By choosing the columns of B ( t ) as reference frame, one gets the following expression of the rate of increase of the mean-square displacement of CM: d d t msd CM ( t ) = tr   ( d d t 〈 r CM ( τ , t ) r CM ( τ , t ) T 〉 all   τ ) = 2   tr   ( D ( t ) ) = 2   tr   ( Λ ( t ) ) = 6 D ( t ) , (25) which leads to the practical formula of the macroscopic diffusion D ( f = 0 ) : D ( f = 0 ) = 1 6 lim t → ∞ d d t ms d CM ( t ) ≈ large   t ms d CM ( t ) 6 t . (26)

The exponent d  ln ( m s d C M ( t ) ) d l n ( t ) plays a central role in identifying the rate of convergence to the hydrodynamic regime, where all displacements are truly uncorrelated (exponent close to unity). Exponents lower than 1 reveal strongly correlated random walks and apparent diffusivities decreasing with simulation time. Such decreases are shown in Figures 3E–I and occur on a broad range of time scales. The hops responsible for the diffusant translation can be visualized with the spectral density of d r CM / d τ , denoted S C M ( f ) , which can be calculated from the Wiener-Khinchin real theorem applied to the random variable r CM ( τ , t ) = ∫ τ τ + t d r CM / d τ | ( θ ) d θ . The pair of transforms is known as MacDonald’s theorem [a demonstration is presented in chapter 18 of Ref. (Van Vliet, 2008)] and reads: { d d t ( ms d CM ( t ) ) = 2 ∫ 0 ∞ S C M ( f ) sin ( 2 π f τ ) d f f S C M ( f ) = t → ∞ 2 f ∫ 0 ∞ d dt ( ms d CM ( t ) ) sin ( 2 π f τ ) d τ . (27)

Since the diffusion of a solute larger than voids in a rubber polymer is strongly related to the relaxation of polymer segments, which could think that the convergence of D ( t ) should obey some linear damping (e.g., Hooke’s law) and consequently that diffusion coefficients could be easily extrapolated, when the associated time constant t 0 is known. This reasoning would lead S C M ( f ) to decrease in 1 f 2 and an exponential convergence to the hydrodynamic regime: { d d t ( ms d CM ( t ) ) = S CM ( f = 0 ) 1 − exp ( − t t 0 ) S CM ( f ) = S CM ( f = 0 ) 1 + 4 π 2 ( f t 0 ) 2 . (28)

Figure 3I shows conversely that S C M ( f ) evolves from a high-frequency white noise (tumbling motions of the diffusant) towards a low-frequency white noise (macroscopic diffusion) by crossing an intermediate regime linear in 1 f (pink noise), separated by two frequencies denoted f 1 and f 2 . This noise corresponds to a random damping process with t 0 distributed between t 2 = 1 / f 2 and t 1 = 1 / f 1 , such reasoning leads to: ∫ t 1 t 2 S CM ( f ) d t t 2 − t 1 = = t 2 > t 1 S CM ( 0 ) 2 π f ( t 2 − t 1 ) ( arctan ⁡ 2 ⁡ π f t 2 − arctan ⁡ 2 ⁡ π f t 1 ) = { S CM ( 0 ) f 2 1 t 1 t 2         for     f ≫ 1 / t 2 S CM ( 0 ) 4 f ( t 2 − t 1 )       for     1 / t 2 ≪ f ≪ 1 / t 1 S CM ( 0 )         for     f ≪ 1 / t 1 (29)

These results have been used to build a solid theory of diffusion in polymers coupling the free-volume theory and the mode-coupled theory (see Section 4.4). They justify the universal cage-escape mechanism to describe random walks in condensed phases.

The frequency f 1 ≤ 1 / t 0 ≤ f 2 and, more particularly, its lower bound control the value of the diffusion coefficient. Ref. (Vitrac and Hayert, 2007) studied the effect of the distribution of sorption energies G i on D systematically. The frequencies f 1 and f 2 are connected with the jumping rates via TST, f = k B T h exp ( − G i → j ‡ − G i k B T ) , with G i the local sorption energy at site i and G i → j ‡ the free energy barrier associated to the translation from i to j . The main results are illustrated in Figure 4.

The strong relationship between sorption free energy and diffusion can be easily understood on a periodic surface of free energy, comprising only two states separated by a fixed distance l (Figure 4A). Although it may seem surprising to simulate a random walk on a two-state line, the distance visited is a random variable because the duration of jumps/transitions is also a random variable. The diffusion coefficient depends on two characteristics, the height of the smallest energy barrier G Δ and the partition coefficient K c o n t r a s t between the two sites. Increasing K c o n t r a s t while keeping the same average sorption leads to a dramatic decrease in D values. This effect was verified on random surface energies in higher dimensions. Beyond a critical dispersion value, preferred sites behave as attractors creating negative correlations in the long term. Simulated assumptions are close to what happens when the size of the diffusant is increased (Neyertz and Brown, 2004). On non-random surfaces (i.e., with channels), the reduction of D values is amplified by additional tortuosity effects. Because the results can significantly deviate depending on whether the micro-reversibility hypothesis is verified, it is preferable to restrict the TST approach to the molecular and supramolecular scales. TST-type modeling lends itself very well to the integration of trajectories by Kinetic Monte Carlo (KMC) methods. Classical molecular dynamics simulation and KMC methods can be combined to reproduce the effect of rare events over the long term (Neyertz and Brown, 2010). The methodology allows to obtain statistically independent but without guaranteeing the final D value if limiting events have been insufficiently sampled.

The previous examples justify the concept of time-homogenization to derive macroscopic transport coefficients. Rare events control the entire dynamics and its activation by temperature. One can intuitively think that such phenomena do not exceed a few hundred nanoseconds and are irrelevant at the phase level of the food. This is erroneous because diffusion is a stochastic process in time and in space at all time and length scales. We used to model the diffusion in the hydrodynamic regime in discrete times as a self-similar skewed trajectory, where the direction of the translation and its span are random normal variates. Close to interfaces and, in particular, curved interfaces, the process ceases to be isotropic. By analogy with optics, the constant K c o n t r a s t and the brusque change of diffusion coefficients at interfaces act as refractive indexes and modify the distribution of molecules on both sides of the interfaces. In dense emulsions with globules close to the percolation threshold, diffusants scatter and reflect towards the medium with the highest residence time. The translation from the medium with the highest diffusion coefficient controls the effective diffusion coefficient. These effects are illustrated in Figure 5 in 1, 2, and 3 dimensions with the following conventions based on the results of Ref. (Vitrac and Hayert, 2020). K is the partition coefficient between the continuous and the dispersed phase. r D = D c D d is the normalized diffusion coefficient. For large K or low r D values, the molecules initially contained in the globule create a halo (i.e., close to a diffraction pattern created by an infinitesimal object) around the globule, with a little back diffusion to the reverse direction. At this stage, it is important to notice that K defines only the ratio of diffusant densities across the interface ∂ Ω ( x = 0 ) and not the probability to cross the interface. For the configuration depicted in Figure 5A with a diffusant located initially at x = x 0 , the ratio of probability to remain the globule after a dimensionless time F o = D d t x 0 2 is given by: pr ( d → d | ∂ Ω , Fo , x = x 0 ) pr ( d → c | ∂ Ω , Fo , x = x 0 ) = 1 K D d D c 1 + K D d D c erf ( 1 2 Fo ) 1 − erf ( D d D c 1 2 Fo ) exp [ − ( D d 2 D c 2 − 1 ) 1 F o ] → F o → ∞ 1 K D d D c . (30)

The time limit of Eq. 30 is counter-intuitive and reflects the ability of the diffuser to make countless round trips across the interface but at different speeds. The local thermodynamic equilibrium condition imposes the density ratio (formally K ) and continuity of the fluxes. Additionally, it is noteworthy that globule interfaces do not represent in this representation a transition state, as the density at the interface is not continuous. However, the global escape rate from each globule could be used as input parameters to simulate diffusion in dispersed media. Compared to Eq. 8, this technique accepts larger fluctuations in surface area, globule size and density, large drops in diffusion coefficients, and finally, partition coefficients very different from unity.

Cage-escape mechanisms can be found in various applications, including the percolation of oil in French-fry or potato chips-type products once the product is removed from the bath. Trapped air bubbles in the potato tissue reduce the accessible cross-section to oil and consequently slow down dramatically the capillary flow through crust defects. These effects are illustrated in Figure 6 and are based on simulations and experimental observations reported in (Vauvre et al., 2015) and in (Achir et al., 2010; Vauvre et al., 2014; Patsioura et al., 2015), respectively. We could think that the effective diffusion coefficient of oil should scale as: D o ( z , T ) = 〈 r 〉 γ ( T ) ⁡ cos ⁡ θ 8 η o ( T ) , (31) where here γ is the surface tension, θ the contact angle, η o the viscosity of the oil, 〈 r 〉 the average radius of cells assuming roughly a cylindrical shape.

The micro-model showed that the successive filling times of two cavities of radius r and length Δ z and connected by a defect of radius r 0 and length l 0 increase considerably with each pass. If the oil-air interface is displaced without creating a bubble (Figure 6A), the recurrence relationship between the filling times Δ t i , f i l l is given by: Δ t i + 1 , f i l l = Δ t i , f i l l − ln ( 1 − ς ) 2 ς D o ( T ) ( Δ z + l 0 ( r r 0 ) 4 ) Δ z , (32) with 0 ≤ ς < 1 a factor that controls the gas phase displacement mechanism [see Ref. (Vauvre et al., 2015)].

Due to the hexagonal shape of the cells and the angular details of the defects, the oil flows faster in the corners (Figure 6B) and fills the defects/restrictions without displacing the gas phase. The destabilization of the liquid films in the defects is responsible for the sudden phenomenon of “snap-off” (collapse of two initially disjointed films), which will lead to the creation and trapping of bubbles in the defects. In the presence of an air bubble trapped within the cavity, filling times diverge and Eq. 32 becomes: Δ t i + 1 , f i l l w i t h   b u b b l e = Δ t i + 1 , l a g − ln ( 1 − ς ) 2 ς D o ( ( 1 − υ ) Δ z + l 0 ( r r 0 ) 4 ) Δ z  with  υ > 0 = Δ t i , f i l l w i t h   b u b b l e − ln ( 1 − ς ) 2 ς D o ( υ Δ z + l 0 ( r ξ r 0 ) 4 ) Δ z − ln ( 1 − ς ) 2 ς D o ( ( 1 − υ ) Δ z + l 0 ( r r 0 ) 4 ) Δ z = Δ t i , f i l l w i t h   b u b b l e − ln ( 1 − ς ) 2 ς D o ( Δ z + l 0 ( r r 0 ) 4 ( 1 + 1 ξ 4 ) ) Δ z (33) with ξ the restriction coefficient of the passage related to the presence of the bubble and υ the fraction of the volume of the cell (i.e., formally the cumulative work of the volume forces acting on the bubble) preceding that must flow before the trapped bubble is released. The comparison of Eqs 32, 33 enables to approximate the lag time as Δ t i + 1 , l a g ≈ 1 ξ 4 Δ t i + 1 , f i l l for low values of ξ et υ . Experimental filling kinetics of connected cells in parenchyma cells confirm this evolution and the distribution of waiting times as shown in Ref. (Patsioura et al., 2015).

This example and the previous one illustrate the importance of time homogenization and the sampling of rare critical events. Indeed, snap-off and partitioning are elusive and equilibrium properties, respectively. Because its occurrence depends on history, its effects cannot be correctly described by a homogenization procedure, assuming a particular steady-state.

Molecular modeling uses forcefields to calculate the pair (distance), angular/planar, quadrupole, and octupole interactions between atoms, as sketched in Figure 8A. Most often, forcefields rely on classical Newtonian mechanics. Forcefields representing chemical bonds are either fitted to experimental data (type-I forcefields such as bead-spring model) or fitted to quantum calculations (type-II forcefield with additional corrections to mimic the vibration spectra of molecules). Non-covalent potentials are associated with van-der-Walls and Coulombic interactions defined by partial charges usually applied at the center of atoms. Some software packages, such as LAMMPS (Sandia National Laboratories) and Materials Studio (Dassault Systèmes BIOVIA), offer flexible formulations for soft and solid-state materials and coarse-grain systems. Others, such as Gromacs and Charmm, offers forcefields optimized/specialized for biological molecules, with an emphasis on proteins and lipids. Integrated environments (Avogadro, MedeA LAMMPS from Materials Design, Materials Studio) facilitate the many operations (setup, post-treatments) needed to build and interpret molecular simulations. Some generic visualization tools such as VMD and Ovito offers versatile visualizations for multiscale simulations. For extremely large datasets, the open source ParaView increases the rendering efficiency by focusing on the levels of details important for the observer.

This section introduces briefly the different principles of molecular mechanics. Additional details on forcefields, integration schemes and inversion procedures can be found in the reference textbook of Frenkel and Smit (Frenkel and Smit, 2002). These details are essential to understand the parameterization of open-source (Gromacs, LAMMPS, Charmm, CP2K) or proprietary/commercial software packages (NAMD, Tinker, Materials Studio). At this stage, it is worth noticing that classical oscillators do not behave as quantum ones, whose more probable state coincides with a minimum of potential energy. The likeliest states associated with classical oscillators in dynamics correspond to kinetic energy minima (speed minima) and, therefore, maxima potential energy (see Figure 8C). According to the fluctuation-dissipation theorem, the states corresponding to each maximum are unlikely, but their long-term average in molecular dynamics converges towards the likeliest state at equilibrium. This result will be obtained for microscopically reversible phenomena integrated in time with a time-reversible scheme. However, the convergence to realistic states may be very slow due to explicit numerical schemes, which are stable only for time steps much smaller than the period of the vibration of the lightest atom, typically 1 picosecond in the presence of hydrogens. When trajectories are not required but only configurations or distributions, Monte-Carlo sampling offer atoms better results. It can explore the surface energy in a non-physical manner and overcome any energy barrier separating two feasible states. When the sampling involves a sequence of trial moves (e.g., translation, rotation, volume change), each move must also be picked randomly to enable the reverse course and preserve the detailed balance. The macroscopic quantities are in both cases extracted from ensemble averages.

Due to the implicit time integration, acceleration can be gained only by reducing the number of freedoms coupled with longer time steps made possible by merging several atoms (or even molecules) into one heavier particle and therefore subjected to lower frequency vibrations. The iterative coarse-graining methodology is shown in Figure 9. The final goal is to preserve the convergence towards the final property by devising a proper mean force potential ϕ C G ( R ) from detailed pair energies ϕ i j ( r i j ) . Even if the potential is not defined uniquely, methodologies providing canonical averages will minimize errors on generated configurations, energies, and entropy comparatively to the full detailed ones. Three methodologies have been usually preferred for coarse-graining.

• Direct potential fitting (Cranford et al., 2010): ϕ C G = ∑ l i s t   a t o m s ω i j ϕ i j where ω i j are summation weights representing all contributions (intramolecular, intermolecular such as hydrophobic and other solvent effects, electrostatic interactions, hydrogen bonding between molecules, and/or entropic effects).

Though the same methodologies could be applied to adjust forcefields to experimental data, they are presented in a logic of bottom-up modeling and simulation built on molecular mechanics. The knowledge gained at a small scale is used to reconstruct the macroscopic behavior. At small scales, forcefields are inherently repulsive, and the packing of particles is very compact. Excluded volumes vanish at a higher level of coarse-graining, and some overlap must be accepted to preserve the bulk density. At this stage, it is noteworthy that described procedures apply only to conservative forces. Dissipative forces required for DPD simulations must be adjusted to dynamic properties such as viscosity or self-diffusion. Finally, and compared with atomic force fields, coarse-grained forcefields depend on the system’s state (solid, liquid, gas), possibly on temperature and local composition. They are, therefore, not recommended for the calculation of state diagrams without independent validation.

Studying the effects of temperature and pressure on phase transition is an easy task with the generalization of efficient thermostats and barostats. Enforcing equilibrium between phases in molecular simulation is comparatively much more difficult. Small systems do not behave as large ones. Water did not freeze at the right temperature as it was confined. Water and oil will not create a nice interface within a small cubic simulation box. Droplets are not spherical, and their shape fluctuates with the number and the arrangement of the molecules inside each droplet. Besides, surface tension in droplets smaller than Tolman’s length deviates from its macroscopic value. The finite size of tested systems complicates the construction of large droplets at high dilution. Setting a dimension longer than the two others forces the creation of one flat interface (or two). In this case, the transport of small solutes can be tracked between the two phases at reasonably low concentrations. However, the statistics are too poor to get reliable partition coefficients, and the strategy cannot be generalized to polymers.

Chemical equilibration is the chief difficulty in reaching thermodynamic equilibrium. It requires an exchange of molecules between the two phases. Liquid-vapor and liquid-liquid equilibria can be efficiently studied in the Gibbs ensemble (Panagiotopoulos, 2002) for moderately dense multicomponent fluids. Figure 10 illustrates the principles. The phase coexistence properties of multicomponent systems are studied from a single Monte-Carlo simulation involving two distinct physical regions with different densities and compositions. The regions (boxes) can be dissociated without any common interface. Three types of perturbations are performed on each box: 1) a random displacement of molecules that ensures equilibrium within each region, 2) an equal and opposite change in the volume of the two regions that results in equality of pressures, and 3) random transfers of molecules that equalize the chemical potentials of each component in the two regions.

The great advantage of the Gibbs method over the alternative techniques for studying phase coexistence is that, in the Gibbs method, densities and compositions of the coexisting phases are resolved simultaneously. The coexistence line between phases does not require the determination of chemical potentials as a function of temperature and pressure. The list of operations must be picked from a repertoire of trial moves in a random order to preserve the microscopic reversibility. However, the Gibbs ensemble method remains impractical for very cohesive phases such as polymers and very ordered systems (crystals, liquid-crystals). Free-energy calculation methods must be preferred (see Figure 11).

G. N. Lewis proposed the concept of the thermodynamic potential of substance i in the phase α , μ i , α ( T , P ) , while verifying that ∂ μ i , α ∂ P | T = V ¯ i with V ¯ i = R T P the molar volume of the gas i . After integration at a constant temperature, the change in chemical potential from pressure P r e f to P is R T ⁡ ln P P r e f for a pure ideal gas. Lewis generalized to phase mixtures (ideal or not) by replacing pressure by a function f i , called fugacity and assessing the capacity of a substance literally to flee: μ i , α − μ i , α r e f = R T ⁡ ln f i , α f i , α r e f . Reference values at the temperature T μ i , α r e f and f i , α r e f are physically related, but the choice of the reference chemical potential or the reference fugacity is arbitrary. Fugacities of substances absorbed on solids can be determined directly by molecular modeling using the grand-canonical ensemble (μVT). Its utility lies in its resemblance with experimental measurements (e.g., sorption microbalance): the system interchanges particles with a reservoir so that the chemical potential can be by changing the number of particles changes accordingly. As for the Gibbs ensemble, the method is prohibitive for condensed phases. Flory-Huggins’ theory offers a scale-invariant to estimate activity coefficients, γ i , α v , relative to volume fractions ϕ i , α in any condensed phase, as summarized in Table 5. The presented relationships expressing γ i , α v in function of ϕ i , α are equivalent to conventional sorption isotherms for binary and ternary mixtures. Indeed, activities and γ i , α v are related as: f i , α ( T ) f i , α r e f ( T ) = γ i , α v ( T ) ϕ i , α . (34)

By considering that partial molar volumes are similar in both phases, solute partition coefficients between the phases α and β are given by: K i , α β ( T ) = ϕ i , α / V ¯ i ϕ i , β / V ¯ i = γ i , β v ( T ) γ i , α v ( T ) . (35)

Defining the activity coefficients respectively to volume fractions instead of molar fractions as in the regular solution theory offers several advantages. Volume fractions enable to consider situations where the solute i is larger than surrounding molecules and also the converse situation where i is infinitely smaller than polymer segments. At this stage, it is important to note that the reference state of solute i is an amorphous state. If its reference state is solid at temperature T , a fugacity correction equals to the ratio of fugacities between solid and liquid states must be added [see Ref. (Nguyen et al., 2017b)]. Flory-Huggins sorption isotherms are equivalent to Henry ones at infinite dilution and describe well swelling effects (Oswin type isotherms) at high concentrations. Sorption is assumed to be reversible (no hysteresis). Generalization to sorption-desorption cycles with hysteresis when T g is crossed during sorption can be found in Ref. (Kadam et al., 2014). For liquids, the Kirkwood-Buff (KB) solution theory (Kirkwood and Buff, 1951) at a microscopic scale can provide exact macroscopic predictions even for interacting liquids and vice-versa, as shown by Arieh Ben-Naim (Ben-Naim, 1977). The pro and cons of KB and FH theories are discussed in (Nguyen et al., 2017b). Corrections derived from KB integrals and pair radial distributions are mainly required for compressible mixtures and large solutes, as shown in (Gillet et al., 2010). The term n k   a r o u n d   i c o m p r e s s i b l e captures the positional entropy reduction due to the specific arrangement between a large solute i and the solvent molecules.

The calculation procedure of excess chemical potentials is summarized in Figure 11, with technical details given in Refs. (Gillet et al., 2009; Gillet et al., 2010; Vitrac and Gillet, 2010; Nguyen et al., 2017a; Nguyen et al., 2017b; Zhu et al., 2019a). The Flory-Huggins isotherm comprises a positional entropy term and an enthalpic term proportional to the extent of { χ i , k } k = P , F . The first term is determined by the number of molecules k , which are displaced by the introduction of i assuming that the mixture is incompressible (first approximation) or compressible (from radial distribution). The second term is essentially enthalpic, but it also contains the non-positional entropic contribution capturing the fluctuations of the contacts between molecules. Both properties can be studied by molecular modeling (packaging of molecules, pair interactions energies). The advantage of the approach is the capacity to extrapolate the results to arbitrary binary, ternary and quaternary composition. Aqueous and hydroalcoholic solutions have been specifically discussed in Ref. (Gillet et al., 2010).

Flory-Huggins coefficients { χ i , k } k = P , F are defined as the dimensionless mixing energy (enthalpy) in excess relative to pure compounds; it is defined for polymers ( k = P ) and liquids ( k = F ): χ i , k ( n k , T ) = 〈 h i + k n k 〉 T + 〈 h k + i n k 〉 T − 〈 h k + k n k 〉 T − 〈 h i + i 〉 T 2 R T   for   k = P , F , (36) where 〈 X 〉 T represents an ensemble-average of X , n P is the length used in the approximation of the polymer and n F = 1 .