Self-consistent molecular dynamics calculation of diffusion in higher n-alkanes

N.D. Kondratyuk1,2 G.E. Norman1,2, V.V. Stegailov1,2* 1Joint Institute for High Temperature of RAS, Moscow, Russia 2Moscow Institute of Physics and Technology, Dolgoprudny, Russia *stegailov@gmail.com Diffusion is one of the key subjects of molecular modeling and simulation studies. However, there is an unresolved lack of consistency between Einstein-Smoluchowski (E-S) and Green-Kubo (G-K) methods for diffusion coefficient calculations in systems of complex molecules. In this work [1], we analyze this problem for the case of liquid n-triacontane. The non-conventional long-time tails of the velocity autocorrelation function (VACF) are found for this system. Temperature dependence of the VACF tail decay exponent is defined. The proper inclusion of the long-time tail contributions to the diffusion coefficient calculation results in the consistency between G-K and E-S methods. Having considered the major factors influencing the precision of the diffusion rate calculations in comparison with experimental data (system size effects and force field parameters), we point to hydrogen nuclear quantum effects as, presumably, the last obstacle to fully consistent n-alkane description. The first results of the path-integral molecular dynamics will be presented. The work was supported by the Russian Science Foundation (grant 14-50-00124).

Decades of the molecular modeling and simulation history have not still provided a desirable accuracy of materials properties calculations in many cases.Diffusion rate calculations in liquids of complex molecules are among them.
There are two equivalent approaches for the calculation of the diffusion coefficient: the Einstein-Smoluchowski (E-S) and Green-Kubo (G-K) methods.The first one uses a mean-square displacement (MSD), the second takes an integral of velocity autocorrelation function (VACF).The E-S formula is used more frequently in molecular dynamics (MD) simulations than the G-K method because MSD is easier to interpret and due to the problems with integration to infinite time in the G-K method.The agreement between the both methods is demonstrated in the case of atomic and simple molecular liquids [1].The disagreement appears in polyatomic molecular systems: the G-K method gives overestimated values in comparison with the E-S approach [2,3].The reason of this discrepancy has not been understood yet.
Here, we use n-triacontane (n-C 30 H 62 ) as an example of a polyatomic molecule.The potential energy of the liquid n-triacontane (Fig. 1.) can be modeled as follows: We use three different force fields to parameterize this equation: TraPPE-UA (united-atom) [4], OPLS-AA (all-atom) [5] and L-OPLS-AA [6] for making sure that obtained results are not artifacts of a particular model.L-OPLS-AA is modified OPLS-AA which is calibrated to reproduce thermodynamic and transport properties of long alkanes.
The diffusivity of liquid n-triacontane is calculated using both E-S and G-K relations.The simulations are performed for the systems of 3375 (OPLS-AA, L-OPLS-AA) and 8000 (TraPPE-UA) molecules.These configurations are obtained by replicating the equilibrated systems of 125 molecules.The averages for MSDs and VACFs are taken from 1 ns equilibrated MD trajectories.The simulations are carried out in the LAMMPS package [7].
The molecule center of mass Δr 2 has a subdiffusive part ( Δr 2 ~tα , α < 1 ), caused by the molecular entanglement at low temperatures.The diffusion coefficients are calculated from the linear asymptotes of Δr 2 .The system size effects on the diffusion coefficients [8] are also considered.The size-corrected value can be found as lim D(N ) at N −1/3 → 0 .The correction is about 6 % for 3375 molecules in the all-atom models and 4 % for 8000 molecules in TraPPE-UA.
The VACF v(0)v(t) of molecules are calculated from the same trajectories as Δr 2 .
The numerical integration of v(0)v(t) gives 1.5 times larger diffusion coefficients than the E-S values.The long-time v(0)v(t) asymptotes are analyzed to find the solution of this issue.The asymptotes are found to have a power form At − β .The values of β are shown in Fig. 2 for the different force fields.They differ from the power of hydrodynamic tail t −3/2 demonstrated for atomic liquids [9].β(T ) decreases as 1 / T with the increase of temperature in the all-atom models.On the other hand, β coincides with the conventional 3/2 hydrodynamic power law for TraPPE-UA and does not depend on the temperature.Therefore, unusual values of β in OPLS-AA and L-OPLS-AA can be connected with the explicit hydrogen atoms.Authors [10] showed that the VACF integration can be subdivided into the numerical contribution and the analytical asymptote one in the binary Lennard-Jones liquid.Such an approach is used to achieve convergence between the E-S and G-K methods in this work as well.
This technique gives a compliance of both E-S and G-K relations within the accuracy of the methods.
This fact verifies the non-conventional values of β obtained in the all-atom force fields.Our results are compared with the diffusion coefficients measured in the experiment in the wide temperature range [11].Both our calculation results and the experimental data match the Arrhenius dependence D = D 0 exp(−ΔE / k B T ) .The OPLS-AA and the L-OPLS-AA force fields predict ΔE to be equal 4.9±0.2kcal/mole and 4.6±0.1 kcal/mole respectively whereas the experimental value is 4.5±0.1 kcal/mole.However, the prefactor D 0 is about 1.8 times lower in the OPLS-AA force field.
Conversely to all-atom models, the calculations in the TraPPE-UA model show that the united atom approach overestimates the diffusion coefficient.The difference can be attributed to the absence of hydrogen atoms which increases the molecule free volume.The value of activation energy ΔE is 3.5±0.1 kcal/mole which is much lower then the experimental value.
The work is supported by the RSF grant 14-50-00124.

Fig. 1 .
Fig. 1.The n-triacontane system in the all-atom model.