Prussian blue analogs (PBAs) are an important material class for aqueous electrochemical separations and energy storage owing to their ability to reversibly intercalate monovalent cations. However, incorporating interstitial H2O molecules in the ab initio study of PBAs is technically challenging, though essential to understanding the interactions between interstitial water, interstitial cations, and the framework lattice that affect intercalation potential and cation intercalation selectivity. Accordingly, we introduce and use a method that combines the efficiency of machine-learning models with the accuracy of ab initio calculations to elucidate mechanisms of (1) lattice expansion upon intercalation of cations of different sizes, (2) selectivity bias toward intercalating hydrophobic cations of large size, and (3) semiconductor–conductor transitions from anhydrous to hydrated lattices. We analyze the PBA nickel hexacyanoferrate [NiFe(CN)6] due to its structural stability and electrochemical activity in aqueous electrolytes. Here, grand potential analysis is used to determine the equilibrium degree of hydration for a given intercalated cation (Na+, K+, or Cs+) and NiFe(CN)6 oxidation state based on pressure-equilibrated structures determined with the aid of machine learning and simulated annealing. The results imply new directions for the rational design of future cation-intercalation electrode materials that optimize performance in various electrochemical applications, and they demonstrate the importance of choosing an appropriate calculation framework to predict the properties of PBA lattices accurately.

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]znH2O, 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 (NH4+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 activity15 and the trend of potential13 that produces their selectivity bias toward large hydrated cations, stimulating the use of nickel hexacyanoferrate [NiFe(CN)6] in Cs+ removal from nuclear waste16,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.13 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.26 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.27 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.28 A recent experiment also used a two-phase synthesis method to eliminate interstitial water in high-crystallinity NaxFe[Fe(CN)6] lattice.29 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.30 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+.33 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.34 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.35 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.36 

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,37 if the potential is too high, charge/discharge of PBAs can lead to oxygen reduction and eventually cause the explosion of sealed battery cells.38 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.41 By using a hybrid exchange-correlation functional in DFT calculations, much better agreement with experiments is achieved at the expense of computational cost.37 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.42 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,37 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.37 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.43 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,44 cement hydrates,45 and liquid-water/metal–oxide interface.46 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.47 While the transferability of classical force field can be improved between the small hydrated system by using machine-learned parameters,48 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–water51,52 and gas–water interactions.53 The resulting interaction energies were shown to follow correct energy ordering for water clusters of varied sizes.54 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.26 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)6nH2O 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)6nH2O 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)6nH2O, 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.26 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δn26 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)64H2O and CsNiFe(CN)62H2O, 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.26 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 features26 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 ESPRESSO55 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 104 Ry/bohr (2.57×103eV/Å). The exchange-correlation interactions were approximated using the PBE functional,56 and non-local van der Waals interactions are captured by the rVV10 functional57 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)61 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.62 Also, the crystal field of non-magnetic cyanide ligands screens the superexchange interactions63 between nearest-neighbor metal centers, causing non-collinear magnetic ordering of paramagnetic atoms.64 Finally, interstitial water can induce conversion between ferromagnetic and paramagnetic phases of the NiFe(CN)6 lattice by distorting the framework structure.65 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.66 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.57 

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.55 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 107 eV.

In Secs. III AIII 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.69 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 studies37,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.

Lattice volume ratios of hydrated lattices relative to anhydrous lattices. The unit-cell volumes of anhydrous lattices are 274.048 Å3 (Naf0), 276.024 Å3 (Nah0), 267.270 Å3 (Kf0),282.461 Å3 (Kh0),286.108 Å3 (Csf0), and 275.230 Å3 (Csh0).

NahnNafnKhnKfnCshnCsfn
n=1 1.000 0.999 0.998 1.008 1.000 1.002 
n=2 0.992 0.993 0.981 1.053 1.018 1.031 
n=3 1.027 1.025 0.996 1.097 1.057 1.297 
n=4 1.039 1.062 1.077 1.176 1.166 1.392 
NahnNafnKhnKfnCshnCsfn
n=1 1.000 0.999 0.998 1.008 1.000 1.002 
n=2 0.992 0.993 0.981 1.053 1.018 1.031 
n=3 1.027 1.025 0.996 1.097 1.057 1.297 
n=4 1.039 1.062 1.077 1.176 1.166 1.392 

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 group71 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.

FIG. 1.

Ground states of Nah0, Nah2, Naf0, and Naf2 lattices, in which the Corey–Pauling–Koltun (CPK)72 coloring scheme is applied to all elements. Ni (green) and Fe (brown) atoms are connected by cyanide ligands at which nitrogen atoms (blue) connect to Ni centers and carbon atoms (gray) connect to Fe centers. In the interstitial space, Na-ions are represented by purple spheres, and each oxygen (red) is connected to two hydrogen (white) atoms.

FIG. 1.

Ground states of Nah0, Nah2, Naf0, and Naf2 lattices, in which the Corey–Pauling–Koltun (CPK)72 coloring scheme is applied to all elements. Ni (green) and Fe (brown) atoms are connected by cyanide ligands at which nitrogen atoms (blue) connect to Ni centers and carbon atoms (gray) connect to Fe centers. In the interstitial space, Na-ions are represented by purple spheres, and each oxygen (red) is connected to two hydrogen (white) atoms.

Close modal

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 n3, 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.

FIG. 2.

Ground states of Kh0, Kh2, Kf0, and Kf2 lattices, in which the CPK coloring scheme is applied to all elements (see Fig. 1 for details). In the interstitial space, K-ions are represented by purple spheres.

FIG. 2.

Ground states of Kh0, Kh2, Kf0, and Kf2 lattices, in which the CPK coloring scheme is applied to all elements (see Fig. 1 for details). In the interstitial space, K-ions are represented by purple spheres.

Close modal

The volume of the Kfn lattice remains nearly unchanged from Kf0 to Kf1 but experiences at least 5% expansion compared to Kf0 when n2. 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,73 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).

FIG. 3.

Ground states of Csh0, Csh4, Csf0, and Csf4 lattices, in which CPK coloring scheme is applied to all elements (see Fig. 1 for details). In the interstitial space, Cs-ions are represented by purple spheres. The hydrogen bonding between oxygen and hydrogen atoms is labeled using black dashed lines.

FIG. 3.

Ground states of Csh0, Csh4, Csf0, and Csf4 lattices, in which CPK coloring scheme is applied to all elements (see Fig. 1 for details). In the interstitial space, Cs-ions are represented by purple spheres. The hydrogen bonding between oxygen and hydrogen atoms is labeled using black dashed lines.

Close modal

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 n2. 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 PBELYP1W74+rVV10 functional for the following reasons. While the PBE functional leads to accurate predictions of bond distances in small molecules75 and layered materials,76 it underestimates cation-molecule binding energies77,78 and overstructures cation–water agglomerates.79 The inaccuracy of the PBE functional in predicting cation–water arrangements can be reduced by using the rVV10 functional57 in conjunction, but the PBE+rVV10 combination can be worse at reproducing measured hydration energies of crystalline materials.80 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,42 

ΦAδn=EAδnDFTn×μH2O.
(1)

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

μH2O=EH2ODFT+Δμ,
(2)

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.81 However, the implementation of Pitzer model requires information of electrolyte density as a function of temperature and ionic concentration,82 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

ANiFe(CN)6n1H2O+AA2NiFe(CN)6n2H2O,
(3)

is calculated as

ϕA=(ΦAhn1+EAΦAfn2)/e,,
(4)

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

n~Ah={Δμ|n:min(ΦAhn)}andn~Af={Δμ|n:min(ΦAfn)}.
(5)

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

ϕeq,A=(ΦAhn~Ah+EAΦAfn~Af)/e+EA/A+SHE,
(6)

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:83 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

Etrans=kBTlnPoΛ3kBT,
(7)

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 function84 

qr=π122(2IAkBT2)12(2IBkBT2)12(2ICkBT2)12
(8)

to calculate the rotational energy as Erot=kBTlnqr, 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

Evib=kBTlnqv=kBTln(i=13eΘv,i/2T1eΘv,i/T),
(9)

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.85 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.86 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 Δμ=0n~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 (24 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.87 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.

FIG. 4.

Equilibrium lattice hydration degree n~ as a function of shift in the chemical potential of water Δμ for Naδn (a), Kδn (b), and Csδn (c) lattices. The variation of equilibrium intercalation potential with respect to SHE vs Δμ is shown in (d) for Na-, K-, and Cs-ions.

FIG. 4.

Equilibrium lattice hydration degree n~ as a function of shift in the chemical potential of water Δμ for Naδn (a), Kδn (b), and Csδn (c) lattices. The variation of equilibrium intercalation potential with respect to SHE vs Δμ is shown in (d) for Na-, K-, and Cs-ions.

Close modal

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.13 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.13 Compared to grand potential analysis that only considers anhydrous lattices,39 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.

1. Electronic structures of Ahn lattices

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.

FIG. 5.

PDOS distributions on the orbitals of Fe(3d), C(2p), and N(2p) in Nah0, Kh0, and Csh0 lattices, as well as the valence orbitals of intercalated cations. The Fermi energy levels EFermi are shifted to zero in the subplots and are labeled using dashed lines. To match the locations of PDOS peaks between different distributions, all the subplots adopt the same range of x-axis from 2 to 2 eV.

FIG. 5.

PDOS distributions on the orbitals of Fe(3d), C(2p), and N(2p) in Nah0, Kh0, and Csh0 lattices, as well as the valence orbitals of intercalated cations. The Fermi energy levels EFermi are shifted to zero in the subplots and are labeled using dashed lines. To match the locations of PDOS peaks between different distributions, all the subplots adopt the same range of x-axis from 2 to 2 eV.

Close modal

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: dxzpx, dxypx, dz2pz, dxypy, and dx2y2py. According to the relative locations of PDOS peaks, the Fe(3d) orbitals split into two groups: (dz2,dx2y2) and (dxy, dxz, dyz) in the Kh0 lattice and (dz2) and (dxy, dxz, dyz,dx2y2) 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,dx2y2) 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.

FIG. 6.

PDOS distributions on the orbitals of Fe(3d), C(2p), and N(2p) in Nah2, Kh2, and Csh2 lattices, as well as the valence orbitals of intercalated cations.

FIG. 6.

PDOS distributions on the orbitals of Fe(3d), C(2p), and N(2p) in Nah2, Kh2, and Csh2 lattices, as well as the valence orbitals of intercalated cations.

Close modal

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 C2, however, orbital coupling between Ni(3d) and other orbitals plays a critical role in the insulator/semiconductor transition of Afn lattices.

2. Electronic structures 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(dx2y2) orbitals contribute to the LUCO, and they split into two groups with their energies following a sequence of EdxzEdxy<EdyzEdx2y2. 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.90 

FIG. 7.

PDOS distributions on the orbitals of Fe(3d), Ni(3d), and N(2p) in Naf0, Kf0, and Csf0 lattices, as well as the valence orbitals of intercalated cations.

FIG. 7.

PDOS distributions on the orbitals of Fe(3d), Ni(3d), and N(2p) in Naf0, Kf0, and Csf0 lattices, as well as the valence orbitals of intercalated cations.

Close modal

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 EdxzEdxy<EdyzEdx2y2, though they are nearly degenerate in the Naf4 lattice.

FIG. 8.

PDOS distributions on the orbitals of Fe(3d), Ni(3d), and N(2p) in Naf4, Kf4, and Csf4 lattices, as well as the valence orbitals of intercalated cations.

FIG. 8.

PDOS distributions on the orbitals of Fe(3d), Ni(3d), and N(2p) in Naf4, Kf4, and Csf4 lattices, as well as the valence orbitals of intercalated cations.

Close modal

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.91 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 distortion90 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.

FIG. 9.

Charge distribution on HOCO in Naf0, Naf4, Csf0, and Csf2 lattices. The isosurface value of EDDs is set to be 104bohr3, and the interstitial species are not shown here to make the plots clear. Diagrams to the right demonstrate the co-planar coupling of 3d and 2p orbitals (top), “z-in” Jahn–Teller distortion in the Csf0 lattice (middle), and “z-out” Jahn–Teller distortion in the Csf2 lattice (bottom).

FIG. 9.

Charge distribution on HOCO in Naf0, Naf4, Csf0, and Csf2 lattices. The isosurface value of EDDs is set to be 104bohr3, and the interstitial species are not shown here to make the plots clear. Diagrams to the right demonstrate the co-planar coupling of 3d and 2p orbitals (top), “z-in” Jahn–Teller distortion in the Csf0 lattice (middle), and “z-out” Jahn–Teller distortion in the Csf2 lattice (bottom).

Close modal

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.91 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,26 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.92 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 (XSEDE93 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.

1.
S.
Klink
,
Y.
Ishige
, and
W.
Schuhmann
, “
Prussian blue analogues: A versatile framework for solid-contact ion-selective electrodes with tunable potentials
,”
ChemElectroChem
4
,
490
494
(
2017
).
2.
W.-J.
Li
,
C.
Han
,
G.
Cheng
,
S.-L.
Chou
,
H.-K.
Liu
, and
S.-X.
Dou
, “
Chemical properties, structural properties, and energy storage applications of prussian blue analogues
,”
Small
15
,
1900470
(
2019
).
3.
R.
Chen
,
Y.
Huang
,
M.
Xie
,
Z.
Wang
,
Y.
Ye
,
L.
Li
, and
F.
Wu
, “
Chemical inhibition method to synthesize highly crystalline prussian blue analogs for sodium-ion battery cathodes
,”
ACS Appl. Mater. Interfaces
8
,
31669
31676
(
2016
).
4.
H. L. B.
Boström
,
I. E.
Collings
,
D.
Daisenberger
,
C. J.
Ridley
,
N. P.
Funnell
, and
A. B.
Cairns
, “
Probing the influence of defects, hydration, and composition on prussian blue analogues with pressure
,”
J. Am. Chem. Soc.
143
,
3544
3554
(
2021
).
5.
B.
Xie
,
P.
Zuo
,
L.
Wang
,
J.
Wang
,
H.
Huo
,
M.
He
,
J.
Shu
,
H.
Li
,
S.
Lou
, and
G.
Yin
, “
Achieving long-life prussian blue analogue cathode for Na-ion batteries via triple-cation lattice substitution and coordinated water capture
,”
Nano Energy
61
,
201
210
(
2019
).
6.
Z.
Li
,
T.
Liu
,
R.
Meng
,
L.
Gao
,
Y.
Zou
,
P.
Peng
,
Y.
Shao
, and
X.
Liang
, “
Insights into the structure stability of prussian blue for aqueous zinc ion batteries
,”
Energy Environ. Mater.
4
,
111
116
(
2021
).
7.
A.
Shrivastava
,
S.
Liu
, and
K. C.
Smith
, “
Linking capacity loss and retention of nickel hexacyanoferrate to a two-site intercalation mechanism for aqueous Mg2+ and Ca2+ ions
,”
Phys. Chem. Chem. Phys.
21
,
20177
20188
(
2019
).
8.
M.
Oliver-Tolentino
,
G.
Ramos-Sánchez
,
G.
Guzmán
,
M.
Avila
,
I.
González
, and
E.
Reguera
, “
Water effect on sodium mobility in zinc hexacyanoferrate during charge/discharge processes in sodium ion-based battery
,”
Solid State Ionics
312
,
67
72
(
2017
).
9.
B.
Xie
,
L.
Wang
,
J.
Shu
,
X.
Zhou
,
Z.
Yu
,
H.
Huo
,
Y.
Ma
,
X.
Cheng
,
G.
Yin
, and
P.
Zuo
, “
Understanding the structural evolution and lattice water movement for rhombohedral nickel hexacyanoferrate upon sodium migration
,”
ACS Appl. Mater. Interfaces
11
,
46705
46713
(
2019
).
10.
M.
Okubo
,
D.
Asakura
,
Y.
Mizuno
,
J.-D.
Kim
,
T.
Mizokawa
,
T.
Kudo
, and
I.
Honma
, “
Switching redox-active sites by valence tautomerism in prussian blue analogues AxMny[Fe(CN)6]⋅ nH2O (A: K, Rb): Robust frameworks for reversible Li storage
,”
J. Phys. Chem. Lett.
1
,
2063
2071
(
2010
).
11.
Y.
Mizuno
,
M.
Okubo
,
E.
Hosono
,
T.
Kudo
,
H.
Zhou
, and
K.
Oh-ishi
, “
Suppressed activation energy for interfacial charge transfer of a prussian blue analog thin film electrode with hydrated ions (Li+, Na+, and Mg+)
,”
J. Phys. Chem. C
117
,
10877
10882
(
2013
).
12.
J.
Wu
,
J.
Song
,
K.
Dai
,
Z.
Zhuo
,
L. A.
Wray
,
G.
Liu
,
Z.-X.
Shen
,
R.
Zeng
,
Y.
Lu
, and
W.
Yang
, “
Modification of transition-metal redox by interstitial water in hexacyanometalate electrodes for sodium-ion batteries
,”
J. Am. Chem. Soc.
139
,
18358
18364
(
2017
).
13.
C. D.
Wessells
,
S. V.
Peddada
,
M. T.
McDowell
,
R. A.
Huggins
, and
Y.
Cui
, “
The effect of insertion species on nanostructured open framework hexacyanoferrate battery electrodes
,”
J. Electrochem. Soc.
159
,
A98
(
2011
).
14.
T.
Ikeshoji
, “
Separation of alkali metal ions by intercalation into a prussian blue electrode
,”
J. Electrochem. Soc.
133
,
2108
(
1986
).
15.
K.
Itaya
,
T.
Ataka
, and
S.
Toshima
, “
Spectroelectrochemistry and electrochemical preparation method of prussian blue modified electrodes
,”
J. Am. Chem. Soc.
104
,
4767
4772
(
1982
).
16.
S.
Rassat
,
J.
Sukamto
,
R.
Orth
,
M.
Lilga
, and
R.
Hallen
, “
Development of an electrically switched ion exchange process for selective ion separations
,”
Sep. Purif. Technol.
15
,
207
222
(
1999
).
17.
B.
Sun
,
X. G.
Hao
,
Z. D.
Wang
,
G. Q.
Guan
,
Z. L.
Zhang
,
Y. B.
Li
, and
S. B.
Liu
, “
Separation of low concentration of cesium ion from wastewater by electrochemically switched ion exchange method: Experimental adsorption kinetics analysis
,”
J. Hazard. Mater.
233-234
,
177
183
(
2012
).
18.
K. C.
Smith
, “
Theoretical evaluation of electrochemical cell architectures using cation intercalation electrodes for desalination
,”
Electrochim. Acta
230
,
333
341
(
2017
).
19.
S.
Porada
,
A.
Shrivastava
,
P.
Bukowska
,
P. M.
Biesheuvel
, and
K. C.
Smith
, “
Nickel hexacyanoferrate electrodes for continuous cation intercalation desalination of brackish water
,”
Electrochim. Acta
255
,
369
378
(
2017
).
20.
J.
Lee
,
S.
Kim
, and
J.
Yoon
, “
Rocking chair desalination battery based on prussian blue electrodes
,”
ACS Omega
2
,
1653
1659
(
2017
).
21.
E. R.
Reale
,
A.
Shrivastava
, and
K. C.
Smith
, “
Effect of conductive additives on the transport properties of porous flow-through electrodes with insulative particles and their optimization for Faradaic deionization
,”
Water Res.
165
,
114995
(
2019
).
22.
K.
Singh
,
Z.
Qian
,
P.
Biesheuvel
,
H.
Zuilhof
,
S.
Porada
, and
L. C.
de Smet
, “
Nickel hexacyanoferrate electrodes for high mono/divalent ion-selectivity in capacitive deionization
,”
Desalination
481
,
114346
(
2020
).
23.
E.
Sebti
,
M. M.
Besli
,
M.
Metzger
,
S.
Hellstrom
,
M. J.
Schultz-Neu
,
J.
Alvarado
,
J.
Christensen
,
M.
Doeff
,
S.
Kuppan
, and
C. V.
Subban
, “
Removal of Na+ and Ca2+ with prussian blue analogue electrodes for brackish water desalination
,”
Desalination
487
,
114479
(
2020
).
24.
V.
Pothanamkandathil
,
J.
Fortunato
, and
C. A.
Gorski
, “
Electrochemical desalination using intercalating electrode materials: A comparison of energy demands
,”
Environ. Sci. Technol.
54
,
3653
3662
(
2020
).
25.
E. R.
Reale
,
L.
Regenwetter
,
A.
Agrawal
,
B.
Dardón
,
N.
Dicola
,
S.
Sanagala
, and
K. C.
Smith
, “
Low porosity, high areal-capacity prussian blue analogue electrodes enhance salt removal and thermodynamic efficiency in symmetric Faradaic deionization with automated fluid control
,”
Water Res. X
13
,
100116
(
2021
).
26.
S.
Liu
and
K. C.
Smith
, “
Linking the polyatomic arrangements of interstitial H2O and cations to bonding within prussian blue analogues ab initio using gradient-boosted machine learning
,”
Phys. Rev. Mater.
5
,
035003
(
2021
).
27.
X.
Wu
,
J. J.
Hong
,
W.
Shin
,
L.
Ma
,
T.
Liu
,
X.
Bi
,
Y.
Yuan
,
Y.
Qi
,
T. W.
Surta
,
W.
Huang
,
J.
Neuefeind
,
T.
Wu
,
P. A.
Greaney
,
J.
Lu
, and
X.
Ji
, “
Diffusion-free Grotthuss topochemistry for high-rate and long-life proton batteries
,”
Nat. Energy
4
,
123
130
(
2019
).
28.
L.
Wang
,
Y.
Lu
,
J.
Liu
,
M.
Xu
,
J.
Cheng
,
D.
Zhang
, and
J. B.
Goodenough
, “
A superior low-cost cathode for a Na-ion battery
,”
Angew. Chem.
125
,
2018
2021
(
2013
).
29.
C.
Duan
,
Y.
Meng
,
Y.
Wang
,
Z.
Zhang
,
Y.
Ge
,
X.
Li
,
Y.
Guo
, and
D.
Xiao
, “
High-crystallinity and high-rate prussian blue analogues synthesized at the oil–water interface
,”
Inorg. Chem. Front.
8
,
2008
2016
(
2021
).
30.
L.
Li
,
Z.
Hu
,
Y.
Lu
,
C.
Wang
,
Q.
Zhang
,
S.
Zhao
,
J.
Peng
,
K.
Zhang
,
S.-L.
Chou
, and
J.
Chen
, “
A low-strain potassium-rich prussian blue analogue cathode for high power potassium-ion batteries
,”
Angew. Chem. Int. Ed.
60
,
13050
(
2021
).
31.
J.
Song
,
L.
Wang
,
Y.
Lu
,
J.
Liu
,
B.
Guo
,
P.
Xiao
,
J.-J.
Lee
,
X.-Q.
Yang
,
G.
Henkelman
, and
J. B.
Goodenough
, “
Removal of interstitial H2O in hexacyanometallates for a superior cathode of a sodium-ion battery
,”
J. Am. Chem. Soc.
137
,
2658
2664
(
2015
).
32.
H. L.
Boström
,
I. E.
Collings
,
A. B.
Cairns
,
C. P.
Romao
, and
A. L.
Goodwin
, “
High-pressure behaviour of prussian blue analogues: Interplay of hydration, Jahn-Teller distortions and vacancies
,”
Dalton Trans.
48
,
1647
1655
(
2019
).
33.
L.
Guadagnini
,
D.
Tonelli
, and
M.
Giorgetti
, “
Improved performances of electrodes based on Cu2+-loaded copper hexacyanoferrate for hydrogen peroxide detection
,”
Electrochim. Acta
55
,
5036
5039
(
2010
).
34.
Y.
Li
,
Q.
Gao
,
D.
Chang
,
P.
Sun
,
J.
Liu
,
Y.
Jia
,
E.
Liang
, and
Q.
Sun
, “
Effect of bond on negative thermal expansion of prussian blue analogues MCo(CN)6 (M= Fe, Ti and Sc): A first-principles study
,”
J. Phys.: Condens. Matter.
32
,
455703
(
2020
).
35.
M.
Jiang
,
X.
Fan
,
S.
Cao
,
Z.
Wang
,
Z.
Yang
, and
W.
Zhang
, “
Thermally activated carbon–nitrogen vacancies in double-shelled NiFe prussian blue analogue nanocages for enhanced electrocatalytic oxygen evolution
,”
J. Mater. Chem. A
9
,
12734
12745
(
2021
).
36.
P.
Nie
,
J.
Yuan
,
J.
Wang
,
Z.
Le
,
G.
Xu
,
L.
Hao
,
G.
Pang
,
Y.
Wu
,
H.
Dou
,
X.
Yan
et al., “
Prussian blue analogue with fast kinetics through electronic coupling for sodium ion batteries
,”
ACS Appl. Mater. Interfaces
9
,
20306
20312
(
2017
).
37.
X.
Guo
,
Z.
Wang
,
Z.
Deng
,
X.
Li
,
B.
Wang
,
X.
Chen
, and
S. P.
Ong
, “
Water contributes to higher energy density and cycling stability of prussian blue analogue cathodes for aqueous sodium-ion batteries
,”
Chem. Mater.
31
,
5933
5942
(
2019
).
38.
L.
Chen
,
L.
Cao
,
X.
Ji
,
S.
Hou
,
Q.
Li
,
J.
Chen
,
C.
Yang
,
N.
Eidson
, and
C.
Wang
, “
Enabling safe aqueous lithium ion open batteries by suppressing oxygen reduction reaction
,”
Nat. Commun.
11
,
2638
(
2020
).
39.
S.
Liu
and
K. C.
Smith
, “
Intercalated cation disorder in prussian blue analogues: First-principles and grand canonical analyses
,”
J. Phys. Chem. C
123
,
10191
10204
(
2019
).
40.
Y.
Tang
,
W.
Li
,
P.
Feng
,
M.
Zhou
,
K.
Wang
, and
K.
Jiang
, “
Investigation of alkali-ion (Li, Na and K) intercalation in manganese hexacyanoferrate KxMnFe(CN)6 as cathode material
,”
Chem. Eng. J.
396
,
125269
(
2020
).
41.
H.
Tang
,
J.
Tao
,
J. P.
Perdew
et al., “
Density functionals combined with van der Waals corrections for graphene adsorbed on layered materials
,”
Phys. Rev. B
101
,
195426
(
2020
).
42.
G.
Sai Gautam
,
P.
Canepa
,
W. D.
Richards
,
R.
Malik
, and
G.
Ceder
, “
Role of structural H2O in intercalation electrodes: The case of Mg in nanocrystalline Xerogel-V2O5
,”
Nano Lett.
16
,
2426
2431
(
2016
), pMID: 26982964.
43.
P.
Xiao
,
J.
Song
,
L.
Wang
,
J. B.
Goodenough
, and
G.
Henkelman
, “
Theoretical study of the structural evolution of a Na2FeMn(CN)6 cathode upon Na intercalation
,”
Chem. Mater.
27
,
3763
3768
(
2015
).
44.
O.
Egorova
,
R.
Hafizi
,
D. C.
Woods
, and
G. M.
Day
, “
Multifidelity statistical machine learning for molecular crystal structure prediction
,”
J. Phys. Chem. A
124
,
8065
8078
(
2020
).
45.
K.
Kobayashi
,
H.
Nakamura
,
A.
Yamaguchi
,
M.
Itakura
,
M.
Machida
, and
M.
Okumura
, “
Machine learning potentials for tobermorite minerals
,”
Comput. Mater. Sci.
188
,
110173
(
2021
).
46.
V.
Quaranta
,
J.
Behler
, and
M.
Hellström
, “
Structure and dynamics of the liquid–water/zinc-oxide interface from machine learning potential simulations
,”
J. Phys. Chem. C
123
,
1293
1304
(
2018
).
47.
S.
Dasetty
,
P. J.
Meza-Morales
,
R. B.
Getman
, and
S.
Sarupria
, “
Simulations of interfacial processes: Recent advances in force field development
,”
Curr. Opin. Chem. Eng.
23
,
138
145
(
2019
).
48.
Y.
Li
,
H.
Li
,
F. C.
Pickard IV
,
B.
Narayanan
,
F. G.
Sen
,
M. K.
Chan
,
S. K.
Sankaranarayanan
,
B. R.
Brooks
, and
B.
Roux
, “
Machine learning force field parameters from ab initio data
,”
J. Chem. Theory Comput.
13
,
4492
4503
(
2017
).
49.
F.
Paesani
, “
Getting the right answers for the right reasons: Toward predictive molecular simulations of water with many-body potential energy functions
,”
Acc. Chem. Res.
49
,
1844
1851
(
2016
).
50.
O.
Sode
and
J. N.
Cherry
, “
Development of a flexible-monomer two-body carbon dioxide potential and its application to clusters up to (CO2)13
,”
J. Comput. Chem.
38
,
2763
2774
(
2017
).
51.
P.
Bajaj
,
A. W.
Gotz
, and
F.
Paesani
, “
Toward chemical accuracy in the description of ion–water interactions through many-body representations. I. Halide–water dimer potential energy surfaces
,”
J. Chem. Theory Comput.
12
,
2698
2705
(
2016
).
52.
M.
Riera
,
N.
Mardirossian
,
P.
Bajaj
,
A. W.
Götz
, and
F.
Paesani
, “
Toward chemical accuracy in the description of ion–water interactions through many-body representations. alkali-water dimer potential energy surfaces
,”
J. Chem. Phys.
147
,
161715
(
2017
).
53.
M.
Riera
,
E. P.
Yeh
, and
F.
Paesani
, “
Data-driven many-body models for molecular fluids: CO2/H2O mixtures as a case study
,”
J. Chem. Theory Comput.
16
,
2246
2257
(
2020
).
54.
E.
Lambros
,
S.
Dasgupta
,
E.
Palos
,
S.
Swee
,
J.
Hu
, and
F.
Paesani
, “
General many-body framework for data-driven potentials with arbitrary quantum mechanical accuracy: Water as a case study
,”
J. Chem. Theory Comput.
17
,
5635
5650
(
2021
).
55.
P.
Giannozzi
,
O.
Baseggio
,
P.
Bonfà
,
D.
Brunato
,
R.
Car
,
I.
Carnimeo
,
C.
Cavazzoni
,
S.
De Gironcoli
,
P.
Delugas
,
F.
Ferrari Ruffino
et al., “
Quantum espresso toward the exascale
,”
J. Chem. Phys.
152
,
154105
(
2020
).
56.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
(
1996
).
57.
R.
Sabatini
,
T.
Gorni
, and
S.
De Gironcoli
, “
Nonlocal van der Waals density functional made simple and efficient
,”
Phys. Rev. B
87
,
041108
(
2013
).
58.
S.
Adhikari
,
H.
Tang
,
B.
Neupane
,
A.
Ruzsinszky
, and
G. I.
Csonka
, “
Molecule-surface interaction from van der Waals-corrected semilocal density functionals: The example of thiophene on transition-metal surfaces
,”
Phys. Rev. Mater.
4
,
025005
(
2020
).
59.
H.
Tang
,
J. P.
Perdew
et al., “
van der Waals corrected density functionals for cylindrical surfaces: Ammonia and nitrogen dioxide adsorbed on a single-walled carbon nanotube
,”
Phys. Rev. B
103
,
195410
(
2021
).
60.
V. L.
Deringer
,
M. A.
Caro
, and
G.
Csányi
, “
A general-purpose machine-learning force field for bulk and nanostructured phosphorus
,”
Nat. Commun.
11
,
1
11
(
2020
).
61.
M.
Verdaguer
and
G. S.
Girolami
, “Magnetic prussian blue analogs,” in Magnetism: Molecules to Materials V (Wiley-VCH, 2005), Vol. 5, pp. 283–346.
62.
M.
Cammarata
,
S.
Zerdane
,
L.
Balducci
,
G.
Azzolina
,
S.
Mazerat
,
C.
Exertier
,
M.
Trabuco
,
M.
Levantino
,
R.
Alonso-Mori
,
J. M.
Glownia
et al., “
Charge transfer driven by ultrafast spin transition in a CoFe prussian blue analogue
,”
Nat. Chem.
13
,
10
14
(
2021
).
63.
J.
Kanamori
, “
Superexchange interaction and symmetry properties of electron orbitals
,”
J. Phys. Chem. Solids
10
,
87
98
(
1959
).
64.
V.
Kavečanskỳ
,
M.
Mihalik
,
Z.
Mitróová
,
M.
Lukáčová
, and
S.
Mat’aš
, “
Neutron diffraction study of crystal and magnetic structure of Dy[Fe(CN)6] 4D2O
,”
Czech. J. Phys.
54
,
571
574
(
2004
).
65.
N.
Yanai
,
W.
Kaneko
,
K.
Yoneda
,
M.
Ohba
, and
S.
Kitagawa
, “
Reversible water-induced magnetic and structural conversion of a flexible microporous Ni(II)Fe(III) ferromagnet
,”
J. Am. Chem. Soc.
129
,
3496
3497
(
2007
).
66.
L.
Talirz
,
S.
Kumbhar
,
E.
Passaro
,
A. V.
Yakutovich
,
V.
Granata
,
F.
Gargiulo
,
M.
Borelli
,
M.
Uhrin
,
S. P.
Huber
,
S.
Zoupanos
et al., “
Materials cloud, a platform for open computational science
,”
Sci. Data
7
,
1
12
(
2020
).
67.
T.
Maxisch
,
F.
Zhou
, and
G.
Ceder
, “
Ab initio study of the migration of small polarons in olivine LixFePo4 and their association with lithium ions and vacancies
,”
Phys. Rev. B
73
,
104301
(
2006
).
68.
W. H.
Sio
,
C.
Verdi
,
S.
Poncé
, and
F.
Giustino
, “
Ab initio theory of polarons: Formalism and applications
,”
Phys. Rev. B
99
,
235139
(
2019
).
69.
K.
Hurlbutt
,
S.
Wheeler
,
I.
Capone
, and
M.
Pasta
, “
Prussian blue analogs as battery materials
,”
Joule
2
,
1950
1960
(
2018
).
70.
A.
Mullaliu
,
P.
Conti
,
G.
Aquilanti
,
J. R.
Plaisier
,
L.
Stievano
, and
M.
Giorgetti
, “
Operando XAFS and XRD study of a prussian blue analogue cathode material: Iron hexacyanocobaltate
,”
Condens. Matter
3
,
36
(
2018
).
71.
Y.
Xu
,
J.
Wan
,
L.
Huang
,
M.
Ou
,
C.
Fan
,
P.
Wei
,
J.
Peng
,
Y.
Liu
,
Y.
Qiu
,
X.
Sun
et al., “
Structure distortion induced monoclinic nickel hexacyanoferrate as high-performance cathode for Na-ion batteries
,”
Adv. Energy Mater.
9
,
1803158
(
2019
).
72.
W. L.
Koltun
, “Space filling atomic units and connectors for molecular models,” U.S. patent 3,170,246 (February 23, 1965).
73.
Q.
Gao
,
N.
Shi
,
A.
Sanson
,
Y.
Sun
,
R.
Milazzo
,
L.
Olivi
,
H.
Zhu
,
S. H.
Lapidus
,
L.
Zheng
,
J.
Chen
, and
X.
Xing
, “
Tunable thermal expansion from negative, zero, to positive in cubic prussian blue analogues of GaFe(CN)6
,”
Inorg. Chem.
57
,
14027
(
2018
).
74.
E. E.
Dahlke
and
D. G.
Truhlar
, “
Improved density functionals for water
,”
J. Phys. Chem. B
109
,
15677
15683
(
2005
).
75.
J. M.
del Campo
,
J. L.
Gázquez
,
S.
Trickey
, and
A.
Vela
, “
Non-empirical improvement of PBE and its hybrid PBE0 for general description of molecular properties
,”
J. Chem. Phys.
136
,
104108
(
2012
).
76.
H.
Peng
and
J. P.
Perdew
, “
Rehabilitation of the Perdew-Burke-Ernzerhof generalized gradient approximation for layered materials
,”
Phys. Rev. B
95
,
081105
(
2017
).
77.
J.
Cha
,
S.
Lim
,
C. H.
Choi
,
M.-H.
Cha
, and
N.
Park
, “
Inaccuracy of density functional theory calculations for dihydrogen binding energetics onto Ca cation centers
,”
Phys. Rev. Lett.
103
,
216102
(
2009
).
78.
Y.
Ohk
,
Y.-H.
Kim
, and
Y.
Jung
, “
Comment on ‘Inaccuracy of density functional theory calculations for dihydrogen binding energetics onto Ca cation centers
,’”
Phys. Rev. Lett.
104
,
179601
(
2010
).
79.
J.
Schmidt
,
J.
VandeVondele
,
I.-F. W.
Kuo
,
D.
Sebastiani
,
J. I.
Siepmann
,
J.
Hutter
, and
C. J.
Mundy
, “
Isobaric-isothermal molecular dynamics simulations utilizing density functional theory: An assessment of the structure and density of water at near-ambient conditions
,”
J. Phys. Chem. B
113
,
11959
11964
(
2009
).
80.
S.
Kiyabu
,
J. S.
Lowe
,
A.
Ahmed
, and
D. J.
Siegel
, “
Computational screening of hydration reactions for thermal energy storage: New materials and design rules
,”
Chem. Mater.
30
,
2006
2017
(
2018
).
81.
K. S.
Pitzer
, “Ion interaction approach: Theory and data correlation,” in Activity Coefficients in Electrolyte Solutions (CRC Press, 2018), pp. 75–153.
82.
M. P.
Humphreys
and
A. J.
Schiller
, “Pytzer: The Pitzer model for chemical activities and equilibria in aqueous solutions in python (beta)” (2021).
83.
P.
Vanysek
, “Electrochemical series,” in CRC Handbook of Chemistry and Physics (CRC Press, 2000), Vol. 8.
84.
J. E.
Mayer
,
Equilibrium Statistical Mechanics
(
Elsevier
,
2013
), Vol. 1.
85.
P.
Peter Atkins
and
J.
De Paula
,
Atkins’ Physical Chemistry
(
OUP Oxford
,
2014
).
86.
J.
Fischer
,
R. R.
Gamache
,
A.
Goldman
,
L. S.
Rothman
, and
A.
Perrin
, “
Total internal partition sums for molecular species in the 2000 edition of the HITRAN database
,”
J. Quant. Spectrosc. Radiat. Transfer
82
,
401
412
(
2003
).
87.
L.
Chen
,
H.
Shao
,
X.
Zhou
,
G.
Liu
,
J.
Jiang
, and
Z.
Liu
, “
Water-mediated cation intercalation of open-framework indium hexacyanoferrate with high voltage and fast kinetics
,”
Nat. Commun.
7
,
11982
(
2016
).
88.
C.
Wei
,
X.-Y.
Fu
,
L.-L.
Zhang
,
J.
Liu
,
P.-P.
Sun
,
L.
Gao
,
K.-J.
Chang
, and
X.-L.
Yang
, “
Structural regulated nickel hexacyanoferrate with superior sodium storage performance by k-doping
,”
Chem. Eng. J.
421
,
127760
(
2021
).
89.
J. C.
Wojdeł
and
S. T.
Bromley
, “
Efficient calculation of the structural and electronic properties of mixed valence materials: Application to prussian blue analogues
,”
Chem. Phys. Lett.
397
,
154
159
(
2004
).
90.
R.
Englman
and
R.
Englman
,
The Jahn-Teller Effect in Molecules and Crystals
(
Wiley-Interscience
,
New York
,
1972
).
91.
D.
Meggiolaro
,
F.
Ambrosio
,
E.
Mosconi
,
A.
Mahata
, and
F.
De Angelis
, “
Polarons in metal halide perovskites
,”
Adv. Energy Mater.
10
,
1902748
(
2020
).
92.
Z.
Fang
,
Y.
Yin
,
X.
Qiu
,
L.
Zhu
,
X.
Dong
,
Y.
Wang
, and
Y.
Xia
, “
Prussian blue cathode with intercalation pseudocapacitive behavior for low-temperature batteries
,”
Adv. Energy Sustain. Res.
2
,
2100105
(
2021
).
93.
J.
Towns
,
T.
Cockerill
,
M.
Dahan
,
I.
Foster
,
K.
Gaither
,
A.
Grimshaw
,
V.
Hazlewood
,
S.
Lathrop
,
D.
Lifka
,
G. D.
Peterson
et al., “
XSEDE: Accelerating scientific discovery
,”
Comput. Sci. Eng.
16
,
62
74
(
2014
).

Supplementary Material