Convergence of Electronic Structure Properties in Ionic Oxides Within a Fragment Approach

Embedded-cluster models of crystalline solids are important to allow accurate wave function methods to be applicable to solids. The ab-initio model potential method, in which the crystal is divided into three different fragments, one quantum fragment, one ab-initio model potential fragment and one point-charge fragment, has historically been shown to be a viable tool for describing the electronic structure in ionic solids. The optimal size of these regions is, of course, individual for each crystal. In this study we analyzed the convergence of the electronic structure properties with respect to an increase of the size of the quantum part and the layer of potentials. MgO crystal and Ni: MgO were used for this purpose as examples of an ideal crystal and a crystal with a point defect. We demonstrated that with an increase of the cluster size, the electron density in the inner part of the cluster becomes very similar to the electron density in the periodic model. Clusters, embedded into a layer of model potential and electrostatic field, are a good alternative to periodic description.

1 Introduction

Metal oxides are important materials with many industrial applications in the areas of catalysis, glass manufacture, extractive metallurgy and many others (Jolivet et al., 2000; Di Cosimo et al., 2014; Arunadevi and Kirubavathy, 2022). Many metal oxides are characterized by highly ionic chemical bonds, which plays an important role in the properties of the materials and chemical reactions.

Alkali and alkaline earth metal oxides are good insulators with a very large optical gap, which makes them an ideal candidates for treatment by means of band theory in combination with density functional theory (DFT) (Schönberger and Aryasetiawan, 1995; Canney et al., 1999; Lany, 2014; Heo et al., 2015; Niedermeier et al., 2016). The interest to small alkaline earth metal oxides clusters is stimulated by a potential use of these materials for hydrogen storage (Srivastava et al., 2018; Mojica-Sánchez et al., 2019) and in biomedical applications (Awasthi et al., 2021). In the presence of dopants, for example, transition or rare metals in metal oxides (these compounds are most interesting from a scientific and industrial point of view), the electronic structure is much more complicated. First of all, in order to model a sparse distribution of the dopants in the host material, the unit cell has to be substantially increased to a so-called supercell. The accuracy of DFT might be insufficient to reproduce the electronic structure of highly correlated d- and f-elements. For these materials, the embedded cluster approach can be considered as an alternative to band structure calculations: 1) there is more flexibility in the selection of the size and the shape of a cluster, 2) the electronic structure calculations are not limited by DFT, but instead the full arsenal of accurate molecular methods from coupled cluster to multiconfigurational theory can be applied. These accurate methods provide a consistent way to include electron correlation effects, but only methods based on a single reference wavefunction were recently generalized for the periodic boundary conditions (Ayala et al., 2001; Hirata et al., 2004; McClain et al., 2017).

At the same time, DFT makes it possible to perform a detailed comparison between periodic models and cluster models of increasing size. Ideally, the inner part of a cluster should have the same electron density as in an infinite crystal. An increase of the cluster size should make the difference between the reference electron density (as in periodic model) and the electron density in the cluster smaller and smaller. In the current study we considered two systems: magnesium oxide (MgO) and magnesium oxide with sparse nickel substitution of magnesium atoms (Ni: MgO). Pure MgO is a wide gap insulator and thus one can expect very small dependence of the results on the amount k-points (or, similarly on the size of the supercell for a calculation with only the Γ point). For the cluster models, assuming non stoichiometric clusters made by the neighboring layers of atoms, there are two possibilities - clusters that are either centered at a magnesium atom or at an oxygen atom. Modeling sparse defects is different from modeling ideal crystals. Substitution of Mg to Ni in the cluster model (centered at magnesium atom) is a rather trivial change, but for periodic models it is necessary to increase the size of the supercell, in order to remove or reduce the interaction between dopant atoms.

The drawback of the cluster approach in electronic structure calculations is the cluster border. Using bare clusters, such as ${[Mg{O}_6]}^{10-}$, without compensating the large formal charge, is insufficient for ionic crystals. Such clusters frequently display poor convergence patterns and might result in unphysical charge distributions. Placing a cluster into an electrostatic field in order to reproduce the Madelung potential from the infinite crystal, is not sufficient because of the presence of partially covalent bonds even in highly ionic crystals and potential “electron-leakage” effects Ref. Shi et al. (2022). Thus a typical solution for an embedding will include both terminating atoms or specially prepared potentials and an electrostatic field. We recently suggested a new protocol for computing model potentials, which can be used at the border of ionic clusters Larsson et al. (2022). In the current work, we will investigate the convergence of properties, related to the electronic structure, with increasing cluster sizes. For reference we will use band structure calculations, made with the same basis set and same DFT functional as in the cluster models. Although we focus on the proof that the electron density in the central part of embedded cluster is close to the electron density in periodic calculations at DFT level, our final goal is to apply accurate ab initio methods can be used to study the ground and excited states.

2 Materials and Methods

Selection of clusters for ionic crystals is not an easy task. Stoichiometric clusters have complicated border structures and are generally not suitable for systems with a point defect. Clusters constructed by adding spheres of neighboring atoms starve from having an uncompensated electronic charge in the quantum region. Taking a cluster such as ${[Mg{O}_6]}^{10-}$ as an example again, adding only point charges will only partially compensate the large formal charge, since the electron density in the quantum region can unphysically spread into the point charge region. In order to compensate for the effect of broken bonds at the border of the cluster, a separate layer can be introduced to the model. For covalent systems, a typical solution will include the use of terminating hydrogen atoms, while for ionic systems the termination can be achieved by placing specially prepared model potentials. These potentials can mimic the absent atoms, so the electron density inside the quantum region will be unperturbed by the cluster border. The presence of these potentials in the calculation only slightly increases the computational costs, since there are not necessarily basis set functions associated with these potentials.

The most reliable and straightforward method for construction of these potentials is AIMP (ab initio model potential) (Seijo and Barandiarán, 1999; Seijo and Barandiarán, 2013). In the Hartree-Fock equation, the Fock operator, $\widehat{f}$, is replaced by ${\widehat{f}}^{AIMP}=\widehat{f}-2\sum _k{ϵ}_k|{\psi }_k\rangle \langle {\psi }_k|$. Originally, the index k was used only to describe the core-states, making AIMPs a type of effective core potentials. Later, it was realised that by design, the AIMP methodology could be use to also freeze valence orbitals, leading to their use as embedding potentials Ref. Barandiarán and Seijo, (1988), Seijo and Barandiarán, (1999). Thus, in an embedded cluster calculation, k runs over all occupied orbitals of the ion described as an AIMP. The one electron energies of the states in the original formulation of AIMPs are defined by Hartree-Fock calculations. In Ref. Larsson et al. (2022) we demonstrated that it is also possible to use one electron energies obtained by DFT calculations, in particular with the PBE functional. The practical limitation of the method is related to the fact that these potentials depend on the crystal structure and should be recomputed for each new crystal structure. We recently suggested a procedure for such a purpose: a code named SCEPIC Larsson et al. (2022) allows one to construct the model potentials based on either HF or DFT. Using this code we constructed model potentials for several ionic crystals and demonstrated that the best quality for the embedding, at least for pure crystals, can be obtained with modified AIMPs, where the construction of the potentials uses PBE-based one electron energies. In the current study we will use the same model potentials, since the structure of the crystal is the same as it was used in Ref. Larsson et al. (2022). In addition to the layer of AIMPs, a set of point charges should be implemented to provide electroneutrality of the whole system. For the point charges, modeling an electrostatic field from the infinite crystal, we used the same strategy as in Sushko and Abarenkov, (2010): a finite set of charges were generated using an algorithm from Ref. Abarenkov, (2007).

The MgO crystal has a rock-salt structure with 6 oxygen atoms surrounding a magnesium atom. We can build a set of clusters, centered on an atom, by adding the layers from the next neighbours. The smallest cluster has a composition MgO₆. This cluster formally has a large electronic charge, however this charge will be compensated by terminating model potentials at the border and by a set of point charges, simulating the Madelung potential of a MgO crystal. Adding the next layer of neighboring atoms, we will get a cluster with composition MgO₆Mg₁₈, MgO₆Mg₁₈O₃₈ and finally MgO₆Mg₁₈O₃₈Mg₆₆. One should not expect monotonic changes in local properties for these incrementally built clusters, since the electronic charge, associated with the quantum part, as well as the atoms on the border of the cluster are alternating between Mg and O depending on cluster size. The nomenclature of these clusters will be on the form ${[Mg{O}_{6}M{g}_{18}]}^{26+}$ etc. to emphasize that the additional ions are, strictly speaking, not equivalent to the inner ions. Occasionally, a more extensive nomenclature on the form ${[M{g}^{(i)}{O}_6^{(i)}M{g}_{18}^{(o)}]}^{26+}$, when necessary to distinguish between inner (X⁽ⁱ⁾), outer (X^(o)) and “extra”(X^(†)) ions. Table 1 introduces a short notation for the clusters, and specifies the details of partitioning: the quantum size, the amount of AIMPs and point charges. The increase of size of these clusters exposes an important problem in quantum chemistry. Accurate computational methods scale poorly with respect to the increase of the basis set size. While DFT-like methods scale as N³ − N⁴, where N is the size of the basis set, multiconfigurational methods scale as N⁷. Table 1 also contains the number of basis functions for the PC-1 basis set (used in the current study in order to reach clusters with as many as possible layers) and for a triple-ζ quality basis set (ANO-RCC-VTZP), which should be used for accurate calculations.

TABLE 1

TABLE 1. Notations for the clusters, used in the study, indicating the quantum part, AIMPs and point charges. The basis set size is presented for PC-1 (used in the current study) and for a triple-ζ quality basis set.

Doping ionic metal oxides with transition metals or rare earth metals creates materials with interesting physicochemical properties. Many of these materials are active in chemical reactions, including catalytic. In practically interesting cases, the concentration of impurities is rather low. Ni: MgO is a compound where a Ni atom substitutes a Mg atom. Since the ionic radii of Mg²⁺ and Ni²⁺ are similar, we can construct corresponding clusters similar to those mentioned above with the center in magnesium atom. The smallest cluster, NiO₆, has no Mg atoms at all, but the AIMP layer, as earlier, contains model potentials for magnesium and oxygen. The larger clusters will contain layer(s) of magnesium. In periodic model it is important to use relatively large supercells to provide a spatial separation between the defects.

Each cluster was surrounded by AIMPs up to a cutoff radius of the longest distance between two real atoms in the QM region plus the length of one unit-cell (4.26 Å) surrounding the central cluster, making the cluster plus AIMP-region spherical in shape. For instance, in ${[Mg{O}_6]}^{10-}$, the longest distance is 2.13 Å (i.e., the Mg–O bond length in MgO), making the radius of the sphere 6.39 Å. This was made primarily for convenience; Ref Larsson et al. (2022) showed that the cluster region should be surrounded by at least one unit-cell layer of AIMPs for local properties to converge. This is further supported by Figure 1, which shows that at a range of beyond c:a 6 Å, AIMPs start to behave like a classical charge with a finite volume. Therefore, having more AIMPs than is necessary does not influence the results of the calculation, only the evaluation of the one-electron Hamiltonian which becomes somewhat slower with increasing the number of AIMPs Ref Larsson et al. (2022). Since they behave primarily as charges beyond a certain range, the total charge of the AIMP region bears no influence on the calculation, and the only important charge to consider is that of the QM region itself. The total structure including the point-charge region was constructed using a radius of 50 Å; the total structure is, however, not a sphere, since with the GPEE method, GPEE-modified unit-cells up to a range of 50 Å were included, rather than ions up to that radius. Multipole moments up to third order were canceled by the GPEE method.

FIGURE 1

FIGURE 1. Comparison of the potential from a set of six point charges (black), six Mg²⁺-AIMPs (red) and six Ca²⁺-AIMPs (green) surrounding a O²⁻-ion. The surrounding charges/AIMPs are placed in an octahedral field and stretched symmetrically.

For the smallest cluster, ${[Mg{O}_6]}^{10-}$, the procedure described in the previous paragraph leads to an AIMP region of ${[M{g}_{54}^{AIMP}{O}_{56}^{AIMP}]}^{4-}$ and point-charge region with a total charge of +14 e in order to provide an electroneutrality of the whole system. The AIMP region increased to ${[M{g}_{122}^{AIMP}{O}_{110}^{AIMP}]}^{24+}$ for the ${[Mg{O}_{6}M{g}_{18}]}^{26+}$-cluster, with the point-charge region having a total charge of -50 e. Next, ${[Mg{O}_{6}M{g}_{18}{O}_{33}]}^{50-}$, an AIMP region of ${[M{g}_{230}^{AIMP}{O}_{192}^{AIMP}]}^{76+}$ and point-charge region of total charge -26 e was used. And finally the largest cluster is ${[Mg{O}_{6}M{g}_{18}{O}_{33}M{g}_{66}]}^{82+}$ with an AIMP region ${[M{g}_{368}^{AIMP}{O}_{422}^{AIMP}]}^{108-}$ and point-charge region +26 e. For these clusters we used nomenclature Mg-I, Mg-II, Mg-III and Mg-IV, where the Roman number specifies the number of Mg spheres, as shown in Table 1. At Figure 2 the structure of Mg-IV cluster is presented in a visual form: the quantum part, the layer of AIMPs and finally a cloud of point charges. For Mg-IV cluster calculations the disk space for temporary storage was about 210Gb, making further increases of the cluster size unfeasible. The O-centered clusters are identical in construction to the Mg-centered, and can be derived by simply inverting the positions of Mg and O and changing the signs of all charges. The nomenclature of these clusters is O-I, O-II and O-III. The largest cluster with a center at an oxygen atom (O-IV) is not included in the results, because in this case we did not manage to obtain a stable solution for the wavefunction calculation. Ni-centered clusters were obtained by simply replacing the central Mg with Ni, without any changes in geometry. The nomenclature for these clusters is from Ni-I to Ni-IV.

FIGURE 2

FIGURE 2. The structure of Mg-IV cluster. From left to right: the quantum part, AIMPs layer (model potentials are shown as pale spheres), point charges surrounding the embedded cluster.

All electronic structure calculations in this work utilized the Perdew–Burke–Ernzerhof (PBE) DFT functional with the pc-1 basis set Jensen, (2002). Cluster calculations were performed with OpenMolcas, with the embedding method described in Ref. Larsson et al. (2022). Periodic calculations were made with CP2K, using all-electron basis sets within the Gaussian and augmented plane wave (GAPW) framework, with a plane wave cutoff of 400 eV. The periodic calculations were based only on the Γ-point, using both a 2 × 2 × 2 supercell and a 3 × 3 × 3 supercell.

2.1 Electron Density Comparison

The electron density is a direct and a straightforward way to compare different cluster and periodic model calculations. Ideally, the electron density in a large enough cluster should be identical or close enough to the electron density obtained from a periodic calculation. Unphysical effects arising from the cluster–embedding border can be very strong and influence the central part of the quantum region. Simply increasing the cluster size might, however, not necessarily improve upon the results, as the number of broken bonds at the border of the cluster increases with cluster size. In principle, a large enough bare cluster would be the most realistic model of an ionic crystal. Unfortunately, the size required to reproduce the Madelung field from the extended crystal will in general prohibit such a route. One way to circumvent these problems is to use specially constructed potentials, designed to behave as similar as possible to the real atoms, at the border between the quantum region and the electrostatic embedding.

There are two options for comparing the electron density. First option includes the comparison of so-called local properties of the electronic structure: charges (Q), bond orders (W), covalency (C) and total valency (V). All these properties are computed from the population matrix (DS), where D is the density matrix and S is the overlap matrix. For molecular and cluster calculations the size of the matrix is determined by basis set size, but for periodic calculations it is determined by the basis set size in the supercell used in the calculation Veryazov et al. (1999). Although these properties depend on the atomic basis set, with a fixed basis set they can be used as an indicator of similarities between two electron densities. Atomic charges defined as a number of electrons, involved in ionic bonding and it is defined as a difference between nuclear charge Z_A and the electron charge localized at the atomic region $(Q_A=Z_A-\sum _{a\in A}{(DS)}_{ii})$, where a is a molecular orbital from atom A. The bond order between atoms A and B $(W_{AB}=\sum _{a\in A}\sum _{b\in B}|{(DS)}_{ij}{|}^2)$ corresponds to the number of electrons involved in the bond between the atoms Wiberg, (1968). Covalency $(C_A=\sum _{B\ne A}{W}_{AB})$ counts the number of electrons involved in covalent bonds Armstrong et al. (1973), and finally, the total valency ($V_A=\frac{1}{2}(C_A+\sqrt{C_A^2+4_{*}|Q_A{|}^2})$ is a property, which counts all electrons from a certain atom, that are involved in the creation of chemical bonds in the system Evarestov and Veryazov, (1991).

To compute local properties from the wavefunctions produced by OpenMolcas and CP2K we used homemade scripts. For ionic crystals atomic charges are similar to formal charges, but never perfectly match them. The contributions to the covalent bonds are small and local: all bond indices between non-neighboring atoms are close to zero. It is important to compare the combination of all local properties to confirm the conclusion about the match between different electron densities.

Another alternative is to use spatial electron density ρ(r). Both OpenMolcas and CP2K can produce spatial electron density computed on a Cartesian grid. Since the electron density is a function of coordinates, the comparison can be performed in some arbitrary volume. For that purpose we used the union of seven spheres centered on the atomic positions of the MgO₆ cluster, with radii of 1.6Å (3/4 of MgO bond length). This volume represents the density in the central part of the studied clusters, however for the smallest clusters (Mg-I, O-I) the border of this volume is too close to the cluster border. Different computational codes can use different approximations in computing electronic density in the regions close to atoms, for instance, the GAPW framework in CP2K, used in all-electron calculations, separate the electron density into a core (hard) density and valency (soft) density Ref Hutter et al. (1999). Only the valence region is dealt with in a similar fashion to standard molecular quantum chemistry codes, such as OpenMolcas. Thus, we removed spherical regions with the radius 0.8Å around each atom, which corresponds to the default value for the core density in CP2K.

Using Cartesian grid in the volume described above, we introduce the root mean square deviation (RMSD) of the density: ${\rho }^{RMSD}=\sqrt{\frac{{\sum }_i^n{\left({\rho }_i^{\text{cluster}}-{\rho }_i^{\text{periodic}}\right)}^2}{n}}$, where the index i runs over all points on a grid, which is about 200000 points.

3 Results and Discussion

3.1 Local Property Descriptors in MgO

In this section we present the results for clusters centered on magnesium (Mg-I, Mg-II, Mg-III and Mg-IV) and clusters centered on oxygen (O-I, O-II, O-III). Calculation of the electronic structure of O-IV cluster diverges, for which reason it was excluded from this work. Periodic calculations were made for with both a 2 × 2 × 2 and a 3 × 3 × 3 MgO supercell.

The local properties of the electronic structure of MgO in periodic and embedded cluster models are presented in Table 2 Comparison of local properties between periodic calculations with different supercell sizes shows discrepancy for the individual descriptors, which one can consider as a threshold. For all local properties the difference is either in the second or in the third decimal point.

TABLE 2

TABLE 2. Electronic structure of clusters with Mg- on O- center. Table showing the convergences of Mulliken charges (Q), Wiberg bond indices (W), covalency (C), total valency (V) and integrated density difference Δρ (multiplied by a factor of 10³) with respect to periodic values. * indicates that the value is undefined for the particular cluster.

The smallest clusters, Mg-I and O-I have only one (central) atom with a proper atomic surrounding, thus, it is rather expected that the electronic structure even for the central region of these clusters is different within Mg-I–Mg-IV and O-I–O-III series. The atomic charge on the central Mg atom in Mg-I is 1.02e, while for larger clusters it is in the range 1.13–1.28 e. For the smallest clusters (Mg-I and O-I) the charges on the central atoms are ‘non symmetric’ for Mg in Mg-I it is 1.02 e, but for O in O-I it is −1.20 e. The difference between atomic charges on the central atoms for the larger clusters in by order of magnitude smaller (1.166 e for Mg-IV vs. −1.175 e for O-III). If we exclude the smallest clusters from consideration, the atomic charges are rather close to the charges in the periodic model. The difference between the absolute values of atomic charges on the central atom and on the nearest neighbor is also reducing with the increase of the cluster size. Increasing the cluster sizes from Mg-I to Mg-IV and from O-I to O-III does not result in a monotonic change in atomic charges. Instead of correlating with whether the central atom is Mg or O, the atomic charges correlate with the formal charge of the cluster. Thus, for a monotonic charge the clusters should be viewed in the order of Mg-I to O-II and finally Mg-III, which are all anionic clusters.

On the other hand, the variations of the atomic charges are relatively small: 0.02 e for oxygen atom in clusters O-II and O-III, and 0.1 e for Mg in clusters Mg-II–Mg-IV. The largest deviation is observed for Mg-III cluster, which in the quantum part has the largest negative charge.

The bond orders, also known as Wiberg indices, between Mg and O show only a small variation between clusters and they are similar to periodic values. Again the largest difference from the general trend is observed for the Mg-III cluster. Considering that the lowest bond index for Mg–O bond and largest atomic charge on Mg atom was obtained for the Mg-III cluster, we can conclude that this cluster is the has most ionic character of the clusters studied here. The covalency of central Mg atom in Mg-IV cluster (1.46) is almost the same as the covalency of Mg in periodic calculations (1.45). The covalency of O in O-II and O-III (1.41) is also close to corresponding value (1.45) in periodic calculations. Total valence of Mg and O in all clusters, except the smallest ones, is close to the corresponding values in periodic calculations.

In contrast, the RMSD of the electron density difference gave a smooth convergence towards the periodic limit for the central XY₆ unit. For Mg-I, Mg-II and O-I clusters, the error is in the range of 0.012–0.015 $e/a_0^3$, but for larger clusters it is in the range of 0.002–0.005 $e/a_0^3$, which is not larger than the corresponding difference between densities obtained for different supercell extensions. In Figure 3, the spatial electron density used in the comparison is visualised, along with density differences plots for the Mg-I, Mg-II and Mg-III clusters.

FIGURE 3

FIGURE 3. Total and difference density plots for MgO. (A) Total density retained from the 3 × 3 × 3 periodic calculation of MgO, only a small part of the total structure visualised for brevity, (B) density difference between the ${[Mg{O}_6]}^{10-}$-cluster and the 3 × 3 × 3 calculation, (C) density difference between the ${[Mg{O}_{6}M{g}_{18}]}^{26+}$-cluster and the 3 × 3 × 3 calculation and (D) density difference between the ${[Mg{O}_{6}M{g}_{18}{O}_{33}]}^{50-}$-cluster and the 3 × 3 × 3 calculation. A wireframe model was used for (D) in order to make the isosurface more visible. Isosurface values of 0.04 e/${\mathrm{a}}_0^3$. Yellow isosurfaces correspond to positive values and cyan to negative values. Blue spheres correspond to Mg and red spheres to O.

3.2 Local Electronic Structure of Ni:MgO

Studies of single-ion dopants has been a major undertaking for AIMP studies. In this section, the convergence of the local electronic structure of Ni and O in Ni: MgO is presented. As in the previous sections, the results will be written for all three clusters on the order of Ni-I, Ni-II and Ni-III. We should note here that the calculations of these clusters were made with UHF Hamiltonian with triplet spin multiplicity. Periodic calculations were made by substituting a single Mg with Ni in both the 2 × 2 × 2 and 3 × 3 × 3 MgO supercells.

The local properties of the electronic structure of Ni: MgO in periodic and embedding cluster model are presented in Table 3. The difference between 2 × 2 × 2 and 3 × 3 × 3 periodic calculations is slightly larger, compared to pure MgO due to interaction between Ni atoms. Cluster models with a single nickel atom is free from this problem. Similar to MgO clusters, the smallest cluster (Ni-I) differs the most from the other clusters. The atomic charge on Ni increases when increasing the cluster size to Ni-II and Ni-III, but drops for the Ni-IV cluster, and becomes similar to the atomic charge in both periodic models. The atomic charge on Ni is smaller than the corresponding charge on Mg. But at the same time the covalent bonds between Ni–O are stronger than Mg–O bond. That results in larger, in comparison to Mg, covalency of nickel atom and the total valency 2.3. The ‘elevated’ total valency for nickel atom is not an unexpected result, since Ni atom can be found in different valence states Evarestov et al. (1994).

TABLE 3

TABLE 3. Electronic structure of clusters Ni:MgO. Table showing the convergences of Mulliken charges (Q), Wiberg bond indices (W), covalency (C), total valency (V) and integrated density difference Δρ (multiplied by a factor of 10³) with respect to periodic values. * indicates that the value is undefined for the particular cluster.

Comparison of electron density shows a similar trend as for MgO clusters. The electron density in Ni-I and Ni-II clusters is slightly different from density in the periodic calculations. But for the Ni-III and Ni-IV clusters, the difference is negligibly small.

In Figure 4, the projected density of states (PDOS) for Ni and O are plotted. The d-states of Ni are independent of the computational model used. In the α spin channel, the double peak closest to the Fermi level corresponds to the split between the t_2g and e_g Ni d-orbitals. Since no e_g orbital is occupied in the β channel, there is only a single t_2g peak near the Fermi level. The difference between the t_2g and e_g α peaks change as 0.72 eV, 0.67–0.69 eV when increasing the model size. From the periodic calculations, a split of 0.70 eV is predicted.

FIGURE 4

FIGURE 4. Projected density of states from clusters (A) ${[{\mathrm{M}\mathrm{g}\mathrm{O}}_6]}^{10-}$, (B) ${[{\mathrm{M}\mathrm{g}\mathrm{O}}_6{\mathrm{M}\mathrm{g}}_{18}]}^{26+}$, (C) ${[{\mathrm{M}\mathrm{g}\mathrm{O}}_6{\mathrm{M}\mathrm{g}}_{18}{\mathrm{O}}_{38}]}^{50-}$ and (D) periodic PBE calculations on Ni: MgO. The α spin channels are plotted with positive amplitudes and β spin channels with negative amplitudes. The Fermi level of the α spin channel is set as zero. D.O.S stands for density of states and is given in arbitrary units. The colour scheme of the PDOS is as follows: black—total, red—s, green—p and blue—d. The PDOS were obtained via Gaussian smearing of the canonical molecular orbital energies, using a smearing width of 1 meV.

In addition to these two major peaks, there are several minor peaks at lower energies originating from the Ni d-orbitals. Comparison with the O PDOS at the same energy region suggests these are primarily caused by the Ni d-functions acting as additional polarizing basis functions for the oxygen p-band. Overall, larger changes are visible in the O PDOS compared to the Ni PDOS; primarily a broadening of the PDOS is observed when increasing the cluster size. This is an expected consequence of including more “real” oxygens in the model.

The advantage of embedded cluster model in comparison to periodic model is that all methods used in molecular quantum chemistry can be applied. It is very important for the properties, involving excited states, in particular excitation spectra. In Larsson et al. (2022) we presented the vibronic d–d transition in Ni: MgO, which was modeled as Ni-II embedded cluster. There, it was demonstrated that when using multiconfigurational theory and appropriate basis sets (ANO-RCC), the experimental spectrum of Ni: MgO could be reproduced to within 1000 cm⁻¹ on average.

Obviously, one advantage of a material such as MgO is that, in comparison to transition metal oxides (e.g., TiO₂), it is a relatively simple material. Adding additional layers of Mg or O does not increase the number of orbitals that might influence, for instance, the spectrum of single ion dopants. In materials like TiO₂, the presence of valence d-orbitals, even if they are formally empty, potentially complicates the picture dramatically. More extensive studies are necessary to verify the applicability of multiconfigurational methods in such systems. At the present, the AIMP methodology is presented as a straightforward and accurate way of modelling ionic materials made from s or p group elements.

3.3 Computational Costs

The computational cost to compute DFT density for these clusters consists of computing the integrals (in the compact form of RI/CD technique implemented in OpenMolcas) and self consistent field iterations. Even with the use of RI/CD technique the size of the temporary disk space is large and can easily reach several hundreds of Gb. The calculation of integrals is approximately one third of the total time. The size of AIMPs and electrostatic field has no effect on the timing. All calculations were performed at Tetralith cluster, with Intel Xeon Gold 6130 CPU and SSD hard drive. The calculations with 8 cores of Mg-I cluster take less than a minute, Mg-II—2 min, Mg-III—14 min and finally Mg-IV—2 h. The total timing also depends on the convergence of the optimal wavefunction. We should note here that a calculation with the periodic boundary conditions (performed with CP2K code with the same basis set and the same computational setup) takes even longer time. The computational resources will increase if a more accurate basis will be used, see Table 1. Using multiconfigurational theory for accurate calculations of the electronic structure will increase the timing in the order of magnitude.

4 Conclusion

Using PBE Hamiltonian and a relatively small basis set we investigated the convergence of the electron density with an increase of the size of embedded clusters. The results shown herein demonstrate that good agreement in both local properties and spatial electron density can be achieved for moderately sized clusters using AIMP embedding. Using the spatial electron density to assess the clusters gives a clearer convergence with increasing cluster size, compared to the local properties. One advantage with the local properties, however, is that they are generally faster to compute and therefore give an important first quality check of the clusters.

The smallest clusters considered accurate here (Mg-II, O-II and Ni-II) all achieve the following two criteria: 1) no single element in the material is described only by AIMPs and 2) in the central region, the co-ordination of all ions are completely saturated. This suggests that in future studies, minimal clusters should be constructed based on such criteria.

While the largest discrepancies were noted for the most anionic clusters (Mg-III and Ni-III, both with a formal charge of −50 e), the values are seemingly not far enough from periodic values to completely invalidate the use of anionic clusters with AIMPs. Even more anionic clusters, however, seem to give poor convergence and are thus discouraged from use in production calculations.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

VV thanks eSSENCE@LU grant 7:3 for financial support.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We acknowledge support from Swedish National Infrastructure for Computing (SNIC). The computations were enabled by resources provided by SNIC at Tetralith partially funded by the Swedish Research Council through grant agreement no. 2018-05973 and LUNARC computer center.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2022.951144/full#supplementary-material

References

Abarenkov, I. V. (2007). Unit Cell for a Lattice Electrostatic Potential. Phys. Rev. B 76, 165127. doi:10.1103/physrevb.76.165127

CrossRef Full Text | Google Scholar

Armstrong, D. R., Perkins, P. G., and Stewart, J. J. P. (1973). Bond Indices and Valency. J. Chem. Soc. Dalton Trans., 838–840. doi:10.1039/dt9730000838

CrossRef Full Text | Google Scholar

Arunadevi, N., and Kirubavathy, S. J. (2022). “Chapter 2 - Metal Oxides: Advanced Inorganic Materials,” in Inorganic Anticorrosive Materials. Editors C. Verma, J. Aslam, and C. M. Hussain (Elsevier), 21–54. doi:10.1016/b978-0-323-90410-0.00002-7

CrossRef Full Text | Google Scholar

Awasthi, S., Pandey, S. K., Arunan, E., and Srivastava, C. (2021). A Review on Hydroxyapatite Coatings for the Biomedical Applications: Experimental and Theoretical Perspectives. J. Mater. Chem. B 9, 228–249. doi:10.1039/d0tb02407d

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayala, P. Y., Kudin, K. N., and Scuseria, G. E. (2001). Atomic Orbital Laplace-Transformed Second-Order Møller-Plesset Theory for Periodic Systems. J. Chem. Phys. 115, 9698–9707. doi:10.1063/1.1414369

CrossRef Full Text | Google Scholar

Barandiarán, Z., and Seijo, L. (1988). The Ab Initio Model Potential Representation of the Crystalline Environment. Theoretical Study of the Local Distortion on NaCl:Cu+. J. Chem. Phys. 89, 5739–5746. doi:10.1063/1.455549

CrossRef Full Text | Google Scholar

Canney, S. A., Sashin, V. A., Ford, M. J., and Kheifets, A. S. (1999). Electronic Band Structure of Magnesium and Magnesium Oxide: Experiment and Theory. J. Phys. Condens. Matter 11, 7507–7522. doi:10.1088/0953-8984/11/39/308

CrossRef Full Text | Google Scholar

Di Cosimo, J. I., Díez, V. K., Ferretti, C., and Apesteguía, C. R. (2014). Chapter 1 Basic Catalysis on MgO: Generation, Characterization and Catalytic Properties of Active Sites. Catalysis 26, 1–28. doi:10.1039/9781782620037-00001

CrossRef Full Text | Google Scholar

Evarestov, R. A., Tupitsin, I. I., and Veryazov, V. A. (1994). Quantum Chemical Calculation of Nickel and Copper Atomic Valencies in Crystalline Oxides. Int. J. Quantum Chem. 52, 295–299. doi:10.1002/qua.560520204

CrossRef Full Text | Google Scholar

Evarestov, R. A., and Veryazov, V. A. (1991). Quantum-chemical Definition of the Atomic Valence in Molecules and Crystals. Theor. Chim. Acta 81, 95–103. doi:10.1007/bf01113380

CrossRef Full Text | Google Scholar

Heo, S., Cho, E., Lee, H.-I., Park, G. S., Kang, H. J., Nagatomi, T., et al. (2015). Band Gap and Defect States of MgO Thin Films Investigated Using Reflection Electron Energy Loss Spectroscopy. AIP Adv. 5, 077167. doi:10.1063/1.4927547

CrossRef Full Text | Google Scholar

Hirata, S., Podeszwa, R., Tobita, M., and Bartlett, R. J. (2004). Coupled-cluster Singles and Doubles for Extended Systems. J. Chem. Phys. 120, 2581–2592. doi:10.1063/1.1637577

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutter, J., Parrinello, M., and Lippert, G. (1999). The Gaussian and Augmented-Plane-Wave Density Functional Method for Ab Initio Molecular Dynamics Simulations. Theor. Chem. Accounts 103, 124–140. doi:10.1007/s002149900042

CrossRef Full Text | Google Scholar

Jensen, F. (2002). Polarization Consistent Basis Sets. II. Estimating the Kohn–Sham Basis set Limit. J. Chem. Phys. 116, 7372–7379. doi:10.1063/1.1465405

CrossRef Full Text | Google Scholar

Jolivet, J.-P., Henry, M., and J, L. (2000). Metal Oxide Chemistry and Synthesis: From Solution to Solid State. Chichester, UK: Wiley.

Google Scholar

Lany, S. (2014). “Polymorphism, Band-Structure, Band-Lineup, and Alloy Energetics of the Group II Oxides and Sulfides MgO, ZnO, CdO, MgS, ZnS, CdS,” in Oxide-based Materials and Devices V. Editors F. H. Teherani, D. C. Look, and D. J. Rogers (International Society for Optics and Photonics SPIE). doi:10.1117/12.2043587

CrossRef Full Text | Google Scholar

Larsson, E. D., Krośnicki, M., and Veryazov, V. (2022). A Program System for Self-Consistent Embedded Potentials for Ionic Crystals. Chem. Phys., 111549. doi:10.1016/j.chemphys.2022.111549

CrossRef Full Text | Google Scholar

McClain, J., Sun, Q., Chan, G. K.-L., and Berkelbach, T. C. (2017). Gaussian-based Coupled-Cluster Theory for the Ground-State and Band Structure of Solids. J. Chem. Theory Comput. 13, 1209–1218. doi:10.1021/acs.jctc.7b00049

PubMed Abstract | CrossRef Full Text | Google Scholar

Mojica-Sánchez, J. P., Zarate-López, T. I., Flores-Álvarez, J. M., Reyes-Gómez, J., Pineda-Urbina, K., and Gómez-Sandoval, Z. (2019). Magnesium Oxide Clusters as Promising Candidates for Hydrogen Storage. Phys. Chem. Chem. Phys. 21, 23102–23110. doi:10.1039/c9cp05075b

PubMed Abstract | CrossRef Full Text | Google Scholar

Niedermeier, C. A., Råsander, M., Rhode, S., Kachkanov, V., Zou, B., Alford, N., et al. (2016). Band Gap Bowing in NixMg1−xO. Sci. Rep. 6, 31230. doi:10.1038/srep31230

PubMed Abstract | CrossRef Full Text | Google Scholar

Schönberger, U., and Aryasetiawan, F. (1995). Bulk and Surface Electronic Structures of MgO. Phys. Rev. B 52, 8788–8793. doi:10.1103/physrevb.52.8788

CrossRef Full Text | Google Scholar

Seijo, L., and Barandiarán, Z. (1999). The Ab Initio Model Potential Method: A Common Strategy for Effective Core Potential and Embedded Cluster Calculations. Comput. Chem. Rev. Curr. Trends 4, 55–152. doi:10.1142/9789812815156_0002

CrossRef Full Text | Google Scholar

Seijo, L., and Barandiarán, Z. (2013). Host Effects on the Optically Active 4f and 5d Levels of Ce3+ in Garnets. Phys. Chem. Chem. Phys. 15, 19221–19231. doi:10.1039/c3cp53465k

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, B. X., Kapil, V., Zen, A., Chen, J., Alavi, A., and Michaelides, A. (2022). General Embedded Cluster Protocol for Accurate Modeling of Oxygen Vacancies in Metal-Oxides. J. Chem. Phys. 156, 124704. doi:10.1063/5.0087031

PubMed Abstract | CrossRef Full Text | Google Scholar

Srivastava, A. K., Misra, N., and Pandey, S. K. (2018). DFT Study on Planar (CaO) N Rings (N = 1-5) and Their Hydrogen Storage Behavior: Ca-O versus Mg-O Clusters. J. Clust. Sci. 29, 57–65. doi:10.1007/s10876-017-1306-x

CrossRef Full Text | Google Scholar

Sushko, P. V., and Abarenkov, I. V. (2010). General Purpose Electrostatic Embedding Potential. J. Chem. Theory Comput. 6, 1323–1333. doi:10.1021/ct900480p

CrossRef Full Text | Google Scholar

Veryazov, V. A., Leko, A. V., and Évarestov, R. A. (1999). Local Characteristics of Crystal Electronic Structure in the Hartree-Fock Method. Phys. Solid State 41, 1286–1290. doi:10.1134/1.1130983

CrossRef Full Text | Google Scholar

Wiberg, K. B. (1968). Application of the Pople-Santry-Segal Cndo Method to the Cyclopropylcarbinyl and Cyclobutyl Cation and to Bicyclobutane. Tetrahedron 24, 1083–1096. doi:10.1016/0040-4020(68)88057-3

CrossRef Full Text | Google Scholar