Theoretical Prediction of the Sublimation Behavior by Combining Ab Initio Calculations with Statistical Mechanics

We develop a theoretical model to predict the sublimation vapor pressure of pure substances. Moreover, we present a simple monoatomic molecule approximation, which reduces the complexity of the vapor pressure expression for polyatomic gaseous molecules at a convincing level of accuracy, with deviations of the Arrhenius prefactor for NaCl and NaF being 5.02% and 7.08%, respectively. The physical model is based on ab initio calculations, statistical mechanics, and thermodynamics. We illustrate the approach for Ni, Cr, Cu (metallic bond), NaCl, NaF, ZrO2 (ionic bond) and SiO2 (covalent bond). The results are compared against thermodynamic databases, which show high accuracy of our theoretical predictions, and the deviations of the predicted sublimation enthalpy are typically below 10%, for Cu even only 0.1%. Furthermore, the partial pressures caused by gas phase reactions are also explored, showing good agreement with experimental results.


Introduction
The determination of the thermodynamic parameters of the solid-to-vapor (sublimation) phase transition is significant for chemical thermodynamics. Even if vapor pressures are often low, they can lead to mass loss and accompanying concentration change and resulting in modification of functional behavior. Additionally, the vapor can condense elsewhere and also affect the material behavior there. Consequently, sublimation processes are considered to be important for various applications in different scientific areas [1]. For example, the analysis of the sublimation behavior of cathode materials such as La 0.58 Sr 0.4 Co 0.2 Fe 0.8 O 3−δ (LSCF) in solid oxide fuel cells (SOFCs) and electrolytes can contribute to the reduction of their degradation, which is of high interest for green energy applications [2]. Furthermore, Liu et al. have determined the thermodynamic activities of SrO, FeO and CoO in SrTi 1−x−y Co x Fe y O 3−δ compounds by measurements of the fugacity, whereby SOFC degradation-related chemical reactions can be predicted. In general, the investigation of the sublimation process and expanding the thermodynamic database have huge potential to be applied in various fields, such as material stability, chemical thermodynamics, energy, and chemical engineering [3].
Among the experimental methods to determine vapor pressures is the Knudsen effusion mass spectrometer (KEMS), which offers the highest accuracy for vaporization studies under near equilibrium conditions [4]. Nevertheless, this method is mainly suitable for high-temperature applications to ensure sufficiently high partial pressures, requiring in particular extended calibration processes, and is not suitable for all materials, hence the information in existing thermodynamic databases is still insufficient [5]. Therefore, to complement experimental methods, a theoretical and modeling perspective is desirable, which also helps to understand the details of the sublimation process.
For many years, electronic structure density functional theory (DFT) has played a great role in computational materials science, by solving the Kohn-Sham equations to obtain self-consistent electronic densities and to find ground-state energy. Initially, DFT calculations were designed for the prediction of ground-state energy (at 0 K), allowing the prediction of electronic, optical, mechanic, and thermodynamic properties of materials such as the density of states (DOS), band structure, second-harmonic generation coefficient, Young's modulus, adsorption energy, chemical reaction activation energy, solvation energy, work function, etc.; see, e.g., [6][7][8]. Based on the DFT calculations, certain thermodynamic properties, phase transitions, and the influences of structural defects can be investigated [9]. Cervinka and Fulem investigated the sublimation enthalpy using DFT calculation [10]. Lopes Jesus et al. studied the sublimation enthalpy of the solid state of organic compounds also with the help of DFT calculations [11]. Furthermore, Halpern and Marzzacco used classical and statistical thermodynamics to calculate thermodynamic properties, such as the enthalpy and entropy of fusion and vaporization for water [12]. Zaby et al. used a modified partition function in quantum cluster equilibrium theory to predict the enthalpy and entropy of vaporization by employing statistical mechanics [13]. Thus, the analysis of the sublimation behavior by the combination of DFT calculation and statistical mechanics is promising and used for thermodynamic databases [14]. However, a pragmatic theoretical model which can predict the entire sublimation function for a solid is still lacking.
In this context, a reliable method combining ab initio calculations, statistical mechanics, and thermodynamics is desired and pursued in this article, which is solely dependent on theoretical calculations without any experimental or adjustable parameters. As such expressions can easily become involved, especially for complex molecules, simplifying expressions are desired, which allow for an easy, but still accurate prediction of vapor pressures. Therefore a pragmatic "monoatomic molecule approximation" is developed to simplify the polyatomic gas molecule calculations. This approximation can improve the calculation efficiency while maintaining the accuracy of the predictions and allows the reduction of the parameter requirements for the model. Although such a reduction may not be applicable for all cases, we demonstrate for a variety of molecules that the vapor pressure predictions are indeed reliable.
In our calculations, the gaseous phases are treated as ideal gases, and pure substances with different bonds are selected for demonstration purposes, i.e., Ni, Cr, Cu with a metallic bond, NaCl, NaF, ZrO 2 with an ionic bond and SiO 2 with a covalent bond. First, in our physical model, the cohesive energy E coh of the solid phase is required as the difference between the solid and gaseous molecules, and it is obtained from DFT calculation using the Vienna Ab Initio Simulation Package (VASP) [15,16]. Second, we consider the different degrees of freedom of the molecules in the solid lattice and in the gas phase and construct their equilibrium conditions in the framework of statistical mechanics.
Atoms in crystalline solid phases vibrate with different frequencies in the lattice, and the vibration behavior is described by the phonon spectra [17]. These spectra are obtained via VASP with help of the Phonopy package, which is an open-source package for phonon calculations at harmonic and quasi-harmonic levels [18]. The K-path for the calculation of phonon spectra is obtained from the open-source package Vaspkit [19]. Instead of using the full spectra, we reduce the acoustic and optical branches to mean vibrational frequencies as central simplification.
Therefore, the article is organized as follows. In the next section, we develop expressions for the equilibrium vapor pressure and consider the possibility of further chemical reactions in the gas phase. Therefore, we first demonstrate the basic calculations given the complete molecule's motions and then present the aforementioned monoatomic molecule approximation by considering only the translational degrees of freedom of the gaseous molecules. Consequently, the diatomic or triatomic gaseous molecules are treated similarly to the monoatomic gaseous molecules. Although this suppression is not appropriate, e.g., for the calculation of heat capacities, it turns out that it is sufficient for vapor pressure predictions. Afterward, we use statistical mechanics to connect the microstates to the thermodynamic properties. During sublimation, the exchange of particles between the solid phase and gaseous phases is possible, hence we use the equality of chemical potentials between the gas phase and the solid as an equilibrium condition. Moreover, further dissociation or following gas phase reactions are investigated. Here, we use Ni 2 (g), Cu 2 (g), Na(g), Cl(g), Zr(g) and O 2 (g) as examples to elucidate such situations. After these theoretical foundations, we present the computational details in the following section. The predictions are validated by comparing them to the commercial databases FactPS from the software package FactSage 8.1 [20]. Finally, we summarize and discuss the results of the physical model and elucidate the advantages and disadvantages in detail.

Theory
The sublimation behavior indicates the phase transition of the substance from the solid state to the gaseous state. During sublimation, the atoms vibrate in the lattice, and when the energy is large enough, the atoms are released from the lattice to become gaseous molecules. Typically, the vapor pressure p has an Arrhenius form, in the following denoted as sublimation function, with A being a constant (prefactor). This relationship assumes that sublimation is a reaction with a single relevant activation barrier. The function shows a linear relationship between ln p and 1/T. Often, the sublimation enthalpy is given per mole of the substance, and then the slope of this line is defined as ∆H sub /R with the gas constant R. As we use here an atomic scale description, we use the equivalent normalization of energies per atom or molecule. Moreover, some molecules become unstable at a certain temperature, which leads to dissociation reactions. In our physical model, further possible chemical reactions are also considered. Thus, the main task for predicting the vapor pressure is the determination of the sublimation enthalpy H sub and the prefactor A.

Sublimation Function
The calculation of the sublimation function is performed using statistical mechanics and DFT calculations. The canonical partition function Z for a system is given by [14,21,22] where E i is the energy of the ith quantum state in the system and k B the Boltzmann constant, with the summation running over the entire phase space. The Helmholtz energy can be written as The chemical potential is obtained as where N is the number of particles in this system. In case the Hamiltonian decomposes additively into independent contributions, the canonical partition function becomes For a crystalline phase, z i can in particular be a vibrational oscillator. In general, the vibrational frequency of each oscillator is different and can be obtained from the phonon spectrum. The crystal structure has a significant influence on the phonon spectrum, consisting of acoustic and optical branches [17]. In our model, we use a simple harmonic approximation and anharmonic effects are ignored. Therefore, the expression for the partition function of single vibrational mode is [21] Here, β is equal to 1/(k B T), = h/2π is the reduced Planck constant, ω is the vibrational frequency of the considered mode, for which we will use in the following averaged values. More precisely, we take the average of the highest and lowest frequencies of both the acoustic and optic branches of the phonon spectrum.
Similarly, the partition function of gaseous phases is given by: where the factor 1/N! accounts for the usual indistinguishability of the atoms. In general, the gaseous molecule has different degrees of freedom, such as translational (z trans ), rotational (z rot ), and vibrational (z vib ) degrees of freedom, and their expressions are [23] Here, m is the mass of the gaseous molecules, V is the volume, σ is the symmetry number of gaseous molecules, c is the light speed, B is the rotational constant, I is the moment of inertia, µ is the reduced mass and r is the interatomic gas distance, and v is the vibrational frequency of the gaseous molecules. The energy for the separation of electrons from the ground state is usually very large, thus, for most situations z el = 1 is a suitable approximation. A significant exception is molecules with electronically degenerate ground states, for instance, the alkali metal atoms with z el = 2 [23].
The construction of the equilibrium state is then based on the chemical potential balance, In this study, we exploit a physical model for the monoatomic molecule and polyatomic molecule first and then propose a simple approximation (monoatomic molecule approximation) that simplifies the expressions for practical applications.

Monoatomic Molecule
A monoatomic gaseous molecule has only three translational degrees of freedom. For the solid phase, we use only the mean vibrational frequency of acoustic branches. Thus, the chemical potential of the molecule in the solid phase and gaseous molecules can be expressed as with From the chemical potential balance we obtain the expression for the vapor pressure in ideal gas approximation, which involves the cohesive energy E coh = E gas − E solid per particle, which is obtained from ab initio simulations [24].

Diatomic Molecule
For a diatomic molecule, we need to consider (i) the potential appearance of optical branches in the phonon spectrum of the solid phase, which are related to additional vibrational degrees of freedom, and (ii) the role of rotational and vibrational contributions in the gas phase. As a result, the vapor pressure can be written as In this paper, we use the mean frequency of acoustic and optic branches separately, thus, the resulting expression becomes where ω acoustic and ω optic are the mean frequencies of the acoustic and optical branches by taking the maximum and minimum value, respectively.

Monoatomic Molecule Approximation
As the complexity of the description rises for larger molecules, a central question of the current paper is whether a simplified "monoatomic molecule approximation", where the internal vibrational and rotational degrees of freedom of a molecule are ignored and it is treated as a point-like particle, is sufficient to deliver reliable results for the prediction of the vapor pressure. Essentially, this monoatomic molecule approximation ignores partially heat capacity contributions in both phases, i.e., the internal relative motion in the unit cell and the vibrational, rotational motion from the gaseous molecules. Therefore, a gaseous molecule is simplified to have only the transitional degrees of freedom. Additionally, for the crystalline phase, which contains different atom types inside its unit cell, only the acoustic branch of the phonon spectrum is retained.
This common simplification on both sides of the chemical potential balance equation offsets part of the deviations. Essentially, this simplification has an impact on the contribution of the heat capacity; however, the main part of the sublimation enthalpy is cohesive energy [10], and therefore we expect this approximation to deliver reasonable results.

Chemical Reaction
We use an A m B n crystal as an example to clarify the complex sublimation behavior and its further chemical reaction. Aspects such as defect formation in the crystal are not considered, and we refer to [25] for further details. The gaseous molecule A m B n (g) evaporates from the A m B n (s) solid crystal, and then it dissociates into A(g) gas and B(g) gas atoms. According to Hess's Law, the complete process is written as In equilibrium, the change in Gibbs free energy ∆G must be equal to 0. According to Van't Hoff's equation (∆G = −RT ln K), K is the equilibrium constant of the reaction. We therefore obtain Hence we obtain According to the balance of the chemical reaction we have mp A(g) = np B(g) , hence Equation (19) can be reformulated as with ∆G = ∆H + ∆(TS ).
Practically, we calculate the sublimation enthalpy at 0 K, and therefore ∆G = ∆H , hence Then Equation (20) can be rewritten as Thus, the enthalpy for the above-mentioned chemical reaction is At this point, a consistent normalization of the enthalpies by converting it to 1 mol leads to and For the determination of the prefactor (standard vapor pressure), also the reaction between the gaseous phases has to be considered, A n B m (g) → nA(g) + mB(g) (27) Therefore, the equilibrium pressures of A(g) and B(g) can be calculated based on the above chemical reaction equation. Consequently, the calculation of a complex system with many different gaseous species becomes involved, and we demonstrate it here by looking at examples such as NaCl(s) to the gas phases Na(g) and Cl(g).

Computational Details
This study uses the most stable structure at 0 K, and the basic structures are from the existing Materials Project database [26]. We apply the projector-augmented wave (PAW) method with a 520 eV cutoff energy to represent electron-ion interactions and the generalized gradient approximation (GGA) with the functional proposed by Perdew, Burke, and Ernzerhof (PBE) [27,28] within the plane wave code VASP. Structural and bonding details are listed in Table 1. For the DFT calculations of the gas phase, we put the single gaseous molecule into a 30 × 30 × 30 Å 3 vacuum box to reduce the interaction between the gaseous molecules to mimic an ideal gas. For the solid structure, we use 2 × 2 × 2 or 3 × 3 × 3 supercell sizes, the k-point settings depend on the substances and system size, and the Monkhorst-Pack and the Γ centered scheme are used for crystal and gaseous molecules, respectively. Table 2 shows the central results needed for the following steps. In the thermodynamic database FactPS, Ni 2 and Cu 2 diatomic gaseous phases exist at the equilibrium state; therefore, their related partial pressures are also explored. Na/Cl(NaCl) and Zr/O 2 (ZrO 2 ) express the decomposition reactions of NaCl and ZrO 2 . Na/Cl(NaCl) and Zr/O 2 (ZrO 2 ) mean Na(g)/Cl(g) from NaCl(g) and Zr(g) and O 2 (g) from ZrO 2 (g), respectively.  Figure 1a shows the results of monoatomic molecules: Ni (red solid line), Cr (green solid line), and Cu (orange solid line). The dashed lines represent the results from the calculations with the FactPS database for comparison. In the Arrhenius plot the slope reflects the sublimation enthalpy and the intercept is the prefactor.

Monoatomic Molecule
The trends (sublimation enthalpy) are nearly the same, and the deviations of sublimation enthalpy are below 10% (Ni: 8.45%, Cr: 1.8%, Cu: 0.1%), and the prefactor accuracy is acceptable (Ni: 6.24%, Cr: 10.45%, Cu: 13.91%). The main reason for the deviation of the prefactor may be the inaccurate calculation of the phonon spectra and the used average frequency.

Monoatomic Molecule Approximation
In this section, we use two different methods to proceed with the theoretical prediction of sublimation functions of NaCl, NaF: (i) Considering the complete description of the molecules and (ii) the monoatomic molecule approximation. Furthermore, we extend the latter to triatomic molecules, ZrO 2 , SiO 2 . The polyatomic molecules' motion is significantly more complex than for monoatomic molecules, as they have translational, vibrational, and rotational degrees of freedom. To verify the accuracy of the monoatomic molecule approximation, we compare the results for NaCl, NaF with the data considering all degrees of freedom to our approximation in Figure 1b. For a diatomic molecule, the vibrational and rotational degrees of freedom of diatomic molecules are included, i.e., the frequencies of optical branches from phonon spectra, the vibrational frequencies of gaseous molecules and their related rotational motion parameters need to be taken into account. The results of both methods match those obtained from the FactPS database very well. Our approximation influences mainly the prefactor, and the predictions (Approximation: 5.02%, complete: 2.71%) of NaCl are in convincing agreement with the database. Similar results for NaF with similar deviations of the prefactor (approximation: 7.08%, complete: 1.71%) are obtained. Comparisons between our predictions and the database show that the deviations mainly result from the sublimation enthalpy (slope). The bond lengths of NaCl and NaF are calculated using T = 0 K DFT calculations, which are found by fitting the energy-volume curves, not considering thermal expansion. Therefore, the accuracy of the bond lengths prediction and the further energy calculations may influence the calculation of the cohesive energy, which in turn affects the further sublimation enthalpy prediction. Nevertheless, our monoatomic molecule approximation leads to reliable predictions.
Furthermore, we have extended the monoatomic molecule approximation to triatomic molecules (SiO 2 , ZrO 2 ) to demonstrate the precision of the predictions. Figure 1c illustrates our approximation application for the triatomic molecules, and the results are consistent with data obtained from FactSage. The deviations (sublimation enthalpy: 7.06% prefactor: 10.43%) of SiO 2 are sufficient for many practical applications. However, in contrast to the database, the sublimation enthalpy accuracy of ZrO 2 is relatively low, and possible reasons are: (i) Due to the instability of these gaseous molecules in reality, the determination of the molecule's configuration and related structural information is carried out using DFT calculations alone. (ii) The inaccurate energy of the gaseous molecules induces sublimation enthalpy deviations. The following calculations of ZrO 2 in the next part further support this interpretation.

Chemical Reaction
We study gas phase reactions for Ni 2 , Cu 2 , NaCl and ZrO 2 as examples, based on the following reactions: According to Equations (24)-(26), the total enthalpy (sublimation enthalpy + chemical reaction enthalpy) of the above three chemical reactions are obtained, The Arrhenius prefactors for these gaseous phases are obtained according to the above chemical reactions. Although other chemical reactions in this process are possible, the unique final equilibrium state between the different gaseous resultants is reached. The source of these gaseous phases is the congruent sublimation process, e.g., the NaCl(s) to NaCl(g) sublimation process, such as Na(g), Cl(g) are all products from NaCl(g). Therefore, two different relevant equilibria play a role, namely between NaCl(s) and NaCl(g) and between NaCl(g) and Na(g), Cl(g). In this paper, subdominant contributions are not considered, such as some elements escaping from the lattice with lower energy and leaving vacancies. Hence, we use p Ni 2 (g) = p Ni(g) /2, and p Cu 2 (g) = p Cu(g) /2, respectively, for the Ni 2 and Cu 2 gas. The dissociation of NaCl and ZrO 2 is different, where we have by mass conservation p Na(g) = p Cl(g) = p NaCl(g) . Figure 1d reveals the precision of our model; especially for ZrO 2 the deviations of the enthalpy (1.96%) and the prefactor (0.7%) are both below 2%. Furthermore, in this calculation, the energy of ZrO 2 (g) is based on Hess's law reduced, because the two distinct steps are treated as a single step, i.e., step 1: ZrO 2 (s) to ZrO 2 (g) and step 2: ZrO 2 (g) to Zr(g) and O 2 (g). We calculate the sublimation enthalpy of Zr(g) and O 2 (g) directly using the energy of Zr(g), O 2 (g) and ZrO 2 (s) without the energy of ZrO 2 (g) so that the negative influences from the inaccurate results of energy of ZrO 2 (g) disappear. Thus, the inaccuracy of the ZrO 2 (g) structural prediction may influence the sublimation enthalpy calculation which leads to a relatively higher deviation. However, the model's accuracy compared to the database is still convincing.

Conclusions
In this paper, the possibilities of coupling ab initio calculations, statistical mechanics, and thermodynamic calculations to theoretically investigate the sublimation process, including gas phase reactions, are explored. The following key points can be concluded: • The comparison of the predictions to thermodynamic databases (FactPS database) has proven that the theoretical prediction of the vapor pressure is promising and precise. • Theoretical predictions of NaCl and NaF with two different methods reveal the feasibility of the monoatomic molecule approximation. The deviations of the prefactor are 5.02% and 7.08%, respectively. The application to triatomic molecules (ZrO 2 and SiO 2 ) indicates the benefit for complex situations, especially SiO 2 , where the deviations of the sublimation enthalpy and the prefactor are 7.06% and 10.43%, respectively. • The additional exploration of the formation of Ni 2 (g), Cu 2 (g), Na(g), Cl(g), Zr(g) and O 2 (g) indicates the rationality of our theoretical calculation to explain further chemical reactions. Except for the predicted prefactor of Na/Cl(g) and Ni 2 , the deviations are below 5%. • The determination of unknown gaseous molecule structures, the approximation of the molecules' motion, and inaccuracies of the thermodynamic database may lead to deviations.