Analyzing Melts and Fluids from Ab Initio Molecular Dynamics Simulations with the UMD Package

1. Analysis of the molecular-dynamics runs

NOTE: The package is available via the GitHub website (https://github.com/rcaracas/UMD_package) and via a dedicated page (http://moonimpact.eu/umd-package/) of the ERC IMPACT project as an open-access package.

1. Extract each specific set of physical properties using one or more dedicated Python scripts from the package. Run all the scripts at the command line; they all employ a series of flags, which are as consistent as possible from one script to another. The flags, their meaning, and default values are all summarized in Table 1.

Table 1: Most common flags used in the UMD package and their most common significance.

2. Start with transforming the output of the MD simulation performed in a first-principles code, like VASP⁸ or QBox⁹, into a UMD file.
   1. If the MD simulations were done in VASP, then at the command line type:
      VaspParser.py -f <OUTCAR_filename> -i <InitialStep>
      where –f flag defines the name of the VASP OUTCAR file, and the –i the thermalization length.
      NOTE: The initial step, defined by –i allows for discarding the first steps of the simulations, which represent the thermalization. In a typical molecular-dynamics run, the first part to the calculation represents the thermalization, i.e., the time it takes the system for all atoms to describe a Gaussian-like distribution of the temperature, and for the entire system to exhibit fluctuations of the temperature, pressure, energy, etc. around equilibrium values. This thermalization part of the simulation should not be taken into account when analyzing the statistical properties of the fluid.
3. Transform the .umd files into .xyz files to facilitate visualization on various other packages, like VMD⁴ or Vesta⁵. At the command line type:
   umd2xyz.py -f <umdfile> -i <InitialStep> -s <Sampling_Frequency>
   where –f defines the name of the .umd file, –i defines the thermalization period to be discarded, and –s the frequency of the sampling of the trajectory stored in the .umd file. Default values are –i 0 –s 1, i.e. considering all the steps of the simulation, without any being discarded.
4. Reverse the umd file into VASP-type POSCAR files using the umd2poscar.py script; snapshots of the simulations can be selected with a predefined frequency. At the command line type:
   umd2poscar.py -f <umdfile> -i <InitialStep> -l <LastStep> -s <Sampling_Frequency>
   where –l represents the last step to be transformed into POSCAR file. Default values are -i 0 -l 10000000 -s 1. This value of –l is big enough to cover a typical entire trajectory.

2. Perform the structural analysis

1. Run the gofrs_umd.py script to compute the pair distribution function (PDF) g_ᴀʙ(r) for all the pairs of atomic types A and B (Figure 3). The output is written in one ASCII file, tab-separated, with the extension gofrs.dat. At the command line type:
   gofrs_umd.py -f <UMD_filename> -s < Sampling_Frequency > -d <DiscretizationInterval> -i <InitialStep>
   NOTE: The defaults are Sampling_Frequency (the frequency for sampling the trajectory) = 1 step; DiscretizationInterval (for plotting the g(r)) = 0.01 Å; InitialStep (number of steps in the beginning of the trajectory that are discarded) = 0. The radial PDF, g_ᴀʙ(r) is the average number of atoms of type B at a distance d_ᴀʙ within a spherical shell of radius r and thickness dr centered on the atoms of type A (Figure 3):
   
   
   with ρ the atomic density, N_A and N_B the number of atoms of type A and B, and δ(r−r_ᴀʙ) the delta function which is equal to 1 if the atoms A and B lie at a distance between r and r+dr. The abscissa of the first maximum of g_ᴀʙ(r) gives the highest probability bond length between the atoms of type A and B, which is the closest to an average bond distance that we can determine. The first minimum delimits the extent of the first coordination sphere. Hence the integral over the PDF up to the first minimum gives the average coordination number. The sum of the Fourier transforms of the g_ᴀʙ(r) for all the pairs of atomic types A and B yields the diffraction pattern of the fluid, as obtained experimentally with a diffractometer. However, in reality, as oftentimes the high order coordination spheres are missing from the g_ᴀʙ(r), the diffraction pattern cannot be obtained in its entirety.

Figure 3: Determination of pair distribution function.
(a) For each atom of one species (for example red), all the atoms of the coordinating species (for example grey and/or red) are counted as a function of distance. (b) The resulting distance distribution graph for each snapshot, which at this stage is only a collection of delta functions, is then averaged over all the atoms and all the snapshots and weighted by the ideal gas distribution to generate (c) the pair distribution function that is continuous. The first minimum of the g(r) is the radius of the first coordination sphere, used later on in the speciation analysis. Please click here to view a larger version of this figure.

2. Extract the average interatomic bond distances as the radii of the first coordination spheres. For this, identify the position of the first maximum of the g_ᴀʙ(r) functions: plot the gofrs.dat file in a spreadsheet application and search for the maxima and minima for each pair of atoms.
3. Identify the radius of the first coordination sphere, as the first minimum of the PDF, g_ᴀʙ(r), using spreadsheet software. This is the basis for the entire structural analysis of the fluid; the PDF yields the average bonding status of the atoms in the fluid.
4. Extract the distances of the first minima, i.e., the abscissa, and write them in a separate file, called, for example, bonds.input. Alternatively, run one of the analyze_gofr scripts of the UMD package to identify the maxima and the minima of the g_ᴀʙ(r) functions. At the command line type:
   analyze_gofr_semi_automatic.py
5. Click on the position of the maximum and the minimum of the g_ᴀʙ(r) function displayed in the graph that is opened by the program. The script automatically scans the current folder, identifies all the gofrs.dat files, and performs the analysis for each one of them. Click again on the maximum and the minimum in the window every time the script needs an educated initial guess.
6. Open and look at the automatically generated file called bonds.input that contains the interatomic bond distances.

3. Perform the speciation analysis

1. Compute the topology of bonding between the atoms, using the concept of connectivity within graph theory: the atoms are the nodes and the interatomic bonds are the paths. The speciation_umd.py script needs the interatomic bond distances defined in the bonds.input file.
   NOTE: The connectivity matrix is constructed at each time step: two atoms that lie at a distance smaller than the radius of their corresponding first coordination sphere are considered to be bonded, i.e., connected. Various atomic networks are built by treating the atoms as nodes in a graph whose connections are defined by this geometric criterion. These networks are the atomic species, and their ensemble defines the atomic speciation in that particular fluid (Figure 4).

Figure 4: Identification of the atomic clusters.
The coordination polyhedra are defined using interatomic distances. All atoms at a distance smaller than a specified radius are considered to be bonded. Here the threshold corresponds to the first coordination sphere (the light red circles), defined in Figure 1. Polymerization and thus chemical species are obtained from the networks of the bonded atoms. Note the central Red₁Grey₂ cluster, which is isolated from the other atoms, which form an infinite polymer. Please click here to view a larger version of this figure.

2. Run the speciation script to obtain the connectivity matrix and obtain the coordination polyhedra or the polymerization. At the command line type:
   speciation_umd.py -f <UMD_filename> -s <Sampling_Frequency> -i <InputFile> -l <MaxLength> -c <Cations> -a <Anions> -m <MinLife> -r <Rings>
   where the -i flag gives the file with the interatomic bond distances, which was produced for example in the previous step. Alternatively, run the script with one single length for all bonds defined by the -l flag.
   NOTE: The -c flag specifies the central atoms, and the -a flag the ligands. Both central atoms and ligands can be of different types; in this case, they must be separated by commas. The -m flag gives the minimum time a species must live to be considered in the analysis. By default this minimum time is zero, all occurrences being counted in the final analysis.
   1. Run the speciation_umd.py script with the flag –r 0, which samples the connectivity graph at the first level to identify the coordination polyhedra. For example, a central atom, denoted as a cation may be surrounded by one or more anions (Figure 4). The speciation script identifies every single one of the coordination polyhedra. The weighted average of all the coordination polyhedra gives the coordination number, identical to the one obtained from the integration of the PDF. At the command line type:
      speciation_umd.py -f <UMD_filename> -i <InputFile> -c <Cations> -a <Anions> -r 0
      NOTE: Average coordination numbers in fluids are fractional numbers. This fractionality comes from the average characteristic of the coordination. The definition based on speciation yields a more intuitive and informative representation of the structure of the fluid, where the relative proportions of the different species, i.e., coordinations, are quantified.
   2. Run the speciation_umd.py script with the flag –r 1, which samples the connectivity graph at all depth levels to obtain the polymerization. The network through the atomic graph has a certain depth, as atoms are bonded further away to other bonds (e.g., in sequences of alternating cations and anions) (Figure 4).
3. Open the two files .popul.dat and .stat.dat consecutively; these constitute the output of the speciation script. Each cluster is written on one line, specifying its chemical formula, the time at which it formed, the time at which it died, its lifetime, a matrix with the list of the atoms forming this cluster. Plot the lifetime of each atomic cluster of all the chemical species found in the simulation as found in the .popul.dat file (Figure 5).
4. Plot the population analysis with the abundance of each species, as found in the .stat.dat file. This analysis, both absolute and relative, corresponds to the actual statistics of the coordination polyhedra for the case -r 0; for the case of polymerization, with -r 1 this needs to be treated carefully as some normalization over the relative number of atoms might need to be applied. The abundance corresponds to the integral over the lifetimes. The .stat.dat file also lists the size of each cluster, i.e., how many atoms form it.

4. Compute diffusion coefficients

1. Extract the mean square displacements (MSD) of the atoms as a function of time to obtain the self-diffusivity. The standard formula of the MSD is:
   
   where the prefactors are renormalizations. With the MSD tool, there are different ways to analyze the dynamical aspects of the fluids.
   NOTE: T is the total time of the simulation and N_α is the number of atoms of type α. The initial time t₀ is arbitrary and spans the first half of the simulation. N_init is the number of initial times. τ is the width of the time interval over which the MSD are computed; its maximum value is half the time length of the simulation. In typical MSD implementations, every window starts at the end of the previous one. But a sparser sampling can speed up the computation of the MSD, without altering the slope of the MSD. For this, the i-th window starts at time t₀(i), but the (i+1)-th window starts at time t₀(i) + τ + v, where the value of v is user-defined. Similarly, the width of the window is increased in discrete steps defined by the user, as such: τ(i) = τ(i-1) + z. The values of z (“horizontal step”) and v (“vertical step”) are positive or zero; the default for both is 20.
2. Compute the MSD using the series of msd_umd scripts. Their output is printed in a .msd.dat file, where the MSD of each atomic type, atom, or cluster is printed on one column as a function of time.
   1. Compute the average MSD of each atomic type. The MSD are computed for each atom and then averaged for each atomic type. The output file contains one column for each atomic type. At the command line type:
      msd_umd.py -f <UMD_filename> -z <HorizontalJump> -v <VerticalJump> -b <Ballistic>
   2. Compute the MSD of each atom. The MSD are computed for each atom and then averaged for each atomic type. The output file contains one column for each atom in the simulation, and then one column for each atomic type. This feature allows for identifying atoms that diffuse in two different environments, like liquid and gas, or two liquids. At the command line type:
      msd_all_umd.py -f <UMD_filename> -z <HorizontalJump> -v <VerticalJump> -b <Ballistic>
   3. Compute the MSD of the chemical species. Use the population of clusters identified with the speciation script, and printed in the .popul.dat file. The MSD are computed for each individual cluster. The output file contains one column for each cluster. To avoid considering large-scale polymers, place a limit on the size of the cluster; its default is 20 atoms. At the command line type:
      msd_cluster_umd.py -f <UMD_filename> -p <POPUL_filename> -s <Sampling_Frequency> -b <Ballistic> -c <ClusterMaxSize>
      NOTE: Defaults values are: –b 100 –s 1 –c 20.
3. Plot the MSD using a spreadsheet-based software (Figure 6). In a log-log representation of the MSD versus time, identify the slope change. Separate the first part, usually short, which represents the ballistic regime, i.e. the conservation of the velocity of atoms after collisions. The second longer part represents the diffusive regime, i.e. scattering of the velocity of atoms after collisions.
4. Compute the diffusion coefficients from the slope of the MSD as:
   
   where Z is the number of degrees of freedom (Z = 2 for diffusion in plane, Z = 3 for diffusion in space), and t is the time step.

5. Time correlation functions

1. Compute the time correlation functions as a measure of the inertia of the system using the general formula:
   
   A can be a variety of time-dependent variables, such as the atomic positions, atomic velocities, stresses, polarization, etc., each yielding—via the Green-Kubo relations^12,13—different physical properties, sometimes after a further transformation.
2. Analyze the atomic velocities to obtain the vibrational spectrum of the liquid and alternative expression of the atomic self-diffusion coefficients.
   1. Run the vibr_spectrum_umd.py script to compute the atomic velocity-velocity auto-correlation (VAC) function for each atomic type and perform its fast-Fourier transform. At the command line type:
      vibr_spectrum_umd.py -f <UMD_filename> -t <temperature>
      where –t is the temperature that must be defined by the user. The script prints two files: the .vels.scf.dat file with the VAC function for each atomic type, and the .vibr.dat file with the vibrational spectrum decomposed on each atomic species and the total value.
   2. Open and read the vels.scf.dat. Plot the VAC function from the vels.scf.dat file using spreadsheet-like software.
   3. Keep the real part of the Fourier VAC. This is what yields the vibrational spectrum, as a function of frequency:
      
      where m are the atomic masses.
   4. Plot the vibrational spectrum from the vibr.dat file using spreadsheet-like software (Figure 7). Identify the finite value at ω=0 that corresponds to the diffusive character of the fluid and the various peaks of the spectrum at finite frequency. Identify the participation of each atomic type to the vibrational spectrum.
      NOTE: The decomposition on atomic types shows that different atoms have different ω=0 contributions, corresponding to their diffusion coefficients. The general shape of the spectrum is much smoother with fewer features than for a corresponding solid.
   5. At the shell, read the integral over the vibrational spectrum, which yields the diffusion coefficients for each atomic species.
      NOTE: Thermodynamic properties can be obtained by integration from the vibrational spectrum, but the results should be used with caution because of two approximations: the integration is valid within the quasi-harmonic approximation, which does not necessarily hold at high temperatures; and the gas-like part of the spectrum corresponding to the diffusion needs to be discarded. The integration should then be done only over the lattice-like part of the spectrum. But this separation usually requires several further post-processing steps and calculations¹⁴, which are not covered by the present UMD package.
3. Run the viscosity_umd.py script to analyze the self-correlation of the components stress tensor to estimate the viscosity of the melt. At the command line type:
   viscosity_umd.py -f <UMD_filename> -i <InitialStep> -s <Sampling_Frequency> -o <frequency_of_origin_shift> -l <max_correlation_timelength>
   NOTE: This feature is exploratory and any results must be taken with caution. Firstly, thoroughly check the convergence of the viscosity with respect to the length of the simulation.
   1. Derive the viscosity of the fluid from the self-correlation of the stress tensor¹⁵ as:
      
      where V and T are the volume and the temperature respectively, κ_B is the Boltzmann constant and σ_ij the ij off-diagonal component of the stress-tensor, expressed in Cartesian coordinates.
   2. Use a more adequate fit to obtain a more robust estimate of the viscosity^15,16 and avoid the noise of the stress-tensor auto-correlation function that might arise from the finite size and the finite duration of the simulations. For the auto-correlation function of the stress tensor, use the following functional form^15,16 that yields good results:
      
      where A, B, τ₁, τ₂, and ω are fit parameters. After integrating, the expression for the viscosity becomes:
      

6. Thermodynamic parameters stemming from the simulations.

1. Run averages.py to extract the average values and the spread (as standard deviation) for pressure, temperature, density, and internal energy from the umd files. At the command line type:
   averages.py -f <UMD_filename> -s <Sampling_Frequency>
   with –s 0 as default.
2. Compute the statistical error of the average, using blocking methods.
   NOTE: There are various flavors of this method. Following the work of Allen and Tildesley², it is common to average over sequences of time blocks, of increasingly longer length, and estimate the standard deviation with respect to the arithmetic average¹⁷. Convergence may be reached in the limit of many- and long-enough block sizes, when the sampling is uncorrelated. Though the actual threshold value for the convergence usually needs to be chosen manually.
   1. Use the halving method¹⁸: starting with the initial data sample, at each step κ, halve the number of the samples by averaging over every two corresponding consecutive samples from the previous step κ−1:
      
   2. Run the fullaverages.py script to perform the complete statistical analysis, including the error of the mean. At the command line type:
      fullaverages.py -s <Sampling_Frequency> -u <units>
      NOTE: The script is automatized to the point of searching for all the .umd.dat files in the current directory and performing the analysis for all of them. Defaults are –s 0 –u 0. For -u 0 the output is minimal, and for -u 1 output is in full, with several alternative units printed out. This script requires graphical support, as it creates a graphical image for checking the convergence for estimating the error on the mean.