Density functional theory and thermodynamics analysis of MAl₁₂ Keggin substitution reactions: Insights into ion incorporation and experimental confirmation

Polyoxometalates (POMs) are a class of metal oxide nanoclusters that have physical and chemical properties that can be influenced by altering the identity of one or more constituent atoms.^1–5 For example, Villa et al. studied the effects of atomic substitution with Ti(IV) into [H_xNb₁₀O₂₈]^6−x to generate [H_xTiNb₉O₂₈]^7−x and [H_xTi₂Nb₈O₂₈]^8−x; they found that composition influenced oxygen-isotope exchange rates by at least four orders of magnitude.^6,7 Of the variety of known POM structure types, the Baker–Figgis–Keggin (Keggin) is one of the most characterized and well-known topologies. Keggin-type nanoclusters have been identified as both cationic and anionic nanoclusters, with a range of applications in heterogeneous catalysts, medicine, sensors, pillaring agents, and water remediation devices.^8–15 

Aluminum Keggin polyoxocations have the general formula of [MO₄Al₁₂(OH)₂₄(H₂O)₁₂]^4+m (MAl₁₂^4+m), where the charge m+ depends on the identity/oxidation state of the central, tetrahedral heteroatom M. The Keggin can form five possible isomers (α, β, γ, δ, and ε) by a 60° rotation of trimeric [Al₃(μ₂−OH)₆(H₂O)₃] units, altering the connectivity between two trimers to be either through shared edges or corners via bridging oxo/hydroxo groups. The orientation of the trimers influences the geometry and electronic structure of the Keggin topology, as previously documented in a study by Tézé et al., comparing bridging oxo angles within the α, β, and γ isomers of [SiW₁₂O₄₀]ⁿ⁻ salts.¹⁶ López et al. compared relative energies of different isomers of [XW₁₂O₄₀]³⁻ (X = P, Si, Al), finding that increasing the number of shared edges resulted in decreased cluster stability for the polyoxoanions.^17,18

There are the five possible Al₁₃ Keggin isomers, but only three have been isolated and characterized to date. The ε-Al₁₃⁷⁺ isomer is considered to be the most stable and is synthesized by the hydrolysis of AlCl₃ with NaOH and either selenate or sulfate salts as crystallization agents.¹⁹ The δ-isomer was first isolated and characterized by Rowsell and Nazar, formed by heating aqueous AlCl₃ · 6H₂O at 80 °C.²⁰ Only the δ-Al₁₃ isomer has been shown to nucleate and form larger derivative species such as [Al₂₆O₈(OH)₅₀(H₂O)₂₀]¹²⁺ (Al₂₆),²¹ [Al₃₀O₈(OH)₅₆(H₂O)₂₆]¹⁸⁺ (Al₃₀), and [Al₃₂O₈(OH)₆₀(H₂O)₂₈(SO₄)₄]¹⁶⁺ (Al₃₂), which are important species in the water purification agent polyaluminum chloride (PAC).^10,11,22 More recently, in 2013, Smart et al. synthesized the γ-isomer by adding Ca(OH)₂ to a AlCl₃ and glycine solution at 95 °C followed by hydrothermal heating.²³ Our previous work has been focused on understanding the reactivity and transformation (isomerization and oligomerization) of Keggin-type Al nanoclusters.^24–29 

Of the three isolable Al₁₃ Keggin isomers, the ε-isomer been shown to prefer substitution at the central, tetrahedral heteroatom M site. The first heteroatom analog identified by Bradley et al. was ε-GaAl₁₂⁷⁺.^30,31 Parker et al. attempted to incorporate other metal species (M = Mg²⁺, Fe³⁺, Co²⁺, Ni²⁺, Cu²⁺, Zn²⁺, In³⁺, La³⁺, and Ce³⁺), but no new substituted analogs were obtained under their synthesis conditions.³² In 2001, Lee et al. discovered a second analog with Ge⁴⁺ as the heteroatom, forming ε-GeAl₁₂⁸⁺.³³ More recently, our experimental-theoretical collaboration has led to the discovery of [δ-GaO₄Al₁₂(OH)₂₄(H₂O)₁₂]⁷⁺ (δ-GaAl₁₂) with Ga³⁺ occupying the tetrahedral site.³⁴ Density functional theory (DFT) calculations reveal a 0.5 eV preference for Ga³⁺ substitution at the tetrahedral site over the octahedral sites, in agreement with the observed experimental structure for both δ-GaAl₁₂ and the related [Ga_2.5O₈Al_28.5Ga_0.5(OH)₅₈(H₂O)₂₇(SO₄)₂]¹⁹⁺ cluster. Calculations indicate preference (greater than 1 eV) for Ga³⁺ to occupy both tetrahedral sites and one of the beltway octahedral sites, in agreement with the experimental structure.³⁴ 

In order to guide experiment to the discovery of other substituted Keggin nanoclusters, Reusser, Casey, and Navrotsky combined experimental calorimetry and DFT calculations to examine the energetics of heteroatom substitution.³⁵ The energetic comparison used M(OH)_4(aq) as the heteroatom starting material and showed that most metal cations would substitute for Ge⁴⁺ in ε-GeAl₁₂⁸⁺. Their structural analysis revealed that the size of the heteroatom influences the stability of the Keggin. After completing a mineral survey of tetrahedrally-coordinated M–O bond lengths, they developed structural criteria for likely heteroatom analogs. Following their criteria, they suggest that ε-ZnAl₁₂⁶⁺ and ε-FeAl₁₂⁷⁺ are the most likely heteroatom analogs to be created. To the best of our knowledge, no experimental work to date has confirmed these predictions. To be able to predict and understand the underlying thermodynamics associated with substitution within the MAl₁₂ Keggin series, it is critical to take into account the solvated cation and cluster properties. A combination of DFT calculations and thermodynamics using experimental information provides an improved description of the aqueous cation chemistry that governs MAl₁₂ nanocluster substitution reactions.

By combining electronic structure calculations with thermodynamics modeling, we aim to adapt and experimentally verify a predictive methodology for M^m+ incorporation into Keggin clusters, thus describing the energetics of cation substitution. Our previous work reported on a DFT + thermodynamics analysis to describe the energetics of substituting Al³⁺ cations with Ga³⁺ in δ-Al₁₃ and Al₃₀ parent models, finding that substitution at the heteroatom site is the most favorable.³⁴ In this current work, we first model exchange of the metal occupying the tetrahedral site at the DFT level and combine this with experimental measurements of metal oxidation and hydration. We analyze the known and predicted ε-MAl₁₂ structural features and energetics to layout design rules and guide synthesis for new nanoclusters. We go on to apply our model to other known substituted Keggins. For example, previous experimental and theoretical work by Wang et al. has characterized Zn²⁺ and Cr³⁺-substituted δ-Keggins [δ-Zn(CrAl)₁₂ and δ-ZnCr₁₂], where in both clusters, Cr was only found in the octahedral positions and Zn was found in the heteroatom position.^36,37 Here, we apply our methodology to the Al³⁺/Zn²⁺/Cr³⁺ system of Wang et al. to understand the site-specific substitution observed in the δ-Keggin structure and compare possible substitution pathways. Chromium is a well-known environmental contaminant that readily substitutes into Fe³⁺ and Al³⁺ mineral phases^38–40 but has yet to be observed as the heteroatom in Keggin-type aluminum polycations. Based upon the theoretical results, we then further evaluate the Al³⁺/Cr³⁺ system experimentally, which led to the synthesis, structure determination, and chemical characterization of two mixed Cr/Al δ-Keggin structures, Na[AlO₄Al_9.6Cr_2.4(OH)₂₄(H₂O)₁₂](2,6-NDS)₄(H₂O)₂₂ (δ-Cr_nAl_13-n-1) and Na[AlO₄Al_9.5Cr_2.5(OH)₂₄(H₂O)₁₂](2,7-NDS)₄(H₂O)_18.5 (δ-Cr_nAl_13-n-2). These results highlight the predictive nature of the computational methodology for the isolation of new cation substitutions within polyaluminum phases.

Crystal structures of ε-Al₁₃, ε-GaAl₁₂, and ε-GeAl₁₂ were obtained from the work of Johansson,¹⁹ Parker,³² and Lee,³³ respectively; these crystal structures served as the initial input for the ε-MAl₁₂ atomic coordinates. For the calculations involving Cr³⁺ and Zn²⁺ substitution into δ-Keggins, the δ-Al₁₃ crystal structure was used as a reference.^21,41 Using previously established methods for studying aluminum nanoclusters,^24,26,28,29 non-periodic DFT calculations were performed at the generalized-gradient-approximation Perdew–Burke–Ernzerhof (GGA-PBE)⁴² level, sampling the energy at the gamma point. Previous benchmarking studies have compared the performance of the PBE exchange-correlation functional and the hybrid Becke, 3-parameter, Lee–Yang–Parr (B3LYP)^43,44 functional in describing the large Al₃₀ nanocluster complexed with Cu²⁺ (238 atoms total: 2 Cu, 30 Al, 94 O, and 112 H). The results showed relative energies that differed by at most 0.1 eV or 0.003 eV per Al atom.²⁴ The double-numeric-plus polarization (DNP), atom-centered basis set was used as implemented in the DMol³ package developed by Delley.^45–47 The DNP basis set has been shown to produce relative energies for Al₁₃ Keggin isomers that differ from those obtained with the triple-numeric (TNP) basis set ∼0.01 eV.³⁴ All geometry optimization calculations were converged to at least 1 × 10⁻⁵ eV/Å with a cutoff radius of 4.5 Å. The conductor-like screening model (COSMO) was used to simulate aqueous solvent effects.⁴⁸ 

A DFT + solvent ion⁴⁹ method has been used to compute the change in free energy for ion dissolution in complex metal oxides^50–53 and is adapted here to describe the energetics of heteroatom substitution (∆G_sub) reactions in Keggin nanoclusters. This method has been applied to understand Ga³⁺ substitution site preferences in δ-GaAl₁₂ and Ga_2.5Al_28.5 nanoclusters.³⁴ In this method, summarized in Fig. 1, there are two terms that comprise ∆G_sub: a DFT-computed term that references atoms in their standard state (∆G₁) and the experimental free energy of solvation (∆G₂).

The overall ε-MAl₁₂ substitution reaction modeled is shown in Eq. (1), where the central Al³⁺ cation is replaced by the incoming M^m+ cation,

Equation (1) can be broken down stepwise into two pieces. The value of the ∆G₁ term is based on DFT total energies [with zero-point vibrational energy (ZVPE) and T∆S correction terms based on vibrational analysis] corresponding to the species shown in Eq. (2), where the incoming metal cation M is referenced to M_(s) and the outgoing Al is referenced to Al_(s). For the substitution reaction studied here, the expression for ∆G₂ for each metal is equal to the respective value of ΔG°_SHE, the experimentally tabulated free energy of solvation for cationic species relative to the standard hydrogen electrode (SHE). The full expression for ∆G₂ is given by Eq. (S1). For the thermodynamic process depicted in Fig. 1 for ion substitution, the value of ∆G₂ is equal to the difference in ΔG°_SHE values of the outgoing Al³⁺ cation dissolution (ΔG°_SHE of −5.03 eV) and incoming M^m+ cation incorporation shown in Eq. (3). For example, the ΔG°_SHE of Ga³⁺ is −1.65 eV; for replacing Al ³⁺ with Ga³⁺, the ΔG₂ term is equal to −3.38 eV. All ΔG values can be found in Table S5 of the supplementary material,

Combining the DFT-calculated ∆G₁ and experimentally tabulated ∆G₂ values, we obtain ∆G_sub, the energy for heteroatom substitution within the Keggin structure according to the overall reaction shown in Eq. (1). The form of the ∆G₂ term in the DFT + thermodynamics method can include correction terms from the standard state that take into account tunable experimental parameters, such as ion concentration and pH [Eq. (S1)]. In this work, we reference aqueous M^m+ and Al³⁺ cations, where the oxidation reactions and corresponding ∆G₂ values can be found in the supplementary material (Table S6). The same approach is used to study Cr³⁺ and Zn²⁺ substitution into the δ-Al₁₃ parent structure, sampling substitution of one Cr³⁺/Zn²⁺ cation to form δ-CrAl₁₂⁷⁺ or δ-ZnAl₁₂⁶⁺, where the position of the substituting cation is varied between the four distinct metal sites in the δ-isomer.

All chemicals, AlCl₃·6H₂O (Fisher Scientific), Cr(NO₃)₃.9H₂O (Alfa Aesar), NaOH (Fisher Scientific), 2,6-naphthalenedisulfonate disodium salt (97%, Sigma Aldrich), and 2,7-naphthalenedisulfonate disodium salt (95%, ACROS Organics), were used as received. Aqueous solutions in this study were prepared with ultrapure water (18.2 MΩ cm, Easypure II). In a typical synthesis, 0.33 g of Cr(NO₃)₃ · 9H₂O was added to 24 ml of a 0.25M NaOH solution and stirred until the solid was completely dissolved in the aqueous phase. A 10 ml aliquot of 0.25M AlCl₃ solution was added to a 50 ml Erlenmeyer flask and heated in a water bath until the temperature of the solution reached 65 °C. The Al³⁺ solution was then titrated with the basic chromium stock solution while maintaining the temperature at ∼65 °C. After cooling to room temperature, 10 ml of this partially hydrolyzed Cr³⁺/Al³⁺ stock solution was loaded into a 23-mL Teflon-lined Parr reaction vessel and placed in a gravimetric oven set at 85 °C to hydrothermally age for seven days. This thermally aged solution was then allowed to slowly cool to room temperature and stored in a sealed container. Different crystallizing agents including disodium 2,6-napthalenedisulfonate (2,6-NDS), disodium 2,7-napthalenedisulfonate (2,7-NDS), sodium selenate, or potassium sulfate, were added in aged and unaged solution to isolate polyaluminum species. The optimized crystallization experiments are provided below.

Na[AlO₄Al_9.6Cr_2.4(OH)₂₄(H₂O)₁₂](2,6-NDS)₄(H₂O)₂₂ (δ-Cr_nAl_13-n -1): A 4.0 ml aliquot of the thermally aged Cr³⁺/Al³⁺ solution was added to a 20 ml glass scintillation vial. Subsequently, a 2.6 ml aliquot of the 2,6-napthalenedisulfonate stock solution (0.1M) was added to the aged solution to crystallize δ-Cr_nAl_13-n-1. After approximately seven days of slow evaporation, purple columnar crystals of δ-Cr_nAl_13-n -1 were formed alongside an amorphous light green flocculant. The crystals were washed using deionized water to remove the amorphous precipitant. The calculated yields for the crystalline material were 20% based upon Al.

Na[AlO₄Al_9.5Cr_2.5(OH)₂₄(H₂O)₁₂](2,7-NDS)₄(H₂O)_18.5 (δ-Cr_nAl_13-n -2): A 4.0 ml aliquot of the thermally aged Cr³⁺/Al³⁺ solution was added to a 20 ml glass scintillation vial, followed by addition of 2.6 ml of 2,7-napthalenedisulfonate stock solution (0.1M). After approximately ten days of slow evaporation, columnar purple crystals of δ-Cr_nAl_13-n-2 were formed alongside amorphous flocculant. The crystals were washed using deionized water and the calculated yields for the crystalline material were 15% based upon Al.

Crystalline materials were separated from mother liquor and immediately coated in mineral oil to preserve crystallinity. The high-quality single crystal of each solid phase was mounted on a MiTeGen micromount and then placed on a Bruker D8 Quest single-crystal diffractometer equipped with a microfocus x-ray beam (Mo Kα; λ = 0.71 073 Å) and a CMOS detector. Diffraction experiments were conducted at 100 K using the Oxford low-temperature cryosystem with the Bruker APEX3 operation software package. Peak intensities were corrected for Lorentz polarization, background effects, and absorption effects using the APEX3 software.⁵⁴ Initial structure solution was determined by intrinsic phasing and refined on the basis of F² for all unique data using the SHELXL 5 program.⁵⁵ Select crystallographic details can be found in Table S1, and additional structural information can be found in Tables S2 and S3 in the supplementary material. The crystallographic information files (CIF) can be found on the Cambridge Structural Database upon requesting numbers 2 045 451 and 2 045 452.

The amounts of Al³⁺ and Cr³⁺ present in the synthesized crystalline material were quantified by analysis on an Agilent 7900 Inductively Coupled Plasma Mass Spectrometer (ICP-MS). For each sample, a small quantity (2 mg–5 mg) of crystalline materials from three different sets of experiments was dissolved into individual test tubes containing a 2% HNO₃ solution. A series of Al³⁺ and Cr³⁺ standard solutions were also prepared in a 2% HNO₃ matrix by dilution of 1000 ppm ICP-MS stocks purchased from Ultra Scientific and SPEX CertiPrep (lot no: L00987 and 24-190CR3M). Calibration curves constructed from standard solutions were used to determine unknown concentration of Al³⁺ and Cr³⁺ in each sample.

Additional chemical characterization for the synthesized material included infrared spectroscopy, thermogravimetric analysis, and scanning electron microscopy. Vibrational spectroscopy of synthesized samples was carried out on a Nicolet Nexus 760 FTIR spectrometer. The solid samples were prepared for FTIR analysis by mixing the sample with solid KBr and pressing into a translucent pellet. Spectra were collected from 400 cm⁻¹ to 4000 cm⁻¹. Thermogravimetric analyses on air dried crystalline materials were performed using a TGA Q500 apparatus (TA Instruments). Approximately 5 mg–10 mg of each sample was placed into an aluminum pan and heated at a ramp rate of 10 °C/min from 20 to 600 °C with high-purity N₂ used as the carrier gas. Scanning Electron Microscopy-Energy Dispersive Spectroscopy (SEM-EDS) was done on Hitachi S-3400N instrument equipped with the silicon drift detector (XFlash, Bruker AXS). The measurements were carried out at 15 kV beam energy, 2.5 min exposure time, and 10 mm of working distance. Additional chemical characterization details for δ-Cr_nAl_13-n-1 and δ-Cr_nAl_13-n-2 are provided within the supplementary material (Sec. II).

For the ε-isomer of the Keggin series, three experimental MAl₁₂ (M = Al³⁺, Ga³⁺, Ge⁴⁺) analogs have been synthesized, characterized, and compared in terms of their reactivity. A relationship has been proposed for the tetrahedral M-µ₄O and octahedral μ₄O-Al_o bonds, where a shorter and stronger tetrahedral M-µ₄O bond corresponds to longer and weaker μ₄O-Al_o bonds, and vice versa.

Comparing the structures of the experimentally obtained ε-MAl₁₂ analogs presented in Table I, there is a structural trend where the μ₄O-M bond lengths are shorter than the μ₄O-Al_o bonds (μ₄O-M < μ₄O-Al_o). A previous study showed that the μ₄O-Al_o bond lengths were influenced by anion outer-sphere adsorption, where elongation and dissociation of these bonds can be tuned by anion pKa.²⁹ By comparing the structures of the three isolable Keggins, we can draw connections between the bond lengths and the properties of the heteroatom M. The DFT-computed μ₄O-M bond lengths are more similar for ε-Al₁₃⁷⁺ and ε-GeAl₁₂⁸⁺ which have identical ionic radii of 0.39 Å compared to ε-GaAl₁₂⁷⁺, suggesting that heteroatom M ionic radius is a more important factor in controlling the tetrahedral bond lengths than charge. The μ₄O-Al_o bonds are within 0.02 Å for the two 7 + polyoxocations (ε-Al₁₃⁷⁺ and ε-GaAl₁₂⁷⁺) than compared to ε-GeAl₁₂⁸⁺, indicating that the heteroatom charge has a greater influence on the octahedral bonds.

The final row of Table I corresponds to the ratio of the μ₄O-M and μ₄O-Al_o bond lengths, or the ratio of the tetrahedral to octahedral bond lengths, denoted as R_t/o. This ratio is similar to Goldschmidt’s tolerance factor, t, which is used to compare the geometric stability and strain of crystal structures as a function of ionic radii.⁵⁶ In recent years, the tolerance factor has been applied to understanding ABX₃ (X = O, S) perovskites, useful in predicting structure–reactivity relationships such as ferroelectric properties, Gibbs free energy of formation ∆G^○, and understanding the migration energy of vacancies.^57,58 Goldschmidt’s tolerance factor has been a useful geometric guide to understanding and rationalizing perovskite geometry and can be applied in this work to predict heteroatom substitution candidates. For the three known ε-MAl₁₂, the R_t/o is <1 due to the fact that the μ₄O-M bonds are shorter than the μ₄O-Al_o. The R_t/o range for the known Keggins is between 0.87 and 0.96, with an average bond length ratio of 0.92 (±0.03) Å. The ratio can be used as a metric for comparison, where desirable candidate ε-MAl₁₂ structures have ratios near the range of the known Keggins.

Altering the heteroatom identity in the three known ε-MAl₁₂ analogs has been shown to influence the structure. By replacing the heteroatom M with metal cation species that vary in their formal charge, ionic radii,^59,60 and coordination preference, we delineate trends in the optimized geometries of the resulting ε-MAl₁₂ analogs. We sample several species spanning the s-, p-, and d-blocks of the periodic table and track changes in the μ₄O-M and μ₄O-Al_o bond lengths, highlighted in Fig. 2. Previous work has shown that these two bonds control ε-MAl₁₂ Keggin reactivity and stability.^29,61 The structural results are shown in Fig. 3, arranged by increasing μ₄O-M bonds. There is a reversal in the observed bond length trends going from left to right in Fig. 3. We systematically compare the changes in bond length trends to cation properties such as coordination preference, ionic radius, and charge. Exact values for ionic radii, bond lengths, and ratios can be found in Table S4 of the supplementary material.

We start by comparing the ε-MAl₁₂ structures for cation species that that do not prefer a tetrahedral coordination environment and are usually found in environments with coordination numbers greater than 6: Ca²⁺, Bi³⁺, Sr²⁺, and Ba²⁺, which are all found on the right-hand side of the plot in Fig. 3. These ε-MAl₁₂ Keggin structures have μ₄O-M bonds that are significantly longer than the μ₄O-Al_o bonds, with R_t/o ratios greater than 1.15. No comparison can be made for the ionic radii for these cations in four-fold coordination, but these species have comparably large atomic radii compared to the other cations in this study. The large size of these heteroatom cations requires a larger cavity, which pushes the μ₄O further away, bringing them closer to the Al_o. The rest of the cations in this study all have been shown to tetrahedrally coordinate and do not have as severe of a difference between the μ₄O-M and μ₄O-Al_o bond lengths, with R_t/o ratios below 1.15.

Next, we compare trends in ionic radius by grouping ε-MAl₁₂ species according to similar ionic radii. Found on the left-hand side of Fig. 3, the first group has a similar ionic radius of 0.4 Å and contains Mn⁴⁺, Ge⁴⁺, Co⁴⁺, and Al³⁺.^59,60 These four ε-MAl₁₂ species have similar μ₄O-M bond lengths of 1.84 (±0.01) Å. These ε-MAl₁₂ Keggin structures differ in their μ₄O-Al_o bond lengths and bond length ratios only for ε-Al₁₃⁷⁺, which has a different charge than the other species. The second group has a similar ionic radius of 0.48 Å and contains Ga³⁺ and Fe³⁺;^59,60 these Keggins have similar μ₄O-M bond lengths of 1.93 (±0.01) Å and identical μ₄O-Al_o bond lengths of 2.00 Å; this results in similar bond length ratios of 0.97 (±0.01). The overall ε-MAl₁₂ structures for Ga³⁺ and Fe³⁺ are very similar due to their identical charges and similar ionic radii. The third group has a similar ionic radius of 0.58 Å and contains Hf⁴⁺, Sn⁴⁺, Mg²⁺, Zr⁴⁺, Zn²⁺, and Li⁺.^59,60 These ε-MAl₁₂ structures have similar μ₄O-M bond lengths of 2.01 (±0.01) Å. The μ₄O-Al_o bond similarity depends on the charge of the cation; the 4+ cations have an average bond length of 2.03 Å, and the 2+ cations have 1.92 Å. Of these six candidate heteroatom species, only Hf⁴⁺, Sn⁴⁺, and Sr⁴⁺ have bond length ratios that are less than 1.00.

Finally, we compare bond length trends in terms of cation charge. There are two cases, not already discussed, of different ionic radii in which bond length similarities can be attributed to identical cation charge. The first is for Al³⁺ (0.39 Å), Ga³⁺ (0.47 Å), and Fe³⁺ (0.49 Å), which are all 3+ cations with similar μ₄O-Al_o bond lengths of 2.01 (±0.01) Å. These three ε-MAl₁₂ Keggins differ in the μ₄O-M bond lengths, ionic radii, and bond length ratio. The second is Li⁺ (0.59 Å) and Na⁺ (0.99 Å), which have a more significant difference in their ionic radii than the first case. The 1+ cations result in similar μ₄O-Al_o of 1.87 (±0.01) Å. Both of these 1+ cations have bond length ratios greater than 1.00. In all cases of comparing cations of similar charge, there is an increase in the bond length ratio that corresponds to an increase in the ionic radius.

From comparing bond length trends, we notice the following relationships: the μ₄O-M bonds depend significantly on the ionic radius, and the μ₄O-Al_o bonds are more sensitive to changes in heteroatom charge. For Keggin nanoclusters with heteroatoms of similar charge and ionic radii (i.e., Mn⁴⁺ Ge⁴⁺, and Co⁴⁺), the resulting structures are nearly identical. The central cavity that the heteroatom occupies depends on the cation ionic radius and charge. These relationships have been previously unidentified for ε-MAl₁₂ and may have important geometric implications in the other four Keggin isomers.

Reusser et al. previously predicted ε-ZnAl₁₂⁶⁺ through a different thermodynamic analysis that used M(OH)_4(aq) as the reference for the heteroatom species.³⁵ They performed a structure analysis by surveying tetrahedral Zn²⁺–O bond lengths in mineral oxides and comparing the values to the optimized bond lengths in ε-ZnAl₁₂⁶⁺. However, in this work, we find that ε-ZnAl₁₂⁶⁺ exhibits a reversal in the bond length trend (μ₄O-M > μ₄O-Al_o) compared to the three known ε-MAl₁₂ analogs and thus, by our findings, would not be predicted to form. A key difference in the two approaches is that our model includes the free energy contribution of the heteroatom going from a metal in standard state to a solvated cation. Therefore, the difference in the model predictions could stem from taking into account the role of solvated cation species.

Guided by this structural analysis, our list of candidate M cations for heteroatom substitution can be grouped as members of the s-block (Be), p-block (Ga, Ge, and Sn), and d-block (Fe, Hf, and Zr). We were able to eliminate cations from our initial set that (i) do not prefer tetrahedral coordination (Ca, Ba, etc.) and (ii) that have bond length ratios (R_t/o) greater than 1.00. To the remaining cations, we apply our thermodynamic model to calculate the energetics of heteroatom substitution.

The list of candidate heteroatom cations that run through the ∆G model includes Be²⁺, Ga³⁺, Fe³⁺, Ge⁴⁺, Hf⁴⁺, Sn⁴⁺, and Zr⁴⁺. These cations vary in their charge (2+ to 4+) and their ionic radii (0.27 Å–0.59 Å), but they have published reduction potential values and all follow the μ₄O-M < μ₄O-Al_o requirement derived from the structural analysis.

The ∆G results are arranged by increasing values of ∆G_sub, corresponding to decreasing favorability for the heteroatom to substitute, found in Fig. 4, and the values can be found in Table S3. The parent ε-Al₁₃⁷⁺ nanocluster is included in the plot as a reference point, where all of its ∆G values are zero. All of the ∆G₁ values are positive except for that of Be²⁺; ε-BeAl₁₂⁶⁺ is the only heteroatom analog that has a lower charge than the reference ε-Al₁₃⁷⁺ cluster, increasing the energetic stability relative to the higher-charged reference.

A negative ∆G₂ value reflects a simultaneously favorable Al³⁺ dissolution and M^m+ incorporation. All ∆G₂ values are negative except for Zr⁴⁺ and Hf⁴⁺, which means that it is more favorable for Al to remain within the Keggin and not form solvated Al³⁺ when being substituted by these two cations. It is interesting that Zr⁴⁺ and Hf⁴⁺ have similar ∆G values across the board; these two cations have similar ionic radii (0.59 Å and 0.58 Å, respectively) and the most similar ΔG°_SHE values (−5.8 eV and −6.2 eV, respectively). Additionally, ε-ZrAl₁₂⁸⁺ and ε-HfAl₁₂⁸⁺ have the most similar geometries, with identical bond length ratios of 0.99, corresponding to μ₄O-M and μ₄O-Al_o bond lengths of ∼2.01 Å and 2.02 Å, respectively.

We now consider the energetic findings for the experimentally obtained ε-GaAl₁₂⁷⁺ and ε-GeAl₁₂⁸⁺ analogs. For ε-GaAl₁₂⁷⁺, the ∆G₁ value is positive, which indicates an unfavorable substitution; however, the ∆G₂ value is very large and favorable such that it recovers the favorable substitution of Ga³⁺ for Al³⁺, reflected by the negative value of ∆G_sub. This is not the case for ε-GeAl₁₂⁸⁺, where the ∆G_sub suggests an unfavorable substitution. The ∆G₁ value is very large and unfavorable, and this is caused by the increase in positive charge going from ε-Al₁₃⁷⁺ to ε-GeAl₁₂⁸⁺. Some of the energy is recovered by the negative (favorable) ∆G₂ value, but not enough to drive a favorable substitution according to this model. Yet, as Ge substitution has been achieved, we go on to consider the shortcomings of the model. The synthesis for ε-GeAl₁₂⁸⁺ is complex, in which the starting material is GeO_2(s), as opposed to the aqueous Ge⁴⁺ cation modeled here. It may be that GeO₂ must dissolve to form Ge⁴⁺ that occupies the heteroatom site, which may be controlled by the reaction kinetics. The exact mechanism for ε-MAl₁₂ analog formation is not known, but the data presented here imply that direct substitution of Al³⁺ for Ge⁴⁺ within a stable Keggin is not thermodynamically favored.

Based on our energetic analysis, the most likely candidate for heteroatom substitution is Fe³⁺. The ∆G_sub for ε-FeAl₁₂⁷⁺ is +0.30 eV, shown in Fig. 4, between the two experimentally obtained Keggins, ε-GaAl₁₂⁷⁺ and ε-GeAl₁₂⁸⁺. Although the ∆G_sub is positive and would indicate substitution to be unlikely, the magnitude is smaller than that for ε-GeAl₁₂⁸⁺. This species also has a bond length ratio less than 1.00, signifying that this heteroatom analog obeys the structural criterion derived in Sec. III B.

The DFT + thermodynamics methodology (detailed in Sec. II B) may also be useful in comparing the energetics of competing pathways to cluster formation. To this end, the model was applied to the δ-Al₁₃ Keggin to assess the thermodynamic stability of replacing specific Al³⁺ atoms with Cr³⁺ and Zn²⁺, sampling the site-specific substitution that results in the experimental crystal structures of the δ-Zn(CrAl)₁₂ and δ-ZnCr₁₂ Keggin clusters reported by Wang et. al.^36,37 Within these heteroatom polycations, the Cr³⁺ cations have only been experimentally identified within the octahedral sites of the δ-isomer (sites 2, 3, and 4, Fig. 5) and have not been found to occupy the tetrahedral heteroatom position (site 1, Fig. 5). The goal of is this section is to apply the DFT + thermodynamics cycle to Cr³⁺ substitution, where this cation has not been observed to occupy the heteroatom position. We will begin by benchmarking the known clusters presented by Wang et al. before moving to experimental studies on Al³⁺/Cr³⁺.

We begin by comparing the energetics of substitution (Table II) for one cation into the δ-Al₁₃ parent structure to form either δ-ZnAl₁₂⁶⁺ or δ-CrAl₁₂⁷⁺, where the position of the Zn²⁺/Cr³⁺ cation is varied between one of the four unique metal cation sites identified in Fig. 5. It is important to note that replacing Al³⁺ with Zn²⁺ causes the overall nanocluster charge to become less positive and closer to achieving charge-neutrality, which increases stability. This increase in stability is reflected by the large, negative ΔG_sub values of the δ-ZnAl₁₂⁶⁺ species, where the most favorable Zn²⁺ substitution occurs for site 1, the heteroatom site. While substitution at the three octahedral sites results in favorable ΔG_sub values, the heteroatom position is preferred by at least 0.5 eV.

Alternatively, substitution of Cr³⁺ is the unfavorable at the heteroatom site 1 with a ΔG_sub value of 1.29 eV. Chromium in the 3+ oxidation state is typically observed in octahedral coordination; thus, forcing the cation into the heteroatom environment requires an energy penalty in the ΔG₁ term. These computed ΔG values agree with experimental findings, supporting the fact that Cr³⁺ has not been observed in the tetrahedral position in any reported Keggin structure to date. The ΔG_sub values for three octahedral sites (2–4) are nearly degenerate of one another (differing by at most 0.13 eV) and are close to 0 eV, indicating that substitution is equivalent at the three octahedral positions. This suggests that despite their different connectivity, there is only a small energetic difference in substitution at these three positions. These values are also much lower in energy than that of a tetrahedral substitution, suggesting that Cr³⁺ would preferentially occupy the outer octahedrally coordinated positions within a δ-CrAl₁₂ species.

To assess the energetics of substituting both Zn²⁺ and Cr³⁺ into the δ-Al₁₃ parent structure to form δ-ZnCrAl₁₁⁶⁺, three different substitution pathways are identified and shown in Fig. 6. All pathways reference δ-Al₁₃ as the starting structure to which the cations are added, using the experimental structure of Wang et al. as a reference point.³⁷ Since substitution with Zn²⁺ was most energetically favorable and is experimentally-observed in the tetrahedral position, in examining different pathways, Zn²⁺ cation substitution always occurs at the heteroatom site 1, as identified in Fig. 5. Because there were minimal energetic differences in Cr³⁺ octahedral substitution at the three sites shown in Table II, Cr³⁺ substitution is reported for site 3, as depicted in Fig. 6. Substitution energetics for the other octahedral sites can be found in the supplementary material, Table S9.

We compare three different substitution pathways, as depicted in Fig. 6. Pathway 1 begins with δ-ZnAl₁₂⁶⁺ to which Cr³⁺ substitutes for an Al³⁺ cation, forming δ-ZnCrAl₁₁⁶⁺. Pathway 2 starts with δ-CrAl₁₂⁷⁺ as to which Zn²⁺ is added, forming δ-ZnCrAl₁₁⁶⁺. Pathway 3 entails the simultaneous substitution of both Zn²⁺ and Cr³⁺ into δ-Al₁₃⁷⁺ to form δ-ZnCrAl₁₁⁶⁺.

The computed ΔG_sub values in Table III are all favorable but show a range greater than 4 eV, implying that the pathway is important for these substitution reactions. Pathway 1 entails the addition of Cr³⁺ to δ-ZnAl₁₂⁶⁺ and is the least favorable formation pathway. Considering pathway 2, where Zn²⁺ is added to δ-CrAl₁₂⁷⁺, this pathway is the most energetically favorable formation route of the three studied here. This information implies that a Cr-substituted Keggin is more likely to accept Zn²⁺ substitution compared to a Zn-substituted Keggin accepting Cr³⁺ substitution. Pathway 3 combines Zn and Cr substitution and produces a ΔG_sub value of −2.46 eV. If one were to combine individual terms from Table II for Zn heteroatom substitution (−5.20 eV) and Cr octahedral site 3 substitution (−0.04 eV), the energy would be −5.24 eV, which is nearly two times the value obtained using pathway 3. It is important to note the starting nanocluster charge difference in the presented pathways; pathway 1 begins with the nanocluster charge already reduced to 6+, while pathways 2 and 3 start with a 7 + nanocluster that is reduced by the addition of Zn. Pathways that involve a charge change have larger, more negative ΔG_sub values. Based on the energetic information in Table III and the fact that no heteroatom-substituted δ-ZnAl₁₂ analog has been identified, Cr-substituted Keggin nanoclusters may be more susceptible to substitution with Zn.

In order to understand the energetics of substitution, we investigate the reported crystal structures of the δ-Al₁₃, δ-Zn(CrAl)₁₂, and δ-ZnCr₁₂ Keggins and compare their bond lengths to the theoretical, optimized geometry of the heteroatom-substituted δ-ZnAl₁₂ nanocluster. The bond lengths highlighted in Fig. 2 are compared δ-Keggin nanoclusters and are reported in Table IV.

As evidenced in the previous work, the geometry of the substituted-Keggin is influenced by cation identity. Tetrahedrally-coordinated Zn²⁺ has an ionic radius of 0.60 Å, which is similar to the octahedrally-coordinated Cr³⁺ radius of 0.61 Å.^59,60 Both cations in different coordination environments are larger than Al³⁺, which is 0.39 Å when tetrahedrally-coordinated and 0.535 Å in an octahedral environment. The difference in Al³⁺ ionic radii causes the μ₄O-M bonds to be shorter and the μ₄O-Al_o bonds to be longer, giving a bond length ratio R_t/o of 0.90. In the experimental structure of δ-Zn(CrAl)₁₂, the bond length ratio is greater than 1.00, with μ₄O-Zn²⁺ bonds that are longer than the octahedral bonds, which contradicts the ratio trend observed for the known ε-MAl₁₂ Keggins. When all octahedral sites contain Cr³⁺ in δ-ZnCr₁₂, the bonds are equivalent, with a ratio of 1.00, which is due to the cations having similar ionic radii.

Based upon the computational results, we further evaluated the Cr³⁺/Al³⁺ system experimentally through the crystallization of heteroatom species. To promote the formation of the δ-isomer, the solutions were aged hydrothermally for seven days. This same procedure was used for the Ga³⁺/Al³⁺ system to promote the formation of δ-(GaAl₁₂)⁷⁺ and the related (Ga_2.5Al_28.5)¹⁹⁺ phase.³⁴ Disulfonate ligands were utilized after the extended hydrolysis reaction to induce crystallization of solid phases. This allows the structural characterization of the crystalline solid to confirm the structural features and cation substitution. Using these methods, we were able to crystalize two δ-isomers with Cr³⁺ substitution: δ-Cr_nAl_13-n-1 and δ-Cr_nAl_13-n-2 [Figs. 7(a) and 8(a)]. The four different octahedral positions that were evaluated using computational methods are also labeled with the same numbering scheme used for the DFT calculations [Figs. 7(b) and 8(b)].

The main structural feature of δ-Cr_nAl_13-n-1 and δ-Cr_nAl_13-n-2 is a δ-Keggin unit, which consists of the central tetrahedral positions surrounded by the 12 octahedrally coordinated metal cations. In the δ-isomer, one of the four trimeric units has a different orientation than the others due to a 60° rotation around the center, highlighted in Figs. 7(b) and 8(b) in pink. The M-O bond length of tetrahedral site (labeled 1) for both forms of the cluster ranged from 1.777(10) Å to 1.835(4) Å, which was similar to the distance observed for pure δ-Al₁₃ [1.795(7) to 1.830(6)]. Bond distances for the octahedrally coordinated metal cation could be found in two regions, where sites 1 and 2 were observed [1.879(11) Å to 1.999(9) Å], whereas site 4 was slightly shorter [1.811(11) Å to 1.954(10) Å].

To initially explore substitution for the δ-Cr_nAl_13-n-1 and δ-Cr_nAl_13-n-2 Keggin clusters, we model Cr³⁺ substitution within the octahedral Al³⁺ sites and allowed the occupancy to free refine on subsequent least-squares cycles. The displacement parameters were carefully evaluated to ensure the values remained within a reasonable range (Table VI). Tetrahedral site 1 and octahedral site 4 were found fully occupied with Al³⁺ for both structures whereas Cr³⁺ occupancies differed slightly. The Cr³⁺ occupancy for δ-Cr_nAl_13-n-1 in unique sites 2 and 3 have an average value 40.1 ± 2.7 and 37.3 ± 1.9, respectively (Table VI). In the case of δ-Cr_nAl_13-n-2, the Cr³⁺ occupancy was observed with an average value 46.2 ± 6.2 and 37.5 ± 6.1 for sites 2 and 3, respectively.

Elemental analyses were carried out to confirm the occupancy and presence of heteroatom substitution in the isolated structures. Scanning electron microscopy with energy dispersion x-ray spectroscopy was initially used to qualitatively confirm the presence of Cr³⁺ (Figs. S5 and S6) and found similar mass percentages of O, Na, Al, S, and Cr at different 3–4 positions throughout the solid crystalline materials. Quantification of the Cr³⁺ content was performed using ICP-MS on acid digested crystalline materials (Table V). The experimental Al/Cr ratio (δ-Cr_nAl_13-n-1 = 4.03 ± 0.02; δ-Cr_nAl_13-n-2 = 3.77± 0.01) was similar to the values obtained by single-crystal x-ray diffraction (Cr_nAl_13-n-1 = 4.51; δ-Cr_nAl_13-n-2 = 4.18). This leads to an overall formula for δ-Cr_nAl_13-n-1 and δ-Cr_nAl_13-n-2 of Na[AlO₄Al_9.6Cr_2.4(OH)₂₄(H₂O)₁₂](2,6-NDS)₄(H₂O)₂₂ and Na[AlO₄Al_9.5Cr_2.5(OH)₂₄(H₂O)₁₂](2,7-NDS)₄(H₂O)_18.5 based upon combined elemental analysis and structural refinement.

Predicted energies from the DFT + thermodynamics approach agreed with the overall substitution within the newly isolated δ-Cr_nAl_13-n-1 and δ-Cr_nAl_13-n-2 phases. Computational results indicated that Cr³⁺ substitution is not favored for the central tetrahedrally-coordinated position (site 1), and this matches the data obtained by single-crystal x-ray diffraction. Furthermore, differences between the energy values were minimal for the theoretical octahedral substitutions, suggesting that all three possible sites could incorporate Cr³⁺. For the structural study, we observed higher electron density at sites 2 and 3, and partial occupancy refinement suggested that these sites contained ∼40% Cr. It is possible that increasing the amount of Cr³⁺ may lead to additional substitution within those positions. We increased the Cr³⁺ concentration within the reaction but were unable to isolate crystalline material from these conditions.

One notable structural feature that also must be mentioned is the position of the Na⁺ cation that is associated with the Keggin cluster. Both structures exhibit a similar bonding environment for Na⁺ with bonds to three bridging hydroxyl groups on the surface of the Keggin molecule and additional linkages to three ligated water molecules with bond distances ranging from 2.316(3)-2.437(9). However, the Na⁺ cation is coordinated to the rotated trimeric group for δ-Cr_nAl_13-n-1, but it is associated with an unrotated trimeric group for the δ-Cr_nAl_13-n-2 version. Previous structures, where δ-MAl₁₂ (M = Al³⁺ or Ga³⁺), indicate that the Na⁺ cation is found only at the rotated trimer, and it has been postulated that its presence stabilizes the δ-isomer from further condensation to Al₃₀ or Ga_2.5Al_28.5 forms.³⁴ Changes in the position of the Na⁺ cation may be indicative of differences in surface reactivity due to Cr³⁺ substitution, but additional investigations are necessary to confirm this hypothesis.

Through a combined structural and energetic analysis using a DFT + thermodynamics methodology, new insights have been obtained on the ε-MAl₁₂ heteroatom analogs. Comparing known properties and trends, we have identified desirable characteristics for heteroatom cation substitution. We have also identified relationships between the heteroatom cation charge and μ₄O-Al_o bond lengths as well as between the cation ionic radius and μ₄O-M bond lengths. For the known ε-isomer analogs, the μ₄O-M:μ₄O-Al_o bond length ratios are less than 1.00, and of the sampled ε-MAl₁₂ species studied here, the best candidate for experimental synthesis is ε-FeAl₁₂⁷⁺. Fe³⁺ has the same charge as Al³⁺ and Ga³⁺, and it also has a similar ionic radius to Ga³⁺. Through the DFT and thermodynamics ∆G methodology, we predict that Ga³⁺ will substitute for Al³⁺ as the tetrahedral center, which agrees with the experimental findings that isolate this cluster species. In agreement with the predictions of Ruesser et al., we propose that of the novel ε-MAl₁₂ species studied here, ε-FeAl₁₂⁷⁺ is the most likely structure to be obtained. At present, the DFT and thermodynamic model used here describes simultaneous cation dissolution (Al³⁺) and incorporation (M^m+) for cation exchange.

Extending the methodology to the δ-Keggin species enabled us to provide a mechanistic understanding of the previously reported compounds by Wang et al. and predicted the Cr³⁺ substitution within the Al³⁺/Cr³⁺ system. This result confirms experimental observations that Cr³⁺ will only occupy the octahedral sites in the Keggin structure; Cr³⁺ does not display tetrahedral coordination, and forcing the cation into this environment comes at an energetic cost greater than 1 eV. This confirmation of the model suggests that it could be used to predict both the identity and location of favorable heteroatom substitutions, which could go on to further inform synthetic strategies to isolate new clusters. Two Cr-substituted δ-Keggins were obtained experimentally with substitution occurring within the octahedral positions, confirming the results predicted from DFT. Thus, the DFT + thermodynamics model is capable of differentiating between octahedral and tetrahedral substitution preferences within both the ε- and δ-Keggin isomers.

As previously mentioned, the form of the ∆G₂ term in the DFT + thermodynamics method can include terms to account tunable experimental parameters, including ion concentration, temperature, and pH. Insights into the elementary steps of cluster formation would further direct the model and allow for more representative processes to be used for ∆G₁. We are currently limited by a lack of understanding on the mechanistic insight for heteroatom analog formation. Present findings suggest that direct cation substitution may be the mechanism to generate the experimentally obtained ε-GaAl₁₂⁷⁺ Keggin nanocluster. Some species do not exist as free M^m+ cations in the pH range 4–6 in which Keggin clusters are synthesized; for example, Ge⁴⁺ is added to the synthesis as GeO_2(s) and Zr may exist as ZrO²⁺ or ZrO_2(s) depending on the concentration.^62,63 Future work will entail modifying the thermodynamic cycle to reference alternative starting species and the incorporation of additional synthesis species such as counter ions in order to assess how their interaction with the Keggin heteroatom analogs influences stability. The combined DFT + thermodynamics methodology detailed here may be applied to study additional systems and processes such as substitution into other polyoxoions (both cationic and anionic), ion exchange at mineral interfaces, and adsorption/desorption reactions on both nanoclusters and surfaces.

See the supplementary material for additional information on structural and chemical characterization of experimental nanoclusters, further details on the DFT + thermodynamics approach and results, and optimized XYZ coordinates of MAl₁₂ models.

This paper is dedicated to Professor Alexandra Navrotsky, a pioneering figure in the field of physical chemistry. Professor Navrotsky is a self-described maverick, though we the authors consider her a true “rock” star. Alex is a supportive and luminary mentor to countless female scientists and engineers. She has imparted on us her creative use of thermodynamics to understand nanocluster formation, which inspired this collaborative project.

J.L.B. and M.S. contributed equally to this work.

We would like to acknowledge Mackenzie Cole for her assistance in running preliminary nanocluster calculations and formatting this manuscript. This work was supported by the National Science Foundation Center (Grant No. NSF-CAREER 1254127). This research was supported, in part, through computational resources provided by the High Performance Computing Team at The University of Iowa, Iowa City, Iowa. This work used the Extreme Science and Engineering Discovery Environment (XSEDE),⁶⁴ which is supported by the National Science Foundation (Grant No. ACI-1053575). J.L.B. would like to acknowledge the College of Liberal Arts and Sciences at the University of Iowa for financial support from the Dissertation Writing Fellowship.

The data that support the findings of this study are available from the corresponding author upon reasonable request.