Testing and validation of a self-diffusion coefficient model based on molecular dynamics simulations

Xia Chen ,Yan Wang ,Lianying Wu, *,Weitao Zhang ,Yangdong Hu

1 College of Chemistry and Chemical Engineering,Ocean University of China,Qingdao 266100,China

2 College of Chemistry and Chemical Engineering,Qingdao University,Qingdao 266071,China

Keywords:Diffusion Molecular simulation Parameter identification Characteristic length Diffusion velocity

ABSTRACTIn our previous work,we endowed a new physical meaning of self-diffusion coefficient in Fick’s law,which proposed that the diffusion coefficient can be described as the product of the characteristic length and the diffusion velocity.To testify this simple theory,in this work,we further investigated the underlying mechanism of the characteristic length and the diffusion velocity at the molecular level.After a complete dynamic run,the statistical average diffusion velocity and the characteristic length of molecules can be obtained by scripts,and subsequently the diffusion coefficient was determined by our proposed theory.The diffusion processes in 35 systems with a wide range of pressure and concentration variations were simulated using this model.From the simulated results,diffusion coefficients from our new model matched well with the experimental results from literatures.The total average relative deviation of predicted values with respect to the experimental results is 8.18%,indicating that the novel model is objective and rational.Compared with the traditional MSD-t model,this novel diffusion coefficient model provides more reliable results,and the theory is simple and straightforward in concept.Additionally,the effect of gas pressure and liquid concentration on the diffusion behavior were discussed,and the microscopic diffusion mechanism was elucidated through the distribution of diffusion velocity and the characteristic length analysis.Moreover,we suggested new distribution functions,providing more reliable data theoretical foundations for the future research about the diffusion coefficient.

1.Introduction

Diffusion is the tendency of particles to spread into an available area.It plays an important role in every natural or industrial process involving mass transfer.Therefore,the determination of diffusion coefficient is of great importance for the calculation of mass transfer processes.There are two kinds of classical diffusion theories used to define diffusion coefficient:Maxwell-Stefan’s (MS)equations and Fick’s law[1–4].The diffusion coefficient of the former has a definite physical meaning,which represents the friction factor between two components [5].The latter is commonly used in mass transfer process with a simple expression,but the physical meaning for its diffusion coefficient is not clear [6].While,in our previous work [7],we endowed a novel physical meaning of the diffusion-coefficient in Fick’s law,in which the diffusion coefficient can be described as the product of the characteristic length and the diffusion velocity.

In many systems,diffusion coefficients are usually obtained through experimental approaches [8,9],such as conduct metric method,Taylor dispersion technique and digital holography method,etc.For example,Rodrigo et al.[10] acquired the binary diffusion coefficients of L-histidine methyl ester dihydrochloride solutions by using the Taylor dispersion technique.Zhou et al.[11]obtained the diffusion coefficient of the ethylene glycol–water system by employing digital holography method.Mattisson et al.[12] measured the diffusion coefficients for lysozyme in gels and liquids by using holographic laser interferometry.Jamnongwong et al.[13] studied the oxygen diffusion coefficients in clean water containing salt,glucose or surfactant based on measurements of volumetric mass transfer coefficients.Wang et al.[14] determined the gas diffusion coefficient in porous shale.Albeit these experimental methods are direct and easy to operate,the measurement of diffusion coefficients in certain systems is tedious,costineffective and even impractical.

To this end,computer simulation method has been introduced to obtain the diffusion coefficients.Nowadays,molecular dynamics(MD) simulation has become an effective tool to calculate the diffusion coefficients [15–20].Various studies on the diffusion properties research have been carried out through MD simulation.For example,Yang et al.[21] investigated the diffusion of methanol/water mixture through MFI-type zeolite (HZSM-5 and silicalite-1)membranes by employing MD simulation.Hu et al.[22]and Zhang et al.[23]computed the diffusion of corrosive particle via MD simulation.Ghaffari et al.[24] explored the self-diffusion coefficients in NaCl aqueous solutions by using MD simulation.In recent study,Zhao et al.[25] calculated the transport diffusion coefficients of CO2and CH4in coal via the molecular simulation.Tsimpanogiannis et al.[26]presented a detailed discussion about self-diffusion coefficient of water by classical molecular simulation.Moultos et al.[27–32] reported a series of papers regrading to the systematic errors in the calculation of self-diffusion caused by the longrange interactions and the imposed periodic boundary conditions.

On the other hand,to study the underlying mechanism of the microscopic phenomena and to predict the diffusion coefficient,various physical and mathematical models have been proposed[33,34].For instance,Kamgar et al.[35] suggested a simple model only with the molecular mobility factor and the thermodynamic correction factor.Considering an analogy between nonideal concentrated solutions and solutions near the consolute point,Cussler et al.[36]suggested to mimic the diffusion through the movement of entire molecules clusters.

Although the novel diffusion coefficient model was given in our previous work,we did not give the values of characteristic length and diffusion velocity.In this work,we further investigated the underlying mechanism of our previously proposed diffusion coefficient model at the molecular level.Through MD simulation,the statistical average diffusion velocity and the characteristic length were obtained,and subsequently the diffusion coefficient was given by our proposed model.Based simulated results of 35 systems,this novel diffusion coefficient model was proved to be objective and rational.Moreover,the effect of gas pressure and liquid concentration on the diffusion behavior were discussed,and the microscopic diffusion mechanism was elucidated through the distribution of diffusion velocity and the characteristic length analysis.

2.MD Simulation

2.1.Theory of the diffusion coefficient model

In our previous work,we developed a new diffusion coefficient model [7] by analyzing the unit of Fick’s Law diffusion coefficient through the dimension analysis method.The proposed model contains two parts:the characteristic length item and the diffusion velocity item.The clear physical meaning of Fick’s Law diffusion coefficient can thus be described as the product of the characteristic length and the diffusion velocity,as follows:

where Videnotes the diffusion velocity;Lidenotes the characteristic length.

Viis the velocity of molecular diffusion,which can be considered as the statistical average velocity of the molecular motion.Liis characteristic length,which is the statistical average of the diffusion distance.And the diffusion distance is the distance which traveled by a moving molecule without changing direction.This is similar to the meaning of the free path of gas molecules.To testify this simple theory,we calculated the characteristic length and the diffusion velocity via MD simulation.

2.2.Simulation method

In this study,all simulations were performed by means of COMPASS force field in Materials studio software.Each diffusion simulation was repeated more than 10 times to obtain good statistic.The initial configurations were built individually using the Amorphous Cell module.The appropriate length of the simulation box was given in the Tables 1–3 according to the concentration and density.The simulations were performed in a cubic box with periodic boundary conditions and different directions.The Nose thermostat and the Berendsen barostat was used to control the system temperature and pressure,respectively.Ewald method was used to simulate electrostatic and van der Waals interactions with a cut-off value of 15.5 ?(1 ?=0.1 nm).For the liquid system,the configurations were firstly subjected to geometry optimization,and then run in NPT(1 ns),NVT(1 ns)and NVE(1 ns)ensemble.A small integration time step of 0.1 fs was selected to obtain more statistic samples.The trajectories of liquid molecules were saved every 10 steps for further analysis.For the gas system,the configurations which were built based on the gas density were firstly subjected to geometry optimization,and then run in NVT (8 ns)and NVE (8 ns)ensemble.For the gas system,the integration time step was 1 fs.The trajectories of gas molecules were sampled every 250 steps for further analysis.

MD results can provide the microscopic details of molecular diffusion to reveal the underlying mechanism of diffusion,such details cannot be directly obtained through experimentation.The diffusion velocity is actual the average velocity of the molecular motion,which can be considered as the statistic average velocity form the available statistic samples.After a complete dynamic run,the instantaneous velocity of each molecule per femtosecond directly obtained by scripts.Based on the available statistic samples,the statistical average diffusion velocity was given.Then Vi(the average velocity of the molecular motion)gained.The characteristic length denotes the statistical average of the diffusion distance which traveled by a moving molecule without changing direction.After a complete dynamic run,the position coordinates of each molecule per femtosecond (such as the path in Fig.1)obtained.The deviation angle of ever frame can be calculated by the arccosine function.For example,as shown in the partial enlarged drawing of the Fig.1b,the deviation angle of frame A(i.e.∠BAC) calculated by the arccosine function (Eq.(2)).We defined that the moving molecule change its direction when the moving direction deviating from a straight line exceeds 2 degrees.Thus,the distance between two adjacent frames that deviation angle exceeded 2 degrees(e.g.∠BAC>2°)is the diffusion distance.Then Li(the statistical average of the diffusion distance) obtained.

Through statistical analysis of the displacement and velocity data,the average diffusion distance and velocity can be obtained.And subsequently the diffusion coefficient was determined by our proposed model (Eq.(1)).

For comparison,we also calculated the diffusion coefficient by mean square displacement (MSD) approach.The MSDis the square of particle displacement during a time interval.Based on the well-known Einstein equation (Eq.(3)),the diffusion coefficient can be obtained.

where Nais the number of diffusion molecules,is the position vector of the molecule at time t,andis the initial position.

3.Results and Discussion

3.1.Computation of MSD

The presence of numerous statistic samples indicates high linearity of the MSD-t curve,resulting in accurate calculation of the diffusion coefficient.Abnormal diffusion,which usually occurs at beginning and ending part of the MSD,is nonlinearly related to time.These irregular data should not be used to calculate the diffusion coefficient.The fitting was conducted only in progresses into normal diffusion.Therefore,the determining time interval(tmin≤t ≤tmax)is necessary for diffusion coefficient calculation.The diffusion of CO2was taken as an example to demonstrate the determination of tminand tmax.The MSD–time curve of CO2is presented in Fig.1(a).

Double logarithmic coordinate plots of MSD are shown in Fig.2(b).The MSD-t curve is straight line when the slope of lg(MSD)-lg(t)curve equal to 1.To observe the fluctuation of the slope of lg(MSD)-lg(t),the derivate of the lg(MSD)-lg(t) is depicted,as shown in Fig.2(a).This work treats the discrete data by using a centered difference formula to approximate the derivative.The derivative function applied to discrete data points can therefore be written as follows:

Fig.1.(a) The diffusion pathway of N2;(b) The diffusion pathway of Alanine (0.59 mol?L-1).

The MSD slope is considered as constant when the fluctuation of the slope of lg(MSD)-lg(t) is ranging from 0.8 to 1.2.Thus,tminand tmaxcan be determined,as shown in Fig.2(a).Consequently,the diffusion coefficients can be calculated from this time interval.The simulated results are listed in Tables 1–3 accordingly.

3.2.Results of the diffusion coefficient,characteristic length,and diffusion velocity

After a complete dynamic run,the diffusion velocity and path of each frame for every molecule obtained by scripts.Based on the available statistic samples,the statistical average diffusion velocity(namely the diffusion velocity)was given.The characteristic length got from the diffusion path.As illustrated in Fig.1,it is obvious that the travel distance before changing direction is much shorter in liquid systems than in gas systems,which means that a moving molecule change direction more easily in the liquid than in the gas.The simulated results of Vi,and Lifor the diffusion of gases,organic vapor,and the liquid system are listed in Tables 1–3 accordingly.

After the characteristic length and the diffusion velocity gained.The diffusion coefficient was determined by our proposed model(Eq.(1)).For liquid system,the self-diffusion coefficients were corrected for finite-size effects by using the Yeh-Hummer correction(see Eqs.(5)–(7)) [30].And for the gas system,we used the simulation box with large size (1000 molecules) to avoid the error caused by system size.

where DYHis the YH correction diffusivity.D?cal-MSDand D?cal-LVare the calculated diffusion coefficient from the MSD-t curve and our new model respectively.Dcal-MSDand Dcal-LVare the corrected D?-cal-MSDand D?cal-LV.kBis the Boltzmann constant,T is the temperature,η is the shear viscosity,L is the length of the cubic simulation box,ξ is a constant which depends on the shape of the simulation box (for a cubic simulation box,ξ=2.837297).

Fig.2.(a) The MSD-t curve of CO2;and the derivate of the lg(MSD),tmin and tmax are marked;(b) The lg(MSD)-lg(t) curve.

The simulated results of Dcal-LVand Dcal-MSDfor the diffusion of gases,organic vapor,and the liquid system are listed in Tables 1–3 accordingly.The calculation of relative deviations (RD,%) and the mean relative deviations (ARD,%) are based on the following equations:

3.3.Gas system analysis

The simulated results for the diffusion of various gas and organic vapor systems are listed in Tables 1 and 2.From Tables 1 and 2,the calculated diffusion coefficients by our new model are in best agreement with those experimental results from literatures.The minimum value of RD-1is 1.89% and the ARD-1is 7.53%,indicating that the new model proposed previously are objective and rational.Besides,the diffusion coefficients were calculated after the high linearity of the MSD-t curve was determined,as listed in Tables 1 and 2.Compared with the traditional MSD-t model,this novel diffusion coefficient model provides more reliable results,and the theory is simple and straightforward in concept.

Besides,in this study,the effects of pressure were evaluated by varying pressures of CO2,N2and O2in MD simulations.As presented in Table 2,the value of diffusion coefficients and the characteristic length decreased dramatically with increasing the gas pressure,while the diffusion velocity appeared less influenced by the pressure.This is quite understandable,since the intermolecular distance will reduce when the pressure is rising,the characteristic length will decrease.Meanwhile,the molecular velocity is little changed.As a result,the value of diffusion coefficient will reduce.

Table2 The simulated results of the gas systems at 300 K for different pressures

Moreover,the diffusion coefficients follow the order of CO2

Table1 The simulated results of various gas and organic vapor systems at 300 K and 1.01 MPa

3.4.Liquid system analysis

As shown in Table 3,the calculated diffusion coefficients are perfectly matched with the experimental results.The minimum value of RD-1is 0.06% and the ARD-1is 8.6%.This further evidence that our new model proposed previously accurately described the diffusion coefficient.

Table3 The simulated results of different liquid systems at 300 K and 1.01 MPa

As shown in Table 3,we tested Alanine diffusion stimulation at different concentrations.It was found that the increase of the Alanine concentration in the aqueous solution caused decreases in both the diffusion coefficient and the characteristic length,but slight changes in the diffusion velocity.The trend of the characteristic length is caused by the reduction of distance between alanine molecules.However,the diffusion velocity was more sensitive to temperature changes.

The other notable rule is that the characteristic length and the diffusion velocity in liquid system are smaller than that in gas system.Those results about characteristic length verified above leads to a conclusion:when the intermolecular distance is reduced,the diffusion coefficients will be decreased.Generally,there are three different ways of molecule movement:translation,rotation,and vibration.In gas system,particles are far apart,with large intermolecular space.In this case,the intermolecular force of particles is relatively small.Translational motion makes a great contribution to this gas molecular motion.As a result,the direction of gas molecular displacement is not easy to change.Whereas in liquid system,the interaction force among particles is relatively large.In this case,molecular rotations and vibrations should be considered.The direction of molecular displacement is much easier to change in liquid system than in gas system,leading to smaller diffusion distance in liquid system.

3.5.Statistical analysis

According to a large number of statistical data from the trajectories,the distribution of diffusion velocity and the characteristic length can be analyzed through the probability density function.

Fig.3.(a) The diffusion velocity distribution of Alanine (0.59 mol.L-1);(b) The diffusion velocity distribution of CO2.

3.5.1.For diffusion velocity

By statistically analyzing of 8 liquid systems and 8 gas systems,we found that the probability density dot for the diffusion velocity firstly increased and then decreased.This pattern is similar to the normal distribution.Hence,we developed a new distribution function for describing the probability distribution of the diffusion velocity based on the normal distribution.The probability distribution function was shown in Eq.(10).As an example,the distribution of diffusion velocity of Alanine (0.59 mol?L-1) and CO2are shown in Fig.3,respectively.From Fig.3,it illustrated that the probability density dot of diffusion velocity obeys our new distribution function very well.

where Viis the velocity;p1,p2,p3are parameters.

After calculation,the values of the parameters and the Pearson correlation coefficient(Pearson’s r)are listed in Table 4.The results show that the new distribution function in this work are satisfactory for the distribution of diffusion velocity.The minimum value of Pearson’s r is 0.99 and the mean value is 0.99.It is worth noting that the parameter p3is approximate to the average diffusion velocity.This indicates that the physical meaning of p3is related to the average diffusion velocity.

Table4 The parameters of the distribution function for diffusion velocity

The calculation of the average relative deviation (Pearson’s r)was based on the following equation:

3.5.2.For characteristic length

Based on the statistically analyzing of 8 liquid systems and 8 gas systems,we found that the tendency of the probability density dot for the characteristic length increase at first,and then decrease,and then go to be horizontal with the characteristic length increasing.Take the characteristic length distribution of Alanine(0.59 mol.L-1) and CO2as an example,as shown in Fig.4.

Fig.4.(a) The characteristic length distribution of Alanine (0.59 mol.L-1);(b) The characteristic length distribution of CO2.

Molecular motion is a random event,and the molecular diffusion process can be considered into two random events:One is the small amplitude diffusion,and the other is the saltatory diffusion.The former one is a high probability event.This pattern is similar to the normal distribution.And the later one is a small probability event,the regularity of distribution is similar to the uniform distribution.Therefore,we developed a new distribution function for describing the probability distribution of the characteristic length.This new distribution function contains two parts:the former part obeys a new distribution (see Eq.(12)).While the later part obeys a uniform distribution (see Eq.(13)).After calculation,the values of the parameters and the Pearson’s r are presented in Table 5.It was showed that the new distribution function could describe the distribution of characteristic length well.The minimum value of Pearson’s r is 0.96 and the mean value is 0.99.

Table5 The parameters of the distribution function for diffusion distance

where Liis the diffusion distance;[a,b] is the value range of Li;a1,a2,a3is parameters.

4.Conclusions

Based on the molecular dynamics simulations,we further investigated the underlying mechanism of the novel diffusion coefficient model at the molecular level.In this work,molecular dynamics simulations were performed for diffusion process in 35 systems(including 12 liquid systems as well as 23 gas and organic vapor systems),to investigate our novel diffusion coefficient model.With the dynamic simulation,the statistical average diffusion velocity and the characteristic length can be obtained from the trajectory file,and subsequently the diffusion coefficient defined by proposed model was determined.From the simulated results,diffusion coefficients from our new model matched perfectly with those experimental results from literatures.The ARD-1of the simulations results was 7.53%and 8.60%in gas and liquid system,respectively,indicating that the novel model proposed previously are objective and rational.Compared with the traditional MSD-t model,this novel diffusion coefficient model provides more reliable results,and the theory is simple and straightforward in concept.Moreover,the effect of gas pressure and liquid concentration on the diffusion behavior were discussed,and the microscopic diffusion mechanism was elucidated through the distribution of diffusion velocity and the characteristic length analysis.Meanwhile,new distribution functions were suggested to describe the distribution of the diffusion velocity and the characteristic length,which provided theoretical foundations for the future research about the diffusion coefficient.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (21776264 and 21376231).My sincere appreciation goes to Prof.Hu and Prof.Wu who provided scientific guidance for me.

Nomenclature

a1,a2,a3represent the parameters of the distribution function

c concentration of component I,mol?L-1

Dcal-MSDcalculated diffusion coefficient from the MSD-t curve,m2?s-1

Dcal-LVcalculated diffusion coefficient from our new model,m2?s-1

Dexpexperimental diffusion coefficient,m2?s-1

KBBoltzmann constant

Lcal,Licharacteristic length,m

L length of the cubic simulation box

M molar mass,g?mol-1

Nanumber of diffusion molecules

Npnumber of samples

p pressure,kPa

p1,p2,p3represent the parameters of the distribution function

T temperature,K

t Time,s

Vcal,Vidiffusion velocity,m?s-1

η shear viscosity

ξ a constant which depends on the shape of the simulation box

σ represent the standard deviation

σystandard deviation of X

σ^ystandard deviation of Y