Molecular dynamics lattice gas equilibrium distribution function for Lennard-Jones particles
aa r X i v : . [ phy s i c s . c o m p - ph ] J a n rsta.royalsocietypublishing.org Research
Article submitted to journal
Subject Areas: computational physics,coarse-graining, fluid mechanics
Keywords: molecular dynamics, lattice gasmethod, lattice Boltzmann method,coarse-graining
Author for correspondence:
Alexander J. Wagnere-mail: [email protected]
Molecular dynamics latticegas equilibrium distributionfunction for Lennard-Jonesparticles
Aleksandra Pachalieva , and Alexander J.Wagner Center for Nonlinear Studies, Los Alamos NationalLaboratory, Los Alamos, NM 87545, USA Department of Mechanical Engineering, TechnicalUniversity of Munich, 85748 Garching, Germany Department of Physics, North Dakota StateUniversity, Fargo, ND 58108, USA
The molecular dynamics lattice gas method mapsa molecular dynamics simulation onto a lattice gasusing a coarse-graining procedure. This is a novelfundamental approach to derive the lattice Boltzmannmethod by taking a Boltzmann average over themolecular dynamics lattice gas. A key property ofthe lattice Boltzmann method is the equilibriumdistribution function, which was originally derivedby assuming that the particle displacements inthe molecular dynamics simulation are Boltzmanndistributed. However, we recently discovered that asingle Gaussian distribution function is not sufficientto describe the particle displacements in a broadtransition regime between free particles and particlesundergoing many collisions in one time step. In arecent publication, we proposed a Poisson weightedsum of Gaussians which shows better agreementwith the molecular dynamics data. We derive anlattice Boltzmann equilibrium distribution functionfrom the Poisson weighted sum of Gaussians modeland compare it to a measured equilibrium distributionfunction from molecular dynamics data and toa equilibrium distribution function from a singleGaussian probability distribution function. © The Authors. Published by the Royal Society under the terms of theCreative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author andsource are credited. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
1. Introduction
The molecular dynamics lattice gas (MDLG) method [15, 16] uses a coarse-graining procedureto establish a direct link between microscopic methods – in particular, molecular dynamics(MD) simulation, and mesoscale methods such as lattice gas (LG) [3, 6] and lattice Boltzmannmethods (LBM) [8, 18]. The MDLG fully relies on MD data and as such it rigorously recovers thehydrodynamics of the underlying physical system, and can be used to verify the behavior andexamine the properties of the LG or the LBM methods directly without using the standard kinetictheory approach. Aspects that can be examined include fluctuating [2, 5, 10, 20], thermal [7, 11],multi-phase and multi component systems [4, 8, 12, 19].A key feature in the LBM is the equilibrium distribution function. The LBM equilibriumdistribution was originally derived by analogy to the continuous Boltzmann equation, wherethe equilibrium distribution for the velocities is a Maxwell Boltzmann distribution. Similarly, theLBM moments of the discrete velocity distribution were matched, to the degree possible, with thevelocity moments of the Maxwell Boltzmann distribution. In the alternative derivation of the LBMfrom MD, it was shown that these previously postulated equilibrium distributions are indeed (atleast approximately) consistent with the MDLG approach for specific discretization combinationsfor lattice and time spacing.In the original MDLG calculation of the equilibrium distribution by Parsa et al. [15] it wasassumed that the particle displacements in the molecular dynamics simulation are also Boltzmanndistributed. This assumption gave an adequate prediction of the global equilibrium distributionfunction of the lattice Boltzmann method. However, later on by examining a non-equilibriumsystem undergoing a shear flow, we noticed small deviations (up to 5%) from the analyticallypredicted equilibrium distribution function. These deviations were traced back to the predictionof the one-particle displacement distribution function. In Pachalieva et al. [13], we proposed acorrection of the displacement distribution function, which shows that a dilute gas with areafraction of φ = 0 . and temperature of 20 LJ is better approximated by a Poisson weighted sumof Gaussians (WSG) probability distribution function. This remains true for a purely ballistic andpurely diffusive regimes (for very small or very large time steps respectively), where the PoissonWSG formulation is reduced to a single Gaussian. In the current publication, we derive the MDLGequilibrium distribution function from the Poisson WSG one-particle displacement function andcompare it to a measured equilibrium distribution function from molecular dynamics (MD)simulation and to the single Gaussian equilibrium distribution function. Our findings showthat the Poisson WSG approximates the measured equilibrium distribution function significantlybetter.The rest of the paper is summarized as follows: In Section 2, we briefly describe theMDLG method and how to derive the equilibrium distribution assuming that the displacementdistribution is given by a single Gaussian. We then show how the equilibrium distributionfunction is altered if the displacements are instead distributed according to a Poisson WSGone-particle displacement function. This is followed by a detailed description of the MDsimulation setup used to obtain the MD data given in Section 3. In Section 4, we compare theequilibrium distribution function obtained from the previously suggested single Gaussian, thenovel Poisson weighted sum of Gaussians and the one measured from MD data. Our analysisshows significant improvement of the equilibrium distribution function analytical predictionwhen the Poisson WSG model is used. Finally, we give a brief conclusion and suggestions forfuture work in Section 5.
2. Derivation of the MDLG equilibrium distribution function
In the MDLG analysis, we impose a lattice onto a MD simulation and track the migration ofthe particles from one lattice position to another with displacement v i after a time step ∆t . Aschematic representation of the lattice is given in Fig. 2b where the numbers 0 to 24 represent the r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. i index of the occupation number of a D2Q25 velocity set. The MDLG occupation number is thusgiven by n i ( x, t ) = X j ∆ x [ x j ( t )] ∆ x − v i [ x j ( t − ∆t )] , (2.1)with the delta function ∆ x [ x j ( t )] = 1 , if particle x is in the lattice cell at time t , and ∆ x [ x j ( t )] =0 , otherwise. Here, the x j ( t ) is the position of the j -th particle at time t . The MDLG evolutionequation is derived as n i ( x + v i , t + ∆t ) = n i ( x, t ) + Ξ i , (2.2)where the lattice gas collision operator Ξ i given by Ξ i = n i ( x + v i , t + ∆t ) − n i ( x, t ) . (2.3)The molecular dynamics lattice Boltzmann (MDLB) distribution function is defined as aBoltzmann ensemble average of the MDLG occupation number n i and it is given by f i = h n i i neq . (2.4)By taking the non-equilibrium ensemble average of Eq. (2.2), we obtain the MDLB evolutionequation f i ( x + v i , t + ∆t ) = f i ( x, t ) + Ω i , with Ω i = h Ξ i i neq , (2.5)where Ω i is the MDLB collision operator. A key element of the LBM is the global equilibriumdistribution function, which in MDLB context is defined as an average of the lattice gas densities n i over the whole MD domain and all iterations of an equilibrium MD simulation. The MDLBequilibrium distribution function writes f eq i = h n i i eq = *X j ∆ x [ x j ( t )] ∆ x − v i [ x j ( t − ∆t )] + eq = M Z dx Z dδx P (1) , eq ( x , δx ) ∆ x [ x ] ∆ x − v i [ x − δx ] , (2.6)where M is the total number of particles and P (1) , eq is the one-particle displacement distributionfunction in equilibrium. This allows us to obtain the equilibrium distribution function f eq i analytically from the one-particle displacements Probability Distribution Function (PDF). In theMDLB formulation, the equilibrium distribution function depends solely on the one-particledisplacement distribution function. Thus, knowing P (1) , eq is crucial in predicting the equilibriumdistribution function. (a) Single Gaussian distribution model Initially, Parsa et al. [15] suggested that the one-particle distribution function is given by a singleGussian in one-dimension ( d = 1 ) P G α ( δx ) = 1[2 π h ( δx α ) i ] d/ exp (cid:20) − ( δx α − u α ∆t ) h ( δx α ) i (cid:21) , (2.7)with displacements δx α and second order moment h ( δx α ) i . The solution factorizes for higherdimensions and it is given by P G ( δx ) = d Y α =1 P G α ( δx ) . (2.8)From Eq. (2.6) the equilibrium distribution function can be written as f eq , G i ρ eq = d Y α =1 f eq , G i,α , (2.9) r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. where f eq , G i,α = N e − ( ui,α − a − e − u i,α a + e − ( ui,α +1)22 a ! + u i,α − (cid:20) erf (cid:18) u i,α − a √ (cid:19) − erf (cid:18) u i,α a √ (cid:19)(cid:21) + u i,α + 12 (cid:20) erf (cid:18) u i,α + 1 a √ (cid:19) − erf (cid:18) u i,α a √ (cid:19)(cid:21) , (2.10)with a = h ( δx α ) i ( ∆x ) , N = a √ π , u i,α = v i,α − u α . (2.11)For details regarding the derivation of the Gaussian equilibrium distribution function, please referto [15].Even though this formulation shows very good agreement with the measured equilibriumdistribution function from MD simulations, under more careful investigation we found that thethere are discrepancies of about 5%. This means that the displacement distribution functioncannot be fully captured by a single Gaussian and a more complex distribution function has to beapplied. (b) Poisson weighted sum of Gaussians model In Pachalieva et al. [13], we have introduced a correction of the displacements PDF proposed byParsa et al. [15] using a Poisson weighted sum of Gaussians (WSG) instead of a single Gaussiandistribution function. The Poisson WSG is given by P WSG ( δx ) = ∞ X c =0 e − λ λ c c ! P c ( δx ) , (2.12)where the P c ( δx ) probability distribution function factorizes for higher dimensions as P c ( δx ) = d Y α =1 P cα ( δx ) , (2.13)with α being a Cartesian coordinate index. P cα ( δx ) is given by P cα ( δx ) = (cid:20) ( λ + 1)2 π ( c + 1) h ( δx α ) i (cid:21) d/ exp (cid:20) − ( λ + 1)( δx α − u α ∆t ) c + 1) h ( δx α ) i (cid:21) , (2.14)where δx α is the displacement in one-dimension, h ( δx α ) i is the second-order moment, u α isthe velocity, c is the number of occurrences, and λ is the average number of collisions. Similar toEq. (2.9), the equilibrium distribution function can be expressed as f c, eq i ρ eq = d Y α =1 f c, eq i,α , (2.15)where ρ eq is the mass density and f c, eq i,α in one-dimension is written as f c, eq i,α = ( N c √ π e − ( ui,α − N c − e − u i,αN c + e − ( ui,α +1)2 N c ! + ( u i,α − (cid:20) erf (cid:18) ( u i,α − N c (cid:19) − erf (cid:18) u i,α N c (cid:19)(cid:21) + ( u i,α + 1)2 (cid:20) erf (cid:18) ( u i,α + 1) N c (cid:19) − erf (cid:18) u i,α N c (cid:19)(cid:21)(cid:27) , (2.16) r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. with N c = r a ( c + 1) λ + 1 , a = h ( δx α ) i ( ∆x ) , u i,α = v i,α − u α , (2.17)where the average number of collisions λ could be approximated using the velocityauto-correlation function. However, to reduce the error coming from the theoretical averagenumber of collisions and to eliminate the second-order and the fourth-order moment errors, wematch these moments to the corresponding ones measured directly from the MD simulations.Thus, we obtain the following two solutions for the average number of collisions λ , = − µ ± q µ − µ µ ] + 2 µ µ − µ ] , (2.18)where µ = h ( δx α ) i and µ are the second- and fourth-order displacement moments,respectively. In Pachalieva et al. [13], we show that λ provides an optimal solution, which weuse to derive the Poisson WSG equilibrium distribution function. f eq , WSG i = ∞ X c =0 e − λ λ c c ! f c, eq i . (2.19)To construct the Poisson WSG, we measured the second-order and fourth-order displacementmoments from the MD simulations. For detailed derivation and discussion of the Poisson WSGdisplacement distribution function, please refer to Pachalieva et al. [13]. i / ∆ x00.0050.010.0150.020.025 P ( X i ) P MD (X i )P G (X i )P WSG λ (X i ) (a) i / ∆ x-0.0004-0.000200.00020.00040.00060.0008 K ( X i ) P MD log ( P MD / P G ) P WSG λ log ( P WSG λ2 / P G ) (b) Figure 1. (Color online) (a) Displacements probability distribution functions. The symbols (black)depict a PDF obtained from an MD simulation of LJ particles in equilibrium. The line (red)illustrates a Gaussian probability distribution function defined in Eq. (2.7) with mean-squareddisplacement fitted directly to the MD data. The dashed line (blue) represents the Poisson WSGobtained from Eq. (2.12). Only the data for positive velocities has been depicted due to symmetry.(b) shows the difference between the distributions per interval X i as defined in Eq. (2.20). Thepresented data is for the standard parameters used in the paper and a coarse-grained time step ∆t = 3 . .In Pachalieva et al. [13], we looked at the probability distribution function of the displacementsand we measured precisely the deviations between the theoretical PDFs and the MD simulationdata. Since these deviations are relatively small, we used the following function to quantify the r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. discrepancies K ( X i ) = K ( R k Q ) = R ( X i ) log (cid:18) R ( X i ) Q ( X i ) (cid:19) , (2.20)where R ( X i ) and Q ( X i ) are probability distributions over an interval X i . By performing a sumover all the bins X i , we obtain the well known Kullback-Leibler (KL) divergence [9] defined as D KL ( R k Q ) = X i R ( X i ) log (cid:18) R ( X i ) Q ( X i ) (cid:19) . (2.21)The KL divergence measures the discrepancies of one probability distribution function toanother. It is always non-negative D KL ( R k Q ) ≥ or equal to zero if and only if the probabilitydistribution functions are identical R ( X i ) = Q ( X i ) [9].In Fig. 1a, we see the true probability distribution function obtained from the MD data P MD ( X i ) , the Gaussian probability distribution function P G ( X i ) , and the Poisson WSGdistribution function P WSG ( X i ) . There is a visible divergence between the Gaussian and theother two distribution functions. We measured K ( X i ) , as defined in Eq. (2.20), for P G ( X i ) compared to the MD data and the Poisson WSG distribution function as shown in Fig. 1b. Theresults suggest that even though the Gaussian and the Poisson WSG probability distributionfunctions have the same second moment, their deviations in the higher moments influencestrongly the form of the distribution function. In Section 4, we show how these deviations effectthe LBM equilibrium distribution function.
3. Simulations setup
The measured equilibrium distribution function depicted as f eq , MD i in Figs. 2a- 4a is obtainedfrom molecular dynamics simulations. To perform the MD simulations we used the open-sourcemolecular dynamics framework LAMMPS [17, noa] developed by Sandia National Laboratories.The MD simulations consist of particles interacting with the standard 6-12 Lennard-Jones (LJ)intermolecular potential given by V LJ = 4 ε (cid:20)(cid:16) σr (cid:17) − (cid:16) σr (cid:17) (cid:21) , (3.1)with σ being the distance at which the inter-particle potential goes to zero, r is the distancebetween two particles, and ε is the potential well depth. The particle mass and the LJ particlediameter are set to m = 1 and σ = 1 , respectively. The number of particles in each simulation hasbeen fixed to N = 99 856 which fills a two-dimensional (2D) square with length L = 1000 LJ units.The area fraction φ of the domain is calculated from the area of the circular LJ particles multipliedby the number of particles divided by the area of the domain, where the radius of the circular LJparticle is given by σ/ . The MD simulations considered in this publication have an area fractionof φ = 0 . . We initialised the simulations using homogeneously distributed particles withkinetic energy corresponding to 20 in LJ units. We focus our attention to MD simulations of afairly dilute gas in equilibrium, since the assumption that the collision times is Poisson distributedis correct only for dilute systems.The LJ timescale is given by the time needed for a particle with kinetic energy of half thepotential energy well ε to traverse one diameter σ of an LJ particle. This can be also expressed as τ LJ = r mσ ε . (3.2)The MD step size is set to . τ LJ which is considerably small to ensure high accuracy of theMD data. We define a dimensionless coarse-grained time step ∆t being the product of the MDstep size and the MD output frequency shown in Table 1. The time step ∆t is chosen such thatthe MD simulations are restricted to the ratio of the mean-squared displacement and the squared r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. Table 1.
MD simulation setup MD step MD output Number of Total MD Total MD ∆t ∆x lx size frequency iterations time steps time( τ LJ ) ( /τ LJ ) ( τ LJ )0.3911 4 250 0.0001 3 911 2 000 7 822 000 782.20.5000 5 200 0.0001 5 000 2 000 10 000 000 1 000.00.5626 5.5 180 0.0001 5 626 2 000 11 252 000 1 125.20.6927 6.6(6) 150 0.0001 6 927 2 000 13 854 000 1 385.40.9009 8.3(3) 120 0.0001 9 009 2 000 18 018 000 1 801.81.1261 10 100 0.0001 11 261 2 000 22 522 000 2 252.21.4994 12.5 80 0.0001 14 994 2 000 29 988 000 2 998.81.6342 13.3(3) 75 0.0001 16 342 2 000 32 684 000 3 268.42.0338 15.625 64 0.0001 20 338 2 000 40 676 000 4 067.62.9280 20 50 0.0001 29 280 2 000 58 560 000 5 856.04.1821 25 40 0.0001 41 821 2 000 83 642 000 8 364.26.1751 31.25 32 0.0001 61 751 2 000 123 502 000 12 350.2lattice size being set to a = h ( δx ) i ( ∆x ) ≈ . , (3.3)which corresponds to the parameter a used in earlier publications [14, 15]. By fixing the valueof a to approximately / , we ensure that most of the LJ particles in equilibrium will travel upto one lattice space which corresponds to an D2Q9 lattice Boltzmann method. To verify that thePoisson WSG equilibrium distribution function f eq , WSG i approximated the MD data better thanthe single Gaussian equilibrium distribution function f eq , G i across the length scale, from ballisticto diffusive regime, we vary the coarse-grained time step ∆t ∈ [0 . , . and the lattice size ∆x ∈ [4 , . of the executed simulations. An overview of the MD simulation setup is given inTable 1. The number of lattice points lx varies from 250 to 32 depending on the coarse-grained timestep ∆t . For each coarse-grained time step ∆t we performed iterations which correspondsto total MD time of . τ LJ to
12 350 . τ LJ for the smallest and largest coarse-grained time step ∆t , respectively. In order to bring the molecular dynamics simulations to equilibrium state beforewe start collecting data, the initial 3 000 000 iterations of each simulation were discarded. Thediscarded iterations are not included in Table 1 for clarity.The MD simulation setup characterizes a hot dilute gas in equilibrium with average velocityfixed to zero Nu α = N X j =1 v j,α = 0 , (3.4)where N is the number of LJ particles.We performed standard molecular dynamics simulations without thermostat. In the LAMMPSframework this is called NVE integration. The microcanonical ensamble NVE is characterized byconstant number of particles (N), constant volume (V) and constant energy (E).
4. Results
In order to obtain a measured equilibrium distribution function, we post-process the collectedMD data using the MDLG analysis tool. The MD domain is overlapped with a lattice and wetrace the migration of the particles over time from one lattice to another. By doing this, we obtainthe MDLG occupation numbers n i ( x, t ) as defined in Eq. (2.1) which after sufficient averagingdeliver the MDLB equilibrium distribution function f eq , MD i as defined in Eq (2.6). r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. The analytical models of the equilibrium distribution function defined in Section 2 dependonly on the choice of the one-particle displacement distribution function. Since we definetwo different one-particle distribution, we expect to see also changes in the respectiveequilibrium distribution function derived from them, even though their second-order momentsare equivalent. However, a non-trivial question remains how the migration of particles from onenode to another changes within a lattice.To gain a better understanding, we calculate the equilibrium distribution function for anextended D2Q25 lattice which corresponds to two neighboring cells in X - and Y -directionsfor a two-dimensional domain. A schematic representation of the D2Q25 lattice is given inFig. 2b. In equilibrium state with zero initial velocity, one distinguishes six sets of equilibriumdistribution function contributions: f eq , ∗ , f eq , ∗ − , f eq , ∗ − , f eq , ∗ − , f eq , ∗ − , and f eq , ∗ − , where each sethas a unique displacement length from the central lattice. When measuring the equilibriumdistribution function f eq , MD i from the MD simulations, we average over the number of latticesfor each set to obtain a symmetric probability distribution function. It is worth mentioning thatthe deviations of the f eq , MD i values within each set are relatively small.The MDLG analysis was introduced for an D2Q49 lattice including a third layer ofneighbouring cells, however, the number of considered neighboring layers depends solely onthe problem at hand. For a simulation in equilibrium with zero velocity, and a parameter a as defined in Eq. (3.3) being set to approximately . , we obtain an equilibrium distributionfunction which is symmetric and has significant contributions up to D2Q25 lattice nodes. ∆ t1e-071e-061e-050.00010.0010.010.11 f i e q , * f eq, MD f eq, MD f eq, MD f eq, MD f eq, MD f eq, MD f eq, G f eq, G f eq, G f eq, G f eq, G f eq, G f eq, WSG f eq, WSG f eq, WSG f eq, WSG f eq, WSG f eq, WSG (a) (b) Figure 2. (Color online) (a) Estimated equilibrium distribution functions f eq , ∗ i obtained either fromMD simulation data ( f eq , MD i ), theoretical solution using a single Gaussian distribution function( f eq , G i ) or theoretical solution using Poisson WSG ( f eq , WSG i ). (b) Schematic representation of theD2Q25 lattice. The equilibrium distribution function f eq , ∗ i values are color coded and each colorrepresents one of the six sets of equilibrium distribution function contributions. Here, the asterisk( ∗ ) corresponds to the variety of methods used to obtain an equilibrium distribution function:measured from MD simulation, single Gaussian analytical solution and Poisson WSG analyticalsolution.The estimated equilibrium distribution function f eq , ∗ i for a variety of coarse-grained timesteps ∆t ∈ [0 . , . is shown in Fig. 2a. The equilibrium distribution function f eq , ∗ i , as r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. mentioned above, is obtained from three different methods: f eq , MD i is measured from an MDsimulation, f eq , G i is theoretically estimated using a single Gaussian distribution function and f eq , WSG i is theoretically estimated from a Poisson WSG distribution function. The theoreticalequilibrium distribution function models are described in detail in Sections 2(a) and 2(b),respectively.In Fig. 2a one can see that the largest equilibrium distribution function contributions arecoming from the first layer neighbours f eq , ∗ − . These nodes are approximated very well byboth theoretical models, please refer to Fig. 3a for a detailed comparison of the measured andthe theoretical f eq , ∗ − . The next equilibrium distribution function groups f eq , ∗ − and f eq , ∗ − aresignificantly smaller than f eq , ∗ − with one to two order of magnitude. For f eq , ∗ − and f eq , ∗ − , wesee that the deviations of the measured and the theoretical single Gaussian model become larger.The Poisson WSG f eq , ∗ − show a very good agreement with the measured equilibrium distributionfunction. The diagonal nodes in the second layer f eq , ∗ − are even smaller and their value could beconsidered negligible. However, the measured equilibrium distribution function f eq , MD i showsa good agreement with the theoretical Poisson WSG f eq , WSG i even for very small contributionssuch as f eq , ∗ − . This suggests that these contributions even though really small are not just noisebut theoretically justified. ∆ t0.980.9911.011.02 f i e q , * / f i e q , G f i eq, G / f i eq, G f eq, MD / f eq, G f eq, MD / f eq, G f eq, MD / f eq, G f eq, WSG / f eq, G f eq, WSG / f eq, G f eq, WSG / f eq, G (a) (b) Figure 3. (Color online) (a) First layer equilibrium distribution functions f eq , ∗ − scaled to theGaussian equilibrium distribution function. The equilibrium distribution functions are obtainedeither from MD simulation data ( f eq , MD0 − ), theoretical solution using a single Gaussiandistribution function ( f eq , G0 − ) or theoretical solution using Poisson WSG ( f eq , WSG0 − ). (b) Schematicrepresentation of the D2Q25 lattice. The equilibrium distribution function f eq , ∗ i values arecolor coded and each color represents one of the six sets of equilibrium distribution functioncontributions. Here, the asterisk ( ∗ ) corresponds to the variety of methods used to obtain anequilibrium distribution function: measured from MD simulation, single Gaussian analyticalsolution and Poisson WSG analytical solution.Figures 3a and 4a depict the equilibrium distribution functions scaled to the single Gaussianequilibrium function. They show how the measured from MD simulation and the novel PoissonWSG equilibrium distribution functions deviate from the single Gaussian. The first layer r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. equilibrium distribution function values are shown in Fig. 3a. These nodes have the largestcontribution to the total equilibrium distribution function.Fig. 3a shows more particles staying at node zero and a depression for the first neighbouringlayer (nodes 1 to 8). This very same feature repeats itself in Fig. 1b. The P WSG λ log(P WSG λ / P G ) values depicted in blue, show that the number of small displacements is enhanced X i /∆x ∈ [0 , . while the number of X i /∆x ∈ [0 . , . displacements is suppressed. ∆ t1234 f i e q , * / f i e q , G f i eq, G / f i eq, G f eq, MD / f eq, G f eq, MD / f eq, G f eq, MD / f eq, G f eq, WSG / f eq, G f eq, WSG / f eq, G f eq, WSG / f eq, G (a) (b) Figure 4. (Color online) (a) Second layer equilibrium distribution functions f eq , ∗ − scaled tothe Gaussian equilibrium distribution function. The equilibrium distribution functions areobtained either from MD simulation data ( f eq , MD9 − ), theoretical solution using a single Gaussiandistribution function ( f eq , G9 − ) or theoretical solution using Poisson WSG ( f eq , WSG9 − ). (b) Schematicrepresentation of the D2Q25 lattice. The equilibrium distribution function f eq , ∗ i values arecolor coded and each color represents one of the six sets of equilibrium distribution functioncontributions. Here, the asterisk ( ∗ ) corresponds to the variety of methods used to obtain anequilibrium distribution function: measured from MD simulation, single Gaussian analyticalsolution and Poisson WSG analytical solution.The second layer equilibrium distribution function values are depicted in Fig. 4a. In Fig. 1b wenotice the enhanced probability of large displacements X i /∆x ∈ [0 . , . which corresponds tothe larger values of f eq , WSG9 − in Fig. 4a. The deviations (up to approx. 4.5%) from the theoreticalsingle Gaussian equilibrium distribution function are also larger compared to the first layer nodes f eq , WSG0 − . Since the f eq , WSG9 − true values are smaller by multiple orders of magnitude than thefirst layer neighbours f eq , WSG0 − these deviations have a smaller impact on the total equilibriumdistribution function, even though they are larger. Nevertheless, Fig. 4a shows clearly that thePoisson WSG equilibrium distribution function captures better the MD data.
5. Conclusion
In the pursuit of deriving a universal LBM collision operator, we investigated a non-equilibriumsystem where it is assumed that collisions keep the distribution functions close to localequilibrium, constrained by the conserved quantities, usually mass and momentum. This means r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A .................................................................. that even small errors in the formulation of the equilibrium distribution function influences thedata obtained for the collision operator.A key ingredient to predict the MDLB equilibrium distribution function correctly is thedefinition of the one-particle displacement distribution function. It is clear that even though thesecond-order moment of the single Gaussian and the Poisson WSG distribution function are thesame, the derived equilibrium distribution functions show non-trivial behaviour which correlateswith the raw probability distribution function. In conclusion, the Poisson WSG equilibriumdistribution function captures better the measured equilibrium distribution function from MDdata.Data Accessibility. This manuscript has no further supporting data.
Authors’ Contributions.
AW supervised the research, contributed to it, and revised the manuscript. APcontributed to the research, embedded the proposed model, set up the test cases, performed the data analysisand wrote the manuscript. All authors read and approved the manuscript.
Competing Interests.
The authors declare that they have no competing interests.
Funding.
AP is partially supported by the Center for Nonlinear Studies (CNLS) and the Laboratory DirectedResearch and Development (LDRD) program at Los Alamos National Laboratory (LANL), and the GermanFederal Ministry of Education and Research (BMBF) in the scope of the project Aerotherm (reference numbers:01IS16016A-B).
References noa LAMMPS Official Website: http://lammps.sandia.gov.2 Adhikari, R., Stratford, K., Cates, M. E., and Wagner, A. J. (2005). Fluctuating lattice boltzmann.
Europhysics Letters (EPL) , 71(3):473–479.3 Blommel, T. and Wagner, A. J. (2018). Integer lattice gas with Monte Carlo collision operatorrecovers the lattice Boltzmann method with Poisson-distributed fluctuations.
Physical Review E ,97(2):023310.4 Briant, A., Wagner, A., and Yeomans, J. (2004). Lattice boltzmann simulations of contact linemotion. i. liquid-gas systems.
Physical Review E , 69(3):031602.5 Dünweg, B., Schiller, U. D., and Ladd, A. J. C. (2007). Statistical mechanics of the fluctuatinglattice boltzmann equation.
Phys. Rev. E , 76:036704.6 Frisch, U., Hasslacher, B., and Pomeau, Y. (1986). Lattice-Gas Automata for the Navier-StokesEquation.
Physical Review Letters , 56(14):1505–1508.7 He, X., Chen, S., and Doolen, G. D. (1998). A novel thermal model for the lattice boltzmannmethod in incompressible limit.
Journal of computational physics , 146(1):282–300.8 He, X. and Luo, L.-S. (1997). Theory of the lattice boltzmann method: From the boltzmannequation to the lattice boltzmann equation.
Physical Review E , 56(6):6811.9 Kullback, S. and Leibler, R. A. (1951). On information and sufficiency.
The annals of mathematicalstatistics , 22(1):79–86.10 Ladd, A. J. C. (1993). Short-time motion of colloidal particles: Numerical simulation via afluctuating lattice-boltzmann equation.
Phys. Rev. Lett. , 70:1339–1342.11 McNamara, G. R., Garcia, A. L., and Alder, B. J. (1995). Stabilization of thermal latticeboltzmann models.
Journal of Statistical Physics , 81(1-2):395–408.12 Osborn, W., Orlandini, E., Swift, M. R., Yeomans, J., and Banavar, J. R. (1995). Lattice boltzmannstudy of hydrodynamic spinodal decomposition.
Physical review letters , 75(22):4031.13 Pachalieva, A. and Wagner, A. J. (2020). Non-Gaussian distribution of displacements forLennard-Jones particles in equilibrium.
Phys. Rev. E , 102:053310.14 Parsa, M. R., Pachalieva, A., and Wagner, A. J. (2019). Validity of the molecular-dynamics-lattice-gas global equilibrium distribution function.
International Journal of Modern Physics C ,30(10):1941007.15 Parsa, M. R. and Wagner, A. J. (2017). Lattice gas with molecular dynamics collision operator.
Physical Review E , 96(1):013314.16 Parsa, M. R. and Wagner, A. J. (2020). Large fluctuations in nonideal coarse-grained systems.
Phys. Rev. Lett. , 124:234501. r s t a . r o y a l s o c i e t y pub li s h i ng . o r g P h il . T r an s . R . S o c . A ..................................................................
17 Plimpton, S. (1995). Fast Parallel Algorithms for Short-Range Molecular Dynamics.
Journal ofComputational Physics , 117(1):1–19.18 Qian, Y. H., D’Humières, D., and Lallemand, P. (1992). Lattice BGK Models for Navier-StokesEquation.
Europhysics Letters (EPL) , 17(6):479–484.19 Shan, X. and Doolen, G. (1995). Multicomponent lattice-boltzmann model with interparticleinteraction.
Journal of Statistical Physics , 81(1-2):379–393.20 Wagner, A. J. and Strand, K. (2016). Fluctuating lattice boltzmann method for the diffusionequation.