Rheology of liquid n-triacontane: Molecular dynamics simulation

Molecular dynamics is applied to calculate diffusion coefficients of n-triacontane C30H62 using Einstein-Smoluchowski and Green-Kubo relations. The displacement 〈Δr2〉(t) has a subdiffusive part 〈Δr2〉 ∼ t α, caused by molecular crowding at low temperatures. Longtime asymptotes of 〈v(0)v(t)〉 are collated with the hydrodynamic tail t-3/2 demonstrated for atomic liquids. The influence of these asymptotes on the compliance of Einstein-Smoluchowski and Green-Kubo methods is analyzed. The effects of the force field parameters on the diffusion process are treated. The results are compared with experimental data.


Introduction
Diffusion processes in materials are of a great interest for researchers during the recent decades. Molecular dynamics (MD) methods are applicable for such problems because they allow to get transport properties from the atomic motion. There are two equivalent methods for calculation of the diffusion coefficient: by Einstein-Smoluchowski (E-S) and Green-Kubo (G-K) [1]. The first one uses a mean-square displacement (MSD) of the particles, the second takes an integral of velocity autocorrelation function (VACF). Both methods are widely used [2][3][4][5][6][7][8] and show a good agreement in the case of atomic and simple molecular liquids [9,10]. The correct calculation of the diffusion in polyatomic molecular systems faces difficulties: the G-K method gives overestimated values in comparison with E-S [11][12][13]. The reason of this discrepancy is not understood yet. Viscosity calculations using the G-K method also face problems with convergence. A time composition method is proposed to avoid such problems in [14] and succesfully applied in [15].
The long-time tail of VACF reflects the physical nature of the system. It decays exponentially in Brownian gases. VACF has two asymptotic exponential regimes in nano-particle liquids [16]. The oscillation behavior of VACF is obtained for an ion movement in a liquid [17]. A hydrodynamic power law t −3/2 is predicted in the case of atomic liquids and dense gases [18]. It is a result of the collective motions. Simulations for Lennard-Jones liquids prove this fact [19,20].
The VACF asymptotes are not discussed enough for systems of complex molecules. The first simulation results on diffusion for n-alkanes are obtained in [21]. The united-atom approach with frozen bond and angle values is used to simulate n-butane liquid. The rotation of the end atoms around the central bond (dihedral interaction) is adjusted to the experimental data that are available on that moment. Authors use both calculation methods, the G-K method gives the overestimated values of diffusion. The more complicated form of the united-atom force field (with angle intermolecular interactions) is used for understanding the nature of diffusion in n-pentane and n-decane in [22,23]. They show that diffusion in n-alkanes depends on the dihedral interactions (by varying CH x -CH 2 -CH 2 -CH x torsion barriers) and hypothesize that the molecule flexibility determines the decay of the VACF.
The comparison with experiment is a part of the discussion about simulation results and demonstrates whether the model used is applicable. Authors [10] show a good agreement of calculated H 2 diffusion coefficient solvated in water with experimental values. Such accurate estimates are made for atomic liquids [9]. However, there are not many works on diffusion in complex molecular liquids, in which comparisons with experimental data are provided [5].
Here, n-triacontane C 30 H 62 is used for studying the diffusion in complex liquids. The structure of the paper is the following. The basic relations, model, equilibration and averaging techniques are considered in the second section. The results of the calculation, the influence of the force field parameters on the diffusion processes are presented in the third section. The VACF asymptotic tails and the convergence of both methods to experimental values are discussed in the fourth section. The results are assumed to be the step toward supplementing for the macroscopic models [24,25].

Simulation technique 2.1. Basic relations
The self-diffusion coefficient D can be obtained from the long-time behavior of the center of mass (COM) mean-squared displacement ∆r 2 using the well-known E-S relation [1]: It is given also by the theoretically equivalent G-K integral formula [1]: where C v is a molecule COM velocity autocorrelation function: The angle brackets in equations (1) and (3) mean the averaging over canonical ensemble.

Model
The interaction in the n-triacontane liquid can be subdivided into intramolecular and intermolecular parts. The first part describes interactions between atoms of the same molecule. It contains bonded, angle and torsion interactions as well as Lennard-Jones and Coulumb (figure 1), The form of these terms and their parameterization depend on the model of the force field. The open potential for liquid simulations-all atom (OPLS-AA) [26] is used in this work. The bonded and angle interactions have harmonic forms: where K a and K b are energy constants, R and θ are a bond-length and angle, R 0 and θ 0 are an equilibrium bond-length and an angle. The torsional interaction is presented via the sum of the cosines of dihedral angle φ with different periodicities (that reflect atom hybridization): where K 1 , K 2 , K 3 , K 4 are energy constants that are obtained from the ab initio molecular orbital calculations [26]. Non-bonded interactions in the molecule (between atoms, that are not included in the bonded and angle interactions) are represented via the Lennard-Jones 6-12 form: where ǫ is an energy constant, σ is a zero-crossing parameter, r is a distance between two atoms. Electrostatic interactions have a form Cq i q j /r, where q i and q j are partial charges of the i-th and j-th atoms, r is a distance between them. These parameters are obtained to reproduce thermodynamic and structural properties of hydrocarbon liquids using Monte Carlo statistical mechanics simulations [26]. The contribution of the long-range Coulumb interactions is calculated using the particle-particle particle-mesh [27] method with a real space cutoff 12Å. The intermolecular part describes forces between atoms that belong to different molecules. It consists of the Lennard-Jones and Coulumb parts.

Equilibration
The starting configuration is a gas of the 125 replicated n-triacontane molecules. The distance between them is larger than the force field cutoff radius. At the first stage, the temperature is set to 600 K. This temperature is needed to perform rapid molecule disorder. The molecules get random orientations in the so called NVE ensemble [28] for 0.1 ns. The second stage is a compression to the experimental density 0.77 g/cm 3 [29] for 0.1 ns. The third stage is a relaxation in the NPT ensemble for 2 ns (P ∼ 1 atm, T = 353 K). The average density at the last 0.5 ns is chosen as the equilibrium value and corresponds to the experimental density at 353 K. The fourth stage is the relaxation in the NVT ensemble at the calculated density for 2 ns. The relative shape anisotropy parameter κ 2 is used in addition to the main equilibration parameters such as temperature, pressure and density. It reflects the information about the average molecule conformation in the system (for details see [30]). The integration timestep is 1 fs. The MD simulations are carried out in the LAMMPS package [31].

Averaging
The diffusivity simulation is performed for the system of 3375 molecules in the NVE ensemble. This configuration is obtained by replicating the equilibrated system. The averages for equations (1) and (3) are obtained from 1 ns equilibrated MD trajectories. The following technique is used. The whole trajectory is divided into statistically independent parts. The duration of these intervals τ should be more than the value of the memory time t m [28]. This time corresponds to the moment when the numerical solution of the equations of motion forgets its initial conditions due to the Lyapunov instability.
At the first step, the average is taken over the molecules: where N is the number of the molecules. The second step is the additional averaging over the intervals that is achieved by a shifting of the zero time point: where M is the number of such intervals. The value of t m for the n-triacontane system is about 3 ps. The characteristic times, when the asymptotic tails of the VACF can be obtained (∼ 5 ps), are a hundred times less than the times, that are needed for the observing linear ∆r 2 (t) asymptote in the case of E-S (∼ 750 ps). Consequently, the number of the intervals M , obtained from the same 1 ns trajectory, is about 60 in the E-S method (60 shifts of the zero time point for the value of 4 ps) and M = 300 in the G-K method. In the section four it is shown that the G-K method requires more statistics for the accurate calculation of the VACF.

Abnormal diffusion
The MSD in reference systems (such as Lennard-Jones gases and liquids) has two regimes: ballistic (at the very beginning particle moves freely, ∆r 2 = v 2 t 2 ) and diffusive ( ∆r 2 = 6Dt, caused by stochastic collisions with neighbors). The ballistic part becomes diffusive continuously in the atomic case.
The C 30 H 62 MSD has an intermediate subdiffusive part (figure 2, the black line), where ∆r 2 ∼ t α , α < 1. The situation, when the intermediate regime appears, is studied in many theoretical and experimental works [2,4,8,[32][33][34]. The effect of the diffusion slowing is obtained using MD for ionic liquids in [33] and for hydrogen gas hydrates in [35], and using ab initio MD for hydrogen molecules solvated in water [10]. It is assumed that such behavior can be explained by the molecular crowding [32].
The torsion interactions play a major role in such crowding because they determine how easily the molecule bends. To prove this fact, MD simulations are performed with the different H-C-C-X energetic barriers in equation (7) for hydrogen atoms (0.5K i and 2K i , i =1-4). The molecule velocities are rescaled to the initial temperature to compensate the changes of the system energy. It is needed for performing the simulation with the same temperature and pressure as in the case with the original parameters. The decrease of these barriers makes the diffusion faster and the subdiffusion region shorter (the blue line in figure 2). The molecules are more flexible and find the path through neighbors easier. The increase of the barriers has an opposite effect (the red line): the motion of the molecule COMs slows down. The reason is that the molecules are more rigid and it becomes difficult to change mutual orientations of the neighbors. The averaging technique used is compared with the LAMMPS method (the solid and dashdot lines in figure 2 respectively). The difference is that LAMMPS uses only one zero time point (the beginning of the MD trajectory) and corresponds to the first step of the technique that is used in this work. It allows to compute the MSD on the fly. The accuracy of the LAMMPS method is very sensitive to the amount of molecules. The second stage of our technique allows to use the whole trajectory for the averaging and covers more microstates. Thus, the subdiffusion part has some differences with the LAMMPS method (due to the better statistics), but the short and long-range asymptotes coincide with each other in both methods.
The diffusion coefficients D E−S are obtained from the long-range linear asymptotes (the grey dashed lines in figure 2). They are compared with the experimental value [36] (table 1). OPLS-AA force field with the original parameters underestimates the diffusion coefficient. Such deviation can be explained by the fact that the calibration of OPLS-AA is not based on the values of diffusion.

VACF asymptotes
The molecules COM VACFs are calculated from exactly the same trajectories as the COM MSDs (figures [3][4][5]. The negative region of the VACF is typical for liquids and displays the fact, that rebounding collisions are more frequent than scattering collisions in dense systems [37,38]. C v (t) decays as t −3/2 at the long-time, that is also called hydrodynamic tail, which is typical for all atomic and simple molecular liquids [18]. The reason of this asymptote is a collective motion of particles.
The accurate calculation of the VACF tail is a non-trivial issue. Firstly, the size of the system should be enough to exclude the influence of the periodic boundary conditions [39]. The second condition is averaging over a sufficiently long trajectory because the values of C v (t)/C v (0) become about 10 −2 and the best accuracy is required.  The power of the tail is different in the original OPLS-AA model ( figure 3). The interesting fact is that it becomes closer to the classical value with the decrease of the torsion energetic constants ( figure 4). On the contrary, the increase of these barriers makes the power of the tail larger ( figure 5). These effects are connected with the flexibility of the molecules. The molecules COMs move more freely in the force field with 0.5K i and the nature of their movement becomes closer to the Brownian. The molecular crowding takes place in the case of the higher H-C-C-X barriers.
The G-K integral (2) is taken numerically till 3 ps and analytically after 3 ps using the long-range tails (table 1): The analytical contribution D a plays a major role in the correct calculation of the diffusion coefficient, because the numerical integration until the moment when the integral becomes a constant gives an overestimated value. The reason is that stochastic fluctuations of the VACF around zero at the long times interfere the real behavior. Thus, the sum over these fluctuations is zero and the diffusion coefficient is much higher than the E-S value. The authors of [11][12][13] face similar problems in ionic liquids.
The small change of the asymptotic tail power (within 0.1) causes a great change of the integral (2). This fact makes the G-K method less applicable than E-S, because the analytical contribution D a depends on this parameter. The errors are evaluated in the following way. The maximum D max a (that corresponds to the maximal suitable power) and minimum D min  The values of diffusion that are calculated using the G-K relation agree with the E-S values within the method accuracy (table 1). The usage of the classical t −3/2 tail gives underestimated value of diffusion in comparison with the E-S results. This convergence validates the nonclassical VACF tails that are discovered in this work for the n-triacontane liquid.
The increase of the MD trajectory duration (or of the molecules number) could improve the accuracy of the G-K method but requires more computational resources.

Conclusions
The n-triacontane C 30 H 62 liquid is used as an example for the molecular dynamics simulation of the diffusivity in complex liquids using Einstein-Smoluchowski and Green-Kubo methods. The dependence of the diffusion processes on the OPLS-AA force field parameters is treated.
(i) The diffusion in n-triacontane liquid is abnormal: the mean-squared displacement (MSD) dependence on time has the subdiffusion region. This fact can be caused by the molecular entanglement due to the molecule length. (ii) The velocity autocorrelation function (VACF) asymptote is found to be t −β . β = 2.1 in the OPLS force field with original parameterization. It differs from the hydrodynamic power law t −3/2 predicted for atomic liquids. (iii) The torsion interactions determine the flexibility of the molecules. The change of torsion H-C-C-X energy barriers reflects on the MSD and VACF. Its increase makes the molecules rigid and slows their movement. Thus, the subdiffusion region becomes longer and the value of β becomes higher, up to 2.4. The decrease of the torsion constants lets molecules move more freely and the subdiffusion region tends to disappear. The β value comes closer to the classical limit 3/2. (iv) The analytical contribution of the VACF tail to the Green-Kubo integral provides the agreement of the diffusion coefficient with the Einstein-Smoluchowski values of the diffusion. This result validates the unusual VACF asymptotes that are discovered for the C 30 H 62 liquid.