Introductions

An important advantage of molecular dynamics simulations is the ability to obtain changes in the motion behavior of molecules. Among the various changes in the behavior of molecules, the diffusion behavior of molecules or atoms is the one we focus on. While the diffusion properties of molecules or atoms can be quantitatively described by diffusion coefficients, although the root mean square displacement can be fitted using the compute msd command that comes with LAMMPS software, the results of the diffusion coefficients fitted by this method are inaccurate for the following reasons:

1. The average number of msd curves outputted by the compute msd command is too small, resulting in insufficient data points required for fitting and large errors in fitting results.

2. If the root mean square displacement is output with the compute msd command, it can be found that the fitted curve fluctuates widely and does not appear as a straight line.

3.The fitting time interval is unreasonable, choosing a reasonable fitting interval is one of the keys to ensure the accuracy of the results; usually when fitting the root mean square displacement, the interception is the result of the middle time period, because at the early stage of the dynamics simulation, the parameters have not yet reached equilibrium and fluctuate greatly, while at the later stage of the simulation, the larger the correlation time, the larger the error will be.

Calculation Tutorial

The usual formula used to calculate the self diffusion coefficient is MSD(t) = n×D×Δt, (n is the number of dimensions of the particle motion space). Generally, after we calculate the MSD curve, we have to draw log-log curve for observation, at which time you will find that the curve is not always the 1st power of t. In general, the relationship between MSD and t exists as follows:

MSD=D_α·t^α

When α < 1, it is called subdiffusive behavior, where the particle is bound by other particles in the system and cannot easily escape from the bound potential well; subdiffusive can be interpreted intuitively as the wandering of the particle can be approximated as moving in the "solvent cage" formed by other particles, The MSD of the long-time motion of the particle can be approximated as the activity of the particle between the solvent cages, so it is linear with time, while the MSD of the short-time motion of the particle can be regarded as the particle has not yet come out of the first solvent cage, so it is a non-linear sub-diffusive behavior. When 1<α<2, it is called super-diffusive behavior, and generally there is an energy input to the system, which causes the particle motion to be enhanced. When α=2, it is called ballistic region, and the particle behavior is similar to uniform motion.

From the above figure, it can be seen that the MSD-t curve does not always appear as a one-time square of t throughout the simulation time period: at the beginning of the simulation, α = 2; while α gradually shifts from 2 to 1 as the simulation time continues to increase. therefore, when calculating the diffusion coefficient, we need to choose the interval for the linear fit.

Determine the linear fitting interval

In general, we have the following methods to determine the linear fitting interval.

1. the minimum covariance determinant method.

2. dynamic diffusion coefficients, i.e., observing the change in MSD/6t.

3. random sampling consistency algorithm. For example, take MSD/6t as an example, when MSD/6t is observed to be at a stable value, it can be considered that the system reaches equilibrium at this time, and this interval can be used as the fitting interval for calculating the diffusion coefficient.

Effect of electric field

1. For the strong electric field model, the particle migration coefficient cannot be calculated by Einstein's diffusion law, and the particle migration velocity is generally defined according to the change of the average coordinates of the particles with time, so it is necessary to write your own script to calculate the particle migration velocity by reading the motion trajectory of the ions.

2. For the weak electric field model, the calculation of the migration coefficient by Einstein's diffusion law is equivalent to the method described in Article 1. This is because under a finite electric field, it is possible to define the migration velocity by calculating the average coordinates of a set of particles that vary linearly with time, which is the most native calculation method; under a weak field, the two methods are equivalent, which is the so-called up-and-down dissipation theorem, but under a strong field, the two calculation methods are not equivalent.