Efficient approach to calculating radial distribution function in bcc Fe lattice

Many properties of condensed matter systems can be described by means of pair correlation functions that makes them an important structural characteristic. The shortest-graph interpolation method allows us to calculate pair correlation functions of classical crystals with pairwise interactions between particles. However, there is still no just so simple and practical approach to predict correlation functions in crystals with many-body interactions that are ubiquitous in nature. In this work, a simple modification of the interpolation method is suggested allowing to describe pair correlations bcc Fe lattice, considered as a classical crystal with many-body interactions of embedded atom model type. It is shown that the radial distribution function of the crystal can be calculated with high accuracy if mean square displacements are known. The obtained results would be useful in various fields of condensed matter physics, materials science, and crystallography.


Introduction
Radial distribution function (RDF) allows to calculate thermodynamic properties, e.g., pressure or energy, of a many-body system with pairwise interaction and contains information about its structure. In recent years, there has been developed the interpolation method (IM) [1][2][3] that allows calculating pair correlation functions of classical crystals with pairwise interactions from first-principles. The IM is based on the existence of two regimes of pair correlations: discreteness of the lattice plays a crucial role at short distances, and the correlation peaks have non-Gaussian form, in contrast, at long-range limits, correlation peaks are represented by Gaussians. The method has been proved both experimentally [4] and by means of molecular dynamic simulations [1][2][3] for various crystals with pairwise interactions. At the same time, many-body interactions are common in nature, including metals, covalent solids, oxides, and water. In this case, RDF are still of great interest, since for instance, it is demanded for proper analysis of extended X-ray absorption fine structure (EXAFS) [5][6][7][8][9] and diffraction data experiments [10][11][12][13]. There are several models taking into account many-body interactions and in case of metals, one of the most common is the embedded-atom model (EAM). This model is based on a combination of a pairwise interaction with the energy of embedding every atom in the electron density produced by its neighbors (other particles). Nowadays EAM models exist for most practically important metals and alloys [14][15][16][17][18] allowing to study various properties and phenomena in such systems by means of molecular dynamic simulations, for instance: thermodynamic [19,20], elastic [18,21], transport [18,22] and spectral [20,23] properties. However, there are still no just so simple IOP Publishing doi:10.1088/1742-6596/1697/1/012074 2 and practical approaches as IM or its modifications to describe correlations functions in crystals with many-body interactions.
In the present work, we aim to test the suitability of IM method for RDF calculation in the case of crystals with EAM-type of interactions. To reveal the role of many-body interactions in pair correlations of crystals, molecular dynamic (MD) simulation has been performed. We demonstrated that the IM allows to estimate RDF for a given temperature under the melting point. The accuracy of the method is illustrated by comparing the theoretical correlation functions with the results of MD simulations for bcc Fe lattice.

Details of simulation
We performed molecular dynamic simulation of bcc three dimensional Fe crystals, where full energy was determined by following expression: where F (ρ) is the embedding energy, ρ(r) is the atomic electron density, φ is a pair potential interaction, α and β are the element types of atoms i, j. The summations in Eq. (1) are performed over all neighbors j of atom i within the cutoff distance r c . Numerical values of F (ρ), ρ(r) and φ(r) functions were taken in accordance with [24]. All calculations were performed at different temperatures in the canonical (NVT) ensemble by using LAMMPS program [25]. The system consisted of N = 32000 particles had periodic boundary conditions. The density was taken under normal conditions (7.874 g/cm 3 ). Simulations were run for 1.5 × 10 6 time steps, both to equilibrate the system and to obtain the equilibrium properties. The numerical time step was chosen as 0.001 ps and cutoff radius was r c = 5Å. At an early stage particles velocities were chosen in accordance with Maxwellian distribution. MD simulations were used to calculate mean square displacement (MSD) and further estimation of the radial distribution function through the interpolation method.

Interpolation method
In accordance with interpolation method RDF is calculated as a summation over contributions from all particles: where p α (r) describes the contribution provided by a particle α with equilibrium spatial position r α , N is the total number of particles in the volume V . We used modified variant of the function p α (r) based on the paper [4]: where ϕ(r) is the pair interaction potential, e α = r α /|r α | is a unit vector in the direction of the particle with radius-vector r α , k B is the Boltzmann constant, T is the temperature. To find the constants C α , b α , a 2 α , and a 2 ⊥α the conditions of normalization should be applied [1]: where D is the spatial dimension, and σ 2 α , σ 2 ⊥α are longitudinal and transverse components of the MSD of the node α by relation to the initial particle and γ is the third moment.
The MSDs of all correlation peaks case can be calculated by means of finite-temperature phonon spectra: where m is the particle mass, ω j (q, T ) is the phonon spectrum of the j-th mode at the temperature T and e j is the phonon polarization's unit vector. Summation was performed over all q and j. Thus, knowledge of phonon dispersion relations allows one to estimate the MSDs. However, IM allows interpolating intermediate MSD values based on values of the closest peak and far one [3]: whereσ 2 α = σ 2 α /σ 2 1 , and parameter φ α originates from three-body correlations and has nearconstant value [3]: Calculation of temperature-dependent phonon spectrum especially of system with manybody interactions is complicated and computationally intensive task and due to this we are using required MSD of the nearest and far peaks obtained from MD.

Pair correlation functions
The RDFs calculated for EAM Fe at two different temperatures T = 300 (a) and T = 1200 (b) are shown on the Fig. 1. Red lines represent theoretical results obtained with IM, whereas blue symbols correspond to MD results. In the insets, there are the first and the second correlation peaks, that are close to each other. With the growth in temperature, the correlation peaks become thermally-smeared, but one can see that the results obtained using the IM and MD simulations agree moderately with each other at both temperatures.
The slices of the first peak in color-code format at the same temperatures as in Fig. 1 (T = 300 (a) and T = 1200 (b)) are presented in Fig. 2. Red lines correspond to the isolines of theoretical fit with IM, which are in a good agreement with the MD results. At low temperatures, the first correlation peak tends to be more gaussian. In contrast, at high temperatures, anharmonicity effects play a crucial role, and significant asymmetry of the peak along the x-direction is observed.

Conclusion
In the present paper, we have considered the problem of constructing the radial distribution function in crystals with EAM interactions. As an example, we considered 3D bcc Fe crystal and studied its RDF based on data of MD simulations. It was found out that the presence of many-body interactions (at least considered one) does not noticeable distortion of RDF form. It is shown that RDF of crystal with EAM interactions can be obtained with high accuracy by IM using only pair part of interactions and renormalized MSDs of the nearest and far peaks.
Note that, in the present work, we have been focused on the accurate reproduction of RDF form for bcc Fe lattice with EAM interactions, whereas the problem of MSD renormalization stands beyond the scope of present paper and should be considered in future works, as well as theoretical basis of the form of pair correlation peak used here. We hope that our results will be useful for further studies of systems with many-body interactions, important for condensed matter physics, crystallography, and materials science.