Lattice simulation method to model diffusion and NMR spectra in porous materials

A coarse-grained simulation method to predict NMR spectra of ions diffusing in porous carbons is proposed. The coarse-grained model uses input from molecular dynamics simulations such as the free-energy profile for ionic adsorption, and density-functional theory calculations are used to predict the NMR chemical shift of the diffusing ions. The approach is used to compute NMR spectra of ions in slit pores with pore widths ranging from 2 to 10 nm. As diffusion inside pores is fast, the NMR spectrum of an ion trapped in a single mesopore will be a sharp peak with a pore size dependent chemical shift. To account for the experimentally observed NMR line shapes, our simulations must model the relatively slow exchange between different pores. We show that the computed NMR line shapes depend on both the pore size distribution and the spatial arrangement of the pores. The technique presented in this work provides a tool to extract information about the spatial distribution of pore sizes from NMR spectra. Such information is diffcult to obtain from other characterisation techniques.


I. INTRODUCTION
Transport in porous materials is relevant for phenomena as diverse as energy conversion and storage, gas storage, heterogeneous catalysis, cell growth and drug delivery 1 . In all these cases, the performance of the porous material depends on its structural properties, the latter determining the geometry and dynamical properties of the confined fluid. The characterisation of porosity and the understanding of its effect on performance are therefore essential to optimise the porous materials for their respective applications. Nuclear Magnetic Resonance (NMR) is an experimental method that is uniquely suited to probe both microscopic structure and the dynamical properties of fluids confined in porous media.
NMR has been widely used to investigate the adsorption of various fluids in porous carbons [2][3][4][5] . Most of these studies report the observation of distinct resonances in the NMR spectra corresponding to adsorbed and freelydiffusing probe species. Irrespective of the nature of the fluid or nucleus studied, the resonances of the species adsorbed on carbon appear at lower frequencies than those of the free species, the difference in the adsorbed and free resonant frequencies depending only weakly on the nature of the fluid. The dominant shift mechanism arises from the magnetic fields caused by the microscopic ring currents of π electrons in the graphitic carbon 6 . Consistent with this, variations of the structure of the porous carbon have a pronounced effect on the observed NMR resonances 2,4,[7][8][9] .
Although NMR is a powerful technique to investigate both the structure of the carbon and the microscopic structure of the fluid at the carbon surface, the interpretation of the NMR spectra is by no means straightforward as the measured spectra depend in a non-trivial way on the pore structure and on the structure and dynamics of the confined fluid. As a consequence, simplifying assumptions are commonly made to interpret the experimental NMR spectra.
A necessary (but not sufficient) prerequisite for quantitative interpretation of the NMR spectra is an accurate description of the effect of aromatic ring currents on the local magnetic fields in the porous medium. It is possible to calculate Nucleus Independent Chemical Shifts (NICS) for aromatic molecules using Density Functional Theory (DFT) calculations 7,10-12 . These studies provide important insights, accounting, for example, for the shift to lower frequencies of NMR resonances of ions confined in a porous graphitic matrix. However, in order to arrive at a model that can describe all features of the experimental NMR spectra, the effects of local fluid structure and the ionic dynamics must be taken into account. A fully atomistic approach to compute such spectra, let alone an ab-initio calculation, would be prohibitively expensive in view of the wide range of relevant time scales.
Here, we propose a coarse-grained method that allows us to bridge the gap between molecular simulations and experiment. It is particularly attractive to use lattice models as these allow us to account for the relevant microscopic effects while being computationally much more efficient than off-lattice models. The key computational advantage of the lattice approach used here is that it allows us to sample, in a single simulation, all possible diffusive trajectories 13 , rather than a single one, as would be obtained from an off-lattice simulation or from a conventional lattice simulation of the diffusion of tracers in porous media [14][15][16] . The (exponential) computational advantage of this approach turns out to be crucial for ex-ploring the effect of multiple factors on the NMR spectra.
In the following section, we describe the implementation of the lattice gas model and show how it can be used to calculate diffusion coefficients and NMR spectra. We then focus on a two-site exchange model to validate our approach. In particular, we show that we can reproduce the coalescence of NMR resonances observed experimentally. Next, we apply the model to single mesoporous slit pores with different pore widths. Slit pore models have been used extensively in the literature to interpret adsorption isotherms in some classes of porous materials, such as disordered activated carbons for energy storage 17,18 , and to simulate ion adsorption and charge storage mechanisms in supercapacitors [19][20][21][22] . Their simplicity and wide utilisation make them a good starting point for the application of our model. We describe how the model is parametrised from molecular simulations and explain how diffusion and chemical exchange can lead to the observation of a single chemical shift corresponding to multiple environments. In the last part of this article, we calculate NMR spectra corresponding to ions diffusing in a carbon particle with a realistic pore size distribution and show that the resulting spectra are affected by the spatial distribution of pore sizes. While some experimental techniques, such as adsorption isotherms analysis, can provide the overall pore size distribution, they are not able to probe the spatial distribution of these various pore sizes. As such, the development of appropriate models to predict NMR spectra of species adsorbed in a realistic porous material will open the door to an appropriate interpretation of the NMR experimental data and give insights into both the carbon and liquid structure of the explored systems. This will have implications in various scientific fields where porous materials are used, such as supercapacitors, which will benefit from a detailed characterisation of the electrode/electrolyte interface and the interconnectivity of the pore structure.

A. Description of the model
The simplest lattice-gas models to simulate diffusion describe the diffusing species as non-interacting particles performing kinetic Monte Carlo moves on a lattice that contains both accessible sites (the fluid) and excluded sites (the porous matrix). A schematic representation of such a model is shown in figure 1. The lattice model is characterised by a lattice parameter a, which is the distance between two lattice nodes, and a timestep ∆t, which corresponds to the typical time it takes the probe particle to diffuse over a distance a. The geometry of the matrix is accounted for by the spatial distribution of excluded lattice sites.
At every timestep, the particles attempt with equal probability to perform a step to any of their nearest neighbours (6 for a three-dimensional simple cubic lattice): if the neighbour site is accessible, the probe particle is moved to the new position, otherwise the particle remains at its original position. We consider threedimensional cubic lattices with periodic boundary conditions. When modelling diffusion in a porous carbon matrix, we must account for the fact that different accessible sites have different energies (strictly speaking: free energies) and are, therefore, not equally populated. A site energy (E i for site i) is required to account for ionic adsorption at the carbon surface. In addition, when computing NMR lineshapes, we have to account for the fact that the local magnetic field and thus resonant frequency at site i, ν i , will depend on its position.
To describe diffusion in a spatially varying potential, we vary our Monte Carlo acceptance rules, such that the correct Boltzmann distribution over lattice sites is obtained in equilibrium. Specifically, the conditional probability to accept a Monte Carlo move from site i to site j, p acc (i → j), which is not necessarily equal to the probability of accepting the reverse move, p acc (j → i), is: where k B is the Boltzmann constant and T is the temperature of the system. As a consequence, the larger E j − E i the less likely a particle is to jump from i to j. With these rules, the probability of a particle to visit site i will be given by the Boltzmann distribution: where ρ i is the average density at site i and < ρ > is the density averaged over all lattice sites.

B. Diffusion coefficients
In MD simulations studies, diffusion coefficients are typically calculated using a Green-Kubo relation that relates the velocity autocorrelation function to the selfdiffusion coefficient: where D is the diffusion coefficient, d is the dimensionality of the system, v(0) is the velocity at time t = 0, v(t) is the velocity at time t and ... denotes an average over all particles and all trajectories. Here we use the discrete equivalent of Eqn. 3 (see Appendix A 1): FIG. 1. Illustration of the lattice gas model. The lattice, characterised by a lattice parameter a, is divided between two types of sites: accessible (fluid) and not accessible (porous matrix). If at a given timestep, a particle's velocity points towards an obstacle, then at the next timestep, that particle remains at its original position. The timestep ∆t relates to the lattice parameter a via the bulk diffusion coefficient D = a 2 /6∆t following Eqn. 6. Note that for the sake of simplicity, the system is represented in two dimensions here but the model is actually three-dimensional.
where ... denotes an average over all considered sites and ∆x i denotes the displacement of a particle in the x direction at timestep i. If we were to exploit the analogy with off-lattice systems, we would estimate ∆x 1 ∆x j by generating random trajectories for a number of particles starting from different sites in the lattice. This approach is highly inefficient as a very large number of trajectories would have to be generated to obtain good statistics for a heterogeneous system. Here we use a different approach, namely the so-called 'moment propagation' method, which was proposed by Frenkel and successfully applied to a number of studies on dynamics of particles in confined systems 13,23,24 . The moment-propagation method is a recursive scheme that allows us to sample all possible trajectories of the diffusing particles, rather than a subset. The computational effort scales as t×M , where t is the simulation time and M the number of lattice sites. The computational effort to sample all trajectories in a non-recursive scheme would scale as z t M where z is the coordination number of the lattice. In the expression for diffusion coefficients, the first term is easily calculated as it is equal to 1/2∆t times the mean-square displacement of a particle in the x direction during one timestep: where a x is the x-component of the vector joining a lattice site to its j-th neighbour and the sum runs over all sites j adjacent i (note that p acc (i → j) = 0 if j is occupied by an obstacle). The general moment-propagation approach to compute the diffusion coefficient is described in A 2.
In a homogeneous fluid, successive jumps are uncorrelated and we have: If the probability to jump to a neighbouring site is equal to one, then D = a 2 /(2d∆t), where a is the lattice spacing. However, we can reduce the probability that a particle carries out a jump. If we denote the probability that a particle stays on the same site by 1 − α, then In systems where the diffusivity varies with position, we can define a factor α(ij) for every link between neighbouring lattice sites i and j. It is important to note that if the probability to jump from i to j is reduced, then the probability to jump from j to i should be reduced by the same factor, otherwise the equilibrium distribution over lattice sites would be changed from the Boltzmann distribution.
C. NMR signal and spectrum NMR spectroscopy probes the nuclear magnetic response of a sample whose magnetism is perturbed from equilibrium by a radio frequency pulse. After this pulse, the transverse magnetisation of the sample decays. The NMR signal measured during this decay is commonly referred to as the Free Induction Decay (FID). The NMR spectrum is obtained by Fourier transforming the FID signal. In a heterogeneous sample, different nuclei will experience different local magnetic fields and will therefore have different Larmor frequencies. As the nuclei diffuse, their environments, and hence their resonant frequencies, change. The FID signal is the superposition of the signals corresponding to all excited nuclei in the sample. For an ensemble of probed spins, the NMR signal is given by: where ν i 0 (t) is the Larmor frequency corresponding to spin i at time t, and ... denotes an average over all spins. The spectrum is then obtained by Fourier transforming this signal following: Using the same approach as for the calculation of the diffusion coefficients (see Appendix A 3), we can discretise these expressions and estimate them by the momentpropagation method.
In practice, it is better to define the resonant frequency ν i as the difference beween the Larmor frequency at site i, ν i 0 , and the Larmor frequency of the bulk liquid (the reference). The resulting frequencies are typically in the kilo-Hertz range. The timestep ∆t should be chosen such that ∆t × ν max 1, where ν max is the largest value of ν i . We note here that the spectral width SW , i.e. the range of frequencies that is studied, is given by the dwell time, dwt, i.e. the time between two acquisition sampling points, following: This spectral width in Hertz can then be converted to a range of chemical shifts in ppm using the Larmor frequency of the studied nucleus (ν 0 ): Note that, the frequency differences (in Hz) between the various resonant frequencies experienced by the spins depend on the applied magnetic field, the frequency increasing linearly with field strength. This is taken into account in the model as will be clear from the study of various magnetic fields in the last section of this article.

D. Model parameters
In practical cases, we must map the system that we wish to simulate on a lattice model with discrete lattice sites and timesteps. In particular, we must specify: -the distribution of obstacles on the lattice to represent the porosity of the material -the distribution of local energy values that account for the adsorption of ions and molecules to the carbon surface -the distribution of local Larmor frequencies that account for the spatial variation of the chemical shifts. The allocation of these quantities to the different lattice sites will typically be based on structural information or use input data from more microscopic descriptions such as Density Functional Theory or Molecular Dynamics simulations (MD).
In addition, we must specify the lattice parameter (a), i.e. the distance between two lattice nodes in the model, the temperature (T ) and the timestep (∆t). These parameters are usually determined by the physical properties of the system to be modelled. The lattice parameter a will usually be defined so that it allows a sufficiently accurate description of the liquid structure. The choice of the timestep is then set by its relation with the diffusion coefficient D and the lattice parameter a, namely D = a 2 /(6∆t) in three dimensional lattices. We then check that the condition ∆t × ν max 1 is respected. If the diffusion coefficient or other quantities (free energies for example) correspond to a given temperature, then this temperature should be used in the lattice simulation.

A. Parametrisation of a two-site exchange model
To validate our numerical approach, we show that the model reproduces the lineshape for a two-site exchange model for which analytical solutions exist 25,26 . We focus on a simple model in which a spin can move between two different chemical environments, S 1 and S 2 , with equal free energies and equal populations. The process is described by the following equilibrium: where k is the exchange rate between the two sites. The NMR lineshape depends on the ratio between the exchange rate k and the difference in resonant frequencies ∆ν = |ν 1 − ν 2 |. Two different regimes can be observed: i) a slow-exchange regime where signals due to the two environments can be distinguished in the spectrum and ii) a fast exchange regime where only a single resonance is observed at a frequency corresponding to an average of ν 1 and ν 2 . The limit between these two regimes is known as the coalescence point where the two peaks merge into one. If we ignore the spin-spin and spin-lattice relaxation effects, coalescence occurs for an exchange rate equal to: where ∆ω = 2π∆ν = |2πν 1 − 2πν 2 | 25 .
FIG. 2. Lattice gas model and two site-exchange between sites with frequencies ν1 (-1 kHz) and ν2 (1 kHz). Here we illustrate the ability of the model to reproduce coalescence for a simple system. The NMR lineshape typically depends on the ratio between the effective exchange rate kexc between sites with different frequencies and the frequency difference ∆ν = |ν1 − ν2|.
We explore the effect of the spatial distribution of sites with different Larmor frequencies for a given correlation time of 0.1 ms (a,b,c). Here, each sphere represents a lattice site and the color distinguishes sites having a frequency ν1 or ν2. The spectra in panel d) illustrate the effect of doubling the number of interfaces, between blocks of distinct frequencies, from two to four: in both cases, two resonances are distinguishable but the resonances are broader when the effective exchange rate increases. e) In the planes case, we are close to the coalescence point and we observe a very broad spectrum where we cannot distinguish anymore the two different resonances. f) In the lines cases, we note a narrowing of the spectrum which is accentuated in g) by decreasing the correlation time τ .
To investigate two-site exchange with our model, we build a lattice of 40×40×40 sites (64,000 lattice sites in total) where a specific frequency (ν 1 = -1 kHz or ν 2 = 1 kHz) is assigned to each lattice site. For all the systems investigated, there are no obstacles so all sites are accessible, and the energies on all sites are equal so that all densities are equal. The model is characterised by a correlation time τ (equivalent to the timestep of the simulation) which is the average time that a diffusing species stays on a site. The inverse of this correlation time thus defines the exchange rate between two sites of the lattice. Not all the jumps between sites will lead to a change of frequency. We thus define an effective exchange rate, k exc , which represents the exchange rate for hops between sites with different frequencies, corresponding to the entire system. This is the rate that is measured in the NMR experiments.
In our model, k exc will depend on two key factors: i) the spatial distribution of the sites with different Larmor frequencies on the lattice, i.e. for a given exchange rate, the more contacts there are between sites with different frequencies, the larger k exc will be, ii) the correlation time, i.e. decreasing the correlation time leads to a global increase of the dynamics of the diffusing species and thus an increase of k exc . We first explore the effect of the spatial distribution of frequencies by keeping the correlation time τ constant (equal to 0.1 ms) and arranging chemical environments in different ways. Four different configurations were built (see figure 2): -blocks with 2 planar interfaces: the frequencies of the different sites are distributed in two blocks with equal number of sites, and thus there are only two planar interfaces (one in the middle of the system and one at the edge where periodic boundary conditions are applied) where diffusion will lead to a frequency change; -blocks with 4 planar interfaces; -planes: the frequencies are distributed such that alternating planes are characterised by ν 1 and ν 2 ; -lines: the frequencies are distributed such that one line in two is characterised by ν 1 and the other by ν 2 .
We then investigate the effect of a reduction of the correlation time by a factor two (τ = 0.05 ms) in the lines configurations. We calculate the NMR signals on a time period of 500 ms which corresponds respectively to 5,000 and 10,000 timesteps when τ equals 0.1 ms and 0.05 ms. This time is longer than the decay time so that there is no truncation of the NMR signal.

B. From slow exchange to fast exchange
The NMR spectra predicted for the different setups are given in figure 2. To relate these results to the outcomes of the analytical models described above, we calculate the effective exchange rate corresponding to the different systems. In all cases, the system is three-dimensional so that a particle will jump in one of the six available directions at each timestep. In the case of the blocks configuration with two interfaces, the number of sites where exchange is possible is equal to the number of sites in one plane (40×40) multiplied by two (because exchange can be in both directions) and by the number of interfaces. The average exchange rate will then be given by the probability of being in one of these sites multiplied by the probability of actually jumping in the direction corresponding to a frequency change (1/6): where τ is equal to 0.1 ms. For the blocks configuration with four interfaces, the effective exchange rate doubles to 0.33 kHz. These two frequencies are lower than k coal , equal to 4.44 kHz in our case, leading to a slow exchange regime where both peaks are distinguishable. For the case of planes and lines with a correlation time of 0.1 ms, the number of directions which lead to frequency changes are respectively 2 (out of 6) and 4 (out of 6), giving effective exchange rates of 3.33 kHz and 6.67 kHz, respectively. The first value is very close to the coalescence point and the corresponding spectrum is very broad. The second value falls in the fast exchange regime where peak narrowing is observed. The effective exchange rate is then increased by a factor two using a smaller correlation time of 0.05 ms leading to a further motional narrowing.
This illustration of a simple two-site exchange model shows that the numerical approach that we propose here can reproduce an albeit simple example of the collapse of the NMR spectrum. But while analytical models can only deal with simple cases such as two-site exchange, the numerical moment-propagation approach can be used to model more complex systems with realistic spatial distributions of resonant frequencies.

A. Parametrisation of the lattice model from molecular simulations
We now focus on the case of slit pores in an attempt to represent what would happen in one class of mesoporous materials. As stated in the description of the lattice gas model proposed here, each simulation goes through a parametrisation step where the obstacles' positions, free energies (E i ) and frequencies (ν i ) have to be assigned to discrete lattice sites. In the case of a slit pore, the positions of the obstacles are defined by two planes corresponding to the two pore walls. We assume that the slit is perpendicular to the z-axis and that there is no variation of free energies or Larmor frequencies with x or y. As a consequence of this and the use of periodic boundary conditions, we only need 1 site in the x and y directions. The lattice parameter a is set to 0.5Å and the number of lattice sites in the z direction, N z , is adapted to the pore width such that N z = (R/a)+1 with R being the pore width. We first focus on a slit pore with a pore width equal to 4 nm, and the number of lattice sites for this system is thus 1×1×81.
The values of the site free energies E i were fixed using existing results obtained from MD simulations. Specifically, we use the simulation results for butylmethylimidazolium tetrafluoroborate dissolved in acetonitrile ([BMI][BF 4 ] in ACN, 1.5 M) at the interface with a graphite surface. As described in ref. 27 , it is possible to calculate densities and consequently free energy profiles from MD simulations. Here, we use the free energy profiles for the BF − 4 anion. At the interface with a planar surface, ions tend to organise themselves in a layered structure which can extend up to several molecular diameters away from the surface. The free energy profiles are thus characterised by peaks at various distances from the interface as can be seen in figure 3. We limit ourselves to 11 B NMR but note that we would expect qualitatively similar results for 19 F and 1 H NMR.
To assign free energies to the lattice sites, we simply map the MD values on the different sites, i.e. for each lattice site, we calculate its distance to the pore wall and give it the corresponding energy value from the realistic free energy profile obtained through MD simulations 27 ( figure 3). There are three comments we can make at this point: i) the accuracy of the results obtained through the lattice method will depend on the accuracy of the initial MD simulations, ii) the accuracy of the final results will also depend on the lattice parameter chosen, i.e., the finer the grid, the better the representation, and iii) while we do not explicitly include correlated motions of the ions in the model, the correlation effects that result in the layered structure of the liquid at the interface are accounted for by our model. To investigate the first point, we also calculate the NMR spectra for a flat profile, where the free energy is constant across the pore. In this case, we chose a closest distance of approach of 0.32 nm, in agreement with the MD results 27 . The frequencies, which correspond to the Larmor frequencies of the nuclei in different local environments for a real system, follow from the DFT calculations of ref. 7 . The frequencies on the lattice sites are simply mapped using the DFT results as was done for the free energies (see figure 3). By doing this, we make a number of assumptions. Firstly, we use the chemical shifts calculated with the Gaussian software 28 with the Nucleus Independent Chemical Shift (NICS) approach, i.e. we assume that the chemical shifts originate from ring current effects and that the nature of the ion and its charge do not influence this quantity. Secondly, the calculations reported were not performed on graphene 7 , which is difficult to treat from an ab-initio point of view, but rather on different aromatic molecules such as coronene, circumcoronene and dicircumcoronene. The molecular size has an impact on the NICS values calculated and here, we will show results for these three aromatic molecules. We note here that the chemical shift also depends on the lateral position of the ions over the surface 2,7 . While this is not investigated here, this spatial dependence is also expected to impact the NMR results.
Finally, we need to determine the timestep ∆t corresponding to the lattice parameter a chosen (0.5Å). The timestep would thus be equal to 4.92 10 −13 s following D = a 2 /(2d∆t). We would like to stress here that if we were to simulate a complete NMR spectrum with the parameters used experimentally, it would be very computationally costly. Indeed, the experimental dwell time (the time between two FID points) is chosen depending on the desired spectral width and is typically of the order of 5 µs so we would have to perform a simulation where each data point from the experimental NMR signal corresponds to more than one million timesteps in our model. Nevertheless, we will show that this is not a problem here as the lattice model allows us to investigate slower dynamics, which provides information that allows to understand experimental spectra. In our model, the diffusion coefficient is increased by changing the timestep while keeping all other parameters constant. We note that the simulation time is kept constant for all the simulations, so when the timestep is decreased by a factor K, the total number of steps is also increased by a factor K and the points are sampled to perform the Fourier transformation on the same number of values. All the simulations are long enough to observe the entire decay of the NMR signal.

B. NMR spectra for ultra-slow diffusion and effect of increasing the ion dynamics
We first calculate NMR spectra for 11 B in a 4 nm slit pore for the hypothetical case of ultra-slow diffusion (D init = 84.7 10 −19 m 2 .s −1 , which corresponds to a timestep of ∆t = 4.92 10 −5 s; a = 0.5Å). The NMR spectra obtained with the lattice model using the realistic free energy profile from MD simulations 27 and the flat free energy profile are shown in figure 4. It is clear that assuming a flat energy profile leads to large differences in the resulting NMR spectrum and the disappearance of many features. In the case of the realistic free energy profile, we can identify different peaks corresponding to the different environments experienced by the anions which tend to organise themselves in layers, as can be seen from the free energy profile (see figure 3). The first adsorbed layer gives rise to the peak at most negative chemical shift while ions in the second and third are observed at smaller shifts. The sharp peak close to 0 ppm corresponds to the central region of the slit pore, which is only weakly affected by the graphitic carbon surfaces. In the case of a flat energy profile, we observe a single resonance, with a shoulder on the negative frequency side.
We now investigate the impact of the NICS profiles on the calculated spectra. We focus on three different molecules, namely coronene, circumcoronene and dicircumcoronene which have molecular areas respectively close to 0.362 nm 2 , 0.981 nm 2 and 1.911 nm 2 . The DFT calculations on these molecules show that the larger the molecular size, the larger the chemical shifts, due to the increased number of rings contributing to ring current effects 7 . The NMR spectra obtained with the lattice model using the different NICS profiles (figure 4) are qualita-tively similar but mirror the effect of the DFT calculations in the sense that larger molecules lead to NMR spectra shifted to lower frequencies.
The effect of increasing the diffusion coefficient of the anions on the resulting NMR spectra is explored by working with a realistic free energy profile for the anions and the NICS profile corresponding to circumcoronene.
We start with a diffusion coefficient D init = 84.7 10 −19 m 2 .s −1 and increase it by several orders of magnitude to reach D = 10,000×D init = 84.7 10 −15 m 2 .s −1 (figure 4). When the diffusion is increased by a factor 10 or 100, the peaks corresponding to the different environments broaden and start to merge. With an increase of a factor 1,000, only one peak appears in the spectrum and it becomes sharper as the diffusion is increased further. As the real diffusion coefficient is much larger than the ones studied here (by a factor of 10 4 ), we should expect that a single very sharp peak would be observed experimentally for a porous car-bon consisting of slit pores with a monodisperse pore size distribution: on the timescale probed by the NMR experiment, for the organic electrolyte investigated here, the anions would explore the entire slit pore and only an average chemical shift would be observed.
A much broader lineshape is observed experimentally than the one predicted by the simple slit-pore model 4,7,8 , and the experimental lineshape depends significantly on the temperature 30 . The current results demonstrate that the experimentally observed lineshape is not related to diffusion in a pore but to the fact that ions explore a distribution of pore sizes. Indeed, the experimental carbons usually show some disorder with a distribution of pore sizes and pore geometries. A realistic model would thus have to combine information about the average shift within pores of different sizes and the distribution of exchange rates between pores. With this information, an even more coarse-grained lattice model can be constructed in which each accessible lattice site is a pore. The average chemical shifts corresponding to the various pore sizes can be determined using the current approach and used as an input in this new system.

C. Variation of the average chemical shift as a function of pore width
In this part, we investigate the variation of the average chemical shift as a function of pore width with the following aims: to i) obtain insight in the global trend of the variation and the impact of the input parameters (free energy and NICS profiles) on the resulting data and ii) produce the information necessary for the parametrisation of a new model representing a carbon particle with a realistic pore size distribution. We will thus calculate the average chemical shift for a range of pore widths ranging from 2 to 10 nm for different setups.
We note here that, for the determination of the average chemical shift, we do not need a lattice model calculation. The average chemical shift is simply given by the integrated NICS profile weighted by the free energy profile: This approach was suggested in a previous study by Xing et al. 2 where the authors apply this methodology to obtain insights into 1 H NMR spectra of water in porous carbon materials. While they consider only a flat energy profile for the hydrogen atoms, we show here that the inclusion of a realistic free energy profile leads to some variations in the obtained average chemical shifts. The plots showing the variation of the average chemical shift as a function of pore width for the various conditions investigated are given in figure 4. All the curves are qualitatively similar showing that the absolute value of the chemical shift increases exponentially with decreasing pore width. As a general result, the graphitic domain size has a larger influence on the results than the use of a realistic or a flat free energy profile. Nevertheless, the influence of the free energy profile should not be neglected as it can lead to relative errors of up to 25 % for the flat profile compared to the realistic profile. We note that the use of a flat energy profile, corresponding to a constant density across the pore, leads to an overestimation of the average chemical shifts, i.e. the values are shifted to larger frequencies compared to the realistic free energy profile values. This is a consequence of the neglect of high anionic density close to the surface that exists in the layered structure.  27 . Secondly, for mesopores, the pore walls should appear microscopically as being close to planar surfaces so that the approximation of pores by slit pores might be acceptable. We see on figure 4 that this experimental result seems to be in good agreement with the trend obtained with the lattice model and the circumcoronene molecule. While this agreement is partly fortuitous as there is no obvious reason at this point suggesting that the circumcoronene molecule should be a better model than the other aromatic molecules, this is very promising for further applications of the lattice model. While Borchardt et al. 8 and other authors 4,5 report NMR results for other carbon materials with average pore sizes in the range of 1 nm, these were not investigated here as the current model for slit pores is parameterised through MD simulations of liquid inside a mesopore 27 .

V. NMR SPECTRA CALCULATION FOR A POROUS CARBON PARTICLE
A. Parametrisation of the lattice model for a carbon particle In the final part of this article, we use an even more coarse-grained lattice model to represent a carbon particle with a realistic pore size distribution and the exchange of ions between pores of different sizes. In this new setup, each lattice site represents an entire pore, with a given pore size. We apply this model to study how the spatial distribution of pore sizes in a carbon particle affects the resulting NMR spectra. In what follows, we will consider a model that resemble the experimental system of Bor-chardt et al. 8 . As before, we have to assign a number of quantities to each lattice site before performing the simulations, namely in this model, i) the pore size, ii) the free energy, iii) the chemical shift and iv) the energy barriers associated with the exchange between each lattice site and its neighbours.
To parametrise the model, we first need to choose a pore size distribution. In their work, Borchardt et al. 8 report an experimental pore size distribution obtained through nitrogen adsorption. The reported experimental pore size distribution presents two maxima: one close to 1 nm and the other close to 4.5 nm. As most of the pore volume in this material corresponds to pores larger than 2 nm (93 % of the total pore volume), we only consider this part of the curve in our model. We fit a log normal pore size distribution to the experimental curve. The mean of this log normal distribution is 1.558 while the standard deviation is 0.193. As shown in figure 5, the fit is quite good and the log normal distribution seems to be a good choice to represent the experimental pore size distribution.
The free energy associated to each site depends on the amount of liquid adsorbed at this site and hence on the pore width (R), the pore surface (S), and the free energy profile inside the pore. As before, we will assume that all the pores are slit pores. For each lattice site, we choose randomly a pore width following the log normal distribution described above. We then need to choose a pore surface to completely define or encapsulate the pore volume associated with the lattice site considered. We choose to work with a distribution of pore surfaces centered around the area of a circumcoronene molecule, since the experimental result 8 was close to the simulation results for this molecule, and intermediate between the coronene and dicircumcoronene areas for which we also have chemical shifts values. The distribution is arbitrarily chosen to be log normal with a mean of -0.1 and a standard deviation of 0.25. Once the pore width and the pore surface have been determined, the free energy of the lattice site is defined so that the density inside the pore, ρ i , is proportional to the integrated density inside the slit pore multiplied by the pore surface: In this way, we take into account both the volume of the pore and the liquid structure inside it to define the free energy of the lattice sites. The next step is to assign a chemical shift to each lattice site. This is again done according to the pore width and pore surface previously chosen. The value of the chemical shift for a lattice site is intrapolated from the values calculated for the three aromatic molecules. For example, for a pore width R, and a pore surface S, intermediate between the coronene area and the circumcoronene area, the chemical shift is chosen as a weighted average of the coronene value at R, δ coron (R) and of the circumcoronene value at R, δ circum (R): where S coron , S circum and S dicircum are taken to be respectively 0.362 nm 2 , 0.981 nm 2 and 1.911 nm 2 . In our lattice model representing a carbon particle, it is not particularly relevant to set a lattice parameter as we do not have a clear idea of what would be a realistic distance between pores of different sizes and we do not know how the diffusion coefficients are affected by the confinement. In contrast, it is essential to be able to vary the exchange rate between pores, as this rate will affect the temperature dependence of the NMR lineshape. We will assume that the presence of (positive) energy barriers Ea(ij) reduce the forward and backward jump rate between sites i and j. Compared to free diffusion, the probability of jumping from a site i to a site j is then reduced by a factor α(ij) equal to: Note that, as mentioned in section II, the probability of jumping from j to i has to be reduced by the same factor α(ij) in order to maintain detailed balance. We do not know the inter-pore jump rates a priori. In what follows, we will assume that the energy barriers obey a Gaussian distribution with a mean value E mean and a standard deviation equal to 0.2, irrespective of their position in the lattice. Following this parametrisation scheme, we build three different models of carbon particles with the same pore size distribution. We use cubic lattices with a dimension of 20×20×20. To check a possible size dependence of the resulting NMR spectra, we also performed the calculations at room temperature with lattices of dimensions 30×30×30 and 40×40×40. The obtained spectra are very similar for the different lattice sizes so that we consider only the smaller lattice in the rest of this article. All calculations are performed with a time step equal to 5 µs and the simulations are run for 50,000 time steps (0.25 s), which is much longer than the decay time of the NMR signal.
We first build a model where the various pore sizes, and consequently the chemical shifts associated, are distributed randomly over the lattice. This situation is represented in figure 5b). For this model, we choose the mean value of the energy barriers E mean such that the calculated linewidth (full width at half maximum peak intensity) at room temperature is equal to the experimental value of 0.45 ppm obtained by Borchardt et al. 8 . This results in a value for E mean equal to 11.8 kJ.mol −1 . Note that the calculation is performed with a spectrometer frequency of 300 MHz as the experiments were done with this condition. This first model will be referred to as Random-11.8. To investigate the effect of the spatial distribution of shifts, we then create a setup where the pore sizes are organised in an ordered way such that there is a gradient of chemical shifts in one direction of the lattice as illustrated in figure 5c). The second model, referred to as Grad-11.8, is a lattice where the chemical shifts are ordered following a gradient and the mean energy barrier is E mean = 11.8 kJ.mol −1 . This model allows us to probe the effect of the gradient on the NMR lineshape while keeping all other parameters constant. Finally, the third model, referred to as Grad-1.7 is a lattice where the chemical shifts are ordered following a gradient and E mean is set equal to 1.7 kJ.mol −1 to recover a linewidth of 0.45 ppm at room temperature. The comparison between Random-11.8 and Grad-1.7 allows us to probe the different behaviours of the random and gradient spatial distributions for two systems giving the same linewidth at room temperature.

B. Investigation of various parameters on the NMR spectra of the model carbon particles
We first calculate the NMR spectra corresponding to the three models Random-11.8, Grad-11.8 and Grad-1.7 at room temperature (6). Comparing the results for Random-11.8 and Grad-11.8, which differ only by the spatial distribution of the shifts, we note that they yield very different NMR spectra. The spectrum is much broader in the case of the gradient than in the case of the random distribution. This is not surprising because, in the gradient model it takes longer to sample all Larmor frequencies than in the random model, where a single jump can dramatically change the Larmor frequency. The NMR spectrum for Grad-11.8 mirrors the asymmetrical distribution of chemical shifts resulting from the log normal distribution of pore sizes. By decreasing the energy barrier from 11.8 to 1.7 kJ.mol −1 in the gradient configuration, we recover the NMR spectrum observed for Random-11.8.
The difference in energy barrier for Random-11.8 and Grad-1.7 leads to different conclusions concerning the exchange rates and diffusion times in the pores. The relatively long exchange times necessary to observe lineshape variations on the experimental time scales could originate from different factors: a slow diffusion in the carbon particle due to confinement or a slow exchange between pores due to either small interpore connections or large distances between pores. From the simulations, we can estimate the exchange rate between pores needed to reproduce the experimental linewidth of 0.45 ppm 8 . For Random-11.8, the exchange rate is equal to 1.75 kHz and corresponds to a correlation time of 570 µs. For Grad-1.7, the exchange rate is much higher, equal to 101 kHz, and corresponds to a correlation time of 9.9 µs. If we assume that the diffusion coefficient in the porosity is the same as the bulk diffusion coefficient 29 , i.e. equal to 84.7 10 −11 m 2 .s −1 , and that the long exchange times are due to large distances between pores, we can estimate the distance between pores. In the case of Random-11.8 and Grad-1.7, we obtain distances between pores of different sizes respectively equal to 1,702 nm and 224 nm. These values seem very high considering that the particle sizes are usually in the micrometer range. If we assume that the diffusion coefficient is two orders of magnitude smaller in the particle compared to the bulk, we estimate distances between pores of around 170 nm for Random-11.8 and 22.4 nm for Grad-1.7. While the distance still seems quite large for the random distribution as the average pore size, 4.5 nm, is much smaller than 170 nm, the estimation for Grad-1.7 seems realistic. The point here is that, while we cannot draw unambiguous conclusions about the origin of the slow exchange, we can see that the estimations are strongly affected by the spatial distribution of shifts.
We now explore the effect of temperature applied magnetic field strength on the NMR lineshapes, to determine if further information on porosity can be extracted (fig-FIG. 6. NMR spectra and linewidths for different setups. a) NMR spectra obtained at room temperature for random/gradient distributions of shifts and barrier heights distributions centered around Emean = 11.8 or 1.7 kJ.mol −1 . The gradient of shifts (dashed red line) leads to a much broader spectrum than the random distribution (black solid line). In the gradient case, when the barrier heights distribution is changed to reproduce the experimental linewidth at room temperature (45 ppm 8 ), the spectrum (blue crosses) is superimposed with the one calculated for the random case. b) Linewidth as a function of temperature for the different setups. c) Linewidth as a function of spectrometer frequency. ure 6). Different magnetic fields are usually labelled according to the Larmor frequency of 1 H in that field: we will refer to this quantity as the spectrometer frequency. We first focus on the temperature effect. For all the systems, when the temperature increases, the diffusion of the ions, or equivalently the exchange between different pores, increases leading to a sharpening of the spectra and thus a decrease of the linewidth. Nevertheless, the temperature dependence is not the same for the different systems. The two gradient cases show a slower decrease of the linewidth with temperature compared to the random case. Moreover, the gradient curves appear to be almost linear over this temperature range while for the random case the temperature dependence of the linewidth appears roughly exponential. Over a wider temperature range, all curves should show an exponential behaviour as a result of the exponential dependence of the transition probability on temperature. The interesting point here is that, while the Random-11.8 and Grad-1.7 have the same linewidth at room temperature, this linewidth does not vary in the same way for the two systems.
The last parameter that we explore is the applied magnetic field by varying the spectrometer frequency from 300 MHz to 1 GHz. The results show that the dependence of the linewidth on the spectrometer frequency is strikingly different for the three systems. The linewidth is almost constant for Grad-11.8, showing that the distribution of shifts in the spectrum corresponds to the actual distribution of shifts in the system. In constrast, for Random-11.8 and Grad-1.7, the linewidth depends dramatically on the applied magnetic field.
Although not investigated here, we can expect that very different ions could lead to shifts in the NMR spectra as the adsorption profiles will not be the same and larger molecules would not be able to access all the pores accessible to small molecules. In conclusion, analysing a range of NMR experiments using the lattice model presented here, should provide a promising tool for charac-terising porous carbon structures.

VI. CONCLUSION
In this work, we have described a coarse-grained simulation method to predict diffusion coefficients and NMR spectra of probe particles diffusing in a porous carbon matrix. The model that we propose allows us to account for relevant microscopic information whilst exploiting the computational efficiency of the 'moment-propagation' approach, a method that allows us to account for all possible trajectories that particles could follow in a discretised model of a porous network.
Our simulation method is able to reproduce the coalescence effect observed experimentally, the results agree with analytical solutions in the case of a two-site exchange model. The model allows NMR spectra for much more complex systems to be predicted as it can deal with any set of frequencies and spatial distributions of these frequencies. In particular, we show that the model can be parametrised realistically using input from molecular simulations such as free energy profiles, to represent the molecular/ionic adsorption, and nucleus independent chemical shifts estimated through DFT calculations, to account for the spatial dependence of the chemical shifts of the anions within the carbon pores.
Parametrisation was performed for an organic electrolyte confined in slit mesopores of various pore sizes and we could show that, while a number of environments would be observed if the diffusion of probed species was ultra-slow, the exchange rates involved in experiments lead to the detection of a single resonance. This peak is observed for an average chemical shift which depends both on the pore size and on the adsorption profile of the studied species. In this work, we also investigate the effect of the molecular size chosen for the chemical shifts calculations on the resulting NMR spectra and show that the graphitic domain size can lead to large variations of the average chemical shifts. This can be seen as a challenge from the parametrisation point of view but it can also be seen as an opportunity to get insights into these domain sizes from NMR while it is not probed from other characterisation techniques.
The last part of this article focused on the parametrisation of an even more coarse-grained lattice model that can be used to represent a carbon particle with a realistic pore size distribution. In this model, each lattice site represents a pore with a given pore size. The model allows us to explore various spatial distributions of the pore sizes and various conditions, such as applying different temperatures and magnetic fields, which can be related to experimental conditions. While some of the parameters in this model are known from microscopic simulations, others can be estimated by comparing computed and experimental spectra for a range of temperatures and magnetic fields. Such a comparison yields novel insights into the structure of porous carbon materials and the structure of the liquid inside the pores. In particular, this new lattice model is expected to provide new insights into in situ NMR experiments performed on supercapacitors. Moreover, because of its versatility, the lattice model is a powerful tool to investigate a full range of materials, for which NMR parameters can be determined, including battery and fuel cells materials.
where D is the self diffusion coefficient of the particles and t is the time. If particles move by a sequence of jumps ∆x i of magnitude a, then where N is the number of jumps such that N ∆t = t. We can rewrite Eqn. A2 as It is convenient to separate the terms with i = j and i = j: In the second term on the right hand side the terms i, j and j, i contribute equally and hence If N δt is much larger than the decay time of ∆x i ∆x j , we can write: This expression is the discrete analog of the Green-Kubo relation To emphasise this analogy we interpret ∆x i /∆t as the 'velocity' at timestep i. Then we obtain: Note, however, that we do not assume that the above expression is valid in the diffusive regime. In particular, for diffusion in a homogeneous medium, we have (for not too short ∆t, so that v x (1)v x (j) = 0):

Moment-propagation approach
The expression for the diffusion constant given in Eqn. A7 suggests that we should average Eqn. A7 over all trajectories that a particle can follow. As the number of trajectories of a walk of length N is equal to z N , where z is the lattice coordination number, it would seem that