Viscosity calculations at molecular dynamics simulations

Viscosity and diffusion are chosen as an example to demonstrate the universality of diagnostics methods in the molecular dynamics method. To emphasize the universality, three diverse systems are investigated, which differ from each other drastically: liquids with embedded atom method and pairwise interatomic interaction potentials and dusty plasma with a unique multiparametric interparticle interaction potential. Both the Einstein-Helfand and Green-Kubo relations are used. Such a particular process as glass transition is analysed at the simulation of the aluminium melt. The effect of the dust particle charge fluctuation is considered. The results are compared with the experimental data.


Introduction
Molecular dynamics (MD) method [1,2] belongs to methods of the computer experiment. As in a real experiment, one is able to distinguish a part, which creates an object of research, and a part, which comprises a set of diagnostics methods to study the object created. Both MD parts are rather independent from each other at programming. Therefore, one and the same MD diagnostics methods can be applied in studies of different objects of the MD modelling, as it takes place in real experiments.
Viscosity and diffusion are chosen in this paper as an example to demonstrate the universality of the MD diagnostics methods. Two basic methods for the transport coefficients calculation are applied. The first one uses the Green-Kubo formulas which relate the transport coefficient to the integrals of the autocorrelation functions [3]. Another method exploits the Einstein-Helfand expressions [3,4] which permit to obtain transport coefficients directly from the particle displacements and velocities.
To emphasize the universality of the MD diagnostics methods, three diverse systems are investigated, which differ from each other drastically. The Lennard-Jones fluid is an example of the system with a pairwise interatomic interaction potential. Aluminium melt is a liquid with a many-particle interatomic interaction potential. Besides, such a particular process as glass transition is analysed. The third system studied is dusty plasma. It is a system with a unique multiparametric interparticle interaction potential.
Green-Kubo and Einstein-Helfand relations are given briefly for transport coefficients in section 2. Motion equations are presented in section 3, which are used to get trajectories of particles. A particular attention is paid to dusty plasmas where motion equations have a unique form [5,6]. Section 4 is devoted to some simulation details of the molecular dynamics method.
The Lennard-Jones fluid is considered in section 5. The system is used to validate the code. Glass transition is treated in section 6. The dependence of shear viscosity on temperature is calculated for the temperature range of glass transition. An example of the aluminium melt is presented. Dusty plasmas viscosity is studied in section 7. Dependencies of the shear viscosity on the charge fluctuation and friction force are investigated. Section 8 contains our conclusions.

Basic equations
Green and Kubo proved that the coefficients describing the transport properties of the system can be represented as integrals of autocorrelation functions [3]. The self-diffusion coefficient D is calculated using the formula: where ξ is a velocity auto-correlation function. There is an expression for the shear viscosity coefficient η: where V is a volume of the particle system, T is a temperature, k B is the Boltzmann constant, ... is averaging over the ensemble, P xy is the off-diagonal element of the stress tensor, which is defined by the formula: where p i x is a momentum component, u ij is an interaction potential, m i is a mass, x i is a component of a radius vector, i and j are numbers of particles.
There is another method of calculating, based on the Einstein relation. In Einstein's theory, the diffusion coefficient can be calculated by using the formula: where N is the number of particles, 0 is the reference time, r i is a radius vector of a particle. The same methodology is proposed by Helfand to calculate the shear viscosity coefficient. The coefficient is obtained from the slope of the function η(t)t: The coefficient can be determined using the following:

Particle trajectories
The dynamics of N particles is described by the Newton's equation: where F i is a force acting on the i-th particle.
The following expression is used to calculate the acceleration: where △t is a time step. The equation of dusty particles motion [5,6] takes into account the following effects: the features of gas discharge near-electrode layer, dust particle charge fluctuations, the dependence of the dust particle charge on distance from electrode, and distance from other dust particles. The trajectories of each particle are obtained using the following equation: where F i inter is an interaction force, F i trap is a force of the trap, F i f r is a frictional force of the neutral gas, F i grav is a gravitational force, F i el is a force of electric field that effects the charged dust particle near the electrode.
In the case of dusty plasma, the interaction with neutral gas particles provides the major contribution to the frictional force. So, the Langevin equation is used. The friction force can be written as follows: where γ is a coefficient of the frictional force, A is a Langevin thermostat parameter, ξ(t) is a random function.
As the equations of motion have a special form, they are not included in LAMMPS (Largescale Atomic/Molecular Massively Parallel Simulator) [7]. The code is written, which is based on the analysis of submitted particle trajectories. That allows performing calculations without a reference to LAMMPS. That means, that the coefficient of shear viscosity is obtained using the coordinates and velocities of every particle.

Simulation
The velocity autocorrelation function is obtained using the following formulas: where n is a number of time steps, τ △t is an argument of the autocorrelation function, k is a number of a time step, v i is a particle velocity. The diffusion coefficient is calculated using the formula: The same method is used to calculate the shear viscosity. The element of the stress-tensor is found using the following formula: where f iy is a component of the force acting on a i-th particle. The element of the stress tensor can be found using the Newton equation (7).
There is a problem in calculating shear viscosity coefficient. It is important to get the ensemble-averaged value. The molecular dynamics trajectory is broken up into segments. The starting point for every segment is chosen in every certain number of steps. Each segment is chosen to be statistically independent [2] and the coefficient is found for each of them. The ensemble averaging is substituted for the time averaging over reference times in the equations (2), (4) - (6).
The periodic boundary conditions are used to avoid surface effects.

Lennard-Jones system
The Lennard-Jones system is used to validate the code. The Lennard-Jones potential is described by the following equation: Liquid argon is used as an example of the Lennard-Jones system. The parameters σ = 3.4Å and ǫ/k B = 120 K are used. Calculations are performed for 256 argon atoms in a constant volume with a time step 0.05 ps.  Einstein-Helfand approach. Circles are MD results. Viscosity coefficient is defined from the slope of the dashed line.
The velocity autocorrelation function of liquid argon is presented in figure 1. The negative region of the function values indicates that it is probable that particles change their directions to the opposite due to collisions with nearest neighbours [8]. The time, when the function decays to zero, approximately equals L/2c, where L is a size of a computational cell, c is the speed of sound. The noise in the tail of the function makes the calculation numerically difficult.
The self-diffusion coefficient and shear viscosity of liquid argon at 94.4 K and 1 atm are obtained to compare with the values [9]. The results are given in table 1. Diffusivity is in a fair agreement with [9], where only Green-Kubo approach is applied. The discrepancy for the shear viscosity can be attributed to the uncertainty of the evaluation of the tail contribution in (2) and possible difference in the cut-off time at the integration.
Our results within the Green-Kubo and Einstein-Helfand approaches are close to each other for both diffusivity and viscosity. The uncertainty of the Einstein-Helfand approach figure 2 should be added to the uncertainty of the Green-Kubo approach discussed above. Figure 2 shows that the limiting slope of the η(t)t function is not defined enough accurately in our calculations.
The uncertainties of the simulation results can be the reason for the discrepancies with the experimental data. Extended simulations are needed to clear the situation. Different liquid argon models can be considered as well.

Glass transition
The molecular dynamics study of the behaviour of shear viscosity of liquid aluminium is carried out. The embedded atom method potential is used at the simulation of the aluminium melt. The expression for the potential is The first term is a sum of pair potentials ϕ over all atom pairs in the system, r ij is the distance between atoms i and j. A non-linear embedding function F (ρ j ) introduces many-body effects, ρ j is an effective electron density that is induced by the neighbouring atoms at the given atom. The state to be studied is a MD computational cell, which is only half filled by atoms. So, the pressure is not changed and remains normal, when volume is changed. The atomic initial configuration is the fcc lattice with parameter a 0 = 4.08Å. Then the system is heated to the temperature 3000 K and is equilibrated to the liquid state. The final step of the initial state preparation is lowering of the temperature and thermalization at 1500 K. Figure 3. The dependence of the shear viscosity coefficient on temperature for two values |dT /dt|: circles for 4 × 10 13 K/s, squares for 2 × 10 11 K/s. The experimental data [10] are depicted by diamonds. To investigate the glass transition, the system is cooled with a constant rate to the temperature 300 K. The rate is low enough in order not to violate the equilibrium at the instant temperature. The viscosity is calculated using the Green-Kubo formula. The code, which is based on LAMMPS [7], is used to get the dependence of the viscosity coefficient on temperature. The results for two cooling rates are presented in figure 3 together with experimental data [10].
The temperature range of the glass transition is found about 1000 K where a steep change of the shear viscosity takes place. It is in a good agreement with the glass transition region found for liquid aluminium in MD simulations [11] by means of other criteria.
The viscosities for two cooling rates coincide with each other, as well as with the experimental data [10] for temperatures above 1000 K, that is for the liquid state.
The difference between the values of viscosity appears at the temperature about 1000 K and grows with the further lowering of the temperature. The dependence of the glass viscosity for one and the same temperature and pressure on the cooling rate is an expectable effect, since an amorphous state is not an equilibrium one. Its parameters depend on a creation path. The higher is the cooling rate the higher is the viscosity.
The relation between simulation and experimental results is of the same character. The MD cooling rates available are remarkably higher than experimental rates, which are about 10 4 K/s. Correspondingly, the experimental viscosities for the amorphous state are lower than the MD ones.

Dusty plasma
Shear viscosity coefficients are calculated for the model of dusty plasma based on the equation (9). Einstein-Helfand relation is used.
As an example, the dependence of the shear viscosity coefficient on the number of particles N is shown in figure 5. The convergence of the coefficient to the limiting value is achieved, when approximately N = 5-10. The fast convergence at relatively small values of N characterizes systems with short-range interparticle interactions, such as dusty plasmas. The scatter of the values obtained in the interval N = 6-9 can be taken as a measure of the accuracy of the calculations performed.   [12].
The subsequent simulations are carried out for the system that consists of 11 particles. The dependence of the shear viscosity on the amplitude of dust particle charge fluctuation is shown in figure 5 for different coefficients of the frictional force. The coefficient affects the shear viscosity. The main reason is that the coefficient of the frictional force grows as the temperature decreases [5,6]. Therefore, due to the decrease of the temperature, the shear viscosity coefficient rises.
The viscosity does not depend on the amplitude δq/Q 0 almost in the whole area of parameters studied. The reason should be studied additionally.
The intervals of γ and δq/Q 0 presented in figure 5 cover the wide range of dusty plasmas parameters studied in experiments, [12,13] included.

Conclusions
Molecular dynamics simulation is applied to calculate shear viscosity coefficients for three diverse systems, which differ from each other drastically. The Einstein-Helfand and Green-Kubo relations are used. The results are compared with experimental and theoretical data.
1. The Lennard-Jones system is used to validate the approach. 2. The glass transition is studied in liquid aluminium. The temperature range is compared with the results, which are obtained using other criteria.
3. The importance of the dust particle charge fluctuation is shown for dusty plasmas.