Diffusion in Homogeneous and in Inhomogeneous Media: A New

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Diffusion in homogeneous and in inhomogeneous media: a new unified approach Luis Fernando Mercier Franco, Marcelo Castier, and Ioannis George Economou J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00653 • Publication Date (Web): 27 Sep 2016 Downloaded from http://pubs.acs.org on October 1, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Diffusion in homogeneous and in inhomogeneous media: a new unified approach Lu´ıs Fernando Mercier Franco, Marcelo Castier, and Ioannis G. Economou∗ Chemical Engineering Program, Texas A&M University at Qatar, PO Box 23874, Doha, Qatar E-mail: [email protected]

Abstract We propose a new method to calculate the diffusion coefficient within molecular dynamics simulations for either homogeneous or inhomogeneous fluids. We formulate such method by solving analytically the Smoluchowski equation for a linear potential of mean force within a thin layer with absorbing boundary conditions. The bulk, or homogeneous, fluid diffusion emerges as a particular case in this approach. We apply this method to bulk liquid water at atmospheric pressure and different temperatures using the SPC/E water force field. We show that our method gives results as accurate as the traditional Einstein-Smoluchowski method, avoiding the fitting procedure required in the traditional method. We also apply this method for molten sodium chloride showing its applicability for multi-component systems. The water vapor-liquid interface is studied as an example of an inhomogeneous system. We calculate all the components of the diffusion tensor at the interface. We observe the same anisotropy between the perpendicular and the parallel components at the interface as it has been noted in the literature. We also calculate the perpendicular self-diffusion coefficient of methane ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

near the calcite surface showing that this coefficient is much lower than the parallel diffusion coefficients. We believe that this new unified approach is a very promising technique for both bulk and confined media.

1

Introduction

Diffusion is a central phenomenon for many natural and industrial processes. Despite its crucial importance, modeling of diffusion remains an enormous challenge, and, even nowadays, phenomenological equations are still used widely. The first attempt to describe mass diffusion mathematically was formulated by the German physician Adolf E. Fick. The phenomenological equation proposed by Fick follows the same formalism as Fourier’s model for heat conduction, Newton’s equation for momentum transport, and Ohm’s law for electric current: the flux is proportional to the driving force. 1,2 The driving force is commonly represented by the gradient of the transport potential. The proportionality constant between the flux and the driving force is an important property of the system. For the case of diffusion, this constant is the diffusion coefficient. Such framework assumes that the gradient of the transport potential is sufficiently small to consider a linear response, i.e., the flux is linearly proportional to the driving force. Whenever the gradient of the transport potential is substantially large, one must consider higher order terms. In the last decades, molecular simulations emerged as a powerful set of techniques capable of solving problems in statistical mechanics numerically. Nevertheless, the ability to calculate a certain macroscopic property from the simulation data is strongly dependent on the proper statistical mechanics formulation. For a uniform fluid, Brownian motion has been used as a model to interpret the simulation data. 3 Einstein-Smoluchowski relation between the mean square displacement and the diffusion coefficient has been extensively used for this purpose. 4,5 This relation comes from the solution to the so-called Fick’s second law, for which, provided that specific boundary and initial conditions are stipulated, the spread of the particles follows a Gaussian evolution. Another approach consists in the equivalence 2


Page 2 of 35

Page 3 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


between a transport property and the integral of the correlation function of the flux. This entails the Green-Kubo framework, which is also derived considering the linear response theory through the fluctuation-dissipation theorem. 6 For an inhomogeneous fluid (a fluid in confined media or with an interface 3,7 ), however, the relation between the mass flux and the driving force is nonlinear, hence the standard methods ought to be replaced. Although one might think inhomogeneity is an exception, inhomogeneity is rather the rule both in nature and in industrial processes. For such systems, the modeled spatial-time distribution of the particles in a defined control volume, Ω, is governed by the Smoluchowski equation: 8 ∂p(r, t) = −∇ · J = ∇ · De−βW (r) · ∇ eβW (r) p(r, t) , ∂t

(1)

where p(r, t) is the probability density function of finding a particle within a thin layer Ω, r is the position vector, t is time, J is the mass flux, D is the diffusion tensor, β = 1/kB T , where kB is the Boltzmann constant and T is the absolute temperature, W (r) is the potential of mean force related to the pair correlation function: e−βW (r) = g(r) = ρ(r)/hρi, where ρ(r) is the particle density distribution and hρi is the average density in Ω. The diffusion tensor is diagonal. In Cartesian coordinates, assuming inhomogeneity in a single direction, one may define a perpendicular component of the tensor (which is orthogonal to the plane of the confinement or the interface) and the parallel components (which are parallel to the plane of the confinement or the interface). The parallel components of the diffusion tensor may be straightforwardly calculated by the method proposed by Liu et al. 9 This method may be seen as a generalization of the Einstein-Smoluchowski relation considering the survival probability of the particles within a specific Ω. Recently, we have shown that even the parallel components of the tensor may differ in the same Ω near an anisotropic crystal surface such as that of calcite. 10 For the perpendicular component of the diffusion tensor, several methods have been

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

proposed. Liu et al. 9 proposed a dual simulation method in which the Langevin dynamics survival probability at long time is fitted to the survival probability calculated through molecular dynamics simulations. They applied this method to the water vapor-liquid interface, showing an anisotropy between the perpendicular and the parallel components of the diffusion tensor at the interface. Mittal et al. 11 developed a new method to calculate the perpendicular component of the diffusion tensor based on the discretization of the Smoluchowki equation and on a Bayesian analysis of transition rates. 8,12 They applied such method to study the dynamics of a hard-sphere fluid confined between parallel walls. Hansen et al. 13 applied the Mean First-Passage Time analysis, disentangling the free energy and diffusivity contributions, 14,15 to the dynamics of water confined by a dipalmitoylphosphatidylcholine bilayer. Carmer et al. 16,17 proposed a new approach, called steady-state color reactioncounterdiffusion method, to determine the diffusion coefficient based on the calculation of the steady-state flux and the gradient of the local mole fraction of the labeled particle. In this work, we propose a new approach for calculating the perpendicular component of the diffusion tensor for inhomogeneous as well as for homogeneous fluids with a single molecular dynamics simulation. Solving the Smoluchowski equation analytically for a linear potential of mean force, we are able to determine the self-diffusion coefficient from the integration of the survival probability. For bulk systems, the same equation may be used considering a null gradient of the potential of mean force. Thus, we propose a single equation that can be used for calculations of the perpendicular component of the diffusion tensor in both homogeneous and inhomogeneous systems. Needless to say, in bulk systems, all the Cartesian components of the diffusion tensor are equal to each other.

4


Page 4 of 35

Page 5 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

Theoretical development

The Smoluchowski equation for one direction in Cartesian coordinates is given by the following expression: βW (r) ∂p(r, t) ∂ −βW (r) ∂ D⊥ e e p(r, t) , = ∂t ∂r ∂r

(2)

where D⊥ is the perpendicular component of the diffusion tensor in relation to the planes defining an arbitrary thin layer Ω. Assuming a thin layer Ω defined by the spatial interval r = [0, L], one may consider an effective diffusion coefficient that does not depend on the position coordinate. Therefore, the Smoluchowski equation may be rewritten as: 2 ∂p(r, t) ∂ ∂ p(r, t) ∂W (r) + = D⊥ p(r, t)β . ∂t ∂r2 ∂r ∂r

(3)

If one assumes that the solution to the time-dependent second order partial differential Equation (3) is separable in terms of time and space, one may write this solution as:

p(r, t) = R(r)T (t).

(4)

The validity of Equation (4) is assessed in the Results and Discussions section by comparison with molecular dynamics simulations. If Equation (4) is valid, one may rewrite Equation (3) as: 2 ∂ 2 W (r) d R(r) dR(r) ∂W (r) dT (t) . + = D⊥ T (t) β + R(r)β R(r) dt dr2 dr ∂r ∂r2

(5)

Rearranging Equation (5), one has: dT (t) 1 d2 R(r) 1 dR(r) ∂W (r) ∂ 2 W (r) 1 = + β + β . D⊥ T (t) dt R(r) dr2 R(r) dr ∂r ∂r2

(6)

Since the right side of Equation (6) depends only on position, r, and the left side only 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 35

on time, each side must be equal to a constant, hence: 1 dT (t) = −γ. D⊥ T (t) dt

(7)

The solution to Equation (7) is an exponential decay:

T (t) = T (0) exp [−γD⊥ t].

(8)

Because the dimensions of γ are 1/L2 , one may arbitrarily replace this parameter as: α . L2

γ=

(9)

Replacing the solution for time evolution on the general solution: αD⊥ p(r, t) = R(r)T (0) exp − 2 t . L

(10)

The survival probability, P (t), is defined as the integral over Ω:

P (t) =

Z

L

p(r, t)dr = 0

Z

L 0

αD⊥ R(r)T (0) exp − 2 t dr. L

(11)

Equation (11) may be further rearranged as:

αD⊥ P (t) = exp − 2 t L

Z

L

R(r)T (0)dr.

(12)

0

One may define a residence time, τr , as the integral of the survival probability over time:

τr =

Z

+∞

P (t)dt.

(13)

0

Therefore, the residence time for the general solution expressed in Equation (10) is given

6


Page 7 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


by: L2 τr = αD⊥

Z

L

R(r)T (0)dr.

(14)

0

Nevertheless, the integral in Equation (14) equals the survival probability at t = 0, which is equal to 1 following Kolmogorov’s axioms; 18 hence the diffusion coefficient may be written as: D⊥ =

L2 . ατr

(15)

Equation (15) is the general expression for the diffusion coefficient. We must stress that such equation is not defined for a specific limit as in Einstein-Smoluchowski approach. The value of L is arbitrarily specified by who is applying the method. The residence time, τr , might be determined from molecular dynamics simulations. Therefore, the only missing parameter in this equation is α. In the next sections, we shall derive analytical expressions for α for two specific cases: a linear potential of mean force; and a null gradient of the potential of mean force, i.e., a bulk system. For multi-component systems, by solving the Smoluchowski equation for each component independently, the diffusion coefficient of each component i is given by: D⊥,i

L2 = . α i τri

(16)

where αi is calculated taking into consideration the potential of mean force given by the component i density profile, and τri is the residence time calculated through the time integral of the component i survival probability.

2.1

Linear potential of mean force

Since Ω must be sufficiently small so that the diffusion coefficient may be considered independent of the position, it seems reasonable that the position dependence of the potential of mean force in such a small layer may, in some cases, be well approximated by a linear relation: −βW (r) = ωr + ξ. 7


(17)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 35

If the potential of mean force is linear, this implies that the natural logarithm of the equilibrium density is also linear, for e−βW (r) = ρ(r)/hρi: ρ(r) = ωr + ξ, ln ρ∗

(18)

where ρ(r) is the mass density distribution in kg·m3 and ρ∗ = 1.0 kg·m−3 is the normalization parameter. The Smoluchowski equation for one direction - Equation (2) - may be rewritten as: 2 ∂ p(r, t) ∂p(r, t) ∂p(r, t) . −ω = D⊥ ∂t ∂r2 ∂r

(19)

Considering the absorbing Dirichlet boundary conditions:

p(0, t) = p(L, t) = 0, ∀t > 0,

(20)

and the initial condition:

p(r, 0) = R L 0

ρ(r) ρ(r)dr

=

ωeωr , r ∈ Ω, eωL − 1

(21)

one has for this set of conditions that the analytical solution to Equation (19), for ω < 2π/L, is given by: 2 2 +∞ nπr X n 1 − eωL cos (nπ) ω2 nπ 2πω exp − D⊥ t . (22) − p(r, t) = 2 ωL sin 2 2 L (e − 1) n=1 L L2 4 ω 2 + nLπ2 The full derivation of Equation (22) is presented in Appendix A. The survival probability, by its turn, is simply: i h 2 2 2 +∞ exp − (2j+1) π − ω D t X ⊥ L2 4 4ω e + 1 P (t) = . 2 2 ωL L (e − 1) j=0 ω 2 + (2j+1) π ωL

L2

8


(23)

Page 9 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The partial residence time, τ (t), as a function of time, is :

τ (t) =

Z

t 0

ii h h 2 2 2 +∞ 1 − exp − (2j+1) π − ω D t X ⊥ 2 L 4 4ω e + 1 . P (t)dt = 2 2 2 2 ωL 2 (2j+1) π D⊥ L (e − 1) j=0 ω 2 + (2j+1) π −ω ωL

L2

4

(24)

L2

The residence time is given by: +∞ 4ω eωL + 1 X τr = lim τ (t) = t→+∞ D⊥ L (eωL − 1) j=0 (2j+1)2 π2 − L2

1 ω2 ω2 + 4

(2j+1)2 π 2 L2

.

(25)

Therefore, the diffusion coefficient is:

D⊥ =

L2 , ατr

(26)

where:

α

2.2

−1

+∞ −1 eωL + 1 X 3ω 2 L2 ω 4 L4 4 4 2 2 (2j + 1) π + (2j + 1) π − . = 4ωL ωL (e − 1) j=0 4 4

(27)

Bulk systems

For bulk systems, the gradient of the potential of mean force is zero, hence the Smoluchowski equation becomes equal to the so-called Fick’s second law: ∂p(r, t) ∂ 2 p(r, t) . = D⊥ ∂t ∂r2

(28)

Considering a uniform distribution of the particles as the initial condition, one has p(r, 0) = 1/L, and considering the absorbing Dirichlet boundary conditions: p(0, t) = p(L, t) = 0, by the application of the separation variable method, 19 one has that the proba-

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 35

bility density distribution of finding a particle at position r at time t is equal to: +∞ 2 4 X sin (2j + 1) π Lr 2 π p(r, t) = exp − (2j + 1) D⊥ t . πL j=0 (2j + 1) L

(29)

The survival probability is defined by the space integral of the probability density distribution: P (t) =

8 π2

+∞ X

h exp − (2j + 1)2

(2j + 1)2

j=0

π 2 L

D⊥ t

i

.

(30)

The residence time is the time integral of the survival probability: +∞ 8L2 X 1 τr = 4 . π D⊥ j=0 (2j + 1)4

(31)

The series in Equation (31) is the sum of the inverse of odd numbers to the fourth power. This may be written as the difference between the sum of the inverse of all natural numbers to the fourth power and the sum of the inverse of the even numbers to the fourth power: +∞ X j=0

+∞

+∞

+∞

+∞

X 1 X 1 X 1 1 1 X 1 = − = − . 4 4 4 4 n (2n) n 16 n (2j + 1)4 n=1 n=1 n=1 n=1

(32)

Using the definition of Riemann zeta function:

ζ(s) =

+∞ X 1 , s n n=1

(33)

one may rewrite Equation (32) as: +∞ X j=0

15 1 ζ(4). 4 = 16 (2j + 1)

10


(34)

Page 11 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Replacing the value of the Riemann zeta function for s = 4, one has: +∞ X j=0

1 π4 = . 96 (2j + 1)4

(35)

Therefore, Equation (31) may be rewritten as:

τr =

L2 . 12D⊥

(36)

In terms of the diffusion coefficient, one has:

D⊥ =

L2 . 12τr

(37)

This means that for a bulk fluid: α = 12.

(38)

Applying the limit for ω → 0 on Equation (27) - derived for a linear potential of mean force -, one arrives at the same value for α as the one shown on Equation (38).

3 3.1

Computational details Molecular Dynamics simulations

We performed molecular dynamics simulations using GROMACS 4.6.5. 20 The Leap-Frog algorithm was applied with 2 fs as timestep to integrate Newton’s equations of motion, and periodic boundary conditions were applied in all directions. The SPC/E force field 21 for water was used with a cut-off radius of 0.9 nm. For molten sodium chloride, the BornHuggins-Mayer-Tosi-Fumi 22,23 potential was used with a cut-off radius of 0.9 nm. A united atom TraPPE force field 24 was used for methane and Xiao et al. 25 force field was applied to calcite crystal with a cut-off radius of 1.0 nm. The electrostatic interactions were calculated

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 35

using the Particle Mesh Ewald method. 26 For equilibration, the Berendsen thermostat and barostat were applied (τT = 1.0 ps and τp = 1.0 ps). 27 For the production stage, we applied the Nosé-Hoover thermostat (τT = 1.0 ps). 28,29

3.2

Bulk systems

For bulk liquid water, an energy minimization of an initial configuration with 1000 water molecules was run before the equilibration. For bulk molten sodium chloride, we used as initial configuration a face centered cubic crystal composed of 1372 ions. For both systems, two equilibration stages of 1 ns each were carried out, the first one in the NPT ensemble to keep the pressure at 1 bar, and the second in the NVT ensemble. All production stages were run for 40 ns. Positions and velocities were stored every 0.2 ps. For the calculation of the diffusivity, the whole production trajectory was divided in 5 blocks of 8 ns each. Tail corrections were applied for both pressure and energy. 3.2.1

Self-diffusion coefficients

For bulk systems, Ω was arbitrarily defined as the whole simulation box for each direction. The survival probability, P (t), was then calculated as the ratio between two correlation functions: N (t0 , t0 + t), which is the number of centers of mass that remain in Ω between t0 and t; and, N (t0 ), which is the number of centers of mass within Ω at t0 . Considering multiple time origins, one has:

P (t) =

τ −1 1 X N (t0 , t0 + t) . τ t =0 N (t0 )

(39)

0

The residence time was then calculated by numerical integration using trapezoidal rule:

τr =

Z

+∞

P (t)dt.

(40)

0

For each Cartesian direction, one component of the diffusion tensor was calculated using 12


Page 13 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Equation (26) with ω = 0, which is completely analogous to Equation (37) in this case. The value of the self-diffusion coefficient of the bulk liquid water was determined by: 1 D = Tr (D) . 3

3.3

(41)

Water vapor-liquid interface

For water vapor-liquid interface, after energy minimization of an initial configuration containing 1000 water molecules, an equilibration stage of 2 ns was carried out in the NVT ensemble. The production stage was run for 40 ns. Positions and velocities were stored every 0.2 ps. For the calculation of the diffusivity, the whole production trajectory was divided in 5 blocks of 8 ns each. No tail corrections were applied, for such corrections assume a homogeneous system. 4

3.4

Methane confined between two calcite crystals

Two calcite crystals, oriented with the plane {10¯14} orthogonal to the z direction and with dimensions of 4.99 nm × 4.856 nm × 1.212 nm, were used to confine 397 methane molecules (corresponding to an average density of 125 kg·m−3 ) in a space of 3.5 nm in the z direction. After energy minimization, 20 ns were run for equilibration in the NVT ensemble. The production stage was run for 2.5 ns, also in the NVT ensemble. Positions and velocities were stored every 10 fs. For the calculation of the diffusivity, the whole production trajectory was divided in 5 blocks of 0.5 ns each. No tail corrections were applied, for such corrections one needs to assume a homogeneous system. 4 3.4.1

Parallel self-diffusion coefficients

Our method is derived for the perpendicular component of the diffusion tensor. For the parallel components, we used the method proposed by Liu et al. 9 In this method, given a layer

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 35

Ω, the parallel self-diffusion coefficient, Dk , is calculated through the following expression: h∆r2 (t)iΩ , t→+∞ 2tP (t)

Dk = lim

(42)

where h∆r2 (t)iΩ is the mean square displacement of the centers of mass that remain in Ω and P (t) is the survival probability, which may be calculated using Equation (39). 3.4.2

Perpendicular self-diffusion coefficients

For the perpendicular component of the diffusion tensor, we calculated the survival probability using Equation (39) and integrated following Equation (40), then we applied Equation (26) using the equilibrium density profile to calculate ω.

4

Results and Discussions

4.1

Bulk liquid water

Since our method is intended for both homogeneous and inhomogeneous systems, we present results for bulk liquid water as an example of a homogeneous system, for such results can be compared to published experimental data as well as to the traditional Einstein-Smoluchowski method. The Einstein-Smoluchowski method may be represented by the following expression: h∆r2 (t)i , t→+∞ 2dt

D = lim

(43)

where h∆r2 (t)i is the mean square displacement, and d is the dimensionality. Since the diffusion coefficient is defined by a limit, one must fit the mean square displacement at long times as a linear function of time, whose half of the slope value is the diffusion coefficient. Figure 1 shows the self-diffusion coefficients of bulk liquid water at atmospheric pressure for different temperatures. The results from molecular dynamics simulations using SPC/E force field are in excellent agreement with the experimental data. 30 Moreover, both tech14


Page 15 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


niques, the Einstein-Smoluchowski method and the unified approach we present here, give essentially the same results within the same standard deviations. Our approach, however, is more than five time faster in computational terms. For avoiding the square operation, albeit using the same multiple time origins, our approach is only a counting procedure to calculate the survival probability.

4.2

Molten sodium chloride

For bulk multi-component systems, the diffusion coefficient of a component depends only on the residence time of such component. Figures 2 and 3 present the results for self-diffusion coefficients of sodium and chloride in molten sodium chloride systems. For the cation, our results agree pretty well with the experimental data, 31 with an average absolute deviation of 2.7%. The self-diffusion coefficient of the chloride ion, however, is overestimated (with an average absolute deviation of 15.2 %). Figure 4 presents a comparison between the ratio of the sodium and chloride self-diffusion coefficients as function of temperature calculated through the Einstein-Smoluchowski method 32 and through our new approach. The high consistency between both approaches indicates that the overestimation observed for the chloride self-diffusion coefficient is not a consequence of the diffusion calculation itself.

4.3

Water vapor-liquid interface

As an example of inhomogeneous system, we chose the water vapor-liquid interface, so we can compare our results with the theoretical calculations done by other authors. 9 Figure 5 presents the results for the components of the water diffusion tensor for vaporliquid equilibrium at 298.15 K. In the liquid region, all components are equal to each other. At the interface, however, there is anisotropy between the parallel and the perpendicular components. This anisotropy, of the same magnitude, was noted by Liu et al. 9 We calculated the perpendicular component of the diffusion tensor at the interface using Equation (26), i.e., considering a linear potential of mean force. In Figure 5, we also present the results 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for the perpendicular component of the diffusion tensor using the discretized Smoluchowski equation (details of the discretization may be found on Appendix B). Both results are very similar to each other. Figure 6 shows the probability density function time evolution for SPC/E water at the vapor-liquid interface at 298.15 K. The probability density function is calculated through Equation (22) and through the discretized Smoluchowski equation. Despite the different initial configurations that result from the different approaches, quite rapidly both surfaces match each other. This means that the gradient of the potential of mean force is more relevant than the detailed initial configuration for the time evolution of the probability density function. Figure 7 and 8 present a comparison between the linear potential of mean force approach and the discretized Smoluchowski equation with the molecular dynamics simulation data at the water vapor-liquid interface at 298.15 K in terms of the survival probability and the partial residence time, respectively. The agreement between our approach and the simulation data is remarkable. This means that Smoluchowski equation is a good model to interpret this phenomenon. Also, this proves that, at least for this case, the solution to the linear potential of mean force Smoluchowski equation is indeed separable in terms of time and space, otherwise a disagreement in both curves would be expected.

4.4

Methane confined between two calcite crystals

The method proposed here has some limitations. Firstly, the analytical solution is restricted to linear potential of mean forces, i.e., the region to calculate the survival probability must be selected in a way that the density profile is well approximated by a straight line. Secondly, the analytical solution holds only for regions where ω < 2π/L. For strongly inhomogeneous fluids such as dense fluids near solid walls, this may be an issue. To further investigate these systems, we present some results for methane confined between two calcite crystals. Figure 9 displays the methane density profile at 375.0 K within the pore. The adsorp16


Page 16 of 35

Page 17 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tion of methane in calcite surface leads to a large peak of methane density near the wall. According to the Figure 9, from z = 1.35 nm to z = 1.45 nm, the methane density goes from practically 0 to a maximum density of ≈ 550 kg·m−3 . We have recently demonstrated that in such system the parallel diffusion coefficients, Dxx and Dyy , present an anisotropic behavior. 10 The perpendicular diffusion coefficient, Dzz , is also expected to have a different value than the parallel diffusion coefficients. In this region (z = 1.35 nm to z = 1.45 nm), we have found that Dxx = 37.8 ± 0.7 m2 ·s−1 and Dyy = 40.3 ± 0.8 m2 ·s−1 . Using the new approach proposed here, we have also found that Dzz = 5.0 ± 0.3 m2 ·s−1 . This means that the perpendicular self-diffusion coefficient of methane near calcite surface is almost ten times lower than the parallel self-diffusion coefficients, which implies that the mobility of methane molecules is much slower in the direction of the confinement near the wall.

5

Conclusions

We have developed a new unified approach to calculate the perpendicular self-diffusion coefficient in both homogeneous and inhomogeneous systems. Using bulk liquid water and molten sodium chloride as examples of homogeneous systems, we have proved that our method is as reliable as the traditional Einstein-Smoluchowski method. We have also applied our method to the water vapor-liquid interface, finding the same anisotropy noted by Liu et al. 9 We have shown that the method is also applicable for strongly inhomogeneous fluids, such as adsorbed methane in calcite crystal surface. Compared to other methods available in the literature, it is easy to implement the method we propose as it only needs the survival probability. We believe this method might be very helpful for calculations of self-diffusion coefficients in confined media or in systems containing interfaces.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Acknowledgement This publication was made possible by NPRP Grant No. 8-1648-2-688 from the Qatar National Research Fund (a member of Qatar Foundation). The statements made herein are solely the responsibility of the authors. We thank the High Performance Computing Center of Texas A&M University at Qatar for generous resource allocation. We are also indebted to Dr. Vasileios Michalis for his help with the Riemann zeta function to simplify the series for the self-diffusion coefficient of bulk systems.

References (1) Truesdell, C. J. Chem. Phys. 1962, 37, 2336–2344. (2) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; John Wiley & Sons Inc.: New York, 2002. (3) Novikov, D. S.; Fieremans, E.; Jensen, J. H.; Helpern, J. A. Nature Phys. 2011, 7, 508–514. (4) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publications: New York, 1987. (5) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2002. (6) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (7) Granick, S. Science 1991, 253, 1374–1379. (8) Hummer, G. New J. Phys. 2005, 7, 34. (9) Liu, P.; Harder, E.; Berne, B. J. J. Phys. Chem. B 2004, 108, 6595–6602.

18


Page 18 of 35

Page 19 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10) Franco, L. F. M.; Castier, M.; Economou, I. G. J. Chem. Phys. 2016, 145, 084702. (11) Mittal, J.; Truskett, T. M.; Errington, J. R.; Hummer, G. Phys. Rev. Lett. 2008, 100, 145901. (12) Sriraman, S.; Kevrekidis, I. G.; Hummer, G. J. Phys. Chem. B 2005, 109, 6479–6484. (13) Hansen, Y. V.; Gekle, S.; Netz, R. R. Phys. Rev. Lett. 2013, 111, 118103. (14) Hinczewski, M.; Hansen, Y. V.; Dzubiella, J.; Netz, R. R. J. Chem. Phys. 2010, 132, 245103. (15) Sedlmeier, F.; Hansen, Y. V.; Mengyu, L.; Horinek, D.; Netz, R. R. J. Stat. Phys. 2011, 145, 240–252. (16) Carmer, J.; van-Swol, F.; Truskett, T. M. J. Chem. Phys. 2014, 141, 046101. (17) Carmer, J.; Jain, A.; Bollinger, J. A.; van-Swol, F.; Truskett, T. M. J. Chem. Phys. 2015, 142, 124501. (18) Kolmogorov, A. N. Foundations of the theory of probability, 2nd ed.; Chelsea Publishing Company: New York, 1956. (19) Kirkwood, J. R. Mathematical Physics with partial differential equations; Academic Press: Amsterdam, 2013. (20) Hess, B.; Kutzner,; ; van-der-Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (21) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269– 6271. (22) Mayer, J. E. J. Chem. Phys. 1933, 1, 270–279. (23) Tosi, M. P.; Fumi, F. G. J. Phys. Chem. Solids 1964, 25, 45–52. 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 35

(24) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569–2577. (25) Xiao, S.; Edwards, S.; Gräter, F. J. Phys. Chem. C 2011, 115, 20067–20075. (26) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092. (27) Berendsen, H. J. C.; Postma, J. P. M.; van-Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (28) Nosé, S. J. Chem. Phys. 1984, 81, 511–519. (29) Hoover, W. G. Phys. Rev. A 1985, 31, 1695–1697. (30) Holz, M.; Heil, S. R.; Sacco, A. Phys. Chem. Chem. Phys. 2000, 2, 4740–4742. (31) Janz, G. J.; Bansal, N. P. J. Phys. Chem. Ref. Data 1982, 11, 505–693. (32) Koishi, T.; Tamaki, S. J. Non-Cryst. Solids 1999, 250-252, 501–505.

6

Appendix A - Analytical solution to Smoluchowski equation with a linear potential of mean force

As shown in the Theoretical Development, the Smoluchowski equation for a linear potential of mean force may be written as: 2 ∂p(r, t) ∂ p(r, t) ∂p(r, t) . −ω = D⊥ ∂t ∂r2 ∂r

(44)

The absorbing Dirichlet boundary conditions are given by Equation (20). The initial condition is given by Equation (21). If the solution to this partial differential equation is separable, i.e., p(r, t) = R(r)T (t), then the Smoluchowski equation may be rewritten as: 2 dT (t) dR(r) 1 d R(r) 1 = −γ. −ω = D⊥ T (t) dt R(r) dr2 dr 20


(45)

Page 21 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For bulk systems, the evolution of the probability of finding a particle in Ω may be written as the multiplication of a function of time and a function of position that results in the Gaussian evolution. 9 For an inhomogeneous fluid, we do not know a priori whether such assumption is valid. Nevertheless, the perturbation introduced in the Smoluchowski equation to account for the potential of mean force depends exclusively on the position, hence one may expect that the spatial-time evolution of the probability density function might also be separable. Following such argument, one has:

T (t) = T (0) exp (−γD⊥ t) ,

(46)

d2 R(r) dR(r) −ω + γR(r) = 0. 2 dr dr

(47)

and:

The general solution to the ordinary differential Equation (47) is:

R(r) = a1 eλ1 r + a2 eλ2 r ,

(48)

where a1 and a2 are constants, and λ1 and λ2 are the eigenvalues for ω < 2π/L:

λ=

ω ± κi, 2

(49)

where: κ=

r

γ−

ω2 . 4

(50)

Using Euler’s formula, one has:

R(r) = c1 cos (κr) + c2 sin (κr) ,

where c1 and c2 are constants.

21


(51)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 35

Applying the absorbing boundary conditions, one has:

c1 = 0,

(52)

c2 = c,

(53)

κ=

nπ , L

(54)

where n is an integer. Therefore, the position-dependent function, R(r), equals:

R(r) = c sin

nπr L

.

(55)

Substituting Equation (50) on Equation (54), one has:

−γn = −

n2 π 2 ω 2 − L2 4

(56)

Therefore, for any positive integer n, one has:

pn (r, t) = Rn (r)Tn (t) = sin

nπr L

2 2 ω2 nπ exp − D⊥ t . − L2 4

(57)

The superposition ansatz is:

p(r, t) =

+∞ X

cn pn (r, t) =

n=1

+∞ X

cn sin

n=1

nπr L

2 2 ω2 nπ exp − − D⊥ t , L2 4

(58)

which is also a solution to Equation (44). The constants cn must be determined using the initial condition:

p(r, 0) =

+∞ X n=1

cn sin

nπr L

22

= RL 0

ρ(r) ρ(r)dr

=


ωeωr . eωL − 1

(59)

Page 23 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Applying the inverse Fourier transform to Equation (59), one has: 2ω cn = L (eωL − 1)

Z

L

eωr sin 0

nπr L

dr.

(60)

Therefore: +∞

X 2ω p(r, t) = L (eωL − 1) n=1

Z

L

e

ωr

sin

0

nπr L

dr sin

nπr L

2 2 ω2 nπ exp − D⊥ t . − L2 4 (61)

Solving the inner integral on Equation (61) by parts, one has: 2 2 +∞ nπr 2πω X n 1 − eωL cos (nπ) ω2 nπ p(r, t) = ωL D⊥ t . exp − − sin (e − 1) n=1 (n2 π 2 + ω 2 L2 ) L L2 4

(62)

The survival probability, P (t), is given by the integral over Ω:

P (t) =

Z

L 0

i h 2 2 2 +∞ exp − (2j+1) π − ω D t X ⊥ e +1 L2 4 p(r, t)dr = 4ωL ωL . 2 2 (e − 1) j=0 (2j + 1) π + ω 2 L2

ωL

(63)

The partial residence time, τ (t), is defined as:

τ (t) =

Z

+∞

3

P (t)dt = 0

ωL

4ωL e + 1 D⊥ (eωL − 1)

+∞ X j=0

h

h 2 2 π − 1 − exp − (2j+1) L2

ω2 4

D⊥ t

(2j + 1)4 π 4 + 43 (2j + 1)2 π 2 −

ii

ω 4 L4 4

.

(64)

Thus, the residence time is: +∞ −1 4ωL3 eωL + 1 X 3 ω 4 L4 4 4 2 2 τr = lim τ (t) = (2j + 1) π + (2j + 1) π − . t→+∞ D⊥ (eωL − 1) j=0 4 4

(65)

And, finally, the perpendicular self-diffusion coefficient is given by: +∞ −1 ω 4 L4 3 4ωL3 eωL + 1 X 2 2 4 4 . (2j + 1) π + (2j + 1) π − D⊥ = τr (eωL − 1) j=0 4 4

23


(66)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

7

Page 24 of 35

Appendix B - Discretization of Smoluchowski equation

Applying the Forward in Time and Centered in Space (FTCS) method, Smoluchowski equation may be discretized as:

pi+1,j = pi,j +

Dδt [φ0,j pi,j+1 + φ1,j pi,j + φ2,j pi,j−1 ] (δr)2

(67)

where i is the index for time evolution, j is the index for spatial discretization, δt is the increment in time, δr is the increment in space, D is the diffusion coefficient, pi,j is the probability distribution function at position j and at time i, and:

φ0,j φ1,j φ2,j

1 ρj+1 = 1 − ln 4 ρj−1 ρj+1 ρj−1 δr ln =− 2+ 2 ρ2j 1 ρj+1 = 1 + ln 4 ρj−1

(68) (69) (70)

where ρj is the discretized density obtained from molecular dynamics simulations. This method is conditionally stable. A very conservative stability condition that we employed is: N Dδt X φ0,j ≤ 1 (δr)2 j=1

(71)

N Dδt X φ1,j ≤ 1 (δr)2 j=1

(72)

N Dδt X φ2,j ≤ 1 (δr)2 j=1

(73)

where N is the number of points in the spatial discretization of Ω. The same absorbing Dirichlet boundary conditions that we used for the analytical solution 24


Page 25 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


were applied for the discretized method:

pi,j=1 = pi,j=N = 0, ∀i ≥ 1

(74)

The initial condition was calculated from the equilibrium density profile using numerical integration for the normalization:

pi=1,j = R L 0

25

ρj ρ(r)dr


(75)


5.0 D × 109 / m2·s−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.0 3.0 2.0 1.0 280.0

300.0 320.0 T /K

340.0

Figure 1: Self-diffusion coefficient of liquid water at 1 bar at different temperatures. Closed squares, experimental data. 30 Open symbols, results from molecular dynamics simulations with the SPC/E water model. Open triangles, Einstein-Smoluchowski method. Open circles, our new approach. The error bars represent the standard deviation.

26


Page 26 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


20.0 DNa+ × 109 / m2·s−1

Page 27 of 35

15.0 10.0 5.0 0.0

1100

1150 1200 T /K

1250

Figure 2: Diffusion coefficient of Na+ in molten sodium chloride at different temperatures. Closed triangles, Arrhenius fitting from experimental data. 31 Open circles, results from molecular dynamics simulations with our new approach. The error bars represent the standard deviation.

27



20.0 DCl− × 109 / m2·s−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

15.0 10.0 5.0 0.0

1100

1150 1200 T /K

1250

Figure 3: Diffusion coefficient of Cl− in molten sodium chloride at different temperatures. Closed triangles, Arrhenius fitting from experimental data. 31 Open circles, results from molecular dynamics simulations with our new approach. The error bars represent the standard deviation.

28


Page 28 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.3 1.2 DNa+ /DCl−

Page 29 of 35

1.1 1.0 0.9 0.8

1100

1200 T /K

1300

Figure 4: Ratio between sodium and chloride diffusion coefficients in molten sodium chloride at different temperatures, calculated via molecular dynamics simulations. Closed circles, Einstein-Smoluchowski method. 32 Open circles, our new approach. The error bars represent the standard deviation.

29



4.5

800

4.0

600

3.5

400

3.0

200

2.5

0 0.0

0.5

1.0 1.5 z / nm

2.0

D × 109 / m2·s−1

5.0 1000 ρ / kg·m−3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.0 2.5

Figure 5: SPC/E water vapor-liquid equilibrium at 298.15 K. Continuous black line, density profile to be read at the left-hand y-axis. Open circles, Dxx calculated using Liu et al. 9 method. Open squares, Dyy calculated using Liu et al. 9 method. Open down-pointing triangles, Dzz calculated through Equation (26). Open up-pointing triangles, Dzz calculated through the discretized Smoluchowski equation. The error bars represent the standard deviation. All diffusivities ought to be read on the right-hand y-axis.

30


Page 30 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


p(z, t) / nm−1

Page 31 of 35

6.0 4.0 2.0 0.0 2.05 0.0

1.0

1.95 2.0

t / ps

3.0

4.0

1.75

1.85 z

/ nm

Figure 6: Time evolution of the probability density function for SPC/E water at the vaporliquid interface at 298.15 K. Green surface, Equation (22). Red surface, results from the discretized Smoluchowski equation.

31



1.0 0.8 P (t)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.6 0.4 0.2 0.0 0.0

10.0

20.0

30.0

t / ps Figure 7: Survival probability for SPC/E water at the vapor-liquid interface at 298.15 K as funcion of time. Open circles, results from Molecular Dynamics simulations. The error bars represent the standard deviation. Green curve, Equation (23). Red curve, results from the discretized Smoluchowski equation.

32


Page 32 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.0 3.0 τ (t) / ps

Page 33 of 35

2.0 1.0 0.0 0.0

10.0

20.0

30.0

t / ps Figure 8: Partial residence time for SPC/E water at the vapor-liquid interface at 298.15 K as a function of time. Open circles, numerically integrated results from Molecular Dynamics simulations. The error bars represent the standard deviation. Green curve, Equation (24). Red curve, results from the discretized Smoluchowski equation.

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Methane density as a function of the position along the z direction at 375.0 K.

34


Page 34 of 35

Page 35 of 35

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

35


Diffusion in Homogeneous and in Inhomogeneous Media: A New

Recommend Documents