2.28. MedeA Diffusion: Reliable Mass Transport Properties from Classical Simulations

Contents

Introduction
The MedeA Diffusion Module
Theoretical Background

download:	`pdf`

2.28.1. Introduction

Diffusion controls a wide variety of processes and properties, including manufacturing of semiconductor devices, environmental degradation of structural materials, corrosion, and permeability in polymers. The MedeA Diffusion module enhances your calculations by automatically computing the diffusivity of selected species using atomistic molecular dynamics techniques, and facilitates observation of the diffusive behavior of the different components.

Key Benefits

Automated plot creation facilitates analysis of calculation results
Visualization and confirmation of the presence of the diffusive dynamics regime through the provision of log(MSD) versus log(t) plots
Calculation of the mean-square displacement (MSD) of selected atom sets
Determination of diffusion coefficients based on the Einstein diffusion equation
Evaluation of calculation uncertainties

Hint

The MedeA Diffusion module works with classical molecular dynamics simulations using MedeA LAMMPS. Ab initio MD trajectories are not supported with the MedeA Diffusion module.

2.28.2. The MedeA Diffusion Module

The MedeA Diffusion module can be accessed in a MedeA LAMMPS Flowchart and you can insert the Diffusion stage into any LAMMPS stage. The MSD of selected atoms is computed from the simulation in the microcanonical (NVE) ensemble. Automated analysis of simulation results completes the workflow.

The following illustrates the MedeA Diffusion stage:

The parameters are:

Time: Duration of the simulation run. The appropriate value for this parameter depends on the speed of the diffusive motion of the components of the system being simulated. This can be determined by examining the MSD-time, and log₁₀(MSD)-log ₁₀(time) plots produced by calculations for the system in question.
Time step: Time step size employed in solving the equations of motion in an NVE ensemble.
Sampling: Number of samples employed in performing averaging. This parameter does not affect dynamics.
Trajectory: Number of trajectory frames saved during the molecular dynamics calculation. This parameter does not affect dynamics.
Subset: One or more subsets of atoms for which MSD information will be accumulated.

2.28.3. Theoretical Background

The mean-square displacement of a system of 216 argon atoms with a density of 0.68043 g/cm³ at 150 K as a function of time is shown below:

This approximately linear plot indicates the presence of a fully random Brownian motion where the Einstein relation is valid for extracting diffusivity from the mean-square displacement:

(1)\[D = \lim_{t\to\infty}\dfrac{\langle (r(t)-r(t_0))^2 \rangle}{2d(t-t_0)}\]

where the numerator is the mean-squared displacement, r is the position, t the elapsed time, d is the dimensionality of the simulation cell, and D is the self-diffusion coefficient.

As illustrated, the computation of D is best achieved by employing the central portion of the mean-square displacement from a given trajectory. This avoids the intrinsic nonlinearity present at short timescales where few collisions have occurred. (In this region, sometimes called the ballistic regime, displacement is linearly proportional to elapsed time because no collisions have yet occurred and the mean-squared displacement is, therefore, more quadratic than linear.)

Also, the final portion of the mean-square displacement plot may be affected by the sample averaging scheme employed (depending on the number of steps included in each sample). Hence, the Diffusion stage employs the central 50% of the time range in its linear regression analysis of the mean-square displacement plot, as illustrated.

The raw data associated with this plot are stored in the text file, {stageid}_msd.gdt, where {stageid} represents the label of the Diffusion stage in this particular MedeA LAMMPS Flowchart, and may be further analyzed using any appropriate plotting or analysis software.

The MedeA environment will also place graphics files (in both gif and png formats) for the mean-square displacement data in the stage directory. These provide a convenient visual indication of the diffusive behavior of the set of atoms specified in this calculation. The components of the mean-square displacement in the X, Y, and Z simulation cell directions are also provided in these output files, allowing for the straightforward analysis of non-isotropic diffusion.

The Diffusion stage reports the self-diffusion coefficient of the selected atoms based on the Einstein relation and thus relies on the presence of fully random Brownian motion where the Einstein relation is applicable. For many systems, achieving the truly random collisions necessary for Brownian diffusion may require substantial simulation times. The smallest diffusion coefficients that can be determined using this technique will depend on the available computational resource, but as a rule of thumb, the smallest measurable diffusivities will be of the order of 10^-7 cm²/s on most computers. Note however that, provided no changes in structure occur with changes in temperature (e.g. first-order or second-order phase transitions), diffusivities often follow Arrhenius behavior. Consequently, the range of diffusivities that can be predicted using simulations can often be extended significantly - e.g. by at least 3 orders of magnitude - by extrapolation of values determined at a series of higher temperatures to lower temperatures of interest.

Given the Einstein relation, which requires that for any diffusing species the MSD be proportional to time and not to some other power, we can confirm the presence of Brownian diffusion through a plot of log₁₀(MSD) versus log₁₀(t), (base 10 logarithms simply being chosen for ease of interpretation). If the resulting plot is linear with a unit gradient, we have simulated long enough to verify that Brownian diffusion has been observed. Further discussion of this topic may be found in popular reference texts [1], [2].

The MedeA environment automatically carries out this analysis, and reports the gradient of the plot and the characteristics of the best fit line to these data in the Job.out and {stageid}_log10_msd_log10_t.gif and png files in the Diffusion stage folder. A log₁₀(MSD) versus log₁₀(t) plot, for the argon system described above, is provided below:

[1]	M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, 1st Edition 1987, Ch 6.5; 2nd Edition 2017, Ch. 8.5)

[2]	D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, 1st Edition 1996; 2nd Edition 2002, Ch. 4.5)

download:	`pdf`