Temperature-dependent potential for the molecular dynamics of the superionic conductor β-PbF2

Molecular dynamic (MD) calculations were performed to investigate the thermodynamic and structural properties of lead fluoride (PbF2) by using a proposed inter-ionic temperature-dependent potential. This potential allows calculating with high precision the linear thermal expansivity and the lattice parameter as a temperature function. In addition, the potential can be represented as a sum of two contributions, a temperature-independent potential added to another temperature-dependent potential, considered last as a correction justified by the one-dimensional Newtonian quantum equation. Two fitting regions were considered, the first region from 300 to 700 K and the other one from 700 to 900 K. These regions arise naturally due to the smooth and continuous transition that PbF2 undergoes until it reaches the superionic state and, allows us to model with high precision the anomaly in the dependence of the lattice parameter with the temperature of this material, a feature that until now under the molecular dynamic method has not been studied. These results are all in good agreement with the experimental measurements.


Introduction
Computational simulations applied to ionic systems are a powerful tool that allows us to study and predict their behavior when they are subjected to external stimuli, which in most cases cannot be achieved in the laboratory. Particularly, molecular dynamics (MD) is a classical method that bases its algorithm on the resolution of Newtonian motion equations [1], which allows the study of large ionic systems at a low computational cost. The quantum effects of each ion can be disregarded (i.e., electronic structure) because the choice of an effective interaction potential, acts manifesting the characteristics as a global average (see Eq. (1)). Therefore, these potentials are derived from quantum or empirical considerations. For this reason, a right choice in the mathematical form of the potentials is crucial to guarantee a behavior like the real one. In molecular dynamics (MD), these interaction potentials are adjusted to describe the actual behavior of physical systems employing constants, and it depends explicitly on the distance between the ions or atoms. However, if a temperature gradient is applied to the system, the effect of the temperature only manifests itself as a contribution to the kinetic energy of each atom [2]. Nevertheless, it is expected that the temperature slightly modifies the interactions between the ions. Due to this, it is necessary to include temperature-dependent corrections to the interaction potential. Currently, few studies employ a methodology that reveals a dependence on interaction potentials with temperature variation. Subramaniyan and Sun [3] developed an alternative to classical molecular dynamics called engineering molecular mechanics (EMM), which allows simulating physical properties in materials at high temperatures (T > 0 K). They employed the coefficient of thermal expansion (CTE), which is a well-studied property, converting the interaction potential into a temperature-dependent potential by allowing the adjustment constants to depend on the CTE. Although the authors manage to optimize the computation time, this method is based on the molecular statics (MS) and not on the MD. Instead, Nasehzadeh et al, [4] reevaluated the adjustment parameters of the Lennard-Jones potential and improve the calculation of the thermodynamic properties of liquids. They used the Chandler-Week-Andersen perturbation theory, which provides a precise analytical expression and allows determining that both σ and ε can be expressed as a linear function of temperature. Other works that respond to more particular problems and that consider this type of dependency can be found in Refs. [5,6,7].
On the other hand, the study of ionic conductors through the use of MD is a current scientific topic, since there is still no precise understanding of the mechanisms involved in the ionic conductivity of the superionic state [8,9,10,11]. One of the most relevant superionic conductors with fluorite structure is β-PbF 2 , which has a structure described as an FCC array of lead cations (Pb 2þ ), where all the tetrahedrally coordinated interstitial sites are filled with fluorine anions (F -), being Fthe mobile species. The lattice parameter, their density and specific heat capacity are 0.593 nm, 7.750 g/cm 3 and 301 J Kg À1 K À1 , respectively. The phase transition to the superionic state, undergone by this compound, is accompanied by a λ-type peak in the specific heat, no reorganization of the crystalline frame structure, and a smooth and continuous change in conductivity as the temperature increases [12]. Likewise, this material has an anomaly in its thermal expansion when it reaches the transition temperature towards the superionic state (711 K), and there are still no studies that can describe this behavior in the lattice parameter by means of MD. In this work, the theory and methodology proposed by Nagornov and Katz [7] is used, this shows a new approach to the formulation of interaction potentials in MD by means of semiclassical considerations it is possible to express the interaction potential as a sum of an independent part of the temperature and temperature-dependent as a perturbation, as is shown in Section 3.1. Its main application is those physical phenomena, where the variation of the temperature takes a crucial role. This work validates a new methodology for obtaining this type of potentials, using crystalline structural variations with temperature for an ionic crystal-like PbF 2 .

Methods
The simulations were made using the DL_POLY Classic 1.9 package developed by Smith et al. [13] at the Daresbury Laboratory. The initial simulation was performed in a cubic simulation box of long size 23.720 Å, made with VESTA ® software, that corresponds to a 4 Â 4 Â 4 supercell, using 768 atoms, 256 lead cations and 512 fluoride anions. Artificial periodic boundary conditions have been used to reproduce the bulk properties. The initial system was equilibrated at environmental conditions (300 K, 1 atm) during 500000 steps of 2 fs that allows to compare with experimental measurements. An Intel ® Core ™ i5-7200U with turbo boost up to 3.1 GHz and 16 GB RAM PC was employed to perform the simulations.
The calculations were carried out using integration times of 2 fs and cutoff radii of 10 Å employing the Verlet method. The sum of Ewald was used to compute the Coulomb term of the potentials presented in Eq. (4). To perform the calculations of the thermodynamic properties, the system was configured under the canonical NVT ensemble, in which the temperature is increased, and subsequently, the properties in the NpT assembly are registered in order to maintained constant temperature and pressure, using an artificial Nos e -Hoover thermostat and barostat respectively which allows making a comparison between the simulated data with experimental data, during 100000 steps of 2 fs, were used to let the system stabilized until the recording process.

Results and discussion
3.1. Construction of a temperature-dependent interatomic potential

Semi-classical considerations
Following a semi-classical treatment from the Ehrenfest theorems [14], it is possible to derive the one-dimensional (1D) Newtonian quantum equation: where μ is the atom wave packet mass, x is the wave packet coordinate and U eff ðx; TÞ is an effective potential dependent on temperature, where ∂U eff ðx; TÞ ∂x ¼ ∂UðxÞ ∂x In Eq. (2) N is the number of atoms in the assembly, e γ is a normalization coefficient, g s is the number of single-particle states and ε s is the energy of the atom in the s state [7]. The effective potential is a sum of two terms, a potential that depends on the temperature and another independent, such that the temperature-dependent part varies small amounts respect to UðxÞ like: Thus, the problem is reduced to find the temperature-dependence of the interionic interaction potentials satisfying the Eq. (3).

The Buckingham potential in the β-PbF 2
The crystalline ionic solids with fluorite-like structure, such as PbF 2 , CaF 2 , BaF 2 , and so on, are often well described by the potential UðxÞ that involves the Born-Mayer-Huggins form: The first two terms on the right side of Eq. (4) are known as Buckingham's potential, being r ij the distance between the i-th and j-th ions. The exponential term represents electronic overlap (energy repulsion), the second term is proportional to À1=r 6 ij being the dispersion term, and the last term is the potential due to the Coulomb interactions.
For the β-PbF 2 , the adjustment parameters most frequently used in the literature for independent-temperature potentials are those proposed by Walker et al. [15], as shown in Table 1.
In this work, the parameter ρ PbÀF was considered as a function of the temperature, with a linear relation: In a similar way to Lennard-Jones coefficients [4] and the pre-exponential factor of the compact representation of Buckingham's potential [7]. By replacing the Eq. (5) in Uðr PbÀF Þ given by Eq. (4) obtains: The result of the right side (Eq. (6)) is obtained by expanding the argument of the exponential in a series of power of temperature and using the fact that it is possible to rewrite this equation, therefore, the temperature-dependent part is: However, Eq. (7) shows that Uðr ij ;TÞ % Uðr ij Þ AE δUðr ij ;TÞ, which is the condition imposed in Eq. (3), that the potential must satisfy.

Determination of the temperature-dependent parameter
In order to determinate the temperature dependence of the ρ PbÀF parameter, this is slightly modified from 0.490 to 0.520 eV for 300, 400, 500 and 600 K, and modified from 0.490 to 0.510 eV for 700-900 K [16], performing a new simulation set with the previous conditions. For each temperature, the ρ PbÀF is obtained to adjust the experimental lattice Table 1. Adjustment constants of the potentials that describe the β-PbF 2 by MD, as proposed by Walker et al. [15]. parameter, from the ρ PbÀF versus lattice parameter plot. Finally, the ρ PbÀF values versus the temperature are plot as shown in Figure 1, where three regions can be seen: the first at low temperatures, that is, before the superionic state, a second region up to 920 K and a third region for temperatures above this value. The calculations in the third region were not considered, due to inconsistencies in the calculations of the transport properties.
It is interesting to note how linearity is lost when passing from the non-superionic to the superionic region, which naturally justifies the need for two adjustment regions, one linear and one of a higher order > Oð1Þ. To find the best numerical adjustment, the continuity of the function ρðTÞ is waived, however, the discontinuity between both sections is so small that it can be neglected for numerical purposes (%2.275 Â 10 À4 Å). The dependence of the fitting parameter with the temperature was found in Ref [16] as: It is found that as the temperature increases, the parameter ρ PbÀF decreases and therefore the potential suffers a decrease, causing the fluoride and lead ions to experience a lower repulsion force with each other. While it might be expected that the repulsion between the ions would be greater because the crystal expands with temperature, this behavior arises as compensation for the classical thermal energy that is imparted to each ion. Likewise, a decrease in potential energy allows fluorine ions to move more easily between available sites in the lattice. Figure 2 shows the behavior of the potential of the Pb-F interaction, in the region of low (Figure 2a) and high (Figure 2b) temperatures. It is inferred that the proposed modification affects the potential in the order of tenths of eV by varying the thermal bath temperature by AE 100 K, corroborating that the Eq. (8) agrees with the Eq. (3).

Structural stability
To verify that the proposed potential does not lead to the study material undergoing a structural phase transition, or an undesired disorder in any of its sublattice with the rise of the temperature, the radial distribution functions g(r) is calculated. This function allows knowing how the ions are distributed from the atoms of the material with respect to a type of ion in a spherical shell (refer to Figure 3). The g(r) has a first welldefined maximum attributed to the first neighbors, and as the temperature increases, it loses its shape, decreasing their intensity and increasing their width, because the ions experience greater mobility as the temperature increases of the system. However, unlike lead ions (Figure 3c), the g(r) of fluorine ions (Figure 3a) changes to a radial distribution function, similar in shape to liquid materials. This behavior partly explains the superionic conductivity in this kind of material, where the sublattice of the mobile ions, has a quasi-liquid behavior due to its high mobility and diffusion, which can be corroborated in Figure 4, in which the ions of lead remain in their network sites, while those of fluorine begin to diffuse gradually through the crystal as the temperature increases. From Figure 3b, a mean distance of 2.42 Å between the lead and fluorine ions at room temperature was calculated. This distance turned out to be slightly greater than that reported by Silva et al. [17], using temperature-independent potentials, but lower than that experimentally     Figure 3b, show how the coordination distance decreases as the temperature increases. Figure 5, shows the variation of the lattice parameter when the temperature of the system, is increased and compared with the results obtained experimentally by Goff et al. [18], (in red circles), against those calculated in this paper (in cyan diamonds); when using the modification proposed in Eq. (8). A high coincidence with that reported experimentally is shown, even fitting to the anomaly present near 700 K, so the choice of two regions shows to be an appropriate process.

Thermal expansion coefficient
Once the dependency of the lattice parameter a 0 ðTÞ, with the temperature, is found, it is possible to find the linear thermal expansion coefficient α (see Eq. (9)) as [19]: A sixth-grade polynomial fit for the data obtained in this work was done, and its gradient was calculated (see Figure 6) where it is compared with experimental data of the linear thermal expansion coefficient obtained by Roberts and White [20]. It is interesting to note that the calculated curve has a maximum close to 670 K while the experimental curve around 690 K. However, both have the same line shape and order of magnitude, with an error slightly lower than 30% in area, in the range of 350-800 K. Beyond 800 K, the data deviates from experimental behavior. Also, it is the first time that this type of graph has been made by MD calculations for β-PbF 2 .

Heat capacity
The heat capacity (C p ) of any material is defined as the amount of heat per unit of mass required so that its temperature raises 1 K [21]. For the molecular dynamic case, C p is calculate by: Being n the number of moles, H the enthalpy computed in the NpT ensemble of the system, and T is the crystal temperature. The C p curve calculated from Eq. (10), from 300 K to 800 K is plotted in Figure 7. The C p increases in a wide region around 700 K, near the superionic transition, predicting a phase transformation to the superionic state due to an increase in internal energy, related to the mobility of mobile ions. In previous study [11], the C p values employing only the Walker's potential showing a higher value (80.49 J K À1 mol À1 ) than the experimental one (76.81 AE 0.51 J K À1 mol À1 at 400 K) and above 700 K a chaotic behavior arise, without exhibiting a phase transition. In this work, a lower value C p , than the experimental is calculated (52.3 J mol À1 K À1 at 500 K against 79.87 J mol À1 K À1 measured by Volodich et al., [22]), and a phase transition is predicted. The difference between the values (experimental and theoretical) is the non-existence of a nominal value in the literature for temperatures above room temperature. However, a similarity in the shape of the lines can be appreciated. From temperatures above 800 K, a difficulty arises in fitting the enthalpy data.

Conclusions
In summary, the potential proposed by Walker et al. [15], was modified, to describe the structural properties of β-PbF 2 as a temperature function, by making this equation parametrically dependent on temperature. The fitting constant ρ PbÀF was replaced by a linear function ρ PbÀF ðTÞ. This parameter was fitted with experimental data of the lattice parameter versus temperature. It was found that two adjustment regions naturally emerged, the first one from 300 to 700 K and another from 700 to 900 K. The first is characterized by having a linear behavior, while the second exhibits a second-order polynomial form. From these functions, calculations by molecular dynamics were made, and it is possible to conclude that there exists a high coincidence between the data of the evolution of the lattice parameter with the temperature and those calculated with our proposed potential, it even fits the anomaly present in the transition temperature towards the superionic state. From the g(r) function obtained, it is possible to conclude that the proposed potential does not involve a structural phase transition, and the ions distance is almost equal to the experimental one. The C p curve is firstly for this material calculated from 300 to 800 K, exhibiting a growing line shape, until a maximum; then, around it, a phase transition towards the superionic state is predicted.   García G.: Conceived and designed the experiments; Contributed materials, analysis tools or data. Mosquera E.: Analyzed and interpreted the data; Wrote the paper.

Funding statement
This work was supported by Ministerio de Ciencia Tecnología e Innovaci on, Colombia(Grant no. C.I. 1165).