Modified Born method for modeling melting temperature using ab initio molecular dynamics

The prediction of a material’s melting point through computational methods is a very difficult problem due to system size requirements, computational efficiency and accuracy within current models. In this work, we have used a newly developed metric to analyze the trends within the elastic tensor elements as a function of temperature to determine the melting point of Au, Na, Ni, SiO2 and Ti within ±20 K. This work uses our previously developed method of calculating the elastic constants at finite temperatures, as well as leveraging those calculations into a modified Born method for predicting melting point. While this method proves to be computationally expensive, the level of accuracy of these predictions is very difficult to reach using other existing computational methods.


Introduction
It is possible to use the elastic and mechanical behavior of materials to extract a range of various material properties including information related to mechanical strength [1][2][3][4][5][6], stability [1,3,4,6,7], and phase transitions [1,5,6,8]; an improved understanding of hardness [4,6], stiffness [2,[5][6][7]9] and ductility [6,7,[9][10][11][12][13] can also be obtained .It is necessary to gather sufficient information regarding the elastic constants of a material to determine these additional quantities.However, procuring these constants experimentally is often difficult over a large temperature range due to the harshness of the environment.This makes it necessary to continue improving theoretical techniques to determine these quantities * Author to whom any correspondence should be addressed.
from first principles calculations.The inclusion of temperature effects when computing material properties is a broad field with an extensive collection of methods available for use.Some of these techniques include: the exact muffin-tin orbitals method [14], transitive correlation methodology [15], first principles quasistatic approach [16][17][18], Debye-Gruneisen theory [19], and Born-Oppenheimer molecular dynamics [20].This is not an exhaustive list of all methodologies, rather a collection of methods that show the wide variety of techniques being used, to emphasize the need for a more concise and universal approach.
A particular issue with many of the free-energy methods mentioned above (including the hysteresis and coexistence methods) is the distinction between the equilibrium melting point and the superheated melting point [50][51][52].The equilibrium melting point (or heterogeneous melting point) is commensurate with bulk melting due to density fluctuations when approaching a phase transition, while the superheated melting point is when a crystal is heated well above the necessary melting temperature and the effects (system disorder, changes in energy, altered local density, etc) can be computed.Superheated melting points are often what computational modeling captures due to the challenge in running trajectories over sufficient timescales for equilibrium melting to occur [51,52].Another difficulty involved in making this determination is system complexity [50].It has been shown that more complex structures (systems aside from single element bulk materials) do not experience melting in a heterogeneous manner due to differing activation barriers for the process, with varying amounts of time required to melt within the system [53][54][55][56].This requires that each temperature be closely monitored for a significant amount of time to determine if the conditions are suitable for melting.
In 1939 Max Born proposed using the temperaturedependent elastic constants as a method for determining the melting point of various materials [57].The underlying hypothesis behind this model is that a solid has resistance to shearing stress due to its elasticity; a liquid does not have this capability.Born developed three inequalities to help identify melting, with C 44 > 0 being the most vital (where C 44 can be calculated from C 11 and C 12 , depending on symmetry [58]).A violation of this inequality signifies that the crystal cannot withstand shearing, which denotes melting within the Born model.These criteria are considered to be a 'triggering mechanism' for homogeneous melting, as opposed to an indicator that melting has occurred [59][60][61][62][63]. Since then, others have postulated the possibility of using the elastic tensor to assess melting [64,65], but this method has not been examined rigorously through modern computational methods.Current methodologies now tend towards the inclusion of vacancies within the crystal structure in melting point simulations [51], indicating need for nucleation sites for complete melting to occur [61-63, 66, 67].However, AIMD simulations allow for crystals to be held at elevated temperatures without providing these interfaces as nucleation sites [61][62][63].
In this work, we use first principles density functional theory (DFT) calculations and AIMD simulations to study the temperature dependence of the elastic constants, intending to further investigate the postulates of Born with modern techniques.Materials investigated include gold (Au), sodium (Na), nickel (Ni), β-quartz (SiO 2 ) and titanium (Ti).We have shown in our previous work that we are able to capture the variation in the elastic tensor elements of other materials using molecular dynamics simulations above and below the known melting point of the material [68] (as well as phase transitions), with this work being focused at the melting point.

Computational methods
For each of the materials studied in this work, their respective unit cells were used to create a supercell of approximately 100 atoms.For Au and Ni, face-centered cubic unit cells (Fm 3m) were used to create supercells of 108-atoms.SiO 2 has a complex hexagonal unit cell (I 42d) which led to a supercell of 96-atoms.A 108-atom supercell of Ti was generated from the hexagonal close-packed unit cell (R 3m).The authors note that, experimentally, Ti undergoes a phase transition to body-centered cubic (BCC) before melting.Finally, the BCC structure of Na (Im 3m) was used to build a 128-atom supercell.These crystal structures are shown in figure 1, and the specific crystallographic information can be found in tables S1-S5 (see supplementary material).
Each of the investigated materials was fully relaxed using DFT [69,70], as implemented in the Vienna ab initio Simulation Package (VASP) [71][72][73].VASP uses a plane wave basis set, and ion-electron interactions were approximated using the projector augmented wave potentials [74].The kinetic energy cutoff was tested with the criterion for convergence being an energy variation of less than 1 meV/atom.For the DFT calculations on the unit cells, 17 × 17 × 17 (Au), 11 × 11 × 11 (Ni), 5 × 5 × 5 (SiO 2 ), 9 × 9 × 9 (Ti) and 14 × 14 × 14 (Na) Monkhorst-Pack k-point grids were found to be sufficient.The Methfessel-Paxton smearing method was used for sampling the Brillouin zone [75], and all structural relaxations and energy calculations were performed using the Perdew-Burke-Ernzerhof functional [76].Plane wave energy cutoffs of 400 eV (Ti, Na and Au) and 450 eV (Ni and SiO 2 ) were found to be sufficient for the materials discussed in this work.
We have elected to use only Γ-point sampling of the reciprocal space [77] for the AIMD simulations of all supercells studied in this work.It has been shown that Γ-point sampling is adequate for systems of this size [68,78].The band curvature is negligible when the dimensions of the supercell are notably larger than the interatomic spacing [78].This method facilitates precise calculation of material properties such as charges, energies, forces and mechanical properties while drastically decreasing computational expense [68,78,79].
AIMD simulations were performed to include temperature effects within our systems.Two steps of equilibration were performed for each temperature considered: volumetric equilibration followed by simultaneous equilibration of the ionic positions and internal pressure.Volumetric equilibration was done using the isothermal-isobaric ensemble (NPT) and allowing the system to run for 5 picoseconds (ps) with time steps of 1 femtosecond (fs) (5000 steps) at each temperature of interest.A smaller time step is used for this part of the equilibration to ensure that noise within the volumetric fluctuations remains well-behaved [68].For the next set of equilibrations, the canonical ensemble (NVT) was used to allow the ionic Crystallographic information for each structure can be found in the supplemental material.positions and internal pressure to equilibrate.This equilibration was done for 20 ps with time steps of 2 fs (10 000 steps).
For calculation of the elastic tensor elements, we take the fully equilibrated system and apply a series of strains (when applying strain, we use values of ±0.5% and ±1.0%).Under the canonical ensemble each strained system is allowed to run for 20 ps with time steps of 2 fs (10 000 steps).The internal stress of the system is calculated at each time step over the duration of these simulations.Since we have a known applied strain, and we directly measure the stress, we can then use the generalized stress-strain variant of Hooke's Law to determine the elements of the elastic tensor (C ij ): where σ i and ϵ j are the tensile stress and longitudinal strain, respectively.For a more detailed description of the AIMD simulations performed, please see our previously published work [68].
The most seemingly straightforward way of computationally determining a material's melting point is by performing molecular dynamics in which a perfect material is heated from a solid phase and finding the temperature at which the crystalline lattice breaks down.However, owing to the energy required to create a solid-liquid interface the material becomes superheated rather than melts, and the predicted melting point is well-above the experimental one.Similarly, the same process could be done except by cooling the material from a liquid phase, although in this case the material becomes supercooled rather than crystallizing.Thus, the melting point predicted from heating and cooling simulations are different, leading to the so-called hysteresis effect.It has been shown that hysteresis is dependent on the heating (or cooling) rate of the system due to sudden changes in volume related to fluctuations in temperature [21,80].This leads to the prediction of the superheated or supercooled melting point of the solid.In this work, however, we postulate that we do not observe this effect because of our two-stage equilibration process.Rather than use a heating or cooling rate, we perform full equilibration of our system (including the volume) at each temperature before calculation of the elastic tensor elements.This two-stage equilibration approach allows for sufficient amorphization of the crystalline lattice at elevated temperatures.We find negligible fluctuations of temperature in these simulations, resulting in almost zero changes in energy and internal pressure throughout the calculation of the tensor elements.

Results and discussion
Because it is known that the elastic constants of Ni are well behaved as a function of temperature, it was chosen as the initial material of interest in this work [81].The premise for these calculations is that elasticity, ductility and plasticity are not meaningful once the temperature exceeds the melting point of a material.Thus, the elements of the elastic tensor at different temperatures may provide insight into the melting point.In his early computations involving melting point, Born provides the trends of the elastic tensor elements as a function of temperature: the constants will decrease linearly as temperature increases, followed by a steep downward curve as melting occurs [57].For the first material discussed in this work (Ni), a blind study was performed where lower and upper temperature boundaries were chosen (300 K and 2100 K, respectively).Using changes in the calculated values of the elastic tensor, we determined if the material was above or below its melting point.Next, a semi-Newtonian method was used to determine the next temperature to be studied, meaning, if it was found that the melting point lies within the temperature range that was chosen, the next set of calculations was performed at the median temperature value.This process was continued until a final temperature range was determined.
The results of these simulations for C 1j (j = 1, 2, 3) are shown in figure 2. Starting at 300 K, there is a slight increase in each of the tensor elements as the temperature reaches 1200 K, followed by a significant decrease in these values as they approach the melting point.While there are fewer data points above the melting temperature, this sudden change in slope can be seen for simulations on both sides of the melting temperature.With these data, we were able to narrow the melting point to a range of 50 K (1700-1750 K).For all other materials in this study the experimental melting point was obtained in advance in order to minimize the number of temperatures needed for determining the melting point, thereby reducing computational cost.
Instead of using the Newtonian approach to investigate the smaller temperature window of Ni, calculations were performed at 10 K increments across the previously predicted range; the results are shown in the smaller image inset in figure 2. In this plot, we see that there is a significant discontinuity for C 11 at 1720 K. Similarly, for C 12 there is a noticeable change in that element at 1730 K.For C 13 there is significant discontinuity at both of these temperatures.The combination of discontinuities at these points leads us to a melting point range of 1720-1740 K.For all other materials investigated, it is this series of discontinuities that we seek to indicate that our material is experiencing a different mechanical response, indicating that melting could become favorable.With the known melting point of Ni being 1728 K, this provides a range well within the margin of error for AIMD simulations and is more accurate than methods used in previous work, which predicted the melting point of Ni to be 1637 ± 100 K [82].
Upon completion of the elastic tensor calculations of Ni, this method was set for testing using four additional materials: β-quartz (SiO 2 ), Ti, Na and Au.For these materials, the semi-Newtonian method was not used to find temperatures of interest.Instead, the material melting points were looked up beforehand and temperatures were picked so that a sufficient number of data points would be collected below the melting point, above the melting point, with a high density of data points (in 10 K increments) chosen around the melting temperature.The resulting elastic tensor elements as a function of temperature are shown in figure 3, where the inset Temperature dependence of the atomic coordination environment of Ti atoms using the common neighbor analysis tool [83,84] in the Open Visualization Tool (OVITO) [85].Atoms are found to be in the hexagonal close packed (HCP), body-centered cubic (BCC) and face-centered cubic (FCC) structures.Any atoms not found to be in one of these crystal environments are grouped in the 'Other' category.
images show the elastic tensor elements over smaller temperature ranges.
For the larger temperature ranges shown in figure 3, the general trends are similar (with one exception in Ti to be discussed below).At points far below the melting point, consecutive data points of all tensor elements show small decreases in the calculated quantities, as is expected [68].The exception to this is the elastic tensor elements of Ti from 800 K to 1200 K, which show more large changes in two of the three tensor elements.Ti undergoes a phase transition to a BCC structure at 1156 K, which would explain the difference in C 1j values at these temperatures.To further assess this, we have investigated the change in atomic coordination environment of the Ti system as a function of temperature using the common neighbor analysis tool [83,84] in the Open Visualization Tool (OVITO) [85].Common neighbor analysis is readily used to classify the structure within crystalline systems by comparing the local structure via atomic pair coordination.
In figure 4 we show that the initial DFT relaxation of the Ti supercell finds all of the atoms in the hexagonal close packed (HCP) crystal structure.All atoms except for one remain in this structure up to 800 K, where one atom is found in the BCC environment.At 1200 K the common neighbor analysis provides less consistent results with a large number of atoms found in the HCP and BCC structures, but the majority of atoms are in positions that do not fit within the crystal environments defined by OVITO.This gives credence to the fact that many of the atoms are moving between the BCC and HCP environments at this temperature.As the temperature further increases, the number of atoms found in the HCP structure rapidly approaches zero.The number of atoms in the BCC structure stays constant until 1800 K, and then decreases sharply.This indicates that phase transitions are also captured using this approach.
As each material approaches their respective melting temperatures the changes in C 1j become larger in magnitude, with exceptionally sporadic behavior within the smaller temperature window investigated for each materials (see inset images within figure 3).These fluctuations do not appear systematic for any of the materials studied within this work, however the seemingly random behavior of these calculated quantities for each material encourage further analysis of the combination of the these results.
With the goal of adopting an objective metric for assigning the melting temperature rather than relying on a qualitative interpretation of the resulting elastic tensor elements, we introduce the quantity Z total to assess all changes in the elastic tensor simultaneously, where each component (Z 11 , Z 12 and Z 13 ) is computed as, Here, T is the current temperature being considered, we define T p as the next lowest temperature to the one being investigated and T min is the lowest temperature investigated for that material.Z is meant to quantify the relative change in tensor element C 1j over a given temperature difference.Z total = Z 11 = Z 12 = Z 13 = 0 at T min of each material since there is no lower temperature to use as reference.As shown in figures 2 and 3, C 11 tends to be significantly larger in magnitude than C 12 and C 13 , Z total equally combines the fluctuations from each constituent to provide a single meaningful quantity.The results for all investigated materials are shown in figure 5.
Figures 5(a)-(e) show the calculated values of Z 11 , Z 12 , Z 13 and Z total for each material studied, with figure 5(f) showing Z total of all materials in a single image.With the exception of Na, we find that all Z values are exceedingly close to zero far from the melting point.We expect that Na does not follow this trend is because the individual elements of the elastic tensor are much lower in magnitude than all other materials studied, therefore the small fluctuations are amplified when comparing among the individual data points.The melting point of Na is also much lower than all other materials studied, so the changes in temperature between the outlying calculations is not as large.However, we do find that for all materials (including Na) as we approach the experimentally known melting point, the values of Z increase substantially with peaks in Z total that are in excellent agreement with the melting temperature.The resulting melting range for each material is shown in table 1, including the experimentally determined melting point of each material.Where the initial analysis of the elastic tensor elements is a bit more obscure, this Z measure which is derived from those elements, provides a material-agnostic way to determine a material's melting point.

Summary and conclusions
In this work, we have presented an updated version of the Born method for determining the melting point of a material using the elements of the elastic tensor.With each of the materials studied, we have shown that the behavior of the elastic tensor elements becomes much more sporadic when approaching melting; in particular, we show that there are drastic fluctuations in neighboring points, and the typical trend of a (relatively) constant decrease in magnitude of the elastic tensor elements is no longer maintained.These discernible fluctuations make it possible to identify the temperature range when the conditions are right for melting to occur.We have also introduced the Z-metric to provide a single quantity that can be used to determine the melting temperature objectively by showing at which point the changes in the elastic constant become more substantial.While these calculations are computationally expensive, this level of accuracy is remarkably difficult to obtain from first principles.We believe that larger system size and analysis of more tensor elements would provide improved accuracy of these results, which may be possible with the continued development of high-performance computing and algorithms.Although this technique was intended to compute the melting point of known materials, it could potentially be used to investigate more complex systems (alloys, complex oxides, etc) with unknown melting points to further improve on the material design process within various industry applications.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Figure 1 .
Figure 1.Unit cells for each of the structures studied in this work.These structures were used to create supercells of (a) Au, (b) Ni, (c) SiO 2 , (d) Ti and (e) Na that were studied with AIMD.Crystallographic information for each structure can be found in the supplemental material.

Figure 2 .
Figure 2. Temperature dependence of the C 1j (j = 1, 2, 3) elements of the elastic tensor of Ni, with a range of 300 K to 2100 K.The red dot-dashed line denotes the experimental melting point of Ni (1728 K).The solid orange lines indicate the zoomed in temperature range which is shown in detail in the inset image.

Figure 3 .
Figure 3. Temperature dependence of the C 1j (j = 1, 2, 3) elements of the elastic tensor of (a) SiO 2 , (b) Ti, (c) Na and (d) Au.The red dot-dashed lines denote the experimental melting points of each material.The solid orange lines indicate the zoomed in temperature ranges shown in detail in the inset images.

Figure 4 .
Figure 4.Temperature dependence of the atomic coordination environment of Ti atoms using the common neighbor analysis tool[83,84] in the Open Visualization Tool (OVITO)[85].Atoms are found to be in the hexagonal close packed (HCP), body-centered cubic (BCC) and face-centered cubic (FCC) structures.Any atoms not found to be in one of these crystal environments are grouped in the 'Other' category.

Figure 5 .
Figure 5.The computed value Z total as a function of temperature, along with each of the components (Z 11 , Z 12 and Z 13 ) for (a) Au, (b) Na, (c) Ni, (d) SiO 2 and (e) Ti.The red dot-dashed line denotes the experimental melting point of each material.Shown in (f) is Z total for all of the materials as a function of temperature, with vertical lines to indicate the melting point of each material.

Table 1 .
Range of melting points for Au, Na, Ni, SiO 2 and Ti determined based on calculated Z values.The experimental melting points are include for reference.