# Nuclear induction lineshape: Non-Markovian diffusion with boundaries

## Data files

## Abstract

The dynamics of viscoelastic fluids are governed by a memory function, which is critical yet computationally intensive to determine, particularly when diffusion is restricted by boundaries. We introduce a computational method that effectively captures the memory effects by analyzing the time-correlation function of the pressure tensor, an indicator of viscosity, through the analytic continuation of the Stokes-Einstein equation to the Laplace domain. This equation is integrated with molecular dynamics (MD) simulations to obtain necessary parameters. Our method computes NMR lineshapes by employing a generalized diffusion coefficient that incorporates the influences of temperature and confinement geometry. This approach establishes a direct link between the memory function and thermal transport parameters, enabling precise computation of NMR signals for non-Markovian fluids in confined geometries.

## README: Nuclear Induction Lineshape: Non-Markovian Diffusion with Boundaries

https://doi.org/10.5061/dryad.m905qfv7n

The simulation was executed on the Hoffman2 cluster at UCLA. The "submit_job.sh" script was used for running a batch of simulations.

A sample of LAMMPS code for viscosity evaluation can be found in each folder, "in.visc_Bulk" and "in.visc_SiO2" for bulk and restricted simulations, respectively.

For data analysis in Jupyter Lab, we provide a Python code in the "Viscosity Analysis.ipynb" notebook.

#### Description of the data and file structure

The various folders are marked to indicate different boundary conditions. The "Xegas.zip" folder contains the dataset for simulating Xe gas particles in bulk. Folders named A11-A75 hold datasets for cylindrical boundary conditions, with the tube radii in Angstroms represented by two digits in each folder's name. Each folder contains a file with the LAMMPS code named "in.visc_SiO2", a file titled "submit_job.sh" for batch submission of simulations to the HPC facility, several log files documenting code performance, and 90 data files. The data file names have a specific structure: the first part indicates the simulation phase and the monitored parameter, the second part denotes the temperature, and the third part specifies the random seed used for the simulation. Each simulation is conducted with 10 different random seeds across 9 temperatures: 240, 260, 280, 300, 320, 340, 360, 380, and 400 K. For instance, a file named "viscg_280_1.dat" represents the results of a viscosity simulation for Xe gas particles at 280 K using the first set of random distributions.

These data are generated using a "ave/correlate" fix within the LAMMPS code. The resulting output includes rows displaying the timestep in the first column, followed by the simulated parameters in subsequent rows. Specifically, the components of the pressure tensor correlations, denoted as "v_pxy," "v_pxz," and "v_pyz," are presented in the last three columns. The Jupyter notebook contains a function demonstrating how to import and work with such a file.

#### Code/Software

The Lammps code has the following sections:

- Variable assignments: This section involves defining crucial variables, such as temperature (T), simulation dimensions (L), the number of particles (npart), and sampling intervals.
- Unit conversion to SI: Here, unit conversions are encoded to ensure consistency with the International System of Units (SI).
- Molecular Dynamics simulation setup: This part involves configuring the molecular dynamics simulation by specifying boundary conditions, defining the simulation region, setting coupling constants, and determining the time steps.
- NVT canonical ensemble fix: A fix is defined to maintain the NVT canonical ensemble during the simulation.
- Equilibration (Prerun): Prior to data collection, a prerun of 1e6 steps is executed to allow the system to reach an equilibrium state.
- Correlation term computation: This section of the code focuses on calculating the correlation terms for the pressure tensor based on the Kubo formulation. The necessary summations are performed.
- Main simulation run: The code proceeds to run for the desired number of steps, during which preliminary viscosity coefficient data is collected and printed for analysis.

The simulation was performed for various temperatures in the range of 200-400 K for xenon particles in gas phase for both bulk and restricted diffusion.

## Methods

MD simulations were conducted to calculate the shear viscosity of gaseous xenon (Xe), as defined by Eq. (7). We utilized the Lennard-Jones (LJ) pair potential, expressed as *U*(*r*) = 4*ϵ*[(*σ*/*r*)^{12} − (*σ*/*r*)^{6}], where interactions between xenon atoms were characterized by *ϵ* = 1.77 kJ/mol, the depth of the potential well, and *σ* = 4.1 Å, the distance at which the potential energy becomes zero. Our simulations of bulk fluid were conducted for isotropic diffusion by placing 2000 xenon atoms within a box defined by periodic boundary conditions. Throughout the simulations, we maintained a consistent particle count, volume, and temperature, adhering to the canonical ensemble (*NVT* ensemble). Each set of simulations was repeated for ten different random seeds for the initial positions and velocities of the particles to ensure robust statistical sampling and accuracy of the results. In the simulations of restricted diffusion (i.e., diffusion limited by the nanotube geometry), nanotubes of a fixed length and various diameters were employed, and the number of particles was adjusted to maintain a constant particle density.

For the simulations, we used xenon (Xe) particles. The interactions between the Xe particles and the cylindrical boundary were modeled using the Lennard-Jones potential, with parameters *ϵ* = 0.3 kJ/mol and *σ* = 4.295 Å, representing a silica tube. To determine the viscosity coefficient, we integrated all three components of the pressure tensor, referred to as *C*_{αβ}(*τ*). Although the *C*_{xy} component showed significantly higher values than the *C*_{yz} and *C*_{xz} components, given the tube’s orientation along the *z-axis*, we opted to integrate all three components together for a thorough analysis.