Intercalated Cation Disorder in Prussian Blue Analogues: First

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

C: Energy Conversion and Storage; Energy and Charge Transport

Intercalated Cation Disorder in Prussian Blue Analogues: First Principles and Grand Canonical Analyses Sizhe Liu, and Kyle Christopher Smith J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b12455 • Publication Date (Web): 22 Mar 2019 Downloaded from http://pubs.acs.org on March 22, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Intercalated Cation Disorder in Prussian Blue Analogues: First Principles and Grand Canonical Analyses

Sizhe Liua and Kyle C. Smitha,b,c,d,* a Department

b Department

of Materials Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

c Computational

d Beckman

of Mechanical Science and Engineering, University of Illinois at UrbanaChampaign, Urbana, IL 61801, USA

Science and Engineering Program, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

Institute for Advanced Science and Technology, University of Illinois at UrbanaChampaign, Urbana, IL 61801, USA *corresponding author: [email protected]

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 52

Abstract Intercalated cation disorder in Prussian Blue analogues (PBAs) results from cation intercalation among the vacant interstitial sites of the PBA host.

In this study we

investigate the influence of intercalated-cation ordering on the electrochemical equilibrium of PBAs. First principles calculations were performed for supercell configurations with two distinct arrangements of cations, vacancy-pair configurations (VPCs) and vacancy-pair free configurations (VPFCs), to analyze the energetic and electronic structure differences between them. VPFCs were found to be energetically favorable. The energetic cost of vacancy pair formation was also attributed to charge sharing between cation and surrounding atoms. Using these findings, a vacancy-pair free (VPF) model was then constructed based on grand canonical ensemble (GCE) theory to predict the composition-dependent equilibrium potential for intercalation into sodium nickel hexacyanoferrate 𝑁𝑎1 + 𝑥𝑁𝑖𝐹𝑒(𝐶𝑁)6. Close agreement with experiment indicates the critical role of cation ordering in PBAs at equilibrium. In this GCE analysis a hybrid Markov method (HMM) was used to sample VPFCs at various 𝑥, showing two distinct patterns of cation/vacancy ordering. Finally, the VPF model was compared to a noninteracting solution (NIS). The failure of the NIS to predict potential at low degree of intercalation further confirms that cation/vacancy ordering is essential to electrochemical equilibria in PBAs.


2

Page 3 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Prussian blue analogues (PBAs) are receiving increasing interest as electroactive materials due to their low cost, high rate capability, long cycle life, and electrochemical activity toward intercalation of various cations.1–5 These compounds have a general formula of 𝐴𝑥𝑀𝑎𝑦[𝑀𝑏(𝐶𝑁)6]𝑧 ∙ 𝑛𝐻2𝑂, with 𝑀𝑎,𝑏 being a transition metal and 𝐴 being an intercalated alkali or alkaline-earth cation.2 Here, we use MaMb-PBA to denote a given PBA. In the PBA lattice, cyanide ligands are bound to transition metals to form a facecentered cubic phase, while interstitial octahedral sites can be vacant or occupied by intercalated cations.5 During the intercalation and deintercalation of cations, diffusion within the solid PBA host results in dynamic sampling of statistical microstates having distribution of intercalated cations among various sites. The large octahedral sites in PBAs make them a competitive intercalation material for applications in energy storage and desalination that require fast cycling and long life. PBAs have been used to reversibly intercalate cations in Li, Na, K and Ca-ion batteries.6–10 These materials are also potential host

compounds

for

electrochemical

desalination

devices

using

intercalation

electrodes.11–16 The site occupancy statistics of intercalated cations in PBAs have received little attention to date because cation intercalation in many PBAs causes very little lattice strain4,17,18 and their kinetics of cation intercalation are known to be facile.9 Thus, the prevailing assumption in the literature is that the specific arrangement of cations and vacancies does not influence the equilibrium and dynamic behavior of PBAs significantly


3


Page 4 of 52

during intercalation/de-intercalation reactions. However, recent magic angle spinning (MAS) NMR spectroscopy experiments suggest that the arrangement of intercalated cations is restricted by the local distribution of charge during intercalation in PBAs.19 These results imply that cation diffusion in PBAs is unlikely to follow a random walk, and atomic-scale understanding of cation/vacancy disorder mechanisms is critical to understanding how diffusion occurs in these systems. Furthermore, intercalated cation ordering phenomena have been investigated in the past for intercalation compounds with layered structures.20–23 During electrochemical cycling of these materials cation ordering is not preserved, and cation disorder induces severe lattice strain20,21 leading to fracture and capacity fade.20 Because the diffusion of guest cations introduces lattice strain in PBAs,4,24,25 capacity fade can be a problem for certain PBAs.26–28 Thus, a detailed investigation of vacancy/cation disorder within PBAs could aid in the design of intercalation materials with improved rate capability and resilience. Furthermore, changes in the magnetism and the lattice constants of PBAs are known to depend on the relative proximity of adjacent vacancies in the host because the re-arrangement of intercalated cations is coupled with electron transfer29 and lattice deformation.4 The CuFe-PBA switches from a ferromagnet to a paramagnet because disorder in the arrangement of intercalated cations alters ferromagnetic interactions between neighboring Fe atoms.30 A similar magnetic transformation was later observed in Prussian Blue (FeFe-PBA) with high Na content.31 The lattice of CuFe-PBA also expands from a cubic phase to a tetragonal phase at a critical degree of intercalation due to the disordered distribution of intercalated Li+ ions.30 In other PBAs, such as FeMnPBA32 and EuFe-PBA,33 lattice deformation induced by cation insertion has been studied


4

Page 5 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


using X-ray diffraction (XRD) experiments and DFT calculations. The results showed that layered cation ordering causes cooperative rotation of Fe-C6 and metal-N6 groups. Given the influence of vacancy/cation disorder on cation diffusion and lattice stability, the properties of PBAs have been tuned for specific applications by varying the concentrations of intercalated cations.34,35 PBAs with tuned cation disorder have shown increased stability36 and better separation capacity34 in these applications. The tunability of PBA properties motivates the development of understanding concerning cation disorder within them. The equilibrium potential for cation intercalation 𝜙𝑒𝑞is an important reference point of the equilibrium and dynamic behavior of electrodes, and its variation with the composition of intercalated ions significantly influences energy efficiency11 and cell voltage12 of electrochemical devices using PBAs as electrodes. At the cell level the equilibrium potential for intercalation is the working electrode potential (in V) relative to a specified reference at thermodynamic equilibrium.

The equilibrium potential is

determined by the chemical potential of cation/electron polarons (in eV/atom) intercalated within a solid host (in addition to the chemical potentials of solution-phase ions). Thus, the equilibrium potential is a measure of the incremental change in Gibbs free energy upon intercalation of cation/electron polarons into the solid host. For PBAs theory and experiment have been combined to build thermodynamic relations for 𝜙𝑒𝑞 in the case of mixed-cation intercalation.37 This work showed that the equilibrium potential during a charge/discharge cycle increases on average with increasing mole fraction of the larger hydrated cation in the electrolyte containing two cations. Closed form equations have also been used to fit experimental data for the composition dependent equilibrium


5


Page 6 of 52

potential of PBAs using both regular solution theory,7,11,13 the Temkin isotherm,16 and the Frumkin isotherm.38 The composition dependent equilibrium potential is an essential input for porous-electrode modeling of batteries and desalination cells using PBAs.7,11,13,38 However, such semi-empirical models lack a rigorous atomistic description of the interactions between intercalated species and their statistical mechanics.

Also, the

temperature dependence of equilibrium potential can be used to harvest low-grade heat (i.e., < 130℃) as electrical energy using thermally regenerative electrochemical cycles (TRECs) with PBA electrodes.39,40 Electrode materials with heavily temperaturedependent 𝜙𝑒𝑞 must be developed to improve the harvesting efficiency of TREC.39 First-principles modeling, such as density functional theory (DFT), can be an effective tool for the analysis of order/disorder phenomena in PBA materials on the basis of their electronic ground states.8,41,42 DFT predictions of the electronic density of states (DOS) have been used to investigate the interatomic interactions in PBAs.8 Phase transitions in PBAs had been predicted based on the variation of DFT-calculated energies around critical intercalation fractions.43 Periodic boundary conditions are used routinely in DFT to simulate infinitely sized systems, but the size of periodic supercells used in practice is limited by computational time and memory. Thus, in isolation DFT calculations are insufficient to investigate the long-range correlations that are essential to describing disordered vacancy/cation arrangements. A given DFT calculation is performed with fixed number of each atomic species, and, as such, the average equilibrium potential at zero temperature can be estimated from DFT by calculating the difference in electronic ground state energy between configurations with different composition.43 While this approach is useful for high-throughput screening of intercalation materials with disparate average


6

Page 7 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


equilibrium potential (e.g., in Ref. 29), for a given intercalation material very large supercells are required to resolve the variation of equilibrium potential over small differences in composition. Furthermore, at finite temperature configurational entropy contributions to the equilibrium potential must be included, in addition to ground state energies.

Thus, complementary analysis is therefore needed in order to study the

disordered arrangement of cations and vacancies in PBAs. A variety of theoretical methods have been used in the past to study the disorder of cations within host compounds, including mean-field theory (MFT),44 Monte Carlo methods (MC),45,46 and grand canonical ensemble theory (GCE).47 MFT approximates many-body correlations as a product of one-body expectation values, from which the expressions for thermodynamic properties, such as entropy and free energy, are derived as functions of the average occupation of intercalated species.48 By fitting equilibrium potential data using phenomenological interactions, MFT has been used to deduce longrange interactions between intercalated atoms.44 However, MFT often fails to predict thermodynamic property scaling near critical points, and it shifts critical points from their true values.48 The thermodynamic properties calculated from MC simulations, on the other hand, agree well with experimental data in the vicinity of critical points.45,48 The ansatz of typical MC simulations is an Ising-like expression for interactions between intercalated species.49 The effective atom-atom interaction parameters used in MC have also been derived from first-principles calculations and cluster expansion techniques50 to precisely capture correlated many-body interactions.

MC methods can suffer from

statistical noise as a result of the stochastic approach used to solve deterministic problems.51 Therefore, many microstates of each intercalated host configuration must be


7


Page 8 of 52

sampled to obtain converged energy. In contrast GCE theory can be used to investigate the distribution of defects in a finite-sized supercell with fixed chemical potential 𝜇, volume 𝑉, and temperature 𝑇.47 The efficiency of GCE theory can be improved further by binning energetically equivalent configurations on the basis of symmetry.52 Accordingly, we use grand canonical ensemble theory presently to predict the equilibrium potential for cation intercalation in PBAs. In this study, we perform first-principles calculations on NiFe-PBA supercells to show that the arrangement of intercalated cations affects the electronic interactions between intercalated Na+ and other atoms in the PBA host. By analyzing the relative stability of two kinds of vacancy ordering, we show that atomic configurations that are free of nearest-neighbor interstitial vacancy pairs are favored thermodynamically.

We

subsequently use this information to propose a vacancy-pair free (VPF) model, in which only certain types of vacancy/cation configurations are considered. The resulting equilibrium potential curve predicted using the VPF model shows agreement with experimental data, in contrast with a non-interacting model.

2. Computational Methods 2.1 Density Functional Theory Simulations The DFT calculations were performed with collinear-polarized spin using the Quantum ESPRESSO (QE) package53 and correlation interactions were modeled using the generalized gradient approximation with the Perdew-Burke-Ernzerhof (PBE) functional and Hubbard corrections (GGA+U).54-56 The 𝑈 values of 3 eV for Ni and 1 eV


8

Page 9 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for Fe were taken from literature8,32,41,42,57 to stabilize the high-spin/low-spin ground state.58 The valence electrons for Ni, Fe, and Na atoms corresponded to 3d94s1, 3d74s1, and 2s12p63s1, respectively, while remaining electrons were kept frozen as core states. Ultrasoft pseudopotentials59 were generated with a scalar-relativistic calculation to describe the interaction between valence and frozen electrons. An energy cutoff of 476 eV was used in all cases for the truncation of plane-wave basis sets. Monkhorst-Pack 𝑘points were used in all calculations with a 4 × 4 × 4 grid for the 2 × 2 × 2 supercell, producing a real space cutoff 60 of 0.3 Å . In all simulations relaxation of atomic positions was performed, while variable-cell relaxation was performed only for lowest energy configurations with certain vacancy number 𝑚.

Cell parameters obtained from those

variable-cell simulations were then used to predict the ground state energy of all other configurations with the same 𝑚. For all SCF calculations ground state energy converged to within 1.36 × 10 ―4 eV. This convergence was achieved with 1.36 × 10 ―2 eV/a.u. for L2norm of the forces on all atoms and maximum lattice stress of 0.7 kbar. For certain initialized configurations that were especially unstable, the atom mass of C and N atoms was doubled to dampen the fast movement of cyanide ligands33 during the relaxation calculations. To calculate the mixing energies of the supercell with various cation arrangements, 12 symmetrically inequivalent configurations were investigated. For each configuration, a set of calculations with different magnetic states were performed performed by initializing the magnetization numbers at Fe and Ni atoms. This approach enabled us to find the lowest energy spin states of each configuration, even if the initial guess for the magnetization was far from the electronic ground state. Ferromagnetic ordering was assumed in all the calculation because it produced lower total energy than


9


Page 10 of 52

the anti-ferromagnetic ordering. The number of unpaired electrons on transition metal atoms were estimated from Bader analysis.61 The optimized lattice parameter 𝑎 of a fully intercalated and a fully de-intercalated 2 × 2 × 2 supercell were 10.16Å and 10.12Å, which was in good agreement with experimental value of 𝑎 = 10.20Å.58 We also found that the intercalated Na+ cations occupied at the body center of certain unit cells after relaxation. While large intercalated cations prefer body-center occupancy,8,42 bare cations with radius smaller than bare Na+ can stably occupy certain off-center locations.42 Because the correlation between the offcenter cation occupancy and the distortion of the PBA framework62 requires deeper investigation, the vacancy/cation ordering phenomena studied presently are particular to the body-center occupancy of cations. The electronic charge density distributions predicted from SCF calculations were visualized using the VESTA program.63 The QE graphical user interface BURAI64 was also used to perform projected density of states (PDOS) calculations with a 4 × 4 × 4 Monkhorst−Pack 𝑘-point mesh.


10

Page 11 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Crystal structures of sodium nickel hexacyanoferrate ([Na1.5NiFe(CN)6]4) with a vacancy pair configuration (a) and vacancy-pair free configuration (b). These configurations are also shown in simplified representations with occupied (orange) and vacant (transparent) sites in (c) and (d), respectively.

2.2 Vacancy Pair Free Model The vacancy pair free (VPF) model was used to construct a grand canonical ensemble (GCE) to predict the equilibrium potential for cation intercalation into a PBA host. The results of DFT calculations suggest that the vacancy pair free configurations (VPFCs), defined as the cation/vacancy arrangement without nearest-neighbor vacancy pairs (or NNVPs, see Figure 1b), are more accessible at room temperature than the vacancy pair configurations (VPCs) having NNVPs (e.g., Figure 2a). Vacancies in the present framework materials are defined as unoccupied interstitial lattice sites. Inspired


11


Page 12 of 52

by this, we assume no VPC would appear in GCE at equilibria and consider a simple interaction model only for VPFCs where the energy of a microstate with certain vacancy number 𝑚 is 𝐸𝑚 = 𝑚 × 𝜖,

(1)

where 𝜖 is the energy required to oxidize the PBA framework during a de-intercalation reaction. The present VPF model neglects volume changes among microstates, such that the energy levels of all VPFCs with the same vacancy number are degenerate. The grand partition function Ξ is thereby given as:47,65 𝑁/2

Ξ=

∑Ω

(

𝑚𝑒𝑥𝑝

𝑚=0

(𝑁 ― 𝑚)𝜇 ― 𝐸𝑚 𝑘𝐵𝑇

)

(2).

The summation here is taken over VPFCs with vacancy number ranging from 𝑚 = 0 to 𝑚 = 𝑁/2 (i.e., 𝑥 = 0 and 𝑥 = 1 in Na1+xNiFe(CN)6), where 𝑁 is the total number of the sites in the PBA framework, Ω𝑚 is the degeneracy of the energy level for 𝑚-vacancy VPFCs, 𝜇 is the chemical potential of Na+/electron polarons in the solid, and 𝑘𝐵𝑇 is the thermal energy scale. The Na+/electron polaron chemical potential 𝜇 is related to the equilibrium reduction potential 𝜙𝑒𝑞 versus an Ag/AgCl reference electrode with negligible liquid junction potential as: 𝜇 = ―𝑒𝜙𝑒𝑞 + 𝑘𝐵𝑇𝑙𝑛(𝑎𝑁𝑎 + 𝑎𝐶𝑙 ― ) + 𝜇0𝑁𝑎 + + 𝜇0𝐴𝑔 + 𝜇0𝐶𝑙 ― ― 𝜇0𝐴𝑔𝐶𝑙

(3).

Here, 𝑎𝑖 and 𝜇0𝑖 are the solution-phase activity and standard-state chemical potential of species 𝑖. 𝑒 is the elementary charge. While our GCE is posed in terms of the partition function Ξ, we note that the GCE can be posed equivalently in terms of the grand potential. While we do not invoke grand


12

Page 13 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


potential explicitly, it can be expressed in terms of the partition function as ― 𝑘𝐵𝑇𝑙𝑛Ξ. A relationship between ensemble-averaged quantities and that of their microstates in the GCE are predicted using the expectation theorem.

Most important to the present

investigation, for certain 𝜇 the expected value of vacancy number 〈𝑚〉 among the various microstates is: 1 〈𝑚〉 = Ξ

𝑁/2

∑ 𝑚Ω

(

𝑚𝑒𝑥𝑝

𝑚=0

(𝑁 ― 𝑚)𝜇 ― 𝐸𝑚 𝑘𝐵𝑇

The expected intercalation fraction 〈𝑥〉 is 〈𝑥〉 = 1 ―

2〈𝑚〉 𝑁

)

(4).

, and, thus, the relationship

between 𝜙𝑒𝑞 and 〈𝑥〉 can be constructed by calculating 〈𝑚〉 for a discrete set of 𝜇 values. For supercells larger than 3 × 3 × 3 the calculation of Ω𝑚 can be done neither using combinatorial methods66 nor using generalized disorder models.67,68 To reduce computational expense degeneracy has been determined in the past by identifying symmetry among disordered configurations with certain site occupancy,52,69 but significant computational effort is still required to enumerate all inequivalent configurations in large supercells. To reduce computational complexity and to avoid combinatorial explosion70,71 in the direct calculation of Ω𝑚, instead we use a recurrence relationship between the degeneracy of VPFCs with successive vacancy numbers, Ω𝑚 and Ω𝑚 + 1, expressed in terms of an average “guard number” 〈〈𝑔𝑚〉〉:

Ω𝑚 + 1 =

Ω𝑚(𝑁 ― 𝑚 × 〈〈𝑔𝑚〉〉) 𝑚+1


(5)

13


Page 14 of 52

In mathematical terms 〈〈𝑔𝑚〉〉 ―1 is equal to the average number of occupied nearest neighbors (ONN) to a vacancy, also averaged over all 𝑚-vacancy VPFCs. In qualitative terms each vacant site in VPFCs has, on average, 〈〈𝑔𝑚〉〉 ―1 nearest-neighbor cations for which incremental vacancy substitution would result in vacancy-pair formation. Thus, the factor of 𝑁 ― 𝑚 × 〈〈𝑔𝑚〉〉 in the numerator is equal to the number of cations for which incremental vacancy substitution can be performed without forming vacancy pairs. The denominator 𝑚 + 1 is used to avoid double counting. Once 〈〈𝑔𝑚〉〉 is known for all integers 𝑚 between 1 and 𝑁/2 ― 1, Eq. 5 can be used as a recurrence relation to determine the corresponding degeneracies for all 𝑚. According to the definition of the guard number, the calculated value of 〈〈𝑔𝑚〉〉 should vary within the range of [2,7], as described subsequently. If vacancies in a VPFC are sparsely distributed and none of their ONNs overlaps, the value of 〈〈𝑔𝑚→1〉〉 = 7 because each vacancy has 6 nearest neighbors (NNs) in the 3D framework. As the vacancy number increases, the ONNs of vacancies begin to overlap and 〈〈𝑔𝑚〉〉 monotonically decreases as a result. At the terminal configuration (TC, 𝑚 = 𝑁/2), it is no longer possible to introduce new vacancies without producing a vacancy pair. This condition requires that 𝑁

𝑁 ― 2 〈〈𝑔𝑚→𝑁/2〉〉 = 0 and consequently that 〈〈𝑔𝑚→𝑁/2〉〉 = 2. For vacancy numbers other than 𝑚 = 1 and 𝑚 = 𝑁/2 we estimate 〈〈𝑔𝑚〉〉 using a hybrid Markov method (HMM, see Sec. 2.3) to sample configurational space. Because different VPFCs with certain 𝑚 may have different guard number, the HMM is used to generate 𝑁𝑠𝑝𝑚𝑙 distinct VPFCs to calculate 〈〈𝑔𝑚〉〉 as:


14

Page 15 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝑁

〈〈𝑔𝑚〉〉 =

∑𝑖 =𝑠𝑚𝑝𝑙 〈𝑔𝑚,𝑖〉 1 𝑁𝑠𝑝𝑚𝑙

(6),

where 〈𝑔𝑚,𝑖〉 is the guard number of the 𝑖𝑡ℎ sample calculated from the site occupancy of the sample of interest and averaged only over that sample’s vacancies. Accurate estimation of 〈〈𝑔𝑚〉〉 requires that the subset of VPFCs sampled for certain 𝑚 be statistically representative of the entire set. Accordingly, we use an HMM to create VPFC samples randomly using a single adjustable parameter, the smearing coefficient 𝑠, chosen accordingly to match the known degeneracy at 𝑚 = 𝑁/2 (Ω𝑚 = 𝑁/2 = 2).

Presently we use 𝑁𝑠𝑝𝑚𝑙 = 500 VPFC samples for each integer value of 𝑚.

Degeneracy values are then recovered from 〈〈𝑔𝑚〉〉 for ensembles with microstate supercell sizes of 4 × 4 × 4, 6 × 6 × 6 , and 8 × 8 × 8 to test the convergence of final results. 2.3 Hybrid Markov Method We use a hybrid Markov method (HMM) to facilitate efficient calculation of the degeneracy of VPFCs for subsequent analysis based on VPF model. The specific method used here combines a Markov process with random vacant-site sampling to sample the configurational space of VPFCs. For a set of 𝑚-vacancy VPFCs that are sampled by HMM, the average “guard number” 〈〈𝑔𝑚〉〉 in Eq. 5 is calculated and the corresponding degeneracy of VPFCs is then estimated. For a certain VPFC realization with 𝑁 sites, 𝑚 × (〈〈𝑔𝑚〉〉 ― 1) is equal to the number of sites around all vacancies in the VPFC that are excluded from incremental vacancy insertion, so as to preserve the vacancy-pair free condition.


15


Page 16 of 52

For the present HMM sampling process we assume that a certain supercell realization starts at a fully intercalated state, where all sites are labeled by using binary sequences consisting of 1 and 2 (Figure 2a). A binary sequence (12121212) has a primitive period of 2, so that there are 2 symmetry types of the sequences obtained by cyclic permutation. As a result, the sites labeled by “1” (or type-1 sites) are surrounded by type-2 sites and vice versa. A VPFC in which 𝑁/2 vacant sites are at either type-1 or type-2 sites is called a “terminal configuration” (TC). We note that, in combinatorial nomenclature, the TC possesses “complete chessboard ordering”.72 Accordingly, the average guard number for TCs 〈〈𝑔𝑁/2〉〉 = 2 because no occupied sites are spared from 𝑁

guarding their vacant neighbors (i.e., 𝑁 ― 〈𝑔𝑁/2〉 × 2 = 0). On the other hand, the guard number in the dilute vacancy limit is 〈〈𝑔1〉〉 = 7 for 3D lattices because only one is present in this limit, and that vacancy possesses six ONNs. For 2D lattices this limit produces

〈〈𝑔1〉〉 = 5.


16

Page 17 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Schematic of a hybrid Markov process that produces a seven-vacancy VPFC on an 8 × 8 2D lattice. A smearing coefficient of 𝑠 = 0.5 and a random number 𝑅 = 0.4 were chosen to illustrate this process.

The process starts from a fully intercalated

supercell in (a), then 𝑛𝑎𝑐𝑡 = ⌈7 × (0.4 × (1 ― 0.5) + 0.5)⌉ = 5 “active” type-1 sites are selected and made vacant (b). HMM then loops over all the ONNs of all the active vacancies to perform Metropolis algorithm.73 Three of the five vacancies move to new locations in (c) after one loop is complete. The new ONNs (shown as grey squares) are then labeled. During the static step 𝑛𝑖𝑛𝑎𝑐𝑡 = 2 “inactive” sites are randomly selected from the remaining orange squares in (d) and (e), and their neighbors are labeled in green. At the end, the number of all the sites other than the orange ones is divided by 𝑚 to produce the value of 〈𝑔7〉 for the generated VPFC shown in (f).


17


Page 18 of 52

After labeling sites, the HMM introduces a “smearing coefficient” 𝑠 to decide the number of active 𝑛𝑎𝑐𝑡 and inactive vacancies 𝑛𝑖𝑛𝑎𝑐𝑡 that will be permitted during a given HMM step: 𝑛𝑎𝑐𝑡 = ⌈𝑚 × (𝑅 × (1 ― 𝑠) + 𝑠)⌉, 𝑛𝑖𝑛𝑎𝑐𝑡 = 𝑚 ― 𝑛𝑎𝑐𝑡

(7)

Here, 𝑠 takes a fixed value between zero and one that is ultimately chosen to satisfy the degeneracy of TCs (i.e., Ω𝑁/2 = 2), and ⌈.⌉ is the ceil operator. 𝑅 is a random number within a range of (0,1). With 𝑠 = 0.5 and 𝑅 = 0.4, we have 𝑛𝑎𝑐𝑡 = 5 and 𝑛𝑖𝑛𝑎𝑐𝑡 = 2 in the example shown in Figure 2. Splitting the vacancies into two active and inactive sets differentiates the Markov process used in the present HMM from random sampling processes, and it enables exhaustive sampling of the VPFC configurational space. All active vacancies are first introduced at either type-1 or type-2 sites (Figure 2b). During the Markov process that follows the insertion of vacancies, the HMM loops over all the ONNs of the vacancies to perform the Metropolis algorithm.73 The transition possibility matrix for each swap between an ONN with its neighboring vacancy is constructed as follows. First, the underlying stochastic matrix 𝒂 for a swap event 𝑠𝑖 of ONN “𝑖” is defined as:

{

1 𝑠 ∈𝑪 𝑎𝑘𝑙 = 2 𝑖 0 𝑠𝑖 ∉ 𝑪

(8)

Where 𝑪 is a collection of swap events taking the supercell system from microstate 𝑘 into any one of its neighboring microstates 𝑙 without forming a NNVP. The value of 1/2 derives from the fact that each swap has two outcomes: ONN “𝑖” either becomes (1) vacant or (2) unchanged. Because the outcomes belong to a finite set of microstates including a


18

Page 19 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


swapped microstate 𝑙 and an un-swapped microstate 𝑘, the transition probability 𝜏𝑘𝑙 from state 𝑘 to state 𝑙 is equal to 𝑎𝑘𝑙exp ( ― ∆𝐸𝑘𝑙/𝑘𝑇). For the present VPF model ∆𝐸𝑘𝑙 = 0 as long as the number of vacancies remains constant, resulting in 𝜏𝑘𝑙 = 𝑎𝑘𝑙 for 𝑘 ≠ 𝑙 and 𝜏𝑘𝑘 = 1 ― 𝑎𝑘𝑙. Therefore, the possibility of producing a new microstate 𝑝(𝑥) 𝑛𝑒𝑤 and that of (𝒙) keeping the old state 𝑝(𝑥) 𝑜𝑙𝑑 at the step 𝑥 comprise the vector distribution 𝒑 . The elements

of the vector 𝒑(𝒙) depend entirely on the result of step 𝑥 ― 1 (i.e., 𝒑(𝒙) = 𝒑(𝒙 ― 𝟏)𝝉), as required for a Markov process. Figure 2c shows a VPFC after a loop over all the ONNs in Figure 2b performing Markov process while the grey sites are the new ONNs. While the Markov process can be implemented with ease, we find that such a process requires a large number of computational loops to obtain a converged value for

〈〈𝑔𝑚〉〉 when VPFCs with large vacancy numbers are sampled. This effect results from the fact that the initial arrangement of vacancies introduced in the HMM already has cha essboard order of cations, making the sampling process less random. To break the preexisting cation ordering, the inactive vacancies is introduced by choosing the vacant sites randomly. In the example of Figure 2d and 2e, two vacancies are introduced without the rearrangement of ONNs. The ONNs of inactive vacancies are labeled by green squares. In Figure 2f, the guard number 〈𝑔7〉 is calculated as

𝑁𝑤ℎ𝑖𝑡𝑒 + 𝑁𝑔𝑟𝑒𝑒𝑛 + 𝑁𝑔𝑟𝑒𝑦 7

with 𝑁𝑤ℎ𝑖𝑡𝑒,

𝑁𝑔𝑟𝑒𝑒𝑛, and 𝑁𝑔𝑟𝑒𝑦 being the number of white, green, and grey squares in Figure 2e. Figure 2f gives the final configuration of the VPFC sample. It is noteworthy that the smearing coefficient 𝑠 influences the ratio of the type-1 and type-2 vacancies in the final sampled configuration, which produces the distinct vacancy/cation arrangement patterns that are described in Sec. 3.3.


19


Page 20 of 52

Figure 3. The quotient (𝑁 ― 𝑚 × 〈𝑔𝑚〉) (𝑚 + 1) calculated based on the HMM with various values of 𝑠 (dots) is compared to the exact solution based on direct enumeration (red curve) for a 4 × 8 lattice plane. This method was validated using a two-dimensional (2D) 4 × 8 lattice. We use a 4 × 8 lattice for validation because (1) the exact degeneracy of VPFCs can be calculated

by brute-force enumeration of all vacancy configurations, and (2) this lattice is large enough to sample a sufficient number of VPFCs for 〈〈𝑔𝑚〉〉 estimation. The result of this validation process is shown in Figure 3. From Eq. 1, we recognized that HMM becomes a simple Monte Carlo method (MCM) when the smearing factor 𝑠 = 1. However, with a limited number of ONN swap, the results produced by HMM (red dots) is closer to the exact solution than MCM (blue dots). Because performing a lengthy Markov process is not computational efficient for 3D supercell, the HMM provides a better 〈〈𝑔𝑚〉〉 estimations without overly repeating the Markov process. In current study, HMM loops over all the ONNs only twice, and the largest deviation of HMM results from the exact solution in


20

Page 21 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3 is 0.08, which implies at most 8.32% error in the calculation of Ω𝑚 by using the HMM. Moreover, the best choice we found for 𝑠 is 0.77 for other values tend to underestimate or overestimate the ratio when the vacancy number 𝑚 > 𝑁/4.

3. Results and Discussion 3.1 Energetic Difference between VPC and VPFC DFT calculations were performed on a vacancy-pair configuration (VPC) (Figure 1a) and on a vacancy-pair free configuration (VPFC) (Figure 1b) with [Na1.5NiFe(CN)6]4 stoichiometry. The supercells simulated here are set to have two vacant sites (i.e., six intercalated sites) to allow at most one vacancy pair (or one vacant chain due to periodic images of vacancies). Six Na+ ions among eight available sites produces a total of 28 distinct arrangements in the supercell with 16 of which being VPFCs and 12 of which being VPCs. The VPC and VPFC shown in Figs. 1a and 1b, respectively, possess the smallest energy difference among all vacancy/cation arrangements, and their vacancy occupancies differ only by one vacancy swapping event. By analyzing the energy difference between these two configurations, we show that the electronic structure of the supercells changes even with a small change in the cation/vacancy arrangement. Table 1. The difference of the total energy and the energy contributions between the particular VPFC (Figure 1a) and VPC (Figure 1b) of interest. Energy Contribution 𝑬𝑽𝑷𝑭𝑪 ― 𝑬𝑽𝑷𝑪 (𝒆𝑽)

One-electron

Electrostatic

157.44

-157.23


Exchange Correlation -0.28

Total 0.08

21


Page 22 of 52

The DFT implementation of QE calculates the contributions of various electronic interactions to the total energy.53 The difference in energy between VPCs and VPFCs was calculated for each contribution to the total energy (Table 1). The electrostatic energy captures the classic Coulombic interactions by using Ewald summation74,75 and the functional of electron density in the Kohn-Sham equations.76 The one-electron energy contribution includes the effects of kinetic energy calculated from the Laplacian of wavefunctions.77 From Table 1, we observe that one vacancy pair in the VPC raises total energy by 40.5meV (i.e., 81 meV/2) relative to the VPFC, because there are two vacancy pairs in Figure 1a when periodic vacancy images are considered. This energy difference is solely attributable to the different arrangement of vacancies in the host. On the one hand, the separation of the two vacancies in the VPFC causes a decrease in electrostatic potential and exchange-correlation potential. On the other hand, the separated vacancies also introduce extra non-interacting kinetic energy of electrons77,78 (i.e., one-electron energy contribution calculated in QE). As a result, the decrease in the classic Coulombic interactions and the approximated exchange-correlation interaction compensates for the increase of the electron kinetic energy, making the VPFC energetically favorable. Beyond the lowest-energy VPCs and VPFCs the relative stability of all remaining configurations in the 2 × 2 × 2 supercell of Na1+xNiFe(CN)6 was also determined. In this analysis the number of Na atoms was varied from four (𝑥 = 0) to eight (𝑥 = 1). A total of 12 configurations were considered that includes all the symmetrically inequivalent52 vacancy/cation arrangements. The corresponding mixing energy per formula unit ∆𝑚𝑖𝑥𝐸 of a given configuration is:


22

Page 23 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∆𝑚𝑖𝑥𝐸 = 𝐸[𝑁𝑎1 + 𝑥𝑁𝑖𝐹𝑒(𝐶𝑁)6] ― (1 ― 𝑥)𝐸[𝑁𝑎𝑁𝑖𝐹𝑒(𝐶𝑁)6] ― 𝑥𝐸[𝑁𝑎2𝑁𝑖𝐹𝑒(𝐶𝑁)6]

(9),

where 𝐸[𝑁𝑎1 + 𝑥𝑁𝑖𝐹𝑒(𝐶𝑁)6] is the energy of a configuration with an intermediate intercalation degree 𝑥 per formula unit. 𝐸[𝑁𝑎𝑁𝑖𝐹𝑒(𝐶𝑁)6] and 𝐸[𝑁𝑎2𝑁𝑖𝐻𝐹𝑒(𝐶𝑁)6] are the energies of the most stable desodiated and sodiated configurations, respectively. Using this definition, negative mixing energy corresponds to stable vacancy/cation configurations with respect to 𝑁𝑎𝑁𝑖𝐹𝑒(𝐶𝑁)6 and 𝑁𝑎2𝑁𝑖𝐹𝑒(𝐶𝑁)6 separation.79 Figure 4 shows all the calculated ∆𝑚𝑖𝑥𝐸𝑖 for each energetically distinct configuration 𝑖 as a function of 𝑥. Energies calculated for VPFCs and the VPCs are labeled by red and blue bars, respectively.

By connecting the most stable configurations in the plot, we construct a

convex hull.43 The energy of the configurations on the convex hull can be interpreted as the free energy of the system at 0 K neglecting configurational entropy.79 A total of 92 configurations have negative formation energy within the convex hull. More importantly, the most stable configuration is found at 𝑥 = 0.50 (i.e., six Na-ion in the cell), which has a ∆𝑚𝑖𝑥𝐸 of -520 meV and is 62.4 meV more stable than the next-lowest energy of VPFC. On the contrary, the energy difference between the VPFC and the VPC is at least 81 meV, suggesting the inaccessibility of the VPCs at equilibrium. Based on the fact that the DFT calculations for the configurations with the same vacancy number converge at the same lattice constants, we conducted a canonical ensemble analysis to calculate the possibility distribution of the configurations at 𝑥=0, 0.25, and 0.5, respectively. The possibility of finding a configuration 𝑖 in a canonical ensemble of four formula-unit PBA configurations with certain 𝑥 is𝑃 𝑖 =

𝛺iexp ( ― 𝑁𝑓𝑢Δ𝑚𝑖𝑥𝐸𝑖/𝑘𝑇) 𝑛 ∑𝑖 𝑡𝑜𝑡𝛺iexp

( ― 𝑁𝑓𝑢Δ𝑚𝑖𝑥𝐸𝑖/𝑘𝑇)

at 𝑇=300K. Here, 𝑛𝑡𝑜𝑡 is the total

number of the inequivalent configurations for certain 𝑥 (e.g., 𝑛𝑡𝑜𝑡 = 3 at 𝑥 = 0.5 in Figure


23


Page 24 of 52

4), and 𝑁𝑓𝑢 is the number of Fe(CN)6 formula units in the supercell (𝑁𝑓𝑢 = 4 here). The degeneracy Ω𝑖 of each energetically equivalent configuration was obtained by identifying symmetrically equivalent structures. We also calculated the maximum number of vacant NNs of a vacancy in all configurations 𝑛𝑁𝑁,𝑚𝑎𝑥 to indicate the features of cation arrangments, and they are listed along with Ω𝑖 and 𝑃 𝑖 in Figure 4. At each 𝑥, 𝑃𝑖 decreases with increasing 𝑛𝑁𝑁,𝑚𝑎𝑥. In other words, a cation arrangement could become more accessible at equilibrium by reducing NNVPs in the configuration. Meanwhile, the configurations with NNVPs (i.e., VPCs) do not have vanishing possibility across all x values. For example, the possibility of finding VPCs in a canonical ensemble at 𝑥 = 0 is 3.5%. Because of the energetic preference toward VPFCs, we later use grand canonical ensemble theory with simplified interactions to predict vacancy/cation occupancy in large supercells (𝑎 > 40Å). However, the exclusion of VPCs in the simplified interactions causes the discrepancy between the theory and the experimental data shown in Sec 3.3.


24

Page 25 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Calculated mixing energies per formula unit as a function of vacancy number 𝑚. The energy levels labeled by the red bars are for VPFCs while the energy levels of VPCs are labeled by the blue bars. The bracket next to each bar lists the value of the degeneracy Ω𝑖, the maximum number of vacant NN of a vacancy 𝑛𝑁𝑁,𝑚𝑎𝑥, and the possibility 𝑃𝑖 of finding each configuration in a canonical ensemble at T=300 K. The starred brackets indicate the configurations reported in Table 1.

The equilibrium potential 𝜙𝑒𝑞 is a measure of the electrochemical work, specifically Gibbs free energy 𝐺,79 done for an infinitesimal degree of charge transfer. Thus, the average equilibrium potential 𝜙𝑒𝑞 versus the standard 𝐴𝑔/𝐴𝑔𝐶𝑙 reference between the compositions 𝑥2 and 𝑥1 is given by:43

𝜙𝑒𝑞 = ―

𝐺[𝑁𝑎1 + 𝑥2𝑁𝑖𝐹𝑒(𝐶𝑁)6] ― 𝐺[𝑁𝑎1 + 𝑥1𝑁𝑖𝐹𝑒(𝐶𝑁)6] ― (𝑥2 ― 𝑥1)𝐺[𝑁𝑎] 𝑥2 ― 𝑥1

+ ∆𝐸𝑜

(10)

with 𝐺 and 𝑥 respectively being the Gibbs free energy and degree of intercalation (DOI). ∆𝐸𝑜 is the potential difference between 𝐴𝑔/𝐴𝑔𝐶𝑙 and 𝑁𝑎 + /𝑁𝑎 standard reference electrodes equal to 4.3𝑉 based on electrochemical series data.80 Following this approach, we determined the average equilibrium potential variation with 𝑥 at 𝑇 = 0 𝐾 using ground state energies for the configurations on the convex hull (Fig. 4) for the Gibbs free energy of partially intercalated phases (Eq. 10). The 0 K free energy of Na metal 𝐺[𝑁𝑎] was determined based on an SCF calculation of a 2 × 1 × 1 body-centered cubic Na supercell. The resulting values of 𝜙𝑒𝑞 were 1.72, 1.32, -0.41, and -0.57 V vs. Ag/AgCl at 𝑥 = 0.25, 0.5, 0.75 and 1.00. These values deviate from experiment7 by as little as 0.92 V


25


Page 26 of 52

(𝑥 = 0.5) to as large as ~1.20 V (𝑥~0). The discrepancy between experiment and the average potentials suggests that the “ground-state only” methodology is inadequate due to the lack of configurational entropy accounted for. To address this issue, we introduce a VPF model in Sec. 2.2 that incorporates all VPFCs accessible at finite temperature. The discrepancy could also be caused by the choice of the 𝑈 parameter81 and the magnetic disorder of Fe and Ni atoms82 at finite temperature. To test the influence of the 𝑈 parameter on the energy difference between VPFC and VPC, we calculated the total energy of the configurations shown in Figure 1a and 1b, with 𝑈 parameters being 4eV for Ni atoms and 2eV for Fe atoms. The new energy difference between the two is 84 meV, similar to the results listed in Table 1, indicating that our calculation is reliable to predict the energy distributions among vacancy configurations at fixed 𝑥. 3.2 Electron Density Distribution and Density of States In order to gain insight into the electronic structure differences between configurations with (VPCs) and without (VPFCs) NNVPs, the electron density distributions (EDDs) were computed for both configurations from Figure 1 on the (400) and (101) crystallographic planes. On these distributions, shown in Figure 5, red regions indicate the core electrons of different nuclei in the supercell. Among all nuclei the EDD around intercalated sodium nuclei is most spherically symmetric. While intercalated cations in PBAs are typically highly ionized,83,84 the weak covalent interactions between Na atoms and surrounding atoms are evidenced by the envelope of the 10 ―3 𝑎.𝑢. ―3 iso-density contour (green) on the (400) plane (Figures 5a and 5c) that shows charge sharing between these atoms. On the (101) plane (Figures 5b and 5d) the same contour shows that the Na-metal interaction is not as obvious as the Na-CN interaction, because ligand


26

Page 27 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cage effects85 suppress the spreading of the electrons from metal atoms outside of 𝐹𝑒𝐶6 and 𝑁𝑖𝑁6 octahedra. The vacant sites in both configurations shown in Figure 5 have electron density less than 10 ―3 𝑎.𝑢. ―3 (blue), suggesting that cation vacancies are practically vacuum, absent of electrons. With the vacancies in the VPFC well separated by electrons shared among Na and its surrounding atoms, no charge sharing is observed between the two vacancies on either crystal plane. As evidenced from Figures 5a and 5b, the large separation between ligands of approximately 5 Å, together with the ligand cage effect, creates a vacancy pair connected by an electron-deficient transition region between the two vacant unit cells. Figure 6 shows the electron density profile along the pink lines labeled in Figures 5a and 5c. From this profile we observe that the sodium ion/vacancy interaction strengthens the electron density distribution by 7.80% at the midpoint (i.e., 𝑧 = 5.08 Å).


27


Page 28 of 52

Figure 5. The calculated iso-density contours on the (400) plane for a VPC (a) show connected vacuum space between vacancies. The electron density distribution around sodium ions in a VPFC (c) interacts with the center cyanide ligand. On the (101) plane, the continuity of vacuum space between vacant sites in a VPC (b) is pronounced. The sodium ions in both a VPC (c) and a VPFC (d) show negligible charge sharing with Fe and Ni atoms. The electron density distribution along the pink lines labeled in (a) and (c) is plotted in Figure 6.


28

Page 29 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. The electron density profile along the pink lines labeled in Figure 5 for the VPFC (red) and the VPC (blue). Bader analysis61 was conducted based on the EDDs of various cation configurations to study the number of unpaired electrons 𝜎 on Fe atoms. The results show that 𝜎 increases from 0 to 0.77 from a fully-intercalated configuration to a fully deintercalated configuration. This finding indicates that Fe-atoms were at their low-spin state as predicted by other research.41 For the VPFC and the VPC presented in Figure 5, 𝜎 has various values on Fe atoms. It was found to be 0.48, 0.48, 0.43 and 0.43 in the VPC, while in the VPFC 𝜎 had values of 0.49, 0.44, 0.44, and 0.42 among different Fe atoms. The Fe atoms with higher 𝜎 (i.e., less reduced) were found near the vacancies while morereduced Fe atoms were distant from the vacancies. As a result, the extent of Fe reduction also depends on the cation/vacancy arrangement.


29


Page 30 of 52

Figure 7a shows the corresponding electronic density of states (DOS) for the VPFC and VPC of interest. By subtracting absolute DOS at a certain energy level relative to the Fermi energy for the VPC from the VPFC, the VPC shows a deficiency of certain electronic states and an abundance of others spread across the entire electron energylevel spectrum. Such a phenomenon is the result of the simultaneous destruction and formation of certain orbital interactions.84,86 The total DOS was also projected onto Na atoms (Figure 7b), cyanide ligands (Figure 7c), and metal atoms (Figure 7d) to investigate which orbital interaction was responsible for the variations shown in Figure 7a. While the DOS profiles projected onto Na atoms for both configurations are nearly identical (Figure 7b), the difference in DOS spectra are primarily observed in their projections onto transition metal atoms and cyanide ligands. This observation further confirms that vacancy pairing induces rearrangement of electrons on the host compound.


30

Page 31 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. The total DOS of the VPFC and the VPC (a) are projected onto sodium atoms (b), cyanide ligands (c), and transition metal atoms (d). The red dots in each subplot represent the DOS difference between the two configurations at various energy levels.

The metallicity of Na1+xNiFe(CN)6 is also observed with considerable electronic states at the Fermi energy87 in the total DOS plot. The calculated Fermi energies 𝜀𝑓 of the


31


Page 32 of 52

VPC and the VPFC are 1.4565, 1.4528eV, respectively. From a quantum-mechanical point of view, the density of electronic states near the Fermi energy increases electronic conduction rates, among other factors.88 Our calculations reveal negligible differences in DOS at 𝜀𝑓 (Figure 7a), and, accordingly, the electronic conductivity of these configurations is likely to be similar in magnitude ceteris paribus. We note however that the arrangement of cations is likely to influence phonon/electron interactions and ultimately the mean-free path of electrons. 3.3 Results Derived from a Vacancy-Pair Free Ensemble at Finite Temperature The results of the preceding first principles calculations suggest that the vacancy pair free configurations (VPFCs) are the preferred type of cation arrangement in PBA lattices at chemical equilibrium. We note that while VPFCs are not possible when more than half of the alkali-cation sites in 𝐴𝑦𝑀𝑎[𝑀𝑏(𝐶𝑁)6] are vacant (where 𝐴 is the alkali cation of interest), cation occupancy less than 50% can be achieved only if the low-potential reduction reaction is accessed, a scenario that is relevant mainly in non-aqueous cycling of PBAs. We therefore focus our attention on stoichiometries for Na1+xNiFe(CN)6 with 0 < 𝑥 < 1. In the following sections, we first use a hybrid Markov method (HMM) to sample VPFCs with various vacancy numbers. Three classes of VPFCs are then found based on their specific cation arrangements. By summarizing the features of each VPFC class, we discuss the potential influence of cation arrangements on intercalation dynamics using chessboard ordering. After we calculate the average guard numbers for VPFCs required for estimating configuration degeneracy, a grand canonical ensemble (GCE) is


32

Page 33 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


constructed for a vacancy-pair free (VPF) model as described in section 2.2. The expectation value of 𝑥 in Na1+xNiFe(CN)6 as a function of equilibrium reduction potential 𝜙𝑒𝑞 is calculated and compared to experimental data. The satisfactory agreement between these two datasets indicates that vacancy/cation ordering plays an important role on electrochemical equilibria in PBAs. The importance of this ordering is further analyzed by comparing the VPF model with a non-interacting solution. The failure of the non-interacting system to predict the equilibrium potential curve is attributed to its neglect of the energy costs associated with nearest-neighbor vacancy pairing. 3.3.1 Classification and Structure of Vacancy-Pair Free Configurations In this section we present the results of HMM simulations used to sample vacancypair free configurations as a function of vacancy number 𝑚. From the sampled configurations we estimate the average guard numbers 〈〈𝑔𝑚〉〉 that are used in Sec. 3.3.3 to predict Ω𝑚 and, ultimately, 𝜙𝑒𝑞 versus degree of intercalation 〈𝑥〉.

We first use

combinatorial methods to classify the different types of VPFCs, based upon which we subsequently interpret the scaling of 〈〈𝑔𝑚〉〉 with vacancy density 𝜌 = 𝑚/𝑁 simulated using the HMM. Using the recurrence relationship of Eq. 5, limits for the guard number 〈𝑔𝑚〉 of a particular configuration can be obtained as a function of 𝜌. The upper limit of 〈𝑔𝑚〉 = 1/𝜌 is derived by requiring that Ω𝑚 + 1 = 0, such that no VPFCs exist for vacancy numbers greater than 𝑚. We classify VPFCs with such a property as saturated configurations (SCs) because no cations can be incrementally substituted by a vacancy without forming a vacancy pair. The envelope of 〈𝑔𝑚〉 values for SCs is shown in Figure 9a, together with


33


Page 34 of 52

an example SC configuration for a 4 × 4 supercell in the upper right inset. While the four vacancies shown in the inset (white) are smaller than the number of vacancies at the TC (𝑁 2 = 8), substitution of a vacancy onto any occupied site (orange) in this particular VPFC produces a NNVP.

Figure 8. Schematics of complete chessboard ordering (a), partial chessboard ordering (b), and anti-chessboard ordering (c) for an 8 × 8 supercell. The partial chessboard ordering in (b) is derived from (a) by filling the vacancies at crimson type-2 sites.

We derive an “under-saturated” limit for 〈𝑔𝑚〉 by considering the principle of chessboard ordering. VPFCs that coincide with this lower bound we classify as undersaturated configurations (USCs). To explain the principle of chessboard ordering we start with a TC, the projection of which is shown for an 8 × 8 × 8 supercell in Figure 8a. By using cyclic permutations of an eight-element binary sequence along the two directions in the plane of projection, all sites on that plane are labeled as either type-1 and type-2. This permutation analysis applies similarly to the third dimension, though not shown


34

Page 35 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


explicitly. By using the labeled site distributions we identify the complete chessboard ordering (CCO) pattern73, where for a given VPFC vacancies occupy one of the two site types. USCs possess partial chessboard ordering (PCO), where cations can substitute with the vacancies of TCs on the crimson type-2 sites in Figure 8b. An example USC for a 4 × 4 supercell is shown in the lower left inset of Figure 9a. For such a configuration, a chain of vacancies is formed that is oriented along the anti-diagonal direction. Based on these principles the USC limit of 〈𝑔𝑚〉 = 1 (2𝜌) +1, shown in Figure 9a, can be determined by requiring the number of sites available for incremental vacancy insertion (𝑁 ― 𝑚 × 〈𝑔𝑚〉 from Eq. 8) to match the number of cations inserted into the TC (𝑁/2 ― 𝑚). We note that the USC limit is not the lower bound for calculated 〈𝑔𝑚〉 curve when the vacancy density is low, as indicated by the calculated guard number curves in Figure 9a that fall below the USC limit. This effect occurs because the vacancy number of such VPFCs is not large enough to establish PCO, and, thus, a mixture of PCO and antichessboard ordering (ACO) is present in the VPFC. We classify configurations having a mixture of PCO and ACO as normal configurations (NCs). An example VPFC with ACO is shown in Figure 9c by repeating the upper right overlay in Figure 9a in 2D space, where all the vacancies are evenly distributed in both type-1 and type-2 lattice sites. The fraction of VPFCs constituted by the three classes identified (SCs, USCs, and NCs) varies with 𝜌 when simulated using the HMM. Also, guard numbers in the SC and USC limits provide guidelines for interpreting the variation of average guard number

〈〈𝑔𝑚〉〉 simulated using the HMM. While we ultimately use a smearing coefficient 𝑠 = 0.77 to reproduce TC degeneracy, here we also examine the results for 𝑠 = 0.10. Irrespective of the choice of the smearing coefficient 𝑠 from 0 to 1, we find that the HMM produces


35


Page 36 of 52

calculated 〈〈𝑔𝑚〉〉 curves that vary between SC and USC limits at high vacancy density. However, an arbitrary choice of 𝑠 can cause the simulated value of Ω𝑁/2 to deviate from its exact value known from crystal symmetry. This deviation occurs because small 𝑠 values tend to increase the frequency of SCs while large 𝑠 values tend to suppress SC frequency. These effects are apparent from the probability distributions of SCs, USCs, and NCs obtained for 𝑠 = 0.77 (Figure 9b) and 𝑠 = 0.10 (Figure 9c). We obtain these distributions by counting the number of sampled VPFCs which have 〈〈𝑔𝑚〉〉 that coincide with either USC or SC limits. For 𝑠 = 0.10, there is a high probability of sampling SCs in a 4 × 4 × 4 supercell when 𝜌 > 0.25 (i.e., 𝑚 > 16). Such a bias toward SCs causes the corresponding 〈〈𝑔𝑚〉〉 curve to shift toward the SC limit, producing a simulated Ω𝑁/2 value of 9.97 instead of 2.

Figure 9. (a) Average guard number 〈〈𝑔𝑚〉〉 as a function of vacancy density 𝜌 simulated with 𝑠 equal to 0.77 and 0.10 for a 4 × 4 × 4 supercell. The inset diagrams in (a) show a saturated configuration (SC, upper right) and an undersaturated configuration (USC, lower left), where white squares are vacancies. The cyan (𝑠 = 0.77) and green (𝑠 = 0.10) lines correspond to the VPFC class distributions shown in (b) and (c), respectively. Each


36

Page 37 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


column in these distributions shows the probability of finding different VPFC types with a certain vacancy number 𝑚.

In contrast, by using an 𝑠 value equal to 0.77 we observe that an increased population of USCs is required to enforce the Ω𝑁/2 = 2 condition for the TC. Figure 9a shows that the 〈〈𝑔𝑚〉〉 curve for 𝑠 = 0.77 and the USC limit collapse on each other when 𝜌 > 0.25, suggesting that the vacancy/cation arrangement approaches CCO as vacancy number increases. In Figure 9b, the corresponding probability distribution indicates that the SC is rarely found due to the abundance of USCs and NCs.

3.3.2 Equilibrium Reduction Potential and Fluctuations In this section we use the results for the degeneracy Ω𝑚 of VPFCs with vacancy number 𝑚 to predict the equilibrium potential 𝜙𝑒𝑞 for cation intercalation into PBAs. We use the grand canonical ensemble (GCE) formalism from Sec. 2.2, and we compare the results of this model to experiment and to a non-interacting solution model. The effect of microstate size on the results is examined by simulating 4 × 4 × 4, 6 × 6 × 6, and 8 × 8 × 8 supercells at 𝑇 = 300𝐾. Figure 10a shows the calculated profiles of Ω𝑚 versus vacancy density 𝜌. The degeneracy value at a given vacancy density increases with the supercell size, but the shape of the profiles is similar among them with all profiles having a maximum at 𝜌 ≈ 0.2 and having Ω𝑁/2 ≈ 2. By substituting the calculated degeneracy values in Figure 10a into Eq. 6 and sweeping over an appropriate range of 𝜙𝑒𝑞, the equilibrium potential curves can be obtained as a function of 〈𝑥〉 for the ensembles with


37


Page 38 of 52

these supercell sizes. Among the three supercell sizes investigated 𝜙𝑒𝑞 at any given 〈𝑥〉 varies by less than 10 ―5 meV, indicating satisfactory convergence.

Figure 10. (a) Degeneracy computed from the HMM for 4 × 4 × 4, 6 × 6 × 6, and 8 × 8 × 8 supercells as a function of vacancy density 𝜌. (b) The calculated equilibrium potential curve at 𝑇 = 300𝐾 based on 4 × 4 × 4 ensemble analysis (red), and the closed form of Ω (black dash line) are compared with experimental measurements7 (blue) and a noninteracting solution model (green) at the same temperature.

Given the negligible influence of supercell size on 𝜙𝑒𝑞 curve, any one of the curves shown in Figure 10a can be used to construct a closed form for Ω as a function of 𝜌. A least-squares fit of the 4 × 4 × 4 data shown in Figure 10a produces the following closedform expression: ln (Ω) = ―2244𝜌4 +2481𝜌3 ―1220𝜌2 +272𝜌 + 0.048.


(10)

38

Page 39 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The degeneracy values predicted by Eq. 10 can be used to approximate 𝜙𝑒𝑞 versus 〈𝑥〉 (Figure 10b) with a maximum deviation from the VPF model of 3 meV. The variation of 𝜙𝑒𝑞 with 〈𝑥〉 is shown in Figure 10b using microstates generated with 4 × 4 × 4 supercells, together with experimental data.7 While we assume an ideal PBA stoichiometry in our model, we note that PBAs tested experimentally often deviate from ideal stoichiometry due to lattice defects that we neglect. Accordingly, we compare our predicted results to the experimental curves of Ref. 7 normalized to the composition range expected based on ideal stoichiometry. The present study is focused on modeling cation/vacancy ordering in PBAs, and thus each microstate in the GCE is a pristine lattice free of defects other than Na+ vacancies. However, the as-synthesized NiFe-PBA tested in Ref. 7 was shown to have excess Ni2+ and a deficiency of Na+, producing the formula 𝑁 𝑎0.57𝑁𝑖1.37𝐹𝑒0.92(𝐶𝑁0.95)6 ∙ 5.53𝐻2𝑂 based on elemental analysis compared with the ideal formula 𝑁𝑎2𝑁𝑖𝐹𝑒(𝐶𝑁)6. While post-synthesis cycling in pure NaCl solution was used to ion-exchange intercalated Ni2+ with Na+, the precise composition of the PBA was not determined. The intercalation of large cations (e.g., Ni2+) distorts lattice space and could affect the diffusion of smaller cations.89 Also [𝐹𝑒(𝐶𝑁)6]3 ― vacancies induce structural instability and can cause the dissolution of metal ions into electrolytes.90,91 Thus, the present model neglects the [𝐹𝑒(𝐶𝑁)6]3 ― vacancies that originate during material synthesis and the occupation of 𝑁𝑖2 + in 𝑁𝑎 + sites in the PBA framework. Furthermore, the current model also excludes the influence of VPCs on the 𝜙𝑒𝑞, which can become statistically possible at low DOI, as shown in Figure 4. Such exclusion may also contribute to the deviation between the VPF model and experimental data when 𝑥 < 0.25.


39


Page 40 of 52

The value of 𝜀 in Eq. 3 is set to 370 𝑚𝑒𝑉 to match the experimental potential at 〈𝑥〉 = 0.50. A maximum deviation of 30 𝑚𝑒𝑉 between the vacancy-pair free model and the experimental data occurs at 𝑥 = 0.15. Considering that this model neglects the effects of volumetric variation of the lattice5 and the role of hydrated ions32 in the interstitial space, the level agreement shown here confirms that the energetic cost of vacancy pairing plays a significant role in the intercalation of cations into PBAs. To further investigate the importance of vacancy pairing in the present system we also compare the VPF model with a non-interacting solution (NIS) model. For the NIS model the energy difference between VPFCs and VPCs is identically zero, and thus all configurations with certain 𝑚 are degenerate irrespective of whether they have vacancy

( )

𝑁! 𝑁 pairs. The value of Ω𝑚 can be simply calculated as 𝑚 = 𝑚!(𝑁 ― 𝑚)!. Using 𝜀 = 401 meV

to match the experimental potential at 〈𝑥〉 = 0.50, we simulated 𝜙𝑒𝑞 versus 〈𝑥〉 for microstates simulated with 4 × 4 × 4 supercells (Figure 10b).

For 〈𝑥〉 > 0.5 the NIS

deviates from experiment by approximately 15 meV, but for 〈𝑥〉 < 0.5 the NIS has large systematic deviation from experiment with a maximum deviation of 120 meV at 𝑥 = 0.08. By comparing 𝜙𝑒𝑞 between the VPF and NIS models we deduce that the persistence of vacancy pairs in the NIS is responsible for these quantitative and qualitative deviations from experiment. To further investigate how the cation/vacancy ordering predicted by the VPF differs from the NIS we analyzed their respective statistical fluctuations in 𝑥. A useful measure of the fluctuations of 𝑥 within the grand canonical ensemble is its variance 𝐹 given as 𝐹 =

〈𝑥2〉 ― 〈𝑥〉2, shown in Figs. 11a and 11b versus 〈𝑥〉. A Dirac delta distribution of


40

Page 41 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


microstates would produce no fluctuations and would possess 𝐹 = 0. Thus, we find that 𝐹 = 0 at 〈𝑥〉 = 1 for both the VPF and NIS models because no vacancies are occupied when the PBA is completed intercalated. The vanishing fluctuation in the NIS system at low 〈𝑥〉 occurs because of the special distribution of 𝑥 of the subsystems. The inset in Figure 11c shows the microstate probability distribution 𝑃(𝑚) =

Ω𝑚 Ξ

(

𝑒𝑥𝑝 ―

𝐸𝑚 + 𝑚𝜇 𝑘𝐵𝑇

) for the

NIS and VPF models at 〈𝑥〉 = 0.036. While the distribution for NIS models has a singletailed shape, the VPF model produces a bell shaped distribution. The ratio of fluctuations 𝐹𝑉𝑃𝐹 𝐹𝑁𝐼𝑆 is shown in Figure 11c. Within the range of 0.10 < 〈𝑥〉 < 0.50 this ratio is approximately 30%, implying that the VPF system generally has smaller fluctuations among microstates than an NIS system of the same size. For 〈𝑥〉 > 0.5 the ratio of 𝐹𝑉𝑃𝐹 𝐹𝑁𝐼𝑆 is independent of system size, and it gradually approaches unity for 〈𝑥〉→1. This limit corresponds to the gradual disappearance of vacancy/cation ordering near the fully-intercalated state. A size dependence of fluctuations is also observed for both systems, where 𝐹 decreases with system size (Figure 11a and 11b). Further, the ratio of fluctuations shifts toward low degree of intercalation as system size increases (Figure 11c). At 〈𝑥〉 = 0, however, the fluctuations calculated from the VPF model 𝐹𝑉𝑃𝐹 are significantly larger than that of the NIS model (Figure 11c). Moreover, the results from the VPF model (Figure 11a) show that the largest fluctuation is obtained at the half-intercalated state, while the NIS model (Figure 11b) has maximum fluctuation at a smaller 〈𝑥〉. This difference can also be attributed to the success of simulating the vacancy/cation ordering in VPF. In the


41


Page 42 of 52

limit of large system size, the peak in 𝐹𝑁𝐼𝑆 is expected to shift to 〈𝑥〉 = 0 or to disappear altogether, such that the value of 𝐹𝑁𝐼𝑆 at 〈𝑥〉 = 0 will be finite.

Figure 11. The statistical fluctuation of 𝑥 as a function of 〈𝑥〉 for (a) the VPF model and (b) the NIS model. (c) The ratio of the fluctuations of these two models is shown as a function of 〈𝑥〉. The inset of (c) shows the occupation probability distribution 𝑃(𝑚) for microstates with certain 𝑥 for the VPF model and NIS model when 〈𝑥〉 = 0.036 for the 4 × 4 × 4 supercell.

4. Conclusions In this study we used DFT calculations to explore the influence of vacancy/Na-ion arrangement on the thermodynamics and electronic structure of the NiFe-PBA material. Analysis of ground-state electronic energies for configurations with (VPC) and without (VPFC) nearest-neighbor vacancy pair revealed that the VPFC is energetically preferred because it reduces electrostatic repulsion between intercalated cations and exchangecorrelation interactions. By visualizing the charge density distribution and analyzing the


42

Page 43 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


electronic density of states projected onto the various atoms comprising Na1+xNiFe(CN)6, electronic charge sharing (i.e., weak covalent interaction) between Na atoms and their surrounding atoms was found to increase charge density around the Na-nuclei (such that the formal valence of Na was weakly reduced), an effect that is ultimately responsible for the energetic preference for vacancy repulsion. Based on our findings from DFT, we then created a vacancy-pair free (VPF) statistical mechanics model to investigate the equilibrium behavior of NiFe-PBA by incorporating the effects of configurational entropy at finite temperature using grand canonical ensemble theory with NiFe-PBA supercells as microstates. The small deviation between the predicted potential curve and experiment7 provides compelling evidence for the substantial influence of vacancy/cation ordering on the equilibrium behavior of PBA materials. The grand canonical ensemble (GCE) method that we have used to predict equilibrium potential 𝜙𝑒𝑞 can be applied to model the equilibrium intercalation of other cations in other host materials, as well. To determine the partition function used in the GCE method, we introduced a recurrence relationship for the degeneracy of VPFC energy levels Ω𝑚, and we used a hybrid Markov method (HMM) to estimate the Ω𝑚 values. Relying on the results from the HMM, we classified three kinds of VPFCs (i.e., under-saturated, saturated, and normal configurations) by identifying the vacancy arrangements of partial-chessboard and antichessboard ordering. Although the dynamic behavior of PBA material beyond the scope of this study, the vacancy arrangements discovered here shed light on the mechanisms of crystal-scale diffusion of vacancies and cations in PBAs. We also constructed a closed form equation for Ω𝑚 as a function of vacancy density to estimate the 𝜙𝑒𝑞 as a function of


43


Page 44 of 52

composition. This approximate degeneracy function produces 𝜙𝑒𝑞 within 3 mV of the result using the HMM, and, thus, it can also be used to calculation of other thermodynamic observables within the GCE. We also compared the VPF model statistics to those of a non-interacting solution (NIS) model to investigate the importance of vacancy/cation ordering. For the NIS energy differences between VPFCs and VPCs were neglected altogether. The predicted 𝜙𝑒𝑞 from NIS at low 〈𝑥〉 in Na1+Fe(CN)6 produces large deviation with experiment, having a maximal deviation of 120 mV at 〈𝑥〉 = 0.08. By analyzing the variance of 𝑥 predicted with the VPF and NIS models, we showed that fluctuations in the NIS system vanish at low 〈𝑥〉, indicating that the discrepancy between NIS predictions and experiment is caused by the abundance of nearest-neighbor vacancy pairs when there is no energetic cost assigned to their formation. There are some limitations of the present VPF model that are worth noting. The present VPF model assumed a negligible dependence of the lattice energy on the local strain induced by cation intercalation. While this assumption is valid for small-ion intercalation in NiFe-PBA,5,6 the intercalation of large ions could cause significant changes in the lattice framework.24,89 The current model also neglects water molecules in the interstitial space. However, the interstitial water in hydrated PBAs could make the lattice structure less stable due to the H2O-enhanced Pauli repulsion within the host.32 In aqueous systems the influence of water molecules on the energetic cost of vacancy pairing deserves further investigation. Finally, the Markov process employed in the HMM ignores the diffusion barrier of cations and vacancies in their movement between sites. For intercalation materials in general diffusion barriers vary with the type of cation.29


44

Page 45 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Therefore, to model dynamic cation diffusion processes using an HMM requires the incorporation of diffusion barriers.

5. Acknowledgments The Department of Mechanical Science and Engineering at the University of Illinois at Urbana-Champaign (UIUC) supported this work. We acknowledge the use of UIUC Campus Computing Cluster resources made possible through the Beckman Institute for Advanced Science and Technology. We are grateful for the constructive feedback given provided by the anonymous reviewers of our manuscript.


45


Page 46 of 52

6. Reference (1)

(2) (3) (4) (5) (6) (7) (8)

(9)

(10) (11) (12) (13) (14)

Nie, P.; Yuan, J.; Wang, J.; Le, Z.; Xu, G.; Hao, L.; Pang, G.; Wu, Y.; Dou, H.; Yan, X.; et al. Prussian Blue Analogue with Fast Kinetics Through Electronic Coupling for Sodium Ion Batteries. ACS Appl. Mater. Interfaces 2017, 9, 20306–20312. https://doi.org/10.1021/acsami.7b05178. Wang, B.; Han, Y.; Wang, X.; Bahlawane, N.; Pan, H.; Yan, M.; Jiang, Y. Prussian Blue Analogs for Rechargeable Batteries. iScience 2018, 3, 110–133. https://doi.org/10.1016/j.isci.2018.04.008. Wessells, C. D.; Peddada, S. V.; McDowell, M. T.; Huggins, R. A.; Cui, Y. The Effect of Insertion Species on Nanostructured Open Framework Hexacyanoferrate Battery Electrodes. J. Electrochem. Soc. 2012, 159, A98. https://doi.org/10.1149/2.060202jes. Wessells, C. D.; Huggins, R. A.; Cui, Y. Copper Hexacyanoferrate Battery Electrodes with Long Cycle Life and High Power. Nat. Commun. 2011, 2, 550. https://doi.org/10.1038/ncomms1563. Wessells, C. D.; Peddada, S. V.; Huggins, R. A.; Cui, Y. Nickel Hexacyanoferrate Nanoparticle Electrodes for Aqueous Sodium and Potassium Ion Batteries. Nano Lett. 2011, 11, 5421–5425. https://doi.org/10.1021/nl203193q. H. Li, C.; Nanba, Y.; Asakura, D.; Okubo, M.; R. Talham, D. Li-Ion and Na-Ion Insertion into Size-Controlled Nickel Hexacyanoferrate Nanoparticles. RSC Adv. 2014, 4, 24955– 24961. https://doi.org/10.1039/C4RA03296A. Shrivastava, A.; Smith, K. C. Electron Conduction in Nanoparticle Agglomerates Limits Apparent Na + Diffusion in Prussian Blue Analogue Porous Electrodes. J. Electrochem. Soc. 2018, 165, A1777–A1787. https://doi.org/10.1149/2.0861809jes. Targholi, E.; Mousavi-Khoshdel, S. M.; Rahmanifara, M.; Yahya, M. Z. A. Cu- and FeHexacyanoferrate as Cathode Materials for Potassium Ion Battery: A First-Principles Study. Chem. Phys. Lett. 2017, 687, 244–249. https://doi.org/10.1016/j.cplett.2017.09.029. Mizuno, Y.; Okubo, M.; Hosono, E.; Kudo, T.; Zhou, H.; Oh-Ishi, K. Suppressed Activation Energy for Interfacial Charge Transfer of a Prussian Blue Analog Thin Film Electrode with Hydrated Ions (Li+, Na +, and Mg2+). J. Phys. Chem. C 2013, 117, 10877–10882. https://doi.org/10.1021/jp311616s. Kuperman, N.; Padigi, P.; Goncher, G.; Evans, D.; Thiebes, J.; Solanki, R. High Performance Prussian Blue Cathode for Nonaqueous Ca-Ion Intercalation Battery. J. Power Sources 2017, 342, 414–418. https://doi.org/10.1016/j.jpowsour.2016.12.074. Liu, S.; Smith, K. C. Quantifying the Trade-Offs between Energy Consumption and Salt Removal Rate in Membrane-Free Cation Intercalation Desalination. Electrochimica Acta 2018, 271, 652–665. https://doi.org/10.1016/j.electacta.2018.03.065. Smith, K. C.; Dmello, R. Na-Ion Desalination (NID) Enabled by Na-Blocking Membranes and Symmetric Na-Intercalation: Porous-Electrode Modeling. J. Electrochem. Soc. 2016, 163, A530–A539. https://doi.org/10.1149/2.0761603jes. Smith, K. C. Theoretical Evaluation of Electrochemical Cell Architectures Using Cation Intercalation Electrodes for Desalination. Electrochimica Acta 2017, 230, 333–341. https://doi.org/10.1016/j.electacta.2017.02.006. Choi, S.; Chang, B.; Kim, S.; Lee, J.; Yoon, J.; Choi, J. W. Battery Electrode Materials with Omnivalent Cation Storage for Fast and Charge-Efficient Ion Removal of Asymmetric Capacitive Deionization. Adv. Funct. Mater. 2018, 1802665. https://doi.org/10.1002/adfm.201802665.


46

Page 47 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) (16)

(17) (18) (19)

(20)

(21) (22)

(23) (24)

(25)

(26) (27) (28)

Kim, T.; Gorski, C. A.; Logan, B. E. Low Energy Desalination Using Battery Electrode Deionization. Environ. Sci. Technol. Lett. 2017, 4, 444–449. https://doi.org/10.1021/acs.estlett.7b00392. Porada, S.; Shrivastava, A.; Bukowska, P.; Biesheuvel, P. M.; Smith, K. C. Nickel Hexacyanoferrate Electrodes for Continuous Cation Intercalation Desalination of Brackish Water. Electrochimica Acta 2017, 255, 369–378. https://doi.org/10.1016/j.electacta.2017.09.137. You, Y.; Wu, X.-L.; Yin, Y.-X.; Guo, Y.-G. A Zero-Strain Insertion Cathode Material of Nickel Ferricyanide for Sodium-Ion Batteries. J. Mater. Chem. A 2013, 1, 14061. https://doi.org/10.1039/c3ta13223d. Margadonna, S.; Prassides, K.; Fitch, A. N. Zero Thermal Expansion in a Prussian Blue Analogue. J. Am. Chem. Soc. 2004, 126, 15390–15391. https://doi.org/10.1021/ja044959o. Flambard, A.; Köhler, F. H.; Lescouëzec, R. Revisiting Prussian Blue Analogues with Solid-State MAS NMR Spectroscopy: Spin Density and Local Structure in [Cd3Fe(CN)62]⋅15 H2O. Angew. Chem. Int. Ed. 2009, 48, 1673–1676. https://doi.org/10.1002/anie.200805415. Wang, H.; Jang, Y.-I.; Huang, B.; Sadoway, D. R.; Chiang, Y.-M. TEM Study of Electrochemical Cycling‐Induced Damage and Disorder in LiCoO2 Cathodes for Rechargeable Lithium Batteries. J. Electrochem. Soc. 1999, 146, 473–480. https://doi.org/10.1149/1.1391631. Zhang, X.; Shyy, W.; Sastry, A. M. Numerical Simulation of Intercalation-Induced Stress in Li-Ion Battery Electrode Particles. J. Electrochem. Soc. 2007, 154, A910–A916. https://doi.org/10.1149/1.2759840. Meethong, N.; Huang, H.-Y. S.; Speakman, S. A.; Carter, W. C.; Chiang, Y.-M. Strain Accommodation during Phase Transformations in Olivine-Based Cathodes as a Materials Selection Criterion for High-Power Rechargeable Batteries. Adv. Funct. Mater. 2007, 17, 1115–1123. https://doi.org/10.1002/adfm.200600938. Liu, W.; Oh, P.; Liu, X.; Lee, M.-J.; Cho, W.; Chae, S.; Kim, Y.; Cho, J. Nickel-Rich Layered Lithium Transition-Metal Oxide for High-Energy Lithium-Ion Batteries. Angew. Chem. Int. Ed. 2015, 54, 4440–4457. https://doi.org/10.1002/anie.201409262. Li, Z.; Xiang, K.; Xing, W.; Carter, W. C.; Chiang, Y.-M. Reversible Aluminum-Ion Intercalation in Prussian Blue Analogs and Demonstration of a High-Power Aluminum-Ion Asymmetric Capacitor. Adv. Energy Mater. 2015, 5, 1401410. https://doi.org/10.1002/aenm.201401410. Felts, A. C.; Andrus, M. J.; Knowles, E. S.; Quintero, P. A.; Ahir, A. R.; Risset, O. N.; Li, C. H.; Maurin, I.; Halder, G. J.; Abboud, K. A.; et al. Evidence for Interface-Induced Strain and Its Influence on Photomagnetism in Prussian Blue Analogue Core–Shell Heterostructures, RbaCob[Fe(CN)6]C·mH2O@KjNik[Cr(CN)6]L·nH2O. J. Phys. Chem. C 2016, 120, 5420–5429. https://doi.org/10.1021/acs.jpcc.5b10761. Okubo, M.; Honma, I. Ternary Metal Prussian Blue Analogue Nanoparticles as Cathode Materials for Li-Ion Batteries. Dalton Trans. 2013, 42, 15881–15884. https://doi.org/10.1039/C3DT51369F. Wu, X.; Deng, W.; Qian, J.; Cao, Y.; Ai, X.; Yang, H. Single-Crystal FeFe(CN)6 Nanoparticles: A High Capacity and High Rate Cathode for Na-Ion Batteries. J. Mater. Chem. A 2013, 1, 10130–10134. https://doi.org/10.1039/C3TA12036H. Asakura, D.; Li, C. H.; Mizuno, Y.; Okubo, M.; Zhou, H.; Talham, D. R. Bimetallic Cyanide-Bridged Coordination Polymers as Lithium Ion Cathode Materials: Core@Shell Nanoparticles with Enhanced Cyclability. J. Am. Chem. Soc. 2013, 135, 2793–2799. https://doi.org/10.1021/ja312160v.


47


(29)

(30)

(31) (32) (33)

(34) (35)

(36) (37) (38) (39) (40) (41) (42) (43)

Page 48 of 52

Ping Ong, S.; L. Chevrier, V.; Hautier, G.; Jain, A.; Moore, C.; Kim, S.; Ma, X.; Ceder, G. Voltage, Stability and Diffusion Barrier Differences between Sodium-Ion and Lithium-Ion Intercalation Materials. Energy Environ. Sci. 2011, 4, 3680–3688. https://doi.org/10.1039/C1EE01782A. Okubo, M.; Asakura, D.; Mizuno, Y.; Kudo, T.; Zhou, H.; Okazawa, A.; Kojima, N.; Ikedo, K.; Mizokawa, T.; Honma, I. Ion-Induced Transformation of Magnetism in a Bimetallic CuFe Prussian Blue Analogue. Angew. Chem. Int. Ed. 2011, 50, 6269–6273. https://doi.org/10.1002/anie.201102048. You, Y.; Yu, X.; Yin, Y.; Nam, K.-W.; Guo, Y.-G. Sodium Iron Hexacyanoferrate with High Na Content as a Na-Rich Cathode Material for Na-Ion Batteries. Nano Res. 2015, 8, 117– 128. https://doi.org/10.1007/s12274-014-0588-7. Xiao, P.; Song, J.; Wang, L.; Goodenough, J. B.; Henkelman, G. Theoretical Study of the Structural Evolution of a Na 2 FeMn(CN) 6 Cathode upon Na Intercalation. Chem. Mater. 2015, 27, 3763–3768. https://doi.org/10.1021/acs.chemmater.5b01132. Kajiyama, S.; Mizuno, Y.; Okubo, M.; Kurono, R.; Nishimura, S.; Yamada, A. Phase Separation of a Hexacyanoferrate-Bridged Coordination Framework under Electrochemical Na-Ion Insertion. Inorg. Chem. 2014, 53, 3141–3147. https://doi.org/10.1021/ic403088r. Boudjema, L.; Mamontova, E.; Long, J.; Larionova, J.; Guari, Y.; Trens, P. Prussian Blue Analogues for the Separation of Hydrocarbons in Humid Conditions. Inorg. Chem. 2017, 56, 7598–7601. https://doi.org/10.1021/acs.inorgchem.7b00563. Gao, Q.; Shi, N.; Sanson, A.; Sun, Y.; Milazzo, R.; Olivi, L.; Zhu, H.; Lapidus, S. H.; Zheng, L.; Chen, J.; et al. Tunable Thermal Expansion from Negative, Zero, to Positive in Cubic Prussian Blue Analogues of GaFe(CN)6. Inorg. Chem. 2018. https://doi.org/10.1021/acs.inorgchem.8b02428. Klink, S.; Ishige, Y.; Schuhmann, W. Prussian Blue Analogues: A Versatile Framework for Solid-Contact Ion-Selective Electrodes with Tunable Potentials. ChemElectroChem 2017, 4, 490–494. https://doi.org/10.1002/celc.201700091. Erinmwingbovo, C.; Palagonia, M. S.; Brogioli, D.; La Mantia, F. Intercalation into a Prussian Blue Derivative from Solutions Containing Two Species of Cations. ChemPhysChem 2017, 18, 917–925. https://doi.org/10.1002/cphc.201700020. Singh, K.; Bouwmeester, H. J. M.; de Smet, L. C. P. M.; Bazant, M. Z.; Biesheuvel, P. M. Theory of Water Desalination with Intercalation Materials. Phys. Rev. Appl. 2018, 9, 064036. https://doi.org/10.1103/PhysRevApplied.9.064036. Gao, C.; Lee, S. W.; Yang, Y. Thermally Regenerative Electrochemical Cycle for LowGrade Heat Harvesting. ACS Energy Lett. 2017, 2, 2326–2334. https://doi.org/10.1021/acsenergylett.7b00568. Lee, S. W.; Yang, Y.; Lee, H.-W.; Ghasemi, H.; Kraemer, D.; Chen, G.; Cui, Y. An Electrochemical System for Efficiently Harvesting Low-Grade Heat Energy. Nat. Commun. 2014, 5, 3942. https://doi.org/10.1038/ncomms4942. C. Wojdeł, J.; R. Moreira, I. de P.; T. Bromley, S.; Illas, F. Prediction of Half-Metallic Conductivity in Prussian Blue Derivatives. J. Mater. Chem. 2009, 19, 2032–2036. https://doi.org/10.1039/B813788A. Ling, C.; Chen, J.; Mizuno, F. First-Principles Study of Alkali and Alkaline Earth Ion Intercalation in Iron Hexacyanoferrate: The Important Role of Ionic Radius. J. Phys. Chem. C 2013, 117, 21158–21165. https://doi.org/10.1021/jp4078689. Zhou, F.; Cococcioni, M.; Kang, K.; Ceder, G. The Li Intercalation Potential of LiMPO4 and LiMSiO4 Olivines with M=Fe, Mn, Co, Ni. Electrochem. Commun. 2004, 6, 1144– 1148. https://doi.org/10.1016/j.elecom.2004.09.007.


48

Page 49 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(44) (45) (46) (47) (48) (49) (50) (51) (52) (53)

(54)

(55) (56) (57)

(58)

Coleman, S. T.; McKinnon, W. R.; Dahn, J. R. Lithium Intercalation in Li x Mo 6 Se 8 : A Model Mean-Field Lattice Gas. Phys. Rev. B 1984, 29, 4147–4149. https://doi.org/10.1103/PhysRevB.29.4147. Chang, D.; Ven, A. V. der. Li Intercalation Mechanisms in CaTi 5 O 11 , a Bronze-B Derived Compound. Phys. Chem. Chem. Phys. 2016, 18, 32042–32049. https://doi.org/10.1039/C6CP05905H. Van der Ven, A.; Thomas, J. C.; Xu, Q.; Swoboda, B.; Morgan, D. Nondilute Diffusion from First Principles: Li Diffusion in Li x TiS 2. Phys. Rev. B 2008, 78. https://doi.org/10.1103/PhysRevB.78.104306. Smith, K. C.; Fisher, T. S.; Waghmare, U. V.; Grau-Crespo, R. Dopant-Vacancy Binding Effects in Li-Doped Magnesium Hydride. Phys. Rev. B 2010, 82. https://doi.org/10.1103/PhysRevB.82.134109. Reimers, J. N.; Dahn, J. R. Application of Ab Initio Methods for Calculations of Voltage as a Function of Composition in Electrochemical Cells. Phys. Rev. B 1993, 47, 2995–3000. https://doi.org/10.1103/PhysRevB.47.2995. Wolverton, C.; Zunger, A. First-Principles Prediction of Vacancy Order-Disorder and Intercalation Battery Voltages in Li x CoO 2. Phys. Rev. Lett. 1998, 81, 606–609. https://doi.org/10.1103/PhysRevLett.81.606. Fontaine, D. D. Cluster Approach to Order-Disorder Transformations in Alloys. In Solid State Physics; Ehrenreich, H., Turnbull, D., Eds.; Academic Press, 1994; Vol. 47, pp 33– 176. https://doi.org/10.1016/S0081-1947(08)60639-6. Drews, T. O.; Braatz, R. D.; Alkire, R. C. Coarse-Grained Kinetic Monte Carlo Simulation of Copper Electrodeposition with Additives. Int. J. Multiscale Comput. Eng. 2004, 2, 313– 327. https://doi.org/10.1615/IntJMultCompEng.v2.i2.90. Grau-Crespo, R.; Hamad, S.; Catlow, C. R. A.; Leeuw, N. H. de. Symmetry-Adapted Configurational Modelling of Fractional Site Occupancy in Solids. J. Phys. Condens. Matter 2007, 19, 256201. https://doi.org/10.1088/0953-8984/19/25/256201. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys. Condens. Matter 2009, 21, 395502. https://doi.org/10.1088/0953-8984/21/39/395502. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671– 6687. https://doi.org/10.1103/PhysRevB.46.6671. Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91. https://doi.org/10.1103/PhysRevLett.91.146401. Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for Mixing Exact Exchange with Density Functional Approximations. J. Chem. Phys. 1996, 105, 9982–9985. https://doi.org/10.1063/1.472933. Ji, Z.; Han, B.; Liang, H.; Zhou, C.; Gao, Q.; Xia, K.; Wu, J. On the Mechanism of the Improved Operation Voltage of Rhombohedral Nickel Hexacyanoferrate as Cathodes for Sodium-Ion Batteries. ACS Appl. Mater. Interfaces 2016, 8, 33619–33625. https://doi.org/10.1021/acsami.6b11070. Hegner, F. S.; Galán-Mascarós, J. R.; López, N. A Database of the Structural and Electronic Properties of Prussian Blue, Prussian White, and Berlin Green Compounds through Density Functional Theory. Inorg. Chem. 2016, 55, 12851–12862. https://doi.org/10.1021/acs.inorgchem.6b02200.


49


(59) (60) (61) (62) (63) (64) (65) (66) (67) (68) (69)

(70)

(71) (72) (73) (74) (75) (76)

Page 50 of 52

Moroni, E. G.; Kresse, G.; Hafner, J.; Furthmüller, J. Ultrasoft Pseudopotentials Applied to Magnetic Fe, Co, and Ni: From Atoms to Solids. Phys. Rev. B 1997, 56, 15629–15646. https://doi.org/10.1103/PhysRevB.56.15629. Yu, Q.; Steen, W. A.; Jeerage, K. M.; Jiang, S.; Schwartz, D. T. Structure-Dependent Solvent and Ion Intercalation in Reduced and Oxidized Nickel Hexacyanoferrates. J. Electrochem. Soc. 2002, 149, E195–E203. https://doi.org/10.1149/1.1474434. Henkelman, G.; Arnaldsson, A.; Jónsson, H. A Fast and Robust Algorithm for Bader Decomposition of Charge Density. Comput. Mater. Sci. 2006, 36, 354–360. https://doi.org/10.1016/j.commatsci.2005.04.010. Matsuda, T.; Kim, J.; Moritomo, Y. Symmetry Switch of Cobalt Ferrocyanide Framework by Alkaline Cation Exchange. J. Am. Chem. Soc. 2010, 132, 12206–12207. https://doi.org/10.1021/ja105482k. Momma, K.; Izumi, F. VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44, 1272–1276. https://doi.org/10.1107/S0021889811038970. Nishihara S. burai | ABOUT http://nisihara.wixsite.com/burai/about (accessed Oct 11, 2018). Reed, T. M.; Gubbins, K. E. Applied Statistical Mechanics. 1973, 270. Young, B.; Bryan, J. Generating Functions for Colored 3D Young Diagrams and the Donaldson-Thomas Invariants of Orbifolds. Duke Math. J. 2010, 152, 115–153. https://doi.org/10.1215/00127094-2010-009. Zbks, B.; Barreto, F. C. S.; Nogueira, R. A. Thermodynamics of an Eight-Site OrderDisorder Model for Ferroelectrics. Ferroelectrics 1977, 16, 93. https://doi.org/10.1080/00150197708237126. Schmechel, R. Gaussian Disorder Model for High Carrier Densities: Theoretical Aspects and Application to Experiments. Phys. Rev. B - Condens. Matter Mater. Phys. 2002, 66, 1–6. https://doi.org/10.1103/PhysRevB.66.235206. Laurencin, D.; Almora-Barrios, N.; de Leeuw, N. H.; Gervais, C.; Bonhomme, C.; Mauri, F.; Chrzanowski, W.; Knowles, J. C.; Newport, R. J.; Wong, A.; et al. Magnesium Incorporation into Hydroxyapatite. Biomaterials 2011, 32, 1826–1837. https://doi.org/10.1016/j.biomaterials.2010.11.017. Henard, C.; Papadakis, M.; Perrouin, G.; Klein, J.; Heymans, P.; Traon, Y. Le. Bypassing the Combinatorial Explosion: Using Similarity to Generate and Prioritize t-Wise Test Configurations for Software Product Lines. IEEE Trans. Softw. Eng. 2014, 40, 650–670. https://doi.org/10.1109/TSE.2014.2327020. Edwards, R.; Glass, L. Combinatorial Explosion in Model Gene Networks. Chaos 2000, 10, 691–704. https://doi.org/10.1063/1.1286997. Gruber, C.; Ueltschi, D. The Falicov-Kimball Model. ArXivmath-Ph0502041 2005. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087–1092. https://doi.org/10.1063/1.1699114. Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N⋅log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. https://doi.org/10.1063/1.464397. Nestler, F.; Pippig, M.; Potts, D. Fast Ewald Summation Based on NFFT with Mixed Periodicity. J. Comput. Phys. 2015, 285, 280–315. https://doi.org/10.1016/j.jcp.2014.12.052. Bickelhaupt, F. M.; Baerends, E. J. Kohn-Sham Density Functional Theory: Predicting and Understanding Chemistry; Wiley-Blackwell, 2007; pp 1–86. https://doi.org/10.1002/9780470125922.ch1.


50

Page 51 of 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(77) (78) (79)

(80) (81) (82) (83) (84) (85) (86) (87) (88) (89) (90)

(91)

Perdew, J. P.; Levy, M. Physical Content of the Exact Kohn-Sham Orbital Energies: Band Gaps and Derivative Discontinuities. Phys. Rev. Lett. 1983, 51, 1884–1887. https://doi.org/10.1103/PhysRevLett.51.1884. Chai, J. Da; Weeks, J. D. Orbital-Free Density Functional Theory: Kinetic Potentials and Ab Initio Local Pseudopotentials. Phys. Rev. B - Condens. Matter Mater. Phys. 2007, 75. https://doi.org/10.1103/PhysRevB.75.205122. Saracibar, A.; Carrasco, J.; Saurel, D.; Galceran, M.; Acebedo, B.; Anne, H.; Lepoitevin, M.; Rojo, T.; Casas Cabanas, M. Investigation of Sodium Insertion–Extraction in Olivine Na x FePO 4 (0 ≤ x ≤ 1) Using First-Principles Calculations. Phys Chem Chem Phys 2016, 18, 13045–13051. https://doi.org/10.1039/C6CP00762G. Trasatti, S. The Absolute Electrode Potential: An Explanatory Note: (Recommendations 1986). https://doi.org/10.1515/iupac.58.0011. Zhou, F.; Cococcioni, M.; Marianetti, C. A.; Morgan, D.; Ceder, G. First-Principles Prediction of Redox Potentials in Transition-Metal Compounds with LDA + U. Phys. Rev. B 2004, 70. https://doi.org/10.1103/PhysRevB.70.235121. Ohkoshi, S.; Iyoda, T.; Fujishima, A.; Hashimoto, K. Magnetic Properties of Mixed FerroFerrimagnets Composed of Prussian Blue Analogs. Phys. Rev. B 1997, 56, 11642– 11652. https://doi.org/10.1103/PhysRevB.56.11642. Clark, S. J.; Robertson, J.; Lany, S.; Zunger, A. Intrinsic Defects in ZnO Calculated by Screened Exchange and Hybrid Density Functionals. Phys. Rev. B 2010, 81. https://doi.org/10.1103/PhysRevB.81.115311. Shi, S.; Liu, L.; Ouyang, C.; Wang, D.; Wang, Z.; Chen, L.; Huang, X. Enhancement of Electronic Conductivity of LiFePO 4 by Cr Doping and Its Identification by First-Principles Calculations. Phys. Rev. B 2003, 68. https://doi.org/10.1103/PhysRevB.68.195108. Shalders, R. D.; Swaddle, T. W. The Ligand Cage Effect in Volumes of Activation for Electron Transfer Reactions of Cobalt(III/II) Complexes in Aqueous Solution. Inorg. Chem. 1995, 34, 4815–4820. https://doi.org/10.1021/ic00123a015. Liu, Z.; Huang, X.; Wang, D. First-Principle Investigations of N Doping in LiFePO4. Solid State Commun. 2008, 147 (11), 505–509. https://doi.org/10.1016/j.ssc.2008.06.013. Jing, Y.; Zhou, Z.; Cabrera, C. R.; Chen, Z. Metallic VS 2 Monolayer: A Promising 2D Anode Material for Lithium Ion Batteries. J. Phys. Chem. C 2013, 117, 25409–25413. https://doi.org/10.1021/jp410969u. Park, M.; Zhang, X.; Chung, M.; Less, G. B.; Sastry, A. M. A Review of Conduction Phenomena in Li-Ion Batteries. J. Power Sources 2010, 195, 7904–7929. https://doi.org/10.1016/j.jpowsour.2010.06.060. Liu, Z.; Pulletikurthi, G.; Endres, F. A Prussian Blue/Zinc Secondary Battery with a BioIonic Liquid–Water Mixture as Electrolyte. ACS Appl. Mater. Interfaces 2016, 8, 12158– 12164. https://doi.org/10.1021/acsami.6b01592. Lee, J.-H.; Ali, G.; Kim, D. H.; Chung, K. Y. Metal-Organic Framework Cathodes Based on a Vanadium Hexacyanoferrate Prussian Blue Analogue for High-Performance Aqueous Rechargeable Batteries. Adv. Energy Mater. 2017, 7, 1601491. https://doi.org/10.1002/aenm.201601491. Wang, R. Y.; Wessells, C. D.; Huggins, R. A.; Cui, Y. Highly Reversible Open Framework Nanoscale Electrodes for Divalent Ion Batteries. Nano Lett. 2013, 13, 5748–5752. https://doi.org/10.1021/nl403669a.


51


Page 52 of 52

Table of Contents Graphic


52

Intercalated Cation Disorder in Prussian Blue Analogues: First

Recommend Documents