Estimation of a Gas Diffusion Coefficient by Fitting Molecular Dynamics Trajectories to Finite-Difference Simulations

TL;DR

This work addresses estimating the diffusion coefficient $D$ of argon diffusing in helium by reconciling particle-based MD trajectories with a continuum FD diffusion model. It fits the MD-derived argon concentration, obtained by binning MD trajectories onto a finite-difference grid, to the FD solution of the diffusion equation using a Levenberg–Marquardt optimization over $D$, implemented with Crank–Nicolson time stepping for numerical stability. The authors demonstrate that the estimated $D$ from 2D simulations can approach the experimental value, with the best agreement observed at an intermediate MD-bin resolution, while discussing factors that contribute to remaining discrepancies such as dimensional reduction, potential model limitations, and statistical noise. The study provides a practical framework for extracting transport coefficients from MD data via continuum-model fitting, with potential extensions to other multiscale transport problems and higher dimensions. The approach combines MD (LAMMPS) with FD solvers and nonlinear least-squares fitting to yield a robust estimate of $D$ from first-principles simulations.

Abstract

A procedure is presented to estimate the diffusion coefficient of a uniform patch of argon gas in a uniform background of helium gas. Molecular Dynamics (MD) simulations of the two gases interacting through the Lennard-Jones potential are carried out using the LAMMPS software package. In addition, finite-difference (FD) calculations are used to solve the continuum diffusion equation for the argon concentration with a given diffusion coefficient. To contain the computational cost and facilitate data visualization, both MD and FD computations were done in two space dimensions. The MD argon trajectories were binned to the FD grid, and the optimal diffusion coefficient was estimated by minimizing the difference between the binned MD data and the FD solution with a nonlinear least squares procedure (Levenberg-Marquardt algorithm). Numerical results show the effect of the MD binning parameter and FD grid spacing. The estimated diffusion coefficient is compared to an experimental measurement.

Figures (5)

• Figure 1: Forward Euler/central difference scheme \ref{['ds:HE_FD']}, the amplification factor $\rho_{\bf m} = 1+kD\lambda_\mathbf{m}$ is plotted versus scaled wavenumber $|\mathbf{m}|/N$ for three time step values near the critical value, $k < k_c, k = k_c, k > k_c$; for simplicity the spatial dimension is $d=1$.
• Figure 1: Lennard-Jones potential $V$ given in \ref{['eq:6-12_pot']}, (a) scaled potential $V/\epsilon$ versus scaled distance $r/\sigma$, (b) dimensional potentials $V$ (kcal/mol) for He-He, Ar-Ar, and mixed Ar-He interactions versus interparticle distance $r$ (Å) with parameters $\epsilon, \sigma$ given in Table \ref{['tb:lj_diff']}heliumargon.
• Figure 1: Numerical results for diffusion of argon in helium, color bar indicates normalized argon concentration, MD binning parameter $N=50$, FD solution of diffusion equation \ref{['eq:HE']} computed on a $50 \times 50$ grid, Columns 1-5 display times $t=0,0.25,1,2.5,5$ ns, Row 1: MD snapshots of argon atoms, Rows 2, 4, 6: binned MD concentration, Rows 3, 5, 7: computed FD concentration.
• Figure 1: The cost function ${\tt cost}(D)$\ref{['eq:cost']} of the diffusion coefficient estimation algorithm is plotted versus the diffusion coefficient $D$ near the global minimum $D_\text{opt}(N)$ with MD binning parameters $N=20, 50, 100$.
• Figure 2: Crank-Nicolson scheme \ref{['ds:HE_CN']}, the amplification factor $\rho_{\bf m} =(1+\frac{1}{2}kD\lambda_\mathbf{m})/(1-\frac{1}{2}kD\lambda_\mathbf{m})$ is plotted versus scaled wavenumber $|\mathbf{m}|/N$ is plotted for the same time steps $k$ as in Figure \ref{['fg:HE_amp']}; for simplicity, the spatial dimension is $d=1$.