Materials design of systems with ionic bonding contributions, i.e., compounds including elements of significantly different electronegativity, necessitates an accurate modeling of their thermodynamic stability [1][2][3][4][5]. The appropriate descriptor is the formation enthalpy [6]-the enthalpy difference between the compound and its elemental references, or its recursive factorization to study multicomponent systems [7]. For metals, high-throughput compatible (semi)local density functional theory (DFT) is known to provide accurate results with errors significantly smaller than the thermal energy at room temperature (≈25 meV/atom) [8,9]. This fueled the construction of large materials databases with millions of entries [10][11][12][13][14][15][16][17][18]. On the contrary, ionic materials pose a much more fundamental challenge for computational approaches.
As outlined in Ref. [1], the formation enthalpy can be subdivided into a total energy difference between the compound and the elements plus a (small) vibrational contribution due to zero-point and thermal effects. As long as all phases involved are chemically similar (in terms of their electronic delocalization character), standard (semi)local DFT's systematic error cancellation allows for a good approximation of the total energy difference [8,9]. This breaks down for ionic systems, such as oxides and nitrides [1][2][3][4]19]: Little error cancellation can be expected between an ionic compound and its metallic/diatomic-gaseous references. Consequently, computing reliable formation enthalpies ab initio would * stefano@duke.edu require accurate total energies for all systems involved. This is generally not possible within a (semi)local approximation.
Significant efforts have been undertaken to investigate the accuracy that can be obtained from a specific level of theory. Many studies demonstrate that compared to experimental formation enthalpies [20][21][22][23], the typical mean absolute error (MAE) for standard functionals-such as LDA [24,25] or PBE [26]-is on the order of several hundred meV/atom [1][2][3][4]19,27,28]. For meta-generalized-gradient approximations, such as the Bayesian error estimation (mBEEF) [29] or the strongly constrained and appropriately normed (SCAN) [30] functionals, an MAE of about 100 meV/atom is obtained [1,27,28,31]. While computationally more demanding, hybrid functionals yield only modest improvements over PBE for transition metal oxides and sulfides [32,33]. Non-self-consistent exact exchange plus random phase approximation (EXX+RPA) and renormalized adiabatic PBE (rAPBE) calculations on PBE orbitals for small sets of about 20 oxides achieved MAEs down to 74-95 meV/atom [32,[34][35][36]. In conclusion, even for the most expensive DFT-based approaches, no satisfactory accuracy (≈ 25 meV/atom) is achieved. Preliminary tests for MgH 2 indicate that quantum Monte Carlo can achieve accurate results with an error of ≈20 meV/atom [37,38], although it remains to be determined whether this applies generally for all materials.
Physically motivated empirical correction schemes parametrizing (semi)local DFT errors with respect to measured values are the only feasible option, achieving accurate formation enthalpies and enabling high-throughput materials design of ionic systems. Initially, a correction for the oxygen reference energy of 1.36 eV per O 2 for PBE was introduced [19]. This scheme was extended to other gases such as H 2 , N 2 , F 2 , and Cl 2 , as well as sulfides for several functionals [39,40]. On top of this, for systems with transition metal ions, an approach for mixing GGA and GGA+U calculations was developed, reducing the MAE to 45 meV/atom for a test set of 49 ternary oxides [3]. Leveraging this method, extensive further parametrization within a local-environment dependent approach was found to lower the MAE to 19 meV/atom [41]. A drawback is, however, that nontransition metals remain uncorrected, which can be particularly problematic for heavy p-block elements [1]. As a complementary approach, the fitted elemental-phase reference energies (FERE) method introduces energy shifts for the elements on an equal footing to minimize the error between measured and calculated results for a large set of binary compounds [2,4]. FERE values for many elements were calculated, yielding an MAE of 48 meV/atom when applied to a test set of 55 ternary compounds. Recently, correction schemes have also been extended to finite temperatures and the Gibbs free energies of solids [42]. It should be noted that the accuracy of schemes fitted to measured values is limited by the experimental error. In the supplementary information of our previous work [1], we investigated the deviation between measured values of different collections for a large set of oxides, indicating that the typical experimental error bar is on the order of 10-20 meV/atom.
While the above correction methods were a major step forward for materials design, their accuracy is limited and the relative stability of polymorphs-sometimes erroneously predicted by DFT [27]-cannot be corrected. Moreover, correction methods based on only composition can lead to incorrect thermodynamic behavior when considering activity vs concentration [1]. To address these shortcomings, we have recently introduced a universal method: coordination corrected enthalpies (CCE) [1]. This advanced correction scheme is the first to leverage structural information by assigning corrections per cation-anion bond, as well as considering the cation oxidation state. CCE achieves an MAE of 27 (24) meV/atom for a test set of 71 (7) ternary oxides (halides), on par with room temperature thermal energy (≈ 25 meV/atom) [1]. It can also correct the relative stability at fixed composition and avoids incorrect thermodynamic behavior by construction.
Here, we present our automated implementation of CCE into the AFLOW framework for computational materials design. It identifies the number of cation-anion bonds, automatically determines oxidation numbers, and includes thermal effects by applying different corrections for designated temperatures. AFLOW-CCE includes a command line interface, a web application, and a Python environment providing useful tools for the scientific community to automatically calculate the CCE correction and formation enthalpies for a given input structure. The paper is organized as follows: After introducing the computational details of the method, a short overview on the CCE functionality in AFLOW is given. Then, the specific analyses within the implementation are described including structural analysis, automatic determination of oxidation numbers, and the inclusion of temperature effects. Available options for the command line interface, CCE corrections for 0 K, and CCE@exp corrections for room temperature are discussed in detail.
The ab-initio calculations for the exchange-correlation functionals LDA [24,25], PBE [26], and SCAN [30] are performed with AFLOW [9,[43][44][45][46][47] and the Vienna ab initio simulation package (VASP) [48] with settings according to Refs. [1,49]. Thermal contributions to the formation enthalpy are calculated using the quasiharmonic Debye model implemented via the AFLOW-Automatic Gibbs Library (AGL) [50][51][52].
Using binary compounds A x 1 Y x 2 as the fit set, the CCE corrections δH T,A +α A-Y per cation-anion A-Y bond and cation oxidation state +α are obtained from the difference between (zero-temperature and zero-pressure) DFT formation enthalpies and experimental standard formation enthalpies at temperature T [1]:
where N A-Y is the number of nearest neighbor A-Y bonds and x i are stoichiometries for the i species. Standard conditions are indicated by the "•" superscript. T can be 298.15 or 0 K, i.e., temperature effects are included in the corrections. A detailed justification of this is presented later.
The corrections can then be applied to any multinary compound
where N i-Y is the number of nearest neighbor bonds between the cation i and anion Y species. For multianion compounds [53], the corrections are summed for all anions separately in Eq. (2).
For the AFLOW-ICSD database, the CCE methodology is applied equivalently with the compound energies partly calculated within DFT+U [49]. In addition, composition dependent energy shifts are applied for the elements for which a U is used to align the related reference energies with the ones calculated from DFT+U .
The room temperature (T r = 298.15 K) formation enthalpy CCE@exp [1] calculated from experimental formation enthalpies per bond δH T r ,i +α i-Y,exp is given by:
These values provide a rough guess with an estimated average accuracy of about 250 meV/atom as obtained from a test for ternary oxides [1].
The automated CCE implementation inside AFLOW enables the correction of an extensive library of ionic materials that are made available via the AFLOW APIs [54,55] and web interfaces [10]. The implementation features three ways of user interaction depicted in Fig. 1: (i) a command line tool, (ii) a web application, and (iii) a Python environment. The The input structure file (here test.POSCAR), precalculated DFT formation enthalpies per cell, and functionals are given via the options --get_cce_corrections < test.POSCAR, --enthalpies_formation_dft=-63.452,-72.084,-72.412, and --functionals=PBE,LDA,SCAN, respectively. The structure can be in any format recognizable by AFLOW, such as VASP POSCAR [48], Quantum Espresso [56], FHI-AIMS [57], ABINIT [58], ELK [59], and CIF [60]. Oxidation numbers for all atoms can be provided as a comma separated list as input using --oxidation_numbers=ox_num_1,ox_num_2,.... (b) When executed, the output includes the CCE corrections and formation enthalpies at both 298.15 and 0 K. If no DFT formation enthalpies are given, an estimate for the formation enthalpy at 298.15 K based on experimental values per bond (CCE@exp, blue) [1] is calculated according to Eq. ( 3). (c),(d) Example command/output when determining oxidation numbers for the structure in test.POSCAR. (e),(f) Example command/output when determining cation coordination numbers for the structure in test.POSCAR. (g) The web interface yields the cation coordination numbers, oxidation numbers, and CCE corrections for the structure provided in the field "Input POSCAR." If DFT formation enthalpies per cell are provided, the output also includes the CCE formation enthalpies. (h) Example Python script using the AFLOW-CCE Python environment. Similar to the command line, functionals, enthalpies_formation_dft, and input oxidation_numbers are optional arguments for the get_corrections method. The results are returned as a dictionary. command line tool [Figs. 1(a)-1(f)] provides the CCE corrections and formation enthalpies, (automatically determined) oxidation numbers, and cation coordination numbers for the given structure file that can be in any format recognizable by AFLOW, such as VASP POSCAR [48], Quantum Espresso [56], FHI-AIMS [57], ABINIT [58], ELK [59], and CIF [60]. Available options are described in section "CCE command line interface." The web interface [Fig. 1(g)] prints the cation coordination numbers, oxidation numbers, and CCE corrections for the selected functionals using the given structure. The output also includes the CCE formation enthalpies when precalculated DFT values are entered in the designated fields. The Python environment is distributed with the AFLOW source and can be generated with the command aflow --cce --print=python. It connects to the command line functionality and imports the results into a CCE class similar to the Python modules of AFLOW-SYM [61] and AFLOW-CHULL [62]. An example script leveraging the functionality is depicted in Fig. 1(h). The CCE object has three built-in methods:
get_corrections(struct_file_path, functionals, enthalpies_formation_dft, oxidation_numbers), get_oxidation_numbers(struct_file_path), and get_cation_coordination_numbers (struct_file_path) corresponding to the command line options mentioned in section "CCE command line interface." Each method requires a path to the input structure file (struct_file_path). For get_corrections, providing functionals, DFT formation enthalpies (enthalpies_formation_dft), and input oxidation_numbers for all atoms in the structure are optional arguments. The results are returned in the form of a Python dictionary.
Prints instructions and example input structure.
Prints the results of the full CCE analysis, i.e., cation coordination numbers, oxidation numbers, and CCE corrections and formation enthalpies, for the given structure. STRUC-TURE_FILE_PATH is the path to the structure file. The file can be in any format supported by AFLOW, e.g., VASP POSCAR, QUANTUM ESPRESSO, AIMS, ABINIT, ELK, and CIF. For VASP, a VASP5 POSCAR is required or, if a VASP4 POSCAR is used, the species must be written on the right side next to the coordinates for each atom just as for the example input structure obtained from --cce.
Determines the CCE corrections and formation enthalpies for the structure in file STRUCTURE_FILE_PATH.
Determines the oxidation numbers for the structure in file STRUCTURE_FILE_PATH.
Determines the number of anion neighbors for each cation for the structure in file STRUCTURE_FILE_PATH.
for --cce=STRUCTURE_FILE_PATH and --get_cce_corrections < STRUCTURE_FILE_PATH:
--enthalpies_formation_dft=enth_1,enth_2,... enth_1,enth_2,... is a comma separated list for precalculated DFT formation enthalpies. They are assumed to be: (i) negative for compounds lower in enthalpy than the elements, (ii) in eV/cell. Currently, corrections are available for PBE, LDA, and SCAN.
--functionals=func_1,func_2,func_3 func_1,func_2,func_3 is a comma separated list of functionals for which corrections should be returned. If used together with --enthalpies_formation_dft, the functionals must be in the same sequence as the DFT formation enthalpies they correspond to. Available functionals are: (i) PBE, (ii) LDA, and (iii) SCAN. Default: PBE (if only one DFT formation enthalpy is provided).
--oxidation_numbers=ox_num_1,ox_num_2,... ox_num_1,ox_num_2,... is a comma separated list of oxidation numbers. It is assumed that: (i) one is provided for each atom of the structure and (ii) they are in the same sequence as the corresponding atoms in the provided structure file.
Obtain output in table format (--print=out) or as JSON (--print=json). Default: out.
For evaluating the number of cation-anion bonds (cation coordination numbers), first the (main) anion species of the system is determined as the one with the highest Allen electronegativity (EN) [63][64][65]. A check is performed whether the material is a multianion system [53], i.e., whether atoms of a type other than the main anion species are only bound to atoms of lower EN or its own type. If such atoms are found, they are designated as additional anions. This is for instance the case for N in HfTaNO 3 , where O is the main anion. Note that in some compounds certain species can occur both as anion and as cation: In ammonium-nitrate (NH 4 NO 3 ) for instance, N occurs both in -3 and +5 oxidation states depending on its neighbors.
Subsequently, the number of anion neighbors for each cation is determined. For this bonding analysis, the nearest neighbor distance is obtained for each species. A bonding cutoff is set by adding a tolerance of 0.5 Å in accordance with Ref. [1]. Then, all anion neighbors between the species selective minimum distance and the cutoff are counted. For the multianion analysis, the tolerance is reduced to 0.4 Å since for larger values, different anions of systems known to be multianion compounds would be detected as being bonded. For instance, in HfTaNO 3 , if the tolerance is not reduced, N and O would be detected as being neighbors and hence nitrogen would not be identified as an anion.
When oxygen is found as an anion, the O-O distances in the system are determined to detect per-and superoxides. The following scenarios can occur: (i) the O-O bond is longer than 1.6 Å indicating an oxide (O 2-ion), (ii) the bond length is between 1.4 and 1.6 Å (peroxide), (iii) the bond length lies between 1.3 and 1.4 Å(superoxide), and (iv) the bond length is shorter than 1.3 Å, i.e., the structure may contain molecular oxygen, the enthalpy of which is not correctable within CCE. For certain special cases such as alkali metal sesquioxides, several of the above scenarios can be fulfilled simultaneously and the implementation will then treat the system as incorporating multiple different oxygen ions. [1]. The number of (su-)peroxide bonds is determined as half the number of (su-)peroxide O atoms.
The default method to automatically determine oxidation numbers is based on Allen ENs [63][64][65]. This choice is in accordance with International Union of Pure and Applied Chemistry (IUPAC) recommendations [68,69] and also conforms with our own tests that this EN scale yields the most reliable oxidation numbers when compared to other scales such as Refs. [67,70]. Table I lists the EN values together with the preferred and all known oxidation numbers for the elements according to Ref. [66], along with a few additions deemed necessary during the test of the implementation. The separation into preferred and all known oxidation numbers is motivated by the finding that in compounds with more than two species, elements tend to occur only in the preferred oxidation states [1]. Missing oxidation states will be added in future releases as needed.
The algorithm (Fig. 2) starts by assigning the anion oxidation numbers. For all anion atoms (main anion species and anions from multianion analysis) the lowest (most negative) oxidation number known for this species is assigned. If the atom was found to belong to a peroxide (superoxide) ion in the structural analysis, the oxidation number is changed to -1 (-0.5).
The set of possible cation oxidation states is first restricted to the preferred values for each species. All cations are then assigned the first (usually most positive) preferred oxidation number for their species. The only exception is Cr for which +3 is the first choice (Table I). After this initial assignment, the sum over all oxidation numbers is evaluated and-if it is zero-the assignment is considered successful. Otherwise, the algorithm proceeds by changing the preferred oxidation numbers according to EN: While checking the sum for each choice, the oxidation states of more electronegative (higher EN) cation species are changed to the next preferred FIG. 2. Oxidation number algorithm. Schematic representation of the algorithm to determine oxidation numbers α of a compound with n species. The part of the algorithm for determining cation oxidation numbers (inside the black dashed box) is first applied making use of the preferred oxidation numbers for all species. If no successful assignment is achieved, it is employed a second time after checking for mixed-valence compounds using all known oxidation states (Table I). The oxidation numbers of more electronegative cation species are iterated faster than the more electropositive ones. Three dots indicate proceeding equivalently for all further cation species. For multianion systems, atoms identified as additional anions during the structural analysis are excluded when assigning cation oxidation numbers. value before the more electropositive (lower EN) ones. It is expected that more electropositive elements occur in higher oxidation states.
If all EN-directed choices of preferred oxidation numbers are exhausted without successful assignment, the system is checked for mixed valence. For these special cases, the oxidation numbers are set explicitly. Currently, this scenario includes Sb 2 O 4 , Pb 3 O 4 , Fe 3 O 4 , Mn 3 O 4 , Co 3 O 4 , Ti-O Magnéli phases, and alkali-metal sesquioxides. If still no successful assignment is achieved, the part of the algorithm for determining cation oxidation numbers (inside the black dashed box in Fig. 2) is repeated with all known oxidation numbers for all cation species. The scheme has been successfully tested on a large number of compounds, including oxides, fluorides, chlorides, and nitrides. The algorithm might not be particularly suited for organic compounds for which the oxidation a Since there are no available Allen electronegativities for La-Yb, the value for Lu is used as these elements are usually very similar. This is confirmed by the Allred and Rochow electronegativities that are very similar for all lanthanides [67].
state of C depends on the functional group. Such materials are presently beyond the scope of AFLOW-CCE. For oxides, the oxidation numbers can also be determined from Bader charges [71], which are compared to the averaged template values of the binary fit set for the respective functional. The formal oxidation number is assigned according to the closest template value. However, this scheme shows systematic difficulties in assigning the correct oxidation numbers for several species in certain oxidation states such as Ti 4+ , V 5+ , Fe 2+ , and Fe 3+ , for which error handling procedures have been implemented. This method is only invoked when specifically requested via the setting DEFAULT_CCE_OX_METHOD=2 in the .aflow.rc setup file of AFLOW. Finally, the user can also provide the oxidation of the CCE method for each cation species A in oxidation states +α for 0 K obtained from binary oxides. The corrections for Si and Ti in oxidation state +4 are obtained from α-quartz (AFLOW label A2B_hP9_152_c_a [74,75]) and rutile (AFLOW label A2B_tP6_136_f_a [74,75]), respectively. The corrections in the last line are for (su-)peroxides according to the approach outlined in Ref. [1]. All corrections are in eV/bond. numbers for all atoms as a comma separated list as input (option --oxidation_numbers= in section "CCE command line interface").
After the determination of the oxidation numbers, the cation-anion and cation oxidation state specific CCE corrections per bond δH T,A +α A-Y [Eq. (1)] can be assigned. As outlined in Ref. [1], vibrational (zero-point and thermal) contributions to the formation enthalpy do not need to be calculated explicitly since they can be parameterized per bond and thus implicitly included into the corrections. Compared to when the vibrational contribution was explicitly included for the fit and test sets, the MAE of the corrected results increased by at most 1 meV/atom. This is negligible considering that the CCE MAE is on the order of 30 meV/atom. Temperature effects are thus included in the corrections according to Fig. 3(a):
The CCE corrections to DFT formation enthalpies are fitted to experimental room temperature formation enthalpies resulting in room temperature corrections. When these are applied according to Eq. ( 2), a direct estimate of the room temperature formation enthalpy is obtained.
For 0 K [Fig. 3(b)], one first subtracts the calculated thermal contribution, deduced from a quasiharmonic Debye model [51] according to Ref. [1], from the experimental formation enthalpy for each functional. This gives a good estimate for the (experimental) 0 K formation enthalpy. Then, the CCE corrections to the DFT formation enthalpies are fitted to these values yielding 0 K corrections (see Table II). The approach does not capture any phase change of the elemental references from 0 K to room temperature. However, this is a rare event that occurs on an energy scale below room temperature, which is smaller than the CCE error.
To test the accuracy of the predicted 0 K formation enthalpies, they are compared to available tabulated values. of the CCE@exp method for each cation species A in oxidation states +α for 298.15 K obtained from binary oxides. The corrections for Si and Ti in oxidation state +4 are obtained from α-quartz (AFLOW label A2B_hP9_152_c_a [74,75]) and rutile (AFLOW label A2B_tP6_136_f_a [74,75]), respectively. The corrections in the last line are for (su-)peroxides according to the approach outlined in Ref. [ meV/atom for PBE, LDA, and SCAN, respectively. When corrected by CCE, the DFT results are drastically improved [Fig. 4(b)] with mean errors reduced to 23, 14, and 13 meV/atom, validating our 0 K approach.
As a future development, the temperature dependence can be implemented as a continuous variable, since it can be parameterized per bond and the thermal contributions at any temperature can be computed from the quasiharmonic Debye model [51] for the fit set. This ansatz, and other approaches directly targeting the Gibbs free energy [42] to include finite temperature effects, pave the way to move beyond stability predictions based only on enthalpies, which are crucial for, e.g., high-entropy materials [72,73].
The corrections are finally used to calculate the total CCE corrections and the CCE formation enthalpies according to Eq. (2) for 298.15 and 0 K for all functionals selected. If needed, corrections for (su-)peroxides are added according to Ref. [1] with the number of respective O-O bonds obtained from the structural analysis. If no precalculated DFT formation enthalpies are provided, an estimate for the formation enthalpy at 298.15 K based on experimental values per bond (CCE@exp) (see Table III) [1] is calculated according to Eq. (3).
We have presented our implementation of the coordination corrected enthalpies (CCE) method into AFLOW for automated correction of DFT formation enthalpies. AFLOW-CCE provides a universal tool to obtain highly accurate formation enthalpies for ionic materials with a typical mean absolute error close to the room temperature thermal energy, i.e., ≈ 25 meV/atom [1]. It interoperates with the existing func-tionality of AFLOW and features a command line tool, a web interface, and a Python environment. Additionally, the AFLOW-CHULL module will be updated with the CCE formation enthalpies where appropriate [62].
The AFLOW-CCE workflow includes a structural analysis to identify the number of cation-anion bonds, an automatic determination of oxidation numbers based on Allen electronegativities, and the inclusion of temperature effects by parametrizing vibrational contributions to the formation enthalpy per bond. With all the required functionality in place, the implementation will be extended to other anion classes beyond oxides such as nitrides, halides, and sulfides by adding the needed corrections in the near future.