Effects of interstitial water and alkali cations on the expansion, intercalation potential, and orbital coupling of nickel hexacyanoferrate from first principles

Prussian blue analog (PBA) materials have found extensive use as Faradaic electrodes for both ion separations and energy storage during the past decade.^1,2 PBAs have a general chemical formula of $AxMy[M′(CN)6]z⋅nH2O$⁠, in which transition metal centers $M$ and $M′$ are linked through cyanide ligands (⁠$CN$⁠). The micro-porous $My[M′(CN)6]z$ framework intercalates and deintercalates monovalent cations while reducing and oxidizing redox-active metal centers on their lattices, respectively. Cations of various kinds can be intercalated into PBAs and inevitably coexist with electrolyte solvent molecules in interstitial space.^3,4 Depending on the composition of intercalated cations and solvent molecules, interstitial species can strengthen/undermine lattice stability,^5–7 improve/reduce interstitial ion mobility,^8,9 and suppress/increase activation energy of intercalation reactions.^10–12 In electrochemical separations using redox-active ion storage materials, the ability to selectivity separate dissolved ionic species requires sufficient gradation in potential. That is, a certain cation is intercalated preferentially over other cations that intercalate at lower potential. The observed correspondence between intercalation potential into PBAs with the hydrated ionic radius of alkali cations in bulk water (NH$4+≈$Cs$+>$Rb$+>$K$+>$Na$+>$Li$+$⁠) implicates interstitial water in the aqueous intercalation of such cations.^13,14 Comparison of the nominal size of PBA interstitial sites (⁠$5$ Å) with the Stokes diameter of alkali cations in bulk water (2.4–4.7 Å from Cs$+$ to Li$+$⁠) further suggests that a steric mechanism is responsible for electrochemical activity¹⁵ and the trend of potential¹³ that produces their selectivity bias toward large hydrated cations, stimulating the use of nickel hexacyanoferrate [$NiFe(CN)6$] in Cs$+$ removal from nuclear waste^16,17 and in the removal of Na$+$ and K$+$ in water desalination.^18–25 However, differences in the measured electrochemical kinetics among cations of varied sizes suggest hydration sheath rotation due to lattice interactions.¹³ Further, our electronic/atomic structure calculations have revealed significant polyatomic interactions between interstitial species in the $NiFe(CN)6$ lattice, including dative bonding between either Lewis acids (interstitial cations or transition metals in the lattice) and Lewis bases (oxygen atoms in water molecules or cyanide ligands in the lattice) with hydrogen bonding between water molecules occurring in concert.²⁶ Links have also been made between the limited expansion of PBA lattices and facile electrochemical kinetics of proton insertion due to extended hydrogen-bonding networks in $CuFe(CN)6$⁠.²⁷ Therefore, a comprehensive study of interstitial species in PBA lattices is essential for the rational design of future redox-active host materials that exhibit high-rate capability, long cycle life, and selective cation intercalation.

Previous studies have established relationships between the electrochemical and mechanical properties of PBA and the composition of interstitial species by varying the intercalation degree of cation, the type of intercalated cation, and lattice hydration degree. The $NaxMn[Fe(CN)6]$ lattice, for instance, transforms from a cubic to a monoclinic phase as intercalation degree of Na-ion $x$ increases from 1.40 to 1.72, and higher intercalation degree causes irreversible phase transition.²⁸ A recent experiment also used a two-phase synthesis method to eliminate interstitial water in high-crystallinity $NaxFe[Fe(CN)6]$ lattice.²⁹ The resulting PBA electrode delivered superior discharge specific capacity even at a current density of 20C. Decreasing lattice hydration degree also increased the intercalation degree of potassium ions in $AxNi[Fe(CN)6]$ electrode as a consequence of improved lattice stability.³⁰ When PBA lattices are under external pressure, increasing interstitial water content modifies low-frequency phonon modes and changes the sequence of lattice phase transitions as pressure increases.^4,31,32 In electrocatalysis applications, the $Cu[Fe(CN)6]$ electrode remains highly sensitive toward hydrogen peroxide at micromolar concentration levels when the framework lattice is intercalated with $Cu2+$⁠.³³ While these experiments show the effects of interstitial species on PBA properties at phenomenological levels, the mechanisms of interstitial species interacting with the framework lattice remain to be elucidated. Therefore, there is a clear need to investigate the impacts of interstitial species on the properties of PBA lattices at the atomic level using theoretical methods.

Density functional theory (DFT) calculation has become a popular method for studying atomic interactions relevant for predicting the properties of PBA electrodes. The DFT-calculated energy difference between perfect and defected lattices allows quantitative study of metal–ligand bond strength, based on which a recent study showed that the negative thermal expansion of PBA lattices is induced by transverse vibrations of cyanide ligands.³⁴ Also, lattice energy varies along the trajectory of cation transport in interstitial space. By studying possible diffusion trajectories of intercalated cations, another study showed that nanocages lower the energy barrier of cation diffusion in rationally designed NiFe-PBA lattices.³⁵ As an orbital-resolved method, DFT calculations are also used to study bonding interactions at interfaces between PBA electrode and solid electrolyte, where orbital coupling functions as an electric “highway” to enhance electron transport.³⁶ 

Among the routinely measured electrochemical properties of intercalation materials, first-principles calculations have shown to be critical for predicting the potentials at which cation can be inserted and removed. While high potentials increase energy density,³⁷ if the potential is too high, charge/discharge of PBAs can lead to oxygen reduction and eventually cause the explosion of sealed battery cells.³⁸ Hence, knowledge of intercalation potential is essential for determining whether PBAs can be used efficiently as electrodes in electrochemical systems. Recent studies evaluate these thermodynamic potentials by calculating Gibbs free-energy change associated with intercalation reactions,^37,39,40 where fully dehydrated framework lattices were considered. However, such calculations are not able to capture the effects of interstitial water on magnetic properties of metal centers, resulting in a large discrepancy between measured and calculated intercalation potential.⁴¹ By using a hybrid exchange-correlation functional in DFT calculations, much better agreement with experiments is achieved at the expense of computational cost.³⁷ Alternatively, another study used grand potentials of a nanocrystalline intercalation material to calculate $Mg2+$ intercalation voltage as an indicator of the host co-intercalating cation with water molecules.⁴² By choosing certain values for the chemical potential of $H2O$ molecules, the calculation managed to reproduce measured intercalation potentials with humble computational cost and found that water co-intercalates with $Mg2+$ depending on the activity of $H2O$ molecules. Similar methods have been used to study hydration phase diagrams of mildly hydrated PBAs,³⁷ but extending the grand potential analysis to highly hydrated PBAs is needed to understand the roles of water molecules during the intercalation of hydrophilic cations that are relevant to aqueous electrochemical separations and energy storage.

To conduct an ab initio study of hydrated PBA lattices requires the determination of the ground-state arrangements of interstitial species. To tackle this, a straightforward strategy is to leverage lattice space groups to find the highly symmetric arrangement of interstitial species as a starting point for geometry optimization.³⁷ By setting strict convergence criteria, such symmetry-adapted strategies have been shown useful for elucidating $H2O$-assisted structural evolution of the $NaxFeMn(CN)6$ lattice.⁴³ However, such DFT-only methods are not economical for finding ground states of highly hydrated PBAs, as the non-covalent interactions in large water clusters can only be captured by expensive density functionals. As an alternative strategy, machine learning (ML) models have been adopted in recent studies to develop force fields for optimizing hydrated systems, such as molecular crystals,⁴⁴ cement hydrates,⁴⁵ and liquid-water/metal–oxide interface.⁴⁶ Given any crystalline structure, these ML-only methods promise fast and accurate predictions of atomic forces when ML models are trained using carefully designed geometric features of atomic clusters. However, the key bottleneck of developing ML-based force field is in validating and extending them to novel hydrated phases.⁴⁷ While the transferability of classical force field can be improved between the small hydrated system by using machine-learned parameters,⁴⁸ reliable force fields for PBAs do not exist to our best knowledge. Hence, it is necessary to develop structural optimization methods for hydrated PBA lattices that produce the accuracy of DFT-level calculations while retaining the efficiency of ML-based methods.

Another challenge of studying hydrated PBAs lies in the difficulty of analyzing many-body interactions induced by the presence of confined interstitial $H2O$ molecules. Two methods for disentangling atomic interactions between water molecules and their surrounding atoms have been used to date. The first involves many-body expansion(MBE) of interaction energy with regard to reference energies obtained from ab initio calculations.^49,50 MBE has been applied to disentangle ion–water^51,52 and gas–water interactions.⁵³ The resulting interaction energies were shown to follow correct energy ordering for water clusters of varied sizes.⁵⁴ The second identifies dominant interactions in atomic clusters by examining correlations between the accuracy of ML models and many-body features. Our previous study examined the accuracy variation of XGBoost ML models using different many-body features of hydrated PBA lattices that are included in training datasets.²⁶ The variation was linked to the interaction strength between distinct polyatomic clusters. The results showed that the nature of water–lattice interactions depends on ionic radius and the hydrophilicity of cations. Both categories of methods use the information of atomic configurations to analyze molecular interactions, so they can be applied to investigate the effects of interstitial species on the structural evolution of PBA lattices. However, complementary methods are still needed to understand the effects of interstitial species on electrochemical and electronic properties of hydrated PBAs.

In response to these challenges, we investigate the expansion, intercalation potential, and orbital coupling of the nickel hexacyanoferrate PBA lattice $AxNiFe(CN)6⋅nH2O$ containing interstitial water and intercalated alkali cations (A). We first combine an ML model, simulated annealing (SA) simulations, and ab initio relaxation calculations (ML+SA+DFT method) to find ground states of $AxNiFe(CN)6⋅nH2O$ lattices. The ordering of interstitial species in ground states is thereby related to the stability of such lattices through analysis of lattice expansion as a function of hydration. Based on the ground-state energies, we then conduct grand potential analysis to determine equilibrium hydration degrees and to predict the potentials of $NiFe(CN)6$ intercalating various cations. Finally, we study orbital coupling on the $NiFe(CN)6$ framework by taking into account non-collinear magnetization in electronic structure calculations. Our results show that (1) the ordering of intercalated cations can be easily frustrated by interstitial water molecules, (2) water-aware grand potential analysis yields accurate predictions of equilibrium intercalation potentials, and (3) the insulator–semiconductor transition of the PBA lattice is caused by symmetry breaking in d-p orbital couplings.

Different methods are used in the current study to find the ground states of hydrated lattices, to predict intercalation potentials, and to examine orbital couplings. In response to the challenge of finding ground states of hydrated $NiFe(CN)6$ lattices, we first introduce a ML+SA+DFT method in Sec. II A where its meta-algorithm is elaborated and its efficiency of relaxing highly hydrated lattices is compared to a DFT-only method. Section II B includes the details for calculating and analyzing electronic structures of the ground states.

The unit-cell configuration that we adopt in all the calculations has a formula of $AxNiFe(CN)6⋅nH2O$⁠, where $A$ is an alkali cation: $Na+$⁠, $K+$⁠, or $Cs+$⁠. Based on their ionic radius and hydrophilicity, these cations are representatives of the light, intermediate, and heavy ion classes, respectively.²⁶ When electrochemically oxidized the $NiFe(CN)6$ lattice is half intercalated with cations (⁠$x=1$⁠), and when reduced it is fully intercalated (⁠$x=2$⁠). The lattice hydration degree varies from anhydrous (⁠$n=0$⁠) to highly hydrated (⁠$n=4$⁠) levels. To make our later discussion concise, we use nomenclature of $Aδn$²⁶ to denote $NiFe(CN)6$ lattices based on the type of intercalated cation $A$⁠, intercalation degree $δ$⁠, and lattice hydration level n. For instance, Naf4 and Csh2 represent lattices of $Na2NiFe(CN)6⋅4H2O$ and $CsNiFe(CN)6⋅2H2O$⁠, respectively.

We search for the ground states of hydrated lattices using the ML+SA+DFT method. The method uses ML-empowered simulated annealing(SA) simulations to conduct coarse search by annealing randomly sampled atomic configurations ${ri}0$ to the configurations near the ground states ${ri}SA$⁠. ${ri}0$ is obtained by placing cations and water molecules randomly at Voronoi points in the $NiFe(CN)6$ framework. The orientation of each water molecule is sampled using unitary quaternions, as described in our previous study.²⁶ The SA simulation uses a XGBoost ML model as the lattice energy predictor, and the detailed simulation workflow is elaborated in Fig. S1 in the supplementary material. Given that the atomic radius and the hydrophilicity vary drastically from the Na-ion to the Cs-ion, we trained distinct XGBoost model for each $Aδn$ lattice using the cation-specific features²⁶ that improve the accuracy of XGBoost models in predicting lattice energy.

After finishing SA simulations, annealed configurations ${ri}SA$ were used as the inputs for high-fidelity relaxation calculations. We used Quantum ESPRESSO⁵⁵ both to relax lattice strains and to minimize atomic forces (i.e., “variable-cell” calculations). Each calculation terminated when the following convergence criteria were satisfied: (1) the lattice pressure changed by less than 0.5 kbar between two consecutive self-consistent field (SCF) steps and (2) all components of atomic forces were smaller than $10−4$ Ry/bohr (⁠$2.57×10−3$eV/Å). The exchange-correlation interactions were approximated using the PBE functional,⁵⁶ and non-local van der Waals interactions are captured by the rVV10 functional⁵⁷ in all calculations. Recent studies show that the PBE+rVV10 functional combination can accurately predict bond lengths between adsorbed molecules and various adsorbents while keeping computational cost manageable.^41,58,59 We refer the interested readers to the supplementary material for descriptions of the energy cutoffs, k-points, and pseudopotentials used in these calculations.

Using this approach, we conducted an extensive search for the ground states of hydrated lattices. For each $Aδn$ lattice, we obtained at least 30 annealed configurations ${ri}SA$ from at least 5 distinct ${ri}0$ by repeating ML-based SA simulations. The annealed configurations were further relaxed by performing high-fidelity DFT relaxation calculations, after which we recognized the relaxed configuration with the lowest energy as the ground state. Based on ground-state configurations, we analyzed the effects of interstitial species on electrochemical and electronic properties of hydrated, cation-intercalated $NiFe(CN)6$ lattices.

We also checked the effectiveness of the ML+SA+DFT method by comparing it to a method that only used the DFT relaxation calculation to find ground states. Figure S2 in the supplementary material shows the results of such a comparison. In so doing, we found that the ML+SA+DFT method accelerates the lattice relaxation calculation by starting from configurations near the ground state. Compared to the DFT-only method, our method also results in relaxed configurations with lower lattice energies. While the search for ground states can be further accelerated by using general-purpose machine-learning force fields,^44,60 to our best knowledge no machine-learning force field has been tested on molecular species in the confined interstitial spaces like those of interest here. These results indicate that hybrid methods that use ML models for lattice relaxation before DFT calculations are effective at finding the ground states of hydrated lattices.

We analyzed the effects of interstitial species on electronic properties of $Aδn$ lattices by examining orbital couplings near the Fermi energy level $EFermi$ after obtaining grounds states using our ML+SA+DFT method. The dominant components of near-Fermi crystal orbitals are the d-orbitals at metal centers, which can be paramagnetic (⁠$NiII$⁠) and ferri-(ferro-)magnetic (⁠$FeII/III$⁠)⁶¹ in this study. By virtue of Pauli’s exclusion rule, the splitting of d-orbitals in crystal fields of cyanide ligands results in two spin orders: high-spin and low-spin states. While the paramagnetic ion $NiII$ remains in a high-spin state on the framework, recent experiments have shown that lattice distortion changes the spin state of Fe centers during electron transfer processes.⁶² Also, the crystal field of non-magnetic cyanide ligands screens the superexchange interactions⁶³ between nearest-neighbor metal centers, causing non-collinear magnetic ordering of paramagnetic atoms.⁶⁴ Finally, interstitial water can induce conversion between ferromagnetic and paramagnetic phases of the $NiFe(CN)6$ lattice by distorting the framework structure.⁶⁵ These results together indicate non-collinear magnetism at metal centers in hydrated, cation-intercalated $NiFe(CN)6$ frameworks and, therefore, we captured such magnetism by enabling magnetization at x, y, and z axes in the subsequent electronic structure calculations that we describe.

We conducted electronic structure calculations by setting the starting percentages of electron spins in the positive z-direction at Ni and Fe centers to be 27.8% and 31.25%, respectively. These values were obtained from the crystal analysis tool provided by Material Cloud platform.⁶⁶ The van der Waals correction to lattice energy was excluded from these calculations because the formalism of the rVV10 kernel is not compatible with non-collinear spin polarization, and it is known to have a limited effect on the electronic structure.⁵⁷ 

To relate near-Fermi orbital couplings to electron transfer in hydrated, cation-intercalated $NiFe(CN)6$ lattices, we performed SCF calculation on $2×2×2$ supercells of the ground-state configurations by adding one extra electron on each framework of interest. A compensating jellium background was used in the calculations to neutralize these supercells.⁵⁵ This allows us to visualize the highest occupied crystal orbitals (HOCOs) within supercells to study the effects of interstitial species on electron density distribution. This “extra-electron” method has found extensive use in studies of polaron transport in crystalline materials where electrons are strongly bounded to defects formed by lattice distortion.^67,68 As we shall see in Sec. III A, lattice distortion also prevails in highly hydrated $NiFe(CN)6$ lattices, and the morphology of HOCOs is affected by the arrangements of interstitial species. To keep computational cost modest, we turn off non-collinear magnetism in HOCO calculations and reduce the SCF convergence threshold to $10−7$ eV.

In Secs. III A–III C, we discuss the effects of interstitial species on the lattice swelling (Sec. III A), cation-intercalation potentials (Sec. III B), and electronic structure (Sec. III C) of hydrated, cation-intercalated $NiFe(CN)6$ lattices. In Sec. III A, we inspect the ordering of interstitial species in ground-state configurations and demonstrate a close relationship between lattice expansion and framework–cation–water interactions. Using the energies of ground states, we conduct grand potential analysis in Sec. III B to predict the potentials for intercalation of alkali cations, from which we obtain agreement with experiment for the observed intercalation selectivity sequence. We also observe a tendency to co-intercalate cations and water molecules from examination of the variation of the equilibrium hydration degree in fully (⁠$Afn$⁠) and half (⁠$Ahn$⁠) cation-intercalated lattices. Finally, Sec. III C reveals crystal orbital coupling near Fermi energy levels, based on which we rationalize the insulator/semiconductor transition for Af0 to Af4 lattices.

The unit-cell volumes of ground-state configurations were calculated after the high-fidelity relaxation calculations obtained from SA+ML+DFT optimization. Table I lists volume ratios of hydrated $Aδn$ (⁠$n>0$⁠) lattices relative to their corresponding anhydrous $Aδ0$ lattices. The volume of each $Af0$ lattice is slightly smaller than the corresponding $Ah0$ lattice because the added electron density at reduced Fe centers causes the contraction of hexacyanometallate complexes.⁶⁹ Note that to obtain the volumes of $Aδ0$ lattices at their ground states, we relaxed the ground-state configurations of $Aδ1$ lattices with interstitial water molecule removed. As a consequence, the body-centered occupancy of cations assumed in recent studies^37,70 is not preserved in the ground states of anhydrous lattices. Figure S3 in the supplementary material shows that the body-centered occupancy not only causes imaginary vibrational modes for Na-ions, but also raises the energy of the Naf0 lattice by 12.6 eV.

Table I indicates that Na-intercalated lattices (⁠$Naδn$⁠) expand less than 3% when hydration degree increases from $n=0$ to $n=3$⁠. To understand this, we visualize the atomic structure for the ground states of Nah0, Nah2, Naf0, and Naf2 lattices in Fig. 1, where Na-ions and oxygen atoms are shown as purple and red spheres, respectively. Most ground states show cubic lattice structures, where the $NiFe(CN)6$ framework is composed of $Ni1/2Fe1/2(CN)3$ subcubes. Hence, we will use Wyckoff notations for the Fm3m space group⁷¹ to refer to symmetric lattice positions to simplify the following discussion. Each $Ni1/2Fe1/2(CN)3$ subcube in the Nah0 lattice contains one Na-ion near the 24d Wyckoff positions. Na-ions in Naf0, on the other hand, are found near horizontal ligands (24d positions) or vertical ligands (24e positions). 24e and 24d occupancy alternates in the subcubes along the direction into the page in Fig. 1. While the near-24d occupancy of Na-ions in Nah0 does not disturb the perpendicular arrangement of $NiN6$ and $FeC6$ octahedra, Na-ion occupancy in Naf0 causes slight tilting of Fe–CN–Ni atomic chains along the basis vector a. From this, we conclude that the change in Na-ion ordering from Nah0 to Naf0 causes mild lattice distortion without the help of interstitial water.

Figure 1 also shows that interstitial water molecules in $Naδ2$ lattices exhibit Na–water coordination that intensifies Na–ligand interactions. Na-ions displace from 24d in Nah0 to 8c positions [body centers of $Ni1/2Fe1/2(CN)3$ subcubes] in Nah2 by bonding Na-ions with neighboring oxygen atoms. Each Na-ion in Nah2 is coordinated with two water molecules, and the proximity of water molecules to ligands deviates the atomic angle $∠NiNC$ from $180°$⁠, causing tilting of metal octahedra. Compared to the Naf0 lattice, the near-24e occupation of Na-ions disappears in the Naf2 lattice, and interstitial water pushes Na-ions toward 24d positions to form nearly co-planar Na–ligand coordination. These adjustments in Na–water arrangement alleviate lattice distortion by restoring $∠NiNC$ to $180°$⁠. Each Na-ion in Naf2 also has two nearest-neighbor $H2O$ molecules of which the dipole moments are parallel and pointing toward the Na-ion. The average distances between Na-ions and their nearest-neighbor oxygen atoms are 2.33 Å in Nah2 and 2.41 Å in Naf2. Because Na-ion ordering is affected by Na–ligand and Na–water interactions, the longer Na–O distance suggests that the Na–ligand interaction prevails over the Na–water interaction in Naf2. At $n=4$⁠, the expansion ratio of $Naδ4$ lattices is larger than $3.9%$⁠, which can be explained by the increased size of Na–water clusters in Nah4 and Naf4 lattices (see Fig. S4 in the supplementary material).

Table I also shows that the expansion of Khn lattices is negligible when $n≤3$⁠, and this insensitivity toward lattice hydration can be related to the effects of interstitial water on the ordering of K-ions. In the Kh0 lattice, K-ions occupy 8c Wyckoff positions without causing octahedral tilting due to the bracing effect of the large K-ion (see Fig. 2). Octahedral tilting is still negligible in Kh2 as the 24e occupancy of K-ions creates enough space for water molecules to occupy the 8c Wykoff positions between K-ions. Also, the average K–N distance in Kh2 is 2.82 Å, which is about 0.2 Å longer than the average Na–N distance in Nah2. Thus, the cation–water and cation–ligand interactions that are responsible for the distortion of Nahn lattices are screened by electrons at K-ions. Figure S4 in the supplementary material shows that K-ions are not able to separate water molecules at n=4, causing the formation of an extended hydrogen bonding system that induces significant tilting and expansion of $Kδ4$ lattices.

The volume of the Kfn lattice remains nearly unchanged from Kf0 to Kf1 but experiences at least 5% expansion compared to Kf0 when $n≥2$⁠. From Fig. 2 the Kf0 lattice has significant tilting of $FeC6$ and $NiN6$ octahedra because of the 24d occupation of K-ions. On the contrary, interstitial water in Kf1 restores the integrity of the $NiFe(CN)6$ framework by moving K-ions to 8c Wyckoff positions. Like the Kh4 lattice, the large expansion ratio of highly hydrated Kfn lattices is also associated with the formation of an extended hydrogen bonding system (see Fig. S4 in the supplementary material).

Compared to $Naδn$ and $Kδn$ lattices, significant lattice expansion is found even in mildly hydrated Cshn and Csfn (⁠$≥2%$ vs $Csδ0$⁠), while no distortion found in $Csδ0$ lattices. The large size of Cs-ions at 8c positions in Csf0 prevents transverse deformation of Ni-CN-Fe atomic chains,⁷³ and the bracing effect of the heavy Cs-ions in Csh0 overcomes attractive Cs-ligand interactions to prevent octahedral tilting. On the contrary, Csh4 and Csf4 expand by 16.6% and 39.2%, respectively, compared to Csh0 and Csf0, indicating large-scale breakdown of the framework structure. Figure 3 shows that the interstitial water molecules in Csh4 and Csf4 also break Ni–N bonds and cause Ni–O coordination. Such bond breakage also creates extra interstitial space for forming an extended hydrogen bonding system (i.e., black dashed bonds in Fig. 3).

We now summarize the effects of interstitial water on lattice distortion and expansion. Water molecules closely coordinate with Na-ions and lock them near 8c Wyckoff positions in Nahn lattices, while the repulsive oxygen–ligand interaction introduces mild lattice distortion. The perpendicular arrangement of metal octahedra is maintained in Nafn (⁠$n<3$⁠) lattices as water molecules fill the interstitial void between Na-ions. The presence of water molecules changes K-ion ordering in Khn lattices without causing significant lattice distortion for two reasons: (1) the bracing effect of K-ions on the framework lattice and (2) the large distance between water molecules and cyanide ligands. Due to the confined interstitial space in Kfn lattices, lattice distortion is curbed for hydration degree $n≤2$⁠. Both Na- and K-intercalated lattices expand significantly at high hydration degrees due to the formation of extended cation–water networks in the interstitial space. The large size of Cs-ions maintains the integrity of the framework in anhydrous lattices (⁠$Csδ0$⁠), but leaves limited space for the occupation of water molecules. As a result, water in confined space tends to attack Ni–N bonds in highly hydrated Cshn and Csfn lattices.

The ordering of interstitial species discussed in Sec. III A causes variations in the ground-state energies of $Aδn$ lattices, based on which we conduct grand potential analysis to predict potentials for $NiFe(CN)6$ to intercalate Na-, K-, and Cs-ions. To do this, we recalculated all ground-state energies using the PBELYP1W⁷⁴+rVV10 functional for the following reasons. While the PBE functional leads to accurate predictions of bond distances in small molecules⁷⁵ and layered materials,⁷⁶ it underestimates cation-molecule binding energies^77,78 and overstructures cation–water agglomerates.⁷⁹ The inaccuracy of the PBE functional in predicting cation–water arrangements can be reduced by using the rVV10 functional⁵⁷ in conjunction, but the PBE+rVV10 combination can be worse at reproducing measured hydration energies of crystalline materials.⁸⁰ Therefore, it is essential to choose a proper functional that captures non-covalent interactions in interstitial bonding systems to explore these effects. Accordingly, Figs. S5 and S6 in the supplementary material display the variations of grand potential vs chemical potential of water molecule for $Naδn$ and $Kδn$ lattices, where curves are obtained based on energies calculated using various exchange-correlation functionals. A comparison of these curves indicates that grand potential analysis based on the PBELYP1W+rVV10 that we subsequently use predicts the correct sequence of intercalation potential and consequently selectivity, as well as the tendency of PBA lattices to remain hydrated during the intercalation of Na- and K-ions.

We now describe the grand potential analysis used to determine the equilibrium hydration of cation-intercalated $NiFe(CN)6$⁠. To investigate the possible co-intercalation of $H2O$ with alkali cations, we first equilibrate $Aδn$ lattices that are open to the exchange between interstitial water and water in a contacting electrolyte as governed by the grand potential,⁴² 

Here, $EAδnDFT$ in Eq. (1) is the DFT-calculated energy of an $Aδn$ lattice per unit cell, and the hydration degree $n$ takes values from 0 to 4. We express the chemical potential of water $μH2O$ in equilibrium with each $Aδn$ lattice as

where $EH2ODFT$ is the DFT-calculated energy of an isolated water molecule in vacuum and $Δμ$ accounts for translational, vibrational, and rotational energies of gaseous water at finite temperature with a pressure corresponding to the electrolyte vapor pressure at the temperature of interest. Other factors, such as ion concentrations in the electrolyte, can also affect the value of $μ$ and can be incorporated into $Δμ$ using Pitzer’s method.⁸¹ However, the implementation of Pitzer model requires information of electrolyte density as a function of temperature and ionic concentration,⁸² which is not readily available for inorganic compounds containing Cs-ion. Therefore, we present the variation of grand potential within a range of $Δμ$ that reflects finite-temperature effects on the chemical potential of water. Given certain $Δμ$⁠, the average intercalation potential for the reaction

is calculated as

where $EA$ is the energy per A atom in the corresponding alkali metal obtained from DFT and “e” is elementary charge. At chemical equilibrium, the grand potentials in Eq. (4) are replaced with $ΦAhn~Ah$ and $ΦAfn~Af$⁠, where the subscripts $n~Ah$ and $n~Af$ are the equilibrium hydration degrees (EHDs) for half- and fully-intercalated $NiFe(CN)6$ lattices, respectively. Given certain $Δμ$⁠, cation type, and intercalation degree, the EHD is the hydration degree of the lowest-energy lattice. Mathematically, $n~Ah$ and $n~Af$ can be written as

While Eq. (4) gives the potential for a certain intercalation reaction versus the corresponding A/$A+$ reference electrode, we add to $⟨ϕA⟩$ the reduction potential of the corresponding alkali metal versus the standard hydrogen electrode in order to compare calculated results to experiments. For instance, the equilibrium potential to intercalate an A-ion versus SHE (standard hydrogen electrode) is calculated as

where $EA/A+SHE$ is the reduction potential of a given A-ion versus SHE. We used $ENa/Na+SHE$⁠, $EK/K+SHE$⁠, and $ECs/Cs+SHE$ values from the standard electrochemical series:⁸³ $−$2.71, $−$2.931, and $−$3.026 V, respectively.

We now identify the relavant range of $Δμ$ for pure liquid water at room temperature, in order to compare calculated $⟨ϕeq,A⟩$ to experiments. To do this, we calculate translational (⁠$Etrans$⁠), vibrational (⁠$Evib$⁠), and rotational (⁠$Erot$⁠) energies of gaseous water as components of $Δμ$⁠. Given certain temperature $T>0$⁠, the translational energy of water vapor per molecule is

where $Po$ and $Λ$ are the vapor pressure of the electrolyte at $T$ and the thermal de Broglie wavelength of a water molecule, respectively. By modeling an individual water molecule as a rigid rotor, we can use the rotational partition function⁸⁴ 

to calculate the rotational energy as $Erot=−kBTln⁡qr$⁠, where $IA$⁠, $IB$⁠, and $IC$ are the rotational moments of inertia for the $H2O$ molecule. Neglecting rotational coupling, the vibrational energy of each $H2O$ molecule is

where $Θv,i$ is a characteristic vibrational temperature. Here, we consider the vibrations of asymmetric and symmetric stretch of H–O bonds, and the scissoring bend, with their corresponding $Θv,i$ values of 5360, 5160, and 2390 K, respectively.⁸⁵ At $T=300$ K, the calculated energies are $Etrans=−0.4755$ eV, $Erot=−0.0984$ eV, and $Evib=0.3382$ eV. We note here that the rotation and vibration motions of water are coupled at moderate temperature, and the product $qv×qr$ is a simplified approximation for the partition function for the internal degrees of freedom of $H2O$ molecules.⁸⁶ In our analysis, we discuss variation of $n~$ and $⟨ϕeq,A⟩$ within the $Δμ$ range of $[Etrans+Erot,0]$ to consider the extremes of $Δμ$ if any of the three finite-temperature effects (rotation, vibration, and translation) are neglected.

Figure 4 shows the variations of $n~$ in Na-, K-, and Cs-intercalated lattices as $Δμ$ varies from $−$1 to 1 eV. Each subplot shows $n~$ vs $Δμ$ for fully and half intercalated conditions, with gaps between the two curves at given $Δμ$ indicating whether a water molecule is co-(de)intercalated with cations. At $Δμ=0$ $n~Ah=n~Af=4$ for Na- and K-intercalated lattices (see the red dashed lines in Fig. 4), a jump from $n~=2$ to $n~=4$ (⁠$2→4$ jump) is found on both $n~Ah$ and $n~Af$ curves in Figs. 4(a) and 4(b), but their positions relative to $Δμ=0$ are different. At room temperature, liquid water deviates from its pure state when mixed with ionic solutes, and the chemical potential of water is smaller in salt solutions. Therefore, a large range of $Δμ$ between the transition from an EHD of 2–4 implies the tendency of $NiFe(CN)6$ to remain highly hydrated against thermal perturbations and salt concentration variations during electrochemical circling. Comparing Figs. 4(a) and 4(b), we deduce that the $Naδn$ lattice is more likely to maintain an EHD of $n~Ah=n~Af=4$ upon Na-ion intercalation than is the $Kδn$ lattice upon K-ion intercalation. In fact, $n~Ah=n~Af=2$ in K-intercalated lattices within a large portion of the range $[Etrans+Erot,0]$⁠, suggesting that K-intercalated lattices maintain $n~=2$ from Khn to Kfn. This observation is in agreement with the fact that water is not transferred with K-ions during intercalation because of the channel space limitation.⁸⁷ As discussed in Sec. III A, the $NiFe(CN)6$ framework loses its chemical identity in highly hydrated $Csδn$ lattices, as water molecules break Ni–N bonds in the limited interstitial space. Therefore, we exclude $ΦCsδ3$ and $ΦCsδ4$ in the calculations of $n~$ and $⟨ϕeq,Cs⟩$⁠. Figure 4(c) shows that the intercalation of Cs-ions is accompanied by de-intercalation of water when $Δμ>Etrans+Erot$⁠. The EHD drops from $n~Ah=4$ in Csh4 to $n~Ah=0$ in Csf0 at $Δμ=0$⁠, and it changes from $n~Ah=1$ in Csh1 to $n~Ah=0$ in Csf0 at $Δμ=Etrans+Evib+Erot$⁠.

Figure 4(d) shows the corresponding variations of $⟨ϕeq,A⟩$ for Na-(blue), K-(orange), and Cs-ions (green) as $Δμ$ varies from $−$1 to 1 eV. For a given cation type, the higher its $⟨ϕeq,A⟩$ is, the more preferably $NiFe(CN)6$ intercalates it over other cations in electrolyte. For $Δμ$ $∈[Etrans+Erot,Etrans+Evib+Erot]$⁠, our analysis predicts $⟨ϕeq,Na⟩<⟨ϕeq,K⟩<⟨ϕeq,Cs⟩$⁠, in agreement with measured sequence of intercalation selectivity.¹³ At $Δμ=Etrans+Erot$⁠, our calculation predicts that $⟨ϕeq,Cs⟩=1.16$ V vs SHE, $⟨ϕeq,K⟩=0.84$ V vs SHE, and $⟨ϕeq,Na⟩=0.79$ V vs SHE. These values deviate from experiment as large as 0.19 V and as small as 0.14 V.¹³ Compared to grand potential analysis that only considers anhydrous lattices,³⁹ our method shows much improved accuracy in predicting intercalation potential and sheds insights into the roles of water in the process of cation intercalation.

We now investigate the effects of interstitial species on the electronic structure of hydrated lattices by performing non-collinear self-consistent field (SCF) calculations based on unit-cell ground states of $Aδn$ lattices obtained by our SA-ML-DFT structural optimization procedure. We focus on the couplings between the 3d-orbitals at metal centers and the valence orbitals of interstitial species to understand the electronic structure of Ahn and Afn lattices separately. By examining the differences in orbital couplings near the Fermi energy level between Ahn and Afn lattices, we elucidate the effects of interstitial cation and water on electron transfer through the framework lattice.

Figure 5 shows the projected density of electronic state (PDOS) distributions for valence orbitals at Fe, C, N, and cations in Nah0, Kh0, and Csh0 lattices. Non-zero PDOS distributions near the Fermi energy level $EFermi$ indicate that Ahn lattices are electronically conductive regardless of the cation type. The Fermi energy levels in Fig. 5 correspond to spin-down states, as Fermi-energy electron at low-spin Fe(III) center occupies a spin-down 3d orbital by the virtue of Pauli’s exclusion rule. Here, we identify orbital coupling between two distinct atoms near $EFermi$ when the PDOS distributions of their valence orbitals have overlap and(or) have peaks located at the same energy level.

Most Fe(3d) orbitals in Fig. 5 are coupled with N(2p) orbitals near $EFermi$⁠, and their PDOS peaks are split into two groups due to the crystal field induced by cyanide ligands. On the contrary, no coupling is found between Fe(3d) and C(2p) orbitals near the Fermi energy level, and thus Fe–C bonds are not involved in electron transfer on the framework lattice. This finding is unexpected because each Fe(III) atom is directly coordinated with six carbon atoms and because the formation of Fe(3d)–N(2p) coupling must overcome energy barriers posed by Fe–C bonds. Non-zero PDOS distribution near $EFermi$ is also found on Na(3s), K(4p), and Cs(5p) orbitals. Thus, we deduce that the highest occupied crystal orbital (HOCO) in Ah0 lattices is a result of Fe–N–cation interaction. By matching the peaks in Fe(3d) and N(2p) distributions, we find five distinct types of orbital coupling between Fe and N atoms: $dxz−px$⁠, $dxy−px$⁠, $dz2−pz$⁠, $dxy−py$⁠, and $dx2−y2−py$⁠. According to the relative locations of PDOS peaks, the Fe(3d) orbitals split into two groups: (⁠$dz2,dx2−y2$⁠) and (⁠$dxy$⁠, $dxz$⁠, $dyz$⁠) in the Kh0 lattice and (⁠$dz2$⁠) and (⁠$dxy$⁠, $dxz$⁠, $dyz,dx2−y2$⁠) in the Nah0 and Csh0 lattices. This difference is at least in part a result of the different cation site occupations in these lattices. As discussed in Sec. III A, the co-planar arrangement of cations and ligands in Nah0 and Csh0 causes anisotropic distortion of metal octahedra, while the body-centered occupancy of K-ions in Kh0 limits lattice distortion.

We now discuss the effects of interstitial water on the electronic structure of Ahn lattices by examining near-$EFermi$ PDOS distributions in the Nah2, Kh2, and Csh2 lattices. Figure 6 shows that Fe–N coupling still persists in Ah2 lattices, as the Fe(3d) and N(2p) orbitals in Nah2 have most of their states below $EFermi$ and the orbital splitting observed for Nah0 is absent for Nah2. We can rationalize the shift in energies of Fe(3d)–N(2p) couplings in Nah2 by noticing that most Na(3s) states are found below $EFermi$ and they couple with O(2p) orbitals. Therefore, the energy of Fe(3d)–N(2p)–Na(3s) coupling is lowered by Na(3s)–O(2p) coupling, and the energies of Na(3s) states are affected by two competing effects: Na–N and Na–O interactions. In the Kh2 lattice, the Fe(3d) orbitals are split into two different groups of (⁠$dz2$⁠) and (⁠$dxy$⁠, $dxz$⁠, $dyz,dx2−y2$⁠) near $EFermi$ with the energy of the Fe(⁠$dz2$⁠) orbital being higher than other Fe(3d) orbitals. Based on the preceding discussion, we deduce that the primary differences between the electronic structure of Kh0 and Kh2 are caused by the ordering of K-ions. Thus, we postulate that the increased energy of the Fe(⁠$dz2$⁠) orbital is caused by strengthened K-ligand interactions. A nonzero PDOS distribution near $EFermi$ is found for O(2p) orbitals in the Csh2 lattice. Given the fact that the Ni–N bond breaks in highly hydrated Cs-intercalated lattices (see Fig. 3), the Fe(3d)–N(2p)–Cs(5p)–O(2p) coupling here indicates that water molecules contribute to electron-conducting orbitals near broken Ni–N bonds.

Figure S8 in the supplementary material shows the PDOS distributions for Ni(3d) orbitals near $EFermi$ in the Ah0 and Ah2 lattices. While non-zero PDOS is found at $EFermi$ in all these distributions, there are much fewer Ni(3d) states near $EFermi$ compared to Fe(3d) states. Thus, Ni(3d) orbitals are not the primary pathways for electron transfer in Ahn lattices. As we shall see in Sec. III C 2, however, orbital coupling between Ni(3d) and other orbitals plays a critical role in the insulator/semiconductor transition of Afn lattices.

Figure 7 shows the PDOS distributions of Fe(3d), Ni(3d), and N(2p) orbitals in Naf0, Kf0, and Csf0 lattices, as well as the valence orbitals of intercalated cations. The bandgap of each Af0 lattice is about 1.3eV regardless of cation types, similar to previous DFT calculations.^88,89 The lowest unoccupied crystal orbital (LUCO) is a result of Ni–N–cation interactions as PDOS peaks are found about 1 eV above $EFermi$ in the PDOS distributions of Ni(3d) and N(2p) orbitals and in those of the valence orbitals of cations. On the contrary, Fe(3d) orbitals do not contribute to the LUCOs in Afn lattices as they have little to no PDOS above $EFermi$⁠. In the Csf0 lattice, Ni(⁠$dxz$⁠), Ni(⁠$dxy$⁠), Ni(⁠$dyz$⁠), and Ni(⁠$dx2−y2$⁠) orbitals contribute to the LUCO, and they split into two groups with their energies following a sequence of $Edxz≃Edxy<Edyz≃Edx2−y2$⁠. As demonstrated in the following discussion, the splitting of the Ni(3d) orbital in the Csf0 lattice indicates Jahn–Teller distortion of $NiN6$ octahedra induced by Cs-ligand interaction.⁹⁰ 

Examination of PDOS for hydrated, fully intercalated configurations in Fig. 8 shows that $EFermi$ is boosted to the lower edge of LUCO at different lattice hydration degrees in different Afn lattices. We, therefore, define a critical hydration degree (CHD) for each Afn lattice as the hydration degree at which PDOS at $EFermi$ becomes non-zero (i.e., at the lower edge of LUCO). CHDs for Nafn, Kfn, and Csfn lattices have values of 4, 4, and 2, respectively. Figure S9 in the supplementary material shows that the relative distance between $EFermi$ and the lower edge of LUCO gradually reduces as the hydration degree of a given Afn lattice increases from n$=0$ to its corresponding CHD. The effect of interstitial water in a given Afn lattice is also reflected in the coupling at LUCO. While Ni(3d) orbitals show no splitting in Af0 lattices (see Fig. S9 in the supplementary material), the energies of Ni(3d) orbitals in the Csf2 follow the sequence of $Edxz≃Edxy<Edyz≃Edx2−y2$⁠, though they are nearly degenerate in the Naf4 lattice.

To understand the difference in CHD between different Afn lattices, we visualize the electron density distribution (EDD) on the HOCOs of Naf0, Naf4, Csf0, and Csf2 lattices using the extra-electron method introduced in Sec. II B. Because the EDDs shown in Fig. 9 are not localized near metal centers or cyanide ligands, we find no formation of any small polaron that would act to reduce the mobility of charge carriers.⁹¹ As labeled in Fig. 9, the EDD of HOCO in Naf0 is only found on the lattice plane of $(111¯)$⁠, but in Csf0 it can be found on $(111¯)$ and $(01¯1)$ planes. The co-planar EDD in Naf0 is caused by the d-$π$ orbital coupling between metal centers and N atoms (see the upper right inlet of Fig. 9), as suggested by the results in Fig. 7. However, the co-planar symmetry of d-$π$ coupling is broken in Csf0, which rearranges the sequence of 3d-orbital energies at metal centers. Because orbital energy rearrangement implies disturbances in weak ligand fields, we calculate diagonal lengths of $NiN6$ octahedra in Csf0 and Naf0 lattices. Each $NiN6$ octahedron in Csf0 has three mutually perpendicular diagonals with their lengths being 4.16, 4.26, and 4.27 Å. The short diagonal of 4.16Å is evidence of “z-in” Jahn–Teller distortion⁹⁰ of the $NiN6$ octahedron, which increases the energies of orbitals by a z factor. On the other hand, the diagonals have an average length of 2.09 Å in Naf0 with a small standard deviation (STD) of 0.002 Å, indicating no distortion at $NiN6$ octahedra in the Naf0 lattice.

Figure 9 also shows that the co-planar d-$π$ coupling is entirely absent from Naf4 and Csf2 lattices. From the perspective of electron transport, Naf4 and Csf2 lattices become electronically conducting because an extended HOCO increases the mobility of electrons (or holes) by reducing energy barriers at lattice defects.⁹¹ From Naf0 to Naf4, the average diagonal length of the $NiN6$ octahedron increases only by 0.05 Å, but the STD among their diagonal lengths increases from 0.002 to 0.05 Å. Given the fact that Na-ions are hydrophilic, the Na–water cluster is not large enough to cause anisotropic distortion of the $NiN6$ octahedron for hydration levels less than the CHD. On the contrary, three diagonals of the $NiN6$ octahedra in Csf2 have lengths of 4.27, 4.27, and 4.49 Å, implying the “z-out” Jahn–Teller distortion (see the lower right inlet in Fig. 9). Therefore, only two water molecules are needed near each Cs-ion to cause anisotropic distortion of metal octahedra.

We now summarize the effects of interstitial species on the electronic structure of $Aδn$ lattices that are deduced from the present first-principles analysis. Cations rearrange the energetic sequence of 3d-orbitals at metal centers through cation–ligand interactions. In Ahn lattices with $n>0$⁠, interstitial water molecules coordinate with hydrophilic cations to lower the energy of metal–ligand–cation orbital coupling. When a given Afn lattice is not sufficiently hydrated, co-planar d-$π$ orbital coupling impedes electron transfer on the framework lattices. To break co-planar coupling, a large distortion at metal octahedra may be caused by the occupation of large cation–water clusters in the interstitial space. Because Na-ions coordinate with more water molecules than Cs-ion and thus form large cation–water clusters, we find a CHD in Nafn lattices that is higher than that in Csfn lattices.

The present ab initio investigation of hydrated, cation-intercalated nickel hexacyanoferrate [$NiFe(CN)6$] establishes relationships between the atomic interactions and electrochemical/electronic properties of framework lattices, identifying principles for the design of high-performance PBA electrodes in aqueous electrochemical systems. In this study, we relate interstitial atomic interactions to the properties of $NiFe(CN)6$ lattices by examining cation ordering, grand potential, and electronic structure of various $Aδn$ lattices. To find the ground states of hydrated lattices, we adopted an ML+SA+DFT structural optimization method that combines the predictive power of XGBoost ML models with the sampling efficiency of simulated annealing processes to explore complicated energy landscapes. The method is found to provide good starting points for high-fidelity relaxation calculations based on DFT, from which we find ground-state configurations that possess symmetries in the arrangements of interstitial species not identified a priori.

We find that the different arrangements of interstitial species observed are related to the tendency of the $NiFe(CN)6$ framework to swell upon lattice hydration or upon cation intercalation. When the $NiFe(CN)6$ lattice is intercalated with small and hydrophilic cations, such as Na-ions, the lattice volume remains unchanged as hydration degree increases, while interstitial water induces lattice distortion by interacting with cyanide ligands. Li-ions are smaller and more hydrophilic than Na-ions,²⁶ so the intercalation of Li-ions is likely to incur more water co-intercalation into PBA lattices, the effect of which has been linked previously to sluggish transport of Li-ions in interstitial space.⁹² On the contrary, interstitial water molecules introduce large expansion in K- and Cs-intercalated lattices with the exception of Khn lattices. In Khn lattices, water molecules rearrange cation ordering by filling voids between K-ions. Also, interstitial water undermines the integrity of framework lattices by attacking Ni–N bonds, indicating the structural instability of highly hydrated $Kδn$ and $Csδn$ lattices.

The energies of ground states are used in grand potential analysis to predict electrochemical potentials for the intercalation of alkali cations. Compared to previous analyses that considered only anhydrous lattices, our predictions based on hydrated lattices agree better with experimental data. Our analysis also predicts the correct sequence of intercalation selectivity being $Na+≤K+<Cs+$ at room temperature. By studying the variation of equilibrium hydration degree, we also find that Na-intercalated lattices remain highly hydrated during the intercalation of cations while Cs-ion intercalation is accomplished by de-intercalation of interstitial water molecules.

The electronic properties of the framework lattice are linked to the orbital couplings between metal centers and other atoms. In half-intercalated lattices, Fe(3d)–N(2p) orbital coupling results in the HOCO for electron transfer, and the energetic order of orbital couplings is found to depend strongly on the interactions between water molecules at tilted metal octahedra. Fully-intercalated $NiFe(CN)6$ lattices are found to transition from insulator to semiconductor at certain hydration degrees that depend on the intercalated cation type. For lattices intercalated with Na-ions, the transition happens when the cation–water cluster is large enough to cause anisotropic distortion of metal octahedra. In Li-intercalated PBAs, the interactions between the framework and Li-water clusters are only possible when the Li-ion is densely covered by water molecules. Hence, such a transition is expected to happen at higher degrees of lattice hydration for Li-intercalated lattices than for Na-intercalated lattices. On the other hand, intercalation of Cs-ions causes the transition to occur at lower hydration degree by intensifying cation–framework interaction that breaks the co-planar symmetry of d-$π$ orbital coupling.

We note that the methods developed in this study can be employed to investigate other classes of intercalation materials. By carefully choosing a particular functional for DFT calculations, we show the possibility of accurately predicting the electrochemical and electronic properties of $NiFe(CN)6$ lattices while keeping computational cost low.

See the supplementary material for the description of the ML+SA+DFT method, the phonon density of states of Naf0 lattices, the ordering of interstitial species in Nah4, Naf4, Kh4, and Kf4 lattices, grand potential analysis using different exchange-correlation functionals, and orbital coupling between Ni centers and cyanide ligands.

This research was funded by the U.S. National Science Foundation (Award No. 1931659) and the Department of Mechanical Science and Engineering at the University of Illinois at Urbana-Champaign. Access to the Blue Waters supercomputer was provided through a Director’s Award from the National Center for Supercomputing Applications (NCSA). Access to the Expanse and Bridges-2 computing platforms, respectively, at the San Diego Supercomputing Center (SDSC) and the Penn State Supercomputing Center (PSC) was supported by the Extreme Science and Engineering Discovery Environment (XSEDE⁹³ Award No. TG-PHY210018), which is supported by the National Science Foundation (NSF) under Grant No. ACI-1548562. We also wish to acknowledge technical support at SDSC and Materials Square (MatSQ) for offering suggestions and help in executing the DFT calculations presented here.

The authors have no conflicts to disclose.

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