UvA-DARE (Digital Academic Repository) Atomic-Scale Design of Anode Materials for Alkali Metal (Li/Na/K)-Ion Batteries Progress and Perspectives

a variety of renewable energy technologies (based on solar, hydro, wind, geothermal power, etc.), 2) provide power sources for electric and hybrid vehicles for low-carbon or zero-carbon emissions of transportation systems, [3,4] and 3) provide power sources for various portable and wearable electronic devices. Alkali metal (AM) ion batteries (AMIBs) including lithium (Li)-ion batteries (LIBs), sodium (Na)-ion batteries (NIBs), and potassium (K)-ion batteries (KIBs) are important rechargeable battery technologies to support the decarboniza-tion of both electricity supply and transportation systems. Currently dominating the rechargeable battery market is LIBs. NIBs and KIBs are developed as more sustainable alternatives and indeed comple-ments, to mitigate challenges associated with the limited and geopolitically isolated Li resources. [1] Post-alkali metal ion batteries are also pursued such as multivalent ion batteries and lithium sulfur batteries. These battery technologies will not be covered in this review, but have been reviewed elsewhere. The three AMIBs the same general cell operation, with in material During charge/discharge the AM ions via an electrolyte between the cathode (positive electrode) the anode rials Computational modeling cannot only provide fundamental insight, but also present input for multiscale models and experimental synthesis, often where quantities cannot readily be obtained by other means. In this review, an overview of three main anode classes; alloying, conversion, and intercalation-type anodes, is provided and how atomic scale modeling is used to understand and optimize these materials for applications in lithium-, sodium-, and potassium-ion batteries. In the last part of this review, a novel type of anode materials that are largely predicted from density functional theory simulations is presented. These 2D materials are currently in their early stages of development and are only expected to gain in importance in the years to come, both within the battery field and beyond, highlighting the ability of atomic scale materials design. cation) compared to the single valent AMIBs, and have also been suggested to have improved safety when coupling the metal anodes with liquid electrolytes. Whilst these battery technologies remain less well studied than their AMIB counterparts, computational studies are excellently poised to translate the same materials design strategy used when adopting LIB anodes to NIBs and KIBs to the multivalent metal-ion battery anodes. Finally, as evident from this review (and not surprising considering the timeline of technology development) KIB anode materials are less studied computationally than NIB and LIB anodes. Hence, there are opportunities for computational studies to screen, validate, optimize, and modify the currently investigate LIB and NIB anode materials for KIBs, and to predict novel ones. The anode materials currently studied computationally for AMIBs should further be evaluated for use in alternative technologies such as divalent batteries, sulfur batteries, air batteries, etc. The current wealth of computational modeling of anode materials have shown to be invaluable for materials design and should continue to be incorporated with experiments based on their success in materials optimization, materials discovery, and materials characterization, and will play an integral part in the future work by also considering nanostruc-tured materials, heterogeneities, composites, and compatibility with electrolytes. Development of a screening methodology with not only anode-specific properties, but also identifying anode properties and materials that will be stable with electrolytes and under different operating conditions will allow the practical advances to accelerate by removing the number of anode-elec-trolyte couples needed to be explored in a lab.


Introduction
The development of sustainable and clean energy sources that are viable alternatives to fossil fuel technologies is one of the most pressing issues in today's society. [1,2] Rechargeable batteries are promising clean energy technologies which can 1) serve as storage sources at local and grid-scale for wide deployment of electrode). [15] A fundamental understanding of the chemical processes governing the AM ion storage and diffusion [16][17][18][19][20] materials is required both for the cathode, electrolyte, and anode materials. Cathode materials and electrolytes have generally been more successfully transferable between LIBs and other AMIBs, than the anode materials. Furthermore, KIBs have not received as much attention as NIBs and LIBs, and there are currently fewer materials that can accommodate the larger K + ions. [21] For the optimization and development of next-generation battery technologies, the development of high capacity anode materials is of utmost importance. [15,22] To this end, computational modeling provides an atomic scale approach to optimize anode materials. For optimum anode performance, the material should have high AM ion storage capacity, rapid ionic conductivity, low voltage, good electronic conductivity, and high phase stability without large volumetric changes during cycling. [22][23][24] Metallic Li and Na anodes were at the beginning of LIB and NIB research considered, but were found to be unsafe for battery applications due to dendrite formation from the reaction between the metallic anode and organic electrolyte solvent molecules. [24] Potassium metal anodes were tested in the early studies of KIBs, but as for Li and Na metal with liquid electrolytes dendrite formation also limits its KIB anode performance. [25][26][27][28][29][30] To design appropriate and high performance anode materials computational modeling has been vital. The challenges and current state of cathode, [31][32][33][34][35][36][37] electrolytes, [38][39][40][41][42][43][44][45] and the solid electrolyte interphase (SEI) [28,[46][47][48] have been covered in detail elsewhere. [27,28,[31][32][33][48][49][50][51][52][53] For further insights on larger scale modeling and battery technologies, the reader is referred to the following review articles. [54][55][56][57] Computational modeling provides unique insight into the atomic scale mechanisms that govern these properties and is an invaluable tool in materials design. The main aim and strength of computational modeling for efficient materials design is its ability to complement and support experimental results and analysis, to unravel important atomic scale properties and features to gain fundamental understanding of properties and mechanisms that are challenging to obtain solely from experimental analysis and measurements. This allows for intelligent materials design via the systematic optimization of anode materials in conjunction with experimental testing, and refinement. In this review, we focus on the atomic scale modeling of AMIB anode materials as detailed in Section 2. AMIB anode materials follow three main mechanisms: intercalation, alloying, and conversion type (Figure 1). [13,[59][60][61] Graphite, an intercalation type anode, will be discussed in the wider context of carbon-based anode materials. The AM ion storage mechanisms in carbon materials such as soft carbon, hard carbon, and graphene are more complex than intercalation, also including adsorption (on defective surfaces), and pore filling (typical hard and soft carbons with amorphous porous regions), in addition to the intercalation in (typically crystalline) interlayer regions. Hence, in this article, we discuss carbon-based (Section 3) and noncarbon based intercalation-type (Section 4) anodes separately, dividing these into two classes due to the large research effort on carbon-based anode materials. After the intercalation-type anodes, the alloying (Section 5) and conversion-type (Section 6) are presented, with 2D materials as novel anode materials for AMIBs being discussed in Section 7.

Computational Methods
Atomic scale computational studies are widely used within the battery modeling community and beyond, forming a synergistic relationship with experimental studies for materials design. [23,62,63] From the simulations of atomic level structures one can obtain both bulk and surface properties. [64,65] A wide range of anode materials properties (e.g., voltage profiles, electronic structures, ionic conduction and migration, defects formation, phase stability, and adsorption of chemical species) can be probed to understand the intercalation, alloying, conversion, and charge/discharge behavior of Li, Na, and K in anode materials. Atomic scale computational methods such as molecular dynamics (MD), and density functional theory (DFT) have been employed for obtaining these materials properties as shown in Figure 2. DFT is a well-utilized and proven electronic structure method that calculates the ground state energy of a many-body Figure 1. Schematic picture of an alkali metal ion (pink spheres) battery, utilizing a layered cathode material, and a generic anode material (represented by a green cuboid). Possible anode materials have been highlighted on the right. Price data for Li, Na, and K refers to the price per ton of their respective ores, with data. [58] system from its electron density. [66][67][68] DFT simulations are typically conducted at 0 K and as snapshots. [32,69] MD simulations uses Newton's equations of motion to model the time evolution of interacting particles by repeatedly integrating many time steps to give a detailed picture of the time evolution, [66][67][68] and are reliant on semiempirical interatomic potentials to describe the interactions between ions in the system, and cannot be used to simulate electronic structures. [68,70,71] Recently, a combination of DFT and MD, ab initio molecular dynamics (AIMD) has gained traction. In AIMD simulations, the interatomic potentials are replaced by forces from DFT to simulate the interatomic interactions, with MD steps to model the time evolution. AIMD is however more computationally expensive than both DFT and MD and is mainly used to model systems at short time scales (≈10 ps as compared to MD where trajectories are normally >1 ns) due to high computational cost. For anode materials, AIMD simulations are becoming increasingly important, and are used to calculate diffusion, thermal expansion, and phase stability. [22,[72][73][74][75][76][77][78][79][80][81] In the next sections, we will present short overviews of what properties can be obtained computationally from these methods, and their relevance to anode materials design. We will also briefly address assumptions and approximations that need to be made when conducting these simulations, such as how the electrons are described, and what implications and errors can arise from this. Both DFT and MD methods have been described in great detail elsewhere, [80,82,83] and hence we will only cover their applications to anode materials here.
DFT simulations are reliant on approximations to the exchange-correlation energy, [82,84] which introduces inherent inaccuracies. One must be aware of such inaccuracies when conducting DFT calculations. The exchange correlation functionals neglect long-range dispersion and van der Waals (vdW) interactions, which are especially important when modeling layered and organic molecular materials. [85][86][87][88][89] For the lack of long-range vdW interactions, a number of corrections can be applied such as vdW functionals (e.g., vdW-optPBE, vdW-DF, and vdW-DF2), and semiempirical corrections such as DFT-D, DFT-D3, and DFT-D3BJ. [87,88,[90][91][92][93][94][95][96][97] The vdW functionals get the dispersion interactions directly from the electron density, whereas the semiempirical corrections include long-range dispersion interactions with coefficients for each atomic species pair. [87,96] The inclusion of these corrections in modeling materials such as graphite, graphene, and phosphorene have shown to improve experimental agreement of structural parameters, interlayer binding energies, formation energies, and adsorption energies. [87,96,97] As described above, the modeling field is diverse, which is further demonstrated in the range of simulation codes that are available. Popular DFT codes include VASP, [98] CP2K, [99] Gaussian, [100] CRYSTAL, [101] CASTEP, [102] DMol3, [103] Wien2K, [104] and Quantum Espresso. [105] Common MD codes include LAMMPS, [106] GROMACS, [107] NAMD, [108] and DL_POLY. [109] The choice of code depends on the material being studied, required information, and code capabilities (such as exchangecorrelation functional, interatomic potential, periodicity, etc.), which can be obtained from each code's website as listed in references. The choice of code should be conducted carefully to avoid implementation errors. [64] It is important to note that not all DFT codes are equally suitable for all applications. In choosing DFT simulation code, one also chooses how electrons are to be treated, which electrons will be explicitly treated (i.e., pseudopotentials or all-electron codes), how the basis sets are to be treated (localized or delocalized), and how to describe the periodicity of the investigated system. [64] Hence, the user must carefully choose the code for their challenges to ensure that the code is best suited for their simulations.

Properties from Simulations
From DFT calculations of the optimized geometrical and electronic structures, properties such as metal intercalation, metal adsorption, voltage, metal migration energies, phase stability, and defect formation energies can be obtained ( Figure 2). A high-performance AMIB anode material should have both high storage capacity of metal ions, high stability during charge/discharge, and allow for rapid metal ion diffusion. [22,24] Metal intercalation-The intercalation energy (E int ), which is related to both the storage capacity and voltage (potential) of an anode material, can be calculated according to Equation (1) [86,110] Metal adsorption-The metal adsorption (or binding) energy (E ads ) is calculated in a similar way to the intercalation energy, comparing the total energies of the system without the adsorbate (E Surface ) and the isolated adsorbate (μ AM ), to the system with the adsorbate (E AM@Surface ) [112] To this end, metal ion adsorption and metal ion intercalation in anode materials have been calculated for a variety of systems to show energetically favorable storage sites, and to evaluate defect and/or dopant engineering in battery materials. [113][114][115][116][117][118] Migration from DFT-Metal-ion diffusion can be calculated and understood from DFT-derived metal ion migration barriers or from MD trajectories. These can then be related to experimentally measured metal-ion migration barriers and diffusion coefficients, obtainable from galvanostatic intermittent titration technique, or muon spin rotation spectroscopy. This can then be used to efficiently screen different materials design strategies (such as dopant concentrations, crystal structure, composition, etc.) to optimize metal diffusion. [22,[119][120][121] Diffusion or metal migration behavior from DFT simulations are commonly obtained via the nudged elastic band (NEB) method. [22,32,59,[122][123][124] NEB simulations require optimized start and end points (typically equivalent lattice sites a set distance apart). Between the start and end points, an initial guess of a migration path is made. The migration path is then divided into a number of steps or reaction coordinates (i.e., intermediate ion positions of the diffusing species) which are optimized (with spring forces acting as constraints to ensure the spacing between the NEB images is retained) within the DFT framework to identify the saddle point or transition state (the individual steps can also be referred to as NEB images). [124] The migration energy barrier (E m ), also referred to as diffusion barrier, is then calculated as the difference in total energy between the saddle point and the start point. [124] The E m obtained for a specific diffusion path from NEB simulations can then be used to obtain the ionic diffusion coefficient (D i ) through an Arrhenius equation (see next paragraph). By repeating the NEB simulation for different diffusion paths, both E m and D i can in theory be mapped for all possible diffusion paths. However, as with the above-mentioned properties, NEB-derived migration barriers are based on DFT total energies, and do not typically take into account the effect of temperature or pressure. Furthermore, whilst metal migration barriers from DFT simulations have been found reliable and in good agreement with experiments, they are computationally expensive and challenging. Especially for complex migration mechanisms, a large number of migration pathways must be sampled, and a large number of NEB images included. For these cases, MD and AIMD simulations are an excellent complement to DFT simulations.
Diffusion coefficients from dynamics simulations-MD and AIMD simulations can be employed to calculate diffusion coefficients and mechanisms by calculating the mean square displacement and tracking the atom trajectory. The mean square displacement (〈r 2 (t)〉) is given by [22,125] where t is time, D i is the self-diffusion coefficient for species i, and B i is a thermal factor associated with atomic vibrations. D i can be obtained by plotting 〈r 2 (t)〉 versus t and can consequently be used to calculate the diffusion activation energy (E a , similar to the migration barrier one obtains from DFT NEB simulations) from Equation (4) exp / ln ln Here, D 0 is a temperature independent pre-exponential term, T is the temperature, and k is Boltzmann's constant. E a is commonly obtained by Arrhenius plots from the gradient of lnD i versus 1/T. However, for applications dependent on the knowledge of the electronic structure of materials, a method such as DFT needs to be employed.
Phase stability-For an anode material to perform well, it also needs to be stable, both on its own and in its lithiated, sodiated, or potassiated phases. [22,126,127] Typically, the structures are generated from structural databases, ab initio random structure searching (AIRSS), or the USPEX code. [128][129][130][131][132][133] Phase stability simulations are vital for the identification and prediction of intermediate stable phases for anodes where the AM ion concentration would be expected to be varied. The stability of different AM site occupations, concentrations, and configurations can be evaluated below [22,128,131] where E s is the formation energy of the investigated AM x Bulk composition, AM Bulk E x is the total energy of the anode material with x number of AM, E Bulk is the total energy of the nonlithiated/sodiated/potassiated anode material, and E AMBulk is the lithiated/sodiated/potassiated anode material. From these relative energies, convex hulls can be constructed by plotting E s as a function of x, to assess the thermodynamic stability of phases, either in a predictive capacity or to confirm experimentally observed phases. [130] Voltage-The charge/discharge profiles relate to different storage and diffusion mechanisms at different voltages. [16] The different metal behaviors in these voltage regions are commonly assigned to the presence of defects, different microstructures, and/or the formation of different alloy compositions. [134][135][136] The electrochemical performance of a battery electrode is commonly measured and presented using a galvanostatic charge/ discharge profile, with voltage (or potential) versus capacity during cycling. These voltage or potential profiles have been simulated with DFT calculations, showing good agreement with experimental anode studies. [16,[137][138][139] These studies have calculated the theoretical capacity (C) from atomic mass, [135,140] or from the metal ion concentration (Equation (6)), [138,141] whereas the equilibrium voltage (V) or potential is calculated from the DFT total energies (Equation (7)) [22,140] AM,max Adv. Energy Mater. 2022, 12,2200662 Graphite, an intercalation-type anode material, is the most commonly used LIB anode material due to its high specific capacity (theoretically 372 mAh g −1 ), [146] good recycling performance, and low operation potential (0.1 V vs Li + /Li). [15,111] Nevertheless, in the development of post-Li battery technologies it is clear that graphite cannot be universally applied. [111,147] Unlike Li, Na does not readily intercalate into graphite. The replacement of Li with larger ions such as K leads to structural changes in graphite. [111] These structural changes in turn affect the storage and metal diffusion behavior, making a direct swap of Li to other battery metals challenging. [111,148] To this end, the disordered carbonaceous anode materials, hard and soft carbon, have shown promising performance, not only for NIBs and KIBs but also for LIBs. [149,150] Hard carbon is a nongraphitizable carbon material with a disordered structure consisting of randomly stacked graphitic domains, with expanded graphitic interlayers, and nanovoids. [151] This leads to hard carbon structurally laying somewhere between crystalline and amorphous, with local order (graphitic segments), but no long-range order. These materials commonly have high defect concentrations, and their structure can be tuned depending on the synthesis technique, precursor, and post-treatment used. [152] Soft carbon can be converted to graphite if annealed at high temperatures, and have less curvature and more ordered alignment of its graphene layers than hard carbon. [16,17,60,153] From experimental evidence from various probing techniques such as Raman spectroscopy, transmission electron microscopy, nuclear magnetic resonance, electron energy-loss spectroscopy, and pair distribution functions from diffraction data, different structural orders for different carbon-based anodes are shown, which adds important insights for guiding the computational models in Figure 3. [63,154]

Computational Models of Carbon Anodes
To study the AM storage, intercalation, adsorption, and migration in carbon-based anodes, computational studies typically utilize monolayer graphene ( Figure 3a) and bilayer (or multilayer) graphene ( Figure 3b) to represent carbon surface and the expanded graphitic layer stacks, respectively. These models apply not only to hard and soft carbons, but also other graphitic anode materials (such as composites), and graphite itself. Graphene is a 2D material with sp 2 -hybridized carbon atoms in a hexagonal lattice. [112,155] The graphite lattice consists of alternate graphene planes with either hexagonal AB, or rhombohedral ABC stacking that are weakly bonded in the vertical direction, but have strong covalent bonds in the graphene planes. [87,110,156] The in-plane carbon atoms are bonded by σ bonds (sp 2 -hybridized carbon) and delocalized π bonds (from the carbon 2p-orbitals), and to accurately model both the graphite structure and metal interaction with graphene and graphite, a variety of dispersion and vdW corrections have been incorporated in the DFT models. [87,110,113,156] To capture the curved morphologies present in hard carbon, carbon nanotubes ( Figure 3c) have been employed. Graphite can be reasonable represented by the monolayer graphene and the bilayer graphene models. Due to the complex nature of the disordered structures in hard and soft carbons, and their structural dependency of precursors, computational models of these carbons focus on the modeling of select motifs, such as using the monolayer graphene, bilayer graphene, and carbon nanotube models. These models are not universal for disordered carbon anodes, and still cannot be modeled at representative length scales without the development of more powerful computers and/or efficient codes.
MD simulations have been used to construct disordered carbon models (Figure 3c-e), typically creating the disordered carbon structures from melt-quenching studies, and annealing predisordered materials. [47,[157][158][159][160][161] Due to the large simulation cells that would be required to capture the complex structure of these disordered and nanoporous carbon anodes, DFT simulations are too computationally expensive to model accurate systems, and MD simulations (which rely on interatomic potentials, often empirically fitted) are unable to describe and capture the different local structure environments and bonding in a hard carbon material. [63] Practically, these anodes do not have the same structure and are sensitive to synthesis techniques and precursors, making the modeling of "one" unified hard carbon model highly challenging.
Glassy carbon (Figure 3c) consists of sp, sp 2 , and sp 3 hybridized carbon, with turbostratic densely packed carbon particles, but does not show the same layering as the previously discussed disordered carbon anodes such as hard and soft carbons. [139,154,162] Depending on external conditions, such as pressure, the sp 3 -character of the glassy carbon can be varied (as shown for glassy carbon models generated using a stochastic quenching method and simulated using DFT). [163] Adv. Energy Mater. 2022, 12,2200662 2200662 (6 of 29) www.advancedsciencenews.com

Metal-Ion Intercalation and Storage in Graphite Anodes
A high-performing anode material must have good storage capacity, and a multitude of storage (also referred to as adsorption at surfaces, or binding or intercalation in bulk) sites available for metal ions. One of the most popular carbon-based anode computational models in the DFT literature is graphene. Graphene is found as a structural motif in a plethora of anode materials, and has also been studied as an anode material in its own right, [135,[164][165][166][167][168][169] making it an important system for battery modelers. Computational studies of metal storage through the calculation of adsorption energies on graphene monolayers bring fundamental insight into the mechanisms of carbonaceous anode materials. Pristine graphene provides three inequivalent sites for metal adsorption; hole site (above the C 6 -ring), bridge site (above a CC bond), and top site (above a C atom). [113,170] DFT studies have found that the majority of elements in the periodic table preferably adsorb at the hole site. [113,114,171] Among the battery metals, Li (−1.22 eV from DFT simulations with D3-BJ dispersion corrections) has the strongest adsorption energy (according to Equation (2)) on graphene, followed by K (−1.05 eV from DFT simulations with D3-BJ dispersion corrections), and Na (−0.62 eV from DFT simulations with D3-BJ dispersion corrections). [113,117,135,170,172] The exact values of these adsorption energies are however dependent on the choice of dispersion correction, making it crucial for computational studies to divulge the use or lack of dispersion correction. The DFT adsorption energies confirm the experimentally observed poorer Na storage on graphene and graphite surfaces, and the computational metal adsorption calculations have also shown that the relative adsorption strength of metals on graphene is governed by the metal ionization potential, which can have implications in the experimental cycling behaviors of these anode materials. [113,114,173] Graphite models (typically multilayer graphene (Figure 3b) have also been employed in DFT simulations to predict the Diagram showing in the center hexagon schematic images of hard carbon (using the falling cards model), graphite (with light and darker gray signaling alternating graphene layers), graphene, and soft carbon (assuming a graphene platelet structure). The surrounding hexagons display different carbon models used in the literature to model AMIB anodes where a) is graphene, b) bilayer graphite, c) carbon nanotube, d) glassy carbon with inserted lithium/sodium (green). Reproduced with permission. [139] Copyright 2015, American Chemical Society. e) A disordered graphitic system obtained by combining machine learning and experimental data. Reproduced with permission. [137] Copyright 2013, Royal Society of Chemistry. f) A disordered carbon model with curved pore structure from DFT-MD annealing and cooling. Reproduced under terms of the CC-BY license. [63] Copyright 2018, The authors, published by Royal Society of Chemistry. (a-d) Brown spheres are carbon, and (e,f) carbon is gray and turquoise spheres, respectively.
Adv. Energy Mater. 2022, 12, 2200662 2200662 (7 of 29) www.advancedsciencenews.com metal storage capacity of graphite anodes, and hard and soft carbon anodes. Graphite forms high stage intercalation compounds with lithium and potassium, but not with sodium (Figure 4b). [111] The theoretical capacity of lithium-graphite intercalation compounds (Li-GICs) is ≈370 mAh g −1 , which corresponds to LiC 6 . [111] The maximum capacity of sodiumgraphite intercalation compounds (Na-GICs) on the other hand is much smaller at 35 mAh g −1 (NaC 64 ). [16,174] This unfavorable sodiation of graphite anodes further leads to sodium plating and consequent dendrite formation, rendering graphite an unsuitable NIB anode material. [175] Potassium GICs (K-GICs) have a theoretical capacity of ≈280 mAh g −1 , corresponding to a concentration of KC 8 . [176,177] The GICs are formed when metal ions intercalate between the graphene planes in the interlayer region ( Figure 4a). [111,178] This leads to the breaking of bonds, and the formation of new ones, with the structure of the GICs being highly dependent on the intercalant (Figure 4b). [111,178] What is common among the different graphite intercalation compounds however is their ionic and electronic conductivity, and their layered structure giving a 2D conductivity as the ions (in defectfree graphite) primarily move within the interlayer, and not through the graphene layers. [111,178] Hence, battery performance issues that lead to anode failure associated with large volumetric changes during cycling in alloy based anodes are not an issue for graphite anodes. [179][180][181] The limitation of graphite as an anode material for post-lithium technologies instead comes from the concentration of metal ions possible to store in graphite. [111] Apart from predicting metal storage capacity, DFT simulations have also been applied to understand the intercalation mechanism, and stability of the carbon structures upon intercalation, using periodic graphite models with two (bilayer model) to four graphene layers. [182][183][184][185][186] Metal intercalation mechanism in graphite (or few-layer graphene model) have been studied for various metal concentrations ( Figure 4). Computational studies of GICs have been widely focused on the difference in intercalation mechanism and stability of Li-GICs, Na-GICs, and K-GICs. [110,111,187] It has been shown that the difference in intercalation behavior between the different AMs is due to the sum of electrostatic and covalent binding contributions for Li-GICs being larger than those for Na-GICs and K-GICs. [85,110,187] The high stability of Li-GICs from formation energy calculations suggests that the small ionic size of lithium leads to small structural deformation in the graphite lattice, giving that only a small structural deformation energy needs to be overcome in Li-GICs. [96,184,188] This small structural deformation in the graphite lattice upon lithium intercalation further leads to the vdW interactions in between the monolayers stabilizing the Li-GICs. [85,110,189] High concentration Na-GICs are from DFT simulations found to be thermodynamically unstable (Figure 4b), which has been attributed to the relatively high redox potential of Na/Na + compared to Li/Li + . [85,190,191] Intercalation of K into graphite leads to larger structural deformations, with the stability of the K-GICs being dependent on these structural deformations. The high capacity is a consequence of the high and exergonic ionization energy and relatively small cohesive energy (Figure 4c,d). [192] The structural deformation energy of the intercalated graphite lattice increases linearly with the size of the AM ion, with the main energy contribution to the thermodynamic stability arising from the expansion of the graphite lattice upon metal intercalation. [85,110,187] To this end, computational studies have been able to directly identify the reason for the seemingly unintuitive thermodynamic stability between . a) Schematic representation of Li graphite intercalation compounds, and b) the corresponding formation energies for Li, Na, and K graphite intercalation compounds from DFT simulations without dispersion corrections. Reproduced with permission. [193] Copyright 2013, Elsevier. c) the different stages of K (blue spheres) intercalation in graphite (yellow) as calculated at DFT-D2 level, with the corresponding calculated potential (voltage) profile for K at different concentrations together with experimentally measured data in d). (c,d) Adapted with permission. [192] Copyright 2015, American Chemical Society.
Na-GICs and K-GICs, proving that the inability of Na to intercalate in graphite is not a result of its larger ionic size than that of Li.

Solvent Co-Intercalation in Graphite
Apart from metal ion storage in graphite, alkali metal atoms can also co-intercalate with solvent molecules. [125,174,[194][195][196][197][198][199] Co-intercalation of metals with solvent molecules in graphite have been shown to improve the intercalation properties not only for Na, [60,[200][201][202][203][204][205] but also for Li [206][207][208] and K. [188] Solvent co-intercalation with Na in graphite (neglecting electrolyte salt and SEI) was shown through dispersion corrected (D3) DFT simulations to be solvent dependent ( Figure 5). [209] Both propylene carbonate (PC) and diglyme form stable initial graphite intercalation compounds with Na, while tetrahydrofuran (THF) and dimethyl carbonate (DMC) (which have weaker solvation energies with Na) do not intercalate in graphite. The high LUMO level (obtained through nonperiodic B3LYP calculations) of diglyme was shown to lead to reversible co-intercalation in graphite (also confirmed by periodic DFT simulations at PBE-D3 level of theory where a complex of two diglyme and one Na ion can intercalate reversibly in expanded graphite with interlayer distance 11.3 Å [202] ), whereas the lower LUMO level of PC would lead to exfoliation of the co-intercalated graphite ( Figure 5). [209] A similar conclusion was made in a later study, where PC was deemed unsuitable for LIBs due to its co-intercalation with Li in graphite leading to exfoliation of the graphite anode, although it had been shown to perform well with hard carbon anodes in NIBs. [52] To form an atomic scale understanding of exfoliation of graphite anodes due to PC co-intercalation (and why this is not an issue for ethylene carbonate (EC)), MD simulations have been employed. [210] These simulations show that the exfoliation of graphene sheets seen with PC cointercalation could be a result of the larger molecular size of PC, as compared to EC, inducing a larger interlayer separation of the graphene sheets. This leads to an 11 times lower barrier to exfoliation for bilayer graphene with intercalated PC than with EC. [210] These studies show the importance of electrolyte in battery performance, and the inclusion of electrolyte in the modeling system for optimizing anode properties.

Effect of Pore Shape and Size on Adsorption and Storage
DFT and MD simulations have shown that the pore structure of carbonaceous anodes directly influences the anode performance. [142,[211][212][213][214][215] From experimental characterization, both planar and curved pore motifs are present in nanoporous hard carbon. The exact contribution to anode performance of different pore shapes and size can be probed from computational studies, allowing for pore engineering to be conducted at an experimental stage. To this end, bilayer graphene and carbon nanotube models (Figure 3b,c) with variable interlayer distances and tube diameters have been used to computationally investigate metal storage. From DFT calculations of Li, Na, and K adsorption energies at both the concave and convex surface of carbon nanotubes (Figure 3c with diameter 1-5 nm), it is shown that these metals have stronger binding in bilayer graphitic pores (using model in Figure 3b) than to the curved carbon nanotubes. [142] Atomic scale studies of different hard carbon motifs have confirmed that metal storage at defect sites in carbon-based anodes can occur, and that the pore shape of nanoporous carbons directly impact metal incorporation. [97,113,121,137,142,216] The suggested Na storage mechanism by pore filling was confirmed through combined experimental Figure 5. Solvent dependency of Na solvent co-intercalation in graphite based on DFT simulations showing the influence of solvation energy and electronic structure on performance, either not intercalating, or leading to intercalation (and the expansion of graphite interlayer distance), with subsequent reversible intercalation or exfoliation. The graphite structures were simulated at DFT-D3 level and the Na solvent interactions to calculate the Na-solvent complex interactions and their HOMO and LUMO levels were simulated at B3LYP/6-311G level. Purple spheres are Na, red O, gray C in solvent molecule, white H, and brown C in graphite. Reproduced with permission. [209] Copyright 2017, Wiley-VCH. Molecule labels have been added with PC being propylene carbonate, DEGDME diethylene glycol dimethyl ether, DMC dimethyl carbonate, THF tetrahydrofuran, and DOL dioxolane. and DFT studies. [142,143] Apart from a pore filling mechanism, surface defects and intercalation in expanded graphite stacks are also responsible for the experimentally observed metal potentials. [121,134,142,156,217] Na storage in small planar graphitic stacks, wide planar graphitic pores, and cylindrical pores (corresponding to the curved morphologies observed in real samples) were investigated using DFT simulations with D3-BJ corrections, showing that both the morphology and size of pores directly influences Na binding energies and migration. The inclusion of Li and K in these studies further shows that the anode performance is not only dependent on morphology, but also on the battery metal ions, as previously observed for metals on graphene. [113,114,137,214,215] Glassy carbon has been employed both as a standalone anode material, and as composites. [146] Glassy carbon structures were simulated using DFT from structures generated from random sampling. [139] From these simulations, it was deduced that Na storage in the two glassy carbon models was more stable than Na storage in graphite, with the carbon amorphitization having a thermodynamically favorable effect on Na insertion. [139] The insertion of Na into glassy carbon was however found to be structure dependent, with some sites returning positive (i.e., energetically unfavorable) formation energies. [139] Contrarily, the defect formation energy for Li insertion in the glassy carbon models remained negative for all concentrations and sites except one, indicating that this is a thermodynamically stable phase. Hence, in addition to providing more energetically favorable Na insertion, glassy carbon could also provide an alternative LIB anode material, as it is expected to increase the anodic voltage and limit electrolyte decomposition in LIBs. Additionally, it would also be an interesting material to study for KIBs, using atomic scale materials design to fine tune the structure for this specific application. Computational modeling opens the possibilities to investigate all three common battery metals simultaneously, and to systematically study not only the effect of different anode material structure on performance, but also how this would differ between LIBs, NIBs, and KIBs.

Effect of Defects on Metal Adsorption and Storage
A materials design strategy particularly well suited for computational studies where one can directly probe the effect of local changes in atomic structure is the introduction of defects. Defects are commonly present in carbonaceous materials including graphene, graphite, hard carbon, and soft carbon, as a result of synthesis methods, battery cycling, and materials engineering. [15,17,20,137,143,[218][219][220][221][222][223][224] Atomic scale computational modeling opens the possibility to directly probe the influence of single defects on several properties, which can be difficult, if not impossible, to discern experimentally. [225,226] Common graphite and graphene defects investigated include single carbon vacancy, also referred to as monovacancy (MV), double carbon vacancies (DV), and Stone-Wales (SW) defects. The defect formation energies (according to Equation (8)) of these vacancy defects in bilayer graphene have been calculated to be 7.9 eV for MV, 9.2 eV for DV, and 6.9 eV for SW, [183] which are similar to the defect formation energies of the same defects in monolayer graphene (7.5, 8, and 5 eV, respectively). [227] The presence of these defects in bilayer graphene improves the Na intercalation as compared to graphite, making Na intercalation in graphite more energetically favorable in graphitic layers with MV, DV, and SW. [183] Similarly, the adsorption of Na on a graphene sheet can be improved by including these defects, [228] and Li and K adsorption have also been shown to be stronger at graphene sheets with vacancy defects. [113,137] Preparing the graphene-based anode materials with appropriate concentrations of these defects could hence provide a promising materials design strategy to increase AM storage.
Jian et al. employed bilayer pristine and defective (MV, DV, and SW) graphene models to simulate sodium intercalation in soft carbon anodes. [16] These DFT-D3-BJ level simulations showed that the higher formation energy defects MV and DV have a larger effect of sodium incorporation and storage than the SW defect. [16] The reason for this was found in the poor electron transfer to the carbon layer. Furthermore, the inclusion of the vacancy defects in the bilayer graphitic model greatly improved the sodium storage, confirming experimental measurements. [16] The sodium potential versus capacity curves (according to Equations (6) and (7)) were also calculated for sodium storage in hard carbon by Reddy et al. [216] The effect of vacancy defects and interlayer distance in bilayer graphene as a proxy for hard carbon was studied through DFT simulations by Tsai et al., [156] Yang et al., [183] and Yang et al. [182] All these studies conclude that the MV and DV defects are beneficial for sodium storage, and that expanded interlayer distances are required for energetically favorable sodium intercalation. [16,137,156,182] Hence, for higher-performance anode materials, defect engineering should be considered as an important route together with control of the interlayer distance.
Nitrogen-containing defects in carbon-based anode materials have also been shown to increase NIB and KIB capacity. [20,229,230] Chen et al. employed expanded graphitic layer models to study the effect of graphitic and pyridinic nitrogen defects on Na and K adsorption energy. [229] As shown previously, high concentration Na GICs (NaC 8 and NaC 18 ) were found to be thermodynamically unstable, resulting in a positive adsorption energy (0.132 eV), whereas K intercalation was shown to be thermodynamically stable (−0.229 eV). Upon the inclusion of a pyridinic nitrogen defect, the Na adsorption energy became negative (−0.267 eV), indicating a thermodynamically stable system. [229] Phosphorous doped hard carbons using graphitic models functionalized with PO, PO, and PC were simulated at GGA+U level without dispersion or vdW corrections and concluded to have a beneficial influence on Na adsorption. [219] Studies of Na adsorption on P-doped and B-doped graphene monolayers showed that the Na adsorption energy was improved at these defects, as compared to pristine graphene. [218] The B substitutional defect strengthened the Na adsorption energy by 1.4 eV, whereas the P substitutional graphene sheet saw a strengthening of the Na adsorption of 0.854 eV. [218] The computational defect studies have helped to explain the experimentally observed battery performance. For example, a revised mechanistic model for Na insertion in hard carbon has been proposed. [151] Based on the experimental charge/discharge curves, it has been postulated that the metal incorporation mechanism is voltage dependent. The measured sodium potential of the sloping region on the charge/discharge curves has Adv. Energy Mater. 2022, 12,2200662 2200662 (10 of 29) www.advancedsciencenews.com been assigned to sodium adsorption at defect sites and sodium intercalation between graphitic layers; whereas the plateau region can be assigned to pore filling. [17,137,143] Apart from defects in the graphitic layers themselves in GICs, defects within the metallic layers can influence anode performance. Wang et al. used dispersion corrected DFT (optpbe-vdW) to study the effect of metal layer defects in Li-, Na-, and K-GICs (Figure 6). [231] In this study, vacancy, self-interstitials, and Frenkel defects were investigated in the Li, Na, and K layers of AMC 6 , AMC 8 , AM 7 C 48 , AM 7 C 64 , AM 7 C 48 , AM 9 C 48 , and AM 9 C 64 (Figure 6a). Metal vacancies were found to have the lowest defect formation energies (Figure 6b), showing that vacancy defects will dominate over both interstitial and Frenkel defects in first stage AM graphite intercalation compounds. [231] It was addressed that both NaC 6 and NaC 8 are thermodynamically unstable, and that the real concentration in Na-GICs will be lower, but was included to show that the sodium concentration in GICs contribute to the diffusion behavior. [231] A similar study was reported in 2014 by Thinius et al. [189] at DFT and DFT-D3-BJ level showing that the Li vacancy in Li-GICs has a lower defect formation energy than the Li interstitial (0.11-0.21 eV vs 0.60-0.75 eV), indicating that this would be the, under equilibrium conditions, most probable defect in the metal layer. [189]

Metal Migration and Diffusion
Apart from storage capacity (which could be derived from the above computational adsorption and intercalation studies), cycling capacity and charge/discharge rate are also key battery performance indicators. Rapid metal diffusion/mobility is required to provide fast charge/discharge, to have high reversible capacity and limit dendritic growth. It is therefore important to also consider the metal migration and diffusion in anode materials. [113,151,177,[232][233][234][235][236] For LIBs, the graphite anode hinders high-rate charge and discharge, limiting the development of electric vehicles and fast-charging consumer technologies. [237,238] Computational studies may screen and predict candidate systems and materials design theories to reduce charging times by calculating migration and diffusion barriers.
The Li migration barrier on the graphene surface (0.31 eV from NEB simulations with D3-BJ corrections) is higher than both the Na and K migration barriers (≈0.1 eV for both simulated using the same methodology as Li). [113,117,171] The difference in migration energy has been attributed to the difference in metal adsorption energy between the different graphene adsorption sites. The difference in lithium adsorption energy between the hole site and the bridge or top site is much larger than that for Na and K. [113,114,168,173] Metal migration through the graphene layer is highly energetically unfavorable, yielding Li migration barriers >10 eV, and Na migration barriers >30 eV. [239][240][241][242] Hence, it is unlikely that metal migration in carbon-based graphitic anode materials would progress through the graphene layers. [189,243] Metal migration in graphite is therefore predominantly 2D, occurring along the graphene planes. The Li migration barrier in graphite (LiC 72 ) is similar (0.34 eV from DFT simulations with D3-BJ corrections) to that of a Li ion on graphene, and has been shown to decrease (indicating more rapid diffusion) with expanding interlayer distance. [142,238] In fully lithiated, sodiated, or potassiated graphite, metal diffusion could also occur through metal vacancy migration ( Figure 6). [231,238,244] Through NEB calculations (Figure 6c), a curved migration path was established for LiC 6 and NaC 6 , whereas a straight migration path was found to have the lowest migration activation energy (Figure 6e) for NaC 8 and KC 8 . [231] It was shown that Li migration following a Frenkel mechanism has a slightly lower Li migration activation energy (E a = 0.42-0.52 eV) than the Li vacancy migration path (E a = 0.42-0.56 eV). [189,231] Both sets of calculated Li migration barriers are in agreement with the experimentally measured Li diffusion barrier of 0.55 eV ( 7 Li NMR). [189] For K-GICs (ranging from KC 8 to KC 24 ), the vacancy migration path was found to have much lower activation barriers (0.11-1.58 eV) than the Frenkel defect migration mechanism (2.42-7.92 eV), indicating much more rapid K diffusion along the vacancy migration path. [177] Hence, the AM ion does not only affect the structural deformation of the GIC lattice, but also the metal diffusion mechanism. [177] The K migration barriers are furthermore lower than the Li migration barriers, suggesting that K mobility is higher than Li mobility in GICs, following the trend found for Li and K migration on graphene. [113,177] DFT AM migration energy barriers in carbon nanotubes showed rapid AM migration, leading to the conclusion that curved morphologies benefit AM diffusion and pore filling. [142] These simulations were directly used by experimental studies to resolve the sodiation/desodiation of hard carbon anodes, showing the power of a synergistic relationship between computational and experimental studies. [121,151,245,246] MD simulations of organic solvents and sodium ions in carbon nanotubes and in slit carbon nanopores show that the dimensionality of the anode directly influences the sodium ion insertion process. The 1D carbon nanotube led to much less sodium ion insertion than in the 2D slit nanopore, in agreement with the DFT simulations on single metals in different carbon morphologies. [212,214,215] As demonstrated above, DFT-based computational simulations have been instrumental in the development and understanding of carbon-based anode materials. DFT simulations have conclusively deciphered the GIC problem for NIB anodes, and shown why graphite can accommodate high concentrations of both Li and K. Furthermore, computational studies have been able to identify both performance limiting and enhancing defects, dopants, pore structures, and local atomic scale structures in highly complex carbon materials such as the hard and glassy carbons, guiding and accelerating experimental anode materials development. The carbon-based anodes are among the most widely investigated materials for AMIBs, and a plethora of other materials have been developed working on the knowledge gained from the carbon-based materials. In the next sections, we will see that the intercalation behavior can also be taken advantage of by noncarbon anodes.

Insertion-Type Anodes
Intercalation-or insertion-type anodes follow a similar mechanism to that described for graphite anodes. These anodes are typically layered materials such as MXenes, spinel titanates, or niobium oxides. [247][248][249][250] As for graphite, the AM intercalation in the insertion-type anode is dependent on the interlayer distance and surface composition, which provides ample opportunity for computational materials design. MXenes (Figure 7) are a relatively novel class of materials, which have applications in energy storage, photocatalysis, gas sensoring, and optoelectronics among others. [250][251][252] MXenes include metal carbides, nitrides, and carbonitrides, with the common formula TM n+1 AX n (n = 1, 2, or 3, TM transition metal, A is a group XIII-XVI elements, and X is C or N). [125] These materials have high conductivity, good mechanical stability, high charging rate, variable interlayer spacing, and provides ample opportunity for surface engineering. [125,[252][253][254] Ti-, Nb-, and V-based MXenes have been identified as the most promising for anode applications. [125,255] The simulated specific capacity for LIB MXene anodes are highest at 276 mAh g −1 for V 2 CO 2 , but remains lower than that of graphite (372 mAh g −1 ). [256] Spinel titanates such as Li 4 Ti 5 O 12 have been investigated as both LIB and NIB anode materials. [13,22,[257][258][259][260][261] No titanates compounds have been studied/reported for KIB anodes. Li 4 Ti 5 O 12 has experimentally been shown to inhibit Li dendrite formation, to have limited SEI formation, and no structural change, making it a so-called zerostrain insertion material. [13,259,260] Its large bandgap (2-3 eV) have led to different doping schemes, carbon layers, and synthesis methods to improve its electronic conductivity. [259,260,262] Na 2+x Ti 3 O 7 , a 2D layered material, and Na 2+x Ti 6 O 13 , which has an open network 3D tunnel structure, have good thermal stability, cycle life, and low AM + intercalation potentials. Niobiumcontaining oxides have been studied both experimentally and computationally as intercalation-type anodes for applications that require fast ionic diffusion, high-voltage, and high-rate anodes. Niobium-doped anatase (TiO 2 ) was recently shown to be a promising anode material, paving the way for other n-type doped anatase anodes to be investigated.
From DFT studies, two Na intercalation sites in Na 2+x Ti 3 O 7 were identified; 7-, and 9-coordinated oxygen sites. [263,264] Upon sodiation (and indeed also lithiation in the iso-structural Li 2+x Ti 3 O 7 ) this leads to Na 4 Ti 3 O 7 being formed during electrochemical cycling, with all Na (and Li) sites being 6-coordianted to oxygen ions. [263][264][265] Na 2+x Ti 6 O 13 has a different crystal structure, which enables Na to intercalate and store in 3D tunnels. The computed voltages of these two NIB anode materials are 0.37 V (for Li 0.90 V) and 0.75 V, respectively, with the maximum x in Na 2+x Ti 6 O 13 being x = 1. [264,266] Na and Li intercalation voltages in Li 2 Ti 3 O 7 are 0.77 and 1.46 V, respectively, at the same level of theory. [264] For both Nb-doped anatase for LIB and NIB anodes, Na intercalation and insertion was found to be improved in these anodes due to the multiple redox couples (Ti 3+/4+ and Nb 5+/4+ ) induced by the Nb and Ti mixing. [267] The Na + occupy the voids that are formed in between the NbO 6 and Adv. Energy Mater. 2022, 12, 2200662  2200662 (12 of 29) www.advancedsciencenews.com TiO 6 octahedra. Experimental preparation of this material utilized K + to open up these voids. [267] Nevertheless these materials have not been investigated for KIB anodes. Orthorhombic Nb 2 O 5 has been investigated for high-rate LIBs and NIBs. Whilst Li can intercalate in Nb 2 O 5 planes, Na cannot (much like previously seen for graphite) and mainly shows surface or near-surface energy storage. [247] To increase performance, functionalization of the anode surfaces with O, OH, and F groups have been employed commonly to tune their properties, in a similar vein to the materials design of the carbon-based intercalation anodes. [268][269][270][271] The performance of both non-, OH-, and O-functionalized titanium carbide (an MXene) were investigated for NIBs and KIBs using DFT calculations. [268] For the OH-functionalized MXene, the strongest Li, Na, and K adsorption is reached at TM n+1 C n (OH) 2 M 0.5 (TM = Ti, V, or Nb). These metal adsorption energies are positive, indicating that Li, Na, and K adsorption is unfavorable at these surfaces (Figure 7a,b). [268] Li and K shows slightly negative adsorption energies (weaker than −0.1 eV) at V 2 C(OH) 2 and Nb 2 C(OH) 2 surfaces (Figure 7b). In all cases, Na was shown to have the strongest repelling energy (Figure 7b), and the OH-functionalized MXenes can thus from computational studies be shown to be unsuitable for AMIBs. O-functionalized MXenes (TM n+1 C n O 2 M 2 ) (Figure 7c) show negative metal adsorption energy and can accommodate a high concentration of Li, Na, and K (Figure 7d). [268] The strongest Li and Na adsorption was observed for TM = Ti and V, whereas K adsorption was equally strong on the MXene surfaces with TM = Nb (Figure 7d). K adsorption on these systems were however found to be much weaker than Li, and Na adsorption, indicating lower KIB capacity for this anode material. [268] The nonfunctionalized MXene and Nb 2 O 5 anodes were however shown to have both the strongest Li, and Na adsorption, and the highest capacities. [268] K did not show any energetically stable interaction with the nonfunctionalized MXene systems, leading to the conclusion that these MXenes are more promising as LIB and NIB anodes than KIB anodes. Fluorination of the Nb 2 O 5 structure (where the fluorine ions substitute the oxygen sites) was shown from both experimental and computational studies to improve the NIB anode performance. [247] Na adsorption energy calculations showed that Na adsorption is highly energetically favorable on both the fluorinated and nonfluorinated surface, [247] although the strong adsorption energies (stronger than 3 eV) could indicate that these sodium ions remain irreversibly bound. The strongest Na adsorption was obtained on the nonfluorinated surface. [247] A separate study compared the electronic structures of non-, F-, and OH-functionalized Ti 3 C 2 and showed that the nonfunctionalized materials are metallic, whereas the functionalization can (depending on the degree of functionalization) open up a narrow bandgap, [272] which might affect the anode performance. Similarly, sodiation can change the electronic properties from insulating to metallic. [273] DFT simulations showed that a Nb concentration of 0.1-1.0% increased the metallic character of TiO 2 due to the induced shift in Fermi level upon doping, whilst retaining the redox properties of the undoped anatase, [274] whereas the semiconducting Na 2 Ti 3 O 7 with bandgap 2.25 eV becomes metallic upon full sodiation. [273] The role of functionalization is not limited to voltage and surface adsorption though, and can be employed as a materials design tool to also modify mobility. Comparing mobility between the functionalized and nonfunctionalized MXene models for Li and Na, the metal migration barriers were found to be lowered on the nonfunctionalized surfaces, indicating more rapid diffusion. [268] For spinel titanates, the Na ions were shown to follow a zigzag diffusion pathway, enabled by vacancy hopping, with E m between 0.186 and 0.226 eV. [273] Na migration barriers in Na 2+x Ti 6 O 13 is 0.24-0.41 eV from NEB calculations. [275] This study also included Na intercalation in K 2 Ti 6 O 13 , which was found to be more sluggish than Na diffusion in Na 2 Ti 6 O 13 due to geometric features. [275] Li diffusion in Na 2 Ti 6 O 13 was found to be energetically unfavorable with E m more than double that of Na. [275] Similar conclusions were made for niobium-based anodes where higher migration barrier for Na migration in the nonfluorinated Nb 2 O 5 (1.33 eV) than that for the fluorinated (≈1 eV) was obtained. [247] Both these migration barriers are however high, and would be expected to lead to limited Na diffusion at room temperature. Niobium tungsten oxide (Nb X W y O z ) anodes' lithium diffusion pathways can be tuned dependent on the oxygen (z) and metal content (x, y) to change the polyhedra structure and organization. These different shear structures were studied using DFT and AIMD, and it was found that the lithium diffusion had migration barriers of 0.1-0.2 eV, making these anodes having rapid lithium diffusion at room temperature. [276] Nb 2 Ti 2 O 9 was found to have a 0.2 eV Na migration barrier, again indicating rapid room-temperature diffusion. [267] To further understand the implication of this on the electrochemical performance of these anodes, joint experimental and computational studies for the functionalized and nonfunctionalized MXenes should be conducted. Particularly it will be Figure 7. Models for Li, Na, and K adsorption on Ti 2 C(OH) 2 nanosheets, and b) their calculated adsorption energies showing that most adsorption sites in a) are energetically unfavorable. c) The model Ti 2 CO 2 nanosheets, and d) the corresponding metal adsorption energies showing energetically favorable Li, Na, and K adsorption. Adapted with permission. [268] Copyright 2014, American Chemical Society.
The niobium-based anodes are so far mainly investigated for LIBs and NIBs, and have thus far not been studied for KIBs. The larger ionic radius of Na + and K + might hinder the intercalation and diffusion of these metals in the polyhedral networks, but would be worth investigating, as would the effect of defects and dopants on anode performance.

Alloying Anode Materials
Alloying materials have shown great promise as anode materials for AMIBs (Figure 8) due to their high energy density and low voltage operation. [180,[277][278][279][280][281][282][283][284][285][286][287] The modeling of the alloying anodes is typically different to the modeling of the insertion/ intercalation-type materials. These anode materials typically undergo an electrochemical alloying process via the reversible reaction AM B e AM B x y x , where AM + is the alkali metal ion, and B the alloying anode material (Figure 8). [288] During electrochemical cycling, the metal ion incorporation and expulsion (i.e., lithiation/delithiation, sodiation/desodiation, potassiation/depotassiation) leads to large volume expansion and compression, which in turn can lead to battery breakdown, failure, and degradation. [76,128,289,290] These large volume changes also influence the SEI formation, [49,61,78,283] making the SEI layer unstable and highly surface dependent. [49,289,291,292] Furthermore, a number of intermediate phases are possible as the AM concentration increases. [288,293,294] To this end, the modeling of these materials often investigates phase stability and crystal structure in terms of structure search, which are then related to battery voltage.
Simulated and measured voltage profiles have shown that the group XIV elements Si, Ge, and Sn form alloying compounds with Li, Na, and K making them suitable as anode materials, but differences in capacities are present when changing the battery metal. [128,60,284,[295][296][297] Alloying anodes based on group XV elements could provide greater specific capacity than group XIV anodes, and might provide better performance for NIBs than LIBs. [180,294] Furthermore, the group XV elements can also form high concentration K-compounds. Phosphorus has received attention as an anode material owing to its low cost, abundance, and large capacity. [126,130,[298][299][300][301][302][303][304] Black phosphorus, which has a layered structure (Figure 9a), has the highest electrical conductivity among the phosphorus allotropes and is the most stable. [298] However, black phosphorous needs to be synthesized under high pressure and temperature, resulting in high costs. As for the other alloying anode materials, phosphorus anodes undergo large volume expansions during the electrochemical reaction, leading to poor cycle performance and low initial Coulombic efficiency. Hence, phosphorous anodes are typically composites with carbon, metals, or chalcogenides. Sb-based metal alloys for LIBs, NIBs, and KIBs are reported to have good cycling stability and performance. [253,285]

Voltage and Phase
Simulated and measured voltage profiles have shown that the group XIV elements form alloying compounds with Li, Na, and K making them suitable as anode materials, but differences in capacities and intermediate phases are present when changing the battery metal. [60,284,[295][296][297] DFT studies on lithium insertion in silicon anodes, leading to the formation of lithium silicide, showed that Li preferably occupy interstitial sites, rather than forming substitutional defects. [305][306][307][308] Electronic structure analysis of the Li x Si further showed that Li x Si softens with increasing x due to the breaking of SiSi bonds. [309] By employing AIRSS, Ge-based LIB anodes were found to be potentially safer than lithium silicide (Figure 9b,c). [128] The calculated voltages for the Li x Ge were higher than the corresponding voltages for Li x Si, which could suggest that the germanides would have a lower energy density (Figure 9d). [128] On the other hand, the lithium germanide anode could provide a safer alternative as higher lithium insertion potentials could reduce lithium plating (which can lead to dendrite formation), and more rapid lithium diffusion due to the lower lithium migration barriers in Ge than Si. [128] Combined experimental and DFT studies identified both crystalline and amorphous Na x Sn phases. Increasing x giving increasing specific capacity in the charge/discharge curve, corresponding well to the experimental voltage/capacity curves (Figure 9e), showing the strength of predictive DFT studies (Figure 9f,g). [133,277,293]  Adv. Energy Mater. 2022, 12, 2200662  2200662 (14 of 29) www.advancedsciencenews.com Figure 9h shows the mixing enthalpies and homogeneous bulk free energies of Na x P and Na 1-x (GeP 3 ), Li x P, and K x P. [126,310] From Figure 9h,i, there are similar phase behaviors between Na 1-x P and Na 1-x (GeP 3 ) in the range of 0.5 ≤ x ≤ 1.0, where both of them have three ground states at x = 0.5, 0.75, 1.0. Whereas there is a ground state at x = 0.25 for Na 1-x (GeP 3 ) but the pseudo ground state for Na 1-x P, implying that the Na 1-x P would be transformed into two phases of NaP Figure 9. a) Crystal structures of black phosphorus and GeP 3 with purple spheres being P, and Ge blue. Reproduced with permission. [310] Copyright 2019, Royal Society of Chemistry. b) Li x Si 1-x and c) Li x Ge 1-x formation enthalpies from DFT simulations of structures (blue and red squares) predicted from crystal databases and AIRSS forming the convex hull (red and blue full lines). d) The simulated voltage profile from the most stable phases in (b) and (c). Reproduced with permission. [128] Copyright 2014, American Physical Society. e) The DFT calculated average voltage relative to sodium metal of Na in tin anode (Na x Sn) and f) the convex hull of Na x Sn incorporating previously known structures and computationally predicted structures. (e) The lowest energy structures in (f) and experimental data. The corresponding models of the phases in (e) are shown in (g). Reproduced under terms of the CC-BY license. [277] Copyright 2017, American Chemical Society. h) Mixing enthalpies (or formation energies) of Na 1−x GeP 3 and formation energies. Reproduced with permission. [310] Copyright 2019, The authors, published by Royal Society of Chemistry. i) Electrode potential versus concentration (x) in AM x P. Reproduced with permission. [126] Copyright 2020, Springer Nature. j) A compilation of amorphous (a-) and crystalline (c-) Si, Ge, and Sn Na alloys and the formation energies of these alloys as a function of Na concentration and whether the alloy is amorphous or crystalline. Adapted with permission. [175] Copyright 2015, American Chemical Society. The amorphous structures were obtained from AIMD simulations, and the formation energies from DFT optimizations.
Adv. Energy Mater. 2022, 12, 2200662  2200662 (15 of 29) www.advancedsciencenews.com and Na 0.5 P in the range of 0≤ x ≤ 0.5, and the Na 1-x (GeP 3 ) would be stable without phase separation in the same range because of the unique stable phase at x = 0.25. As a result, Na(GeP 3 ) favors faster Na + kinetics and diffusion, leading to superior rate capability, attributed to the absence of lattice mismatch caused by phase separation. The evidence of different phase behaviors between Na 1-x P and Na 1-x (GeP 3 ) can be further confirmed by the free energies (Figure 9b,c). The phase of Na 1-x P has a broad spinodal region at the range of 0.221 ≤ x ≤ 0.777, suggesting the continuous phase separation of NaP and Na 0.5 P because of spinodal instability. On the contrary, Na 1-x (GeP 3 ) have two spinodal regions (0.109 ≤ x ≤ 0.392 and 0.608 ≤ x ≤ 0.891) compared with Na 1-x P. [126] However, there is a single-phase region (0.392 ≤ x ≤ 0.608) for Na 1-x (GeP 3 ) at the meanwhile, implying the existence of single phase Na(GeP 3 ), which facilitate Na + kinetics by reducing the disconnection and impedance in the internal Na(GeP 3 ). [310,311] Figure 9i further shows that different intermediate compositions are also part of the voltage curves (in agreement with experimental studies) for the phosphourous anodes between the different AM (LiP 7 → Li 3 P 7 → LiP → Li 3 P, Na 3 P 11 → Na 3 P 7 → NaP → Na 5 P 4 → Na 3 P and K 3 P 11 → K 2 P 3 → KP → K 4 P 3 → K 3 P), although the fully charged anode is found to be AM 3 P for all. [126] This difference in alloying mechanisms as a function on AM have also been predicted from DFT simulations for other alloying anodes such as bismuth, silicon, and antimony. [21,294,[312][313][314][315] In an amorphous silicon (a-Si) anode, the Na intercalation energy becomes thermodynamically favorable. [175,307] This leads to an increase in the average a-Si anode voltage for Na, with further benefit of reduced anode volume expansion. [307] Li incorporation in a-Si than c-Si structure was also found to be more energetically stable. From DFT optimization of MD generated amorphous models, it was shown that amorphous Li x Si crystallizes into Li 15 Si 4 (x = 3.75) upon lithiation. [316] This crystallization does not lead to phase separation, nor large-scale atomic motion. Rather the changes in the local bonding configuration are responsible for the transformation from amorphous to crystalline. [316] It was however noted, that due to the possible Na cluster formation in a-Si, limiting the practical metal intercalation, the charge and discharge process needs to be carefully controlled in order for a-Si to be a viable anode material for NIBs. [81,307] However, the volume expansion in these amorphous anodes is significant; 160% for a-LiSi, 230% for a-NaSi, 140% for a-LiGe, 200% for a-NaGe, 260% for a-Li 3.75 Sn, and 480% for a-Na 3.75 Sn. [175] This difference in volume expansion, and also in mechanical properties related to lattice expansion and metal concentration, [317] was assigned to the disintegration of the sodiated materials at Na/metal ratios above 1:1, whereas the lithiated materials maintained the connectivity even at higher Li concentrations. [175] Other silicon allotropes such as Si 24 , which has an open cage structure, have been considered for NIBs, but sluggish sodium diffusion (due to high Na migration barriers) is currently limiting its use even though they show less volume expansion (2.3%) than the a-NaSi. [290] K x Sb, Na x Sb, and Li x Sb show different alloying behavior. [315,[318][319][320] DFT studies helped to reveal that this arises from the difference in thermodynamically stable alloy structures. [318] Experimental research show that after cycling, the Na-Sb alloy forms a hexagonal Na 3 Sb polymorph, whereas Li 3 Sb is found in a cubic structure. [315,318] DFT simulations were used to calculate phase diagrams of Na-Sb and Li-Sb, which show that the hexagonal phase is significantly more favorable for Na 3 Sb than the cubic phase, and that the cubic phase could be stabilized under high pressure. [318] A combined experimental and DFT study further investigated the intrinsic thermodynamic and kinetic properties of Sb lithiation and sodiation and found that the sodiation of Sb also resulted in amorphous phases. [315] Furthermore, the intermediate states observed during experimental studies of Sb anodes for LIBs (Sb→Li 2 Sb→Li 3 Sb (cubic)) and NIBs (competitive reaction paths between Sb→NaSb→Na 3 Sb (hexagonal), and Sb→Na 3 Sb (hexagonal)) were confirmed, demonstrating the predictive and guiding strengths of theoretical calculations for battery anode materials. [294,318,321] For KIBs, four intermediate phases exists; KSb 2 , KSb, K 5 Sb 4 , K 3 Sb. [320] The calculated equilibrium potentials for the potassiation of Sb from DFT simulations was found to be 0.89 V for KSb 2 , 0.85 eV for KSb, 0.44 eV for K 5 Sb 4 , and finally 0.40 eV for K 3 Sb (hexagonal), which corresponds to the highest potassiation. [320] Hence, the final Sb alloying anode product with Li, Na, and K are all M 3 Sb, but the intermediate phases and the final crystal structures differs. Similar observations were made for bismuth alloying anodes. [294] The binary sodium bismuth alloy has been shown from structure exploration calculations (using both pre-existing crystal structures from the Crystal Structure Database, [294] and employing structure prediction with the USPEX code [131] ) to form the alloying compounds Na 3 Bi, and NaBi. [294,322] This indicates that Na in Bi metal anodes follows a similar alloying mechanism to Na in Sn anodes. [294,322] The Na alloying mechanism in P and As anodes on the other hand was shown to include a larger number of Na concentrations, and undergoes significant elastic softening, whereas the Bi anodes show less loss of elastic strength when sodiated. [126,294,312] Similar elastic softening issues were also previously shown to be present in group XIV anodes. [323] This might lead to large irreversible capacities and battery breakdown, leading to the need for further theoretical and experimental studies to stabilize these anode materials. [294]

Li, Na, and K Migration in Alloying Anodes
The calculated Li migration barrier in Si anodes following an interstitial migration path along a hexagonal transition site is 0.54 eV. [306,324] At higher Li concentrations (Li 0.125 Si 0.875 ) the Li migration barrier increased to 0.69 eV. These higher Li concentrations results in a distorted structure and clustering of Li atoms within the anode structure. [306] Chou et al. performed a comparative DFT study of Li and Na in Si, Ge, and Sn. [175] NEB calculations showed that the migration barriers are lower for Li in these crystals than for Na. [175] The higher Na migration barriers were attributed to the larger ionic size of Na. [175] However, a separate study of Se suggested that the difference in sodiation and lithiation speed was more to do with the kinetics of the lithiation and sodiation reactions than the ionic radii. [279] Following an interstitial migration path, the Na migration barriers in crystalline Si, Ge, and Sn were calculated to be 1.08, 0.78, and 0.53 eV, respectively. [175] The same migration path gave migration barriers of 0.62, 0.44, and 0.39 eV for Li in crystalline Si, Adv. Energy Mater. 2022, 12, 2200662  2200662 (16 of 29) www.advancedsciencenews.com Ge, and Sn, respectively. [175] The reason for the higher migration barriers of Na in relation to Li could be due to Na insertion in crystalline silicon (c-Si) anodes being thermodynamically unfavorable (Na intercalation energies 1.90 eV observed from DFT studies) (Figure 9j). [175,307,308] The introduction of a second Na at the interstitial site next to the first Na interstitial in bulk Si was suggested to reduce the Na diffusion barrier in c-Si, arguing that an increasing Na concentration in c-Si anode material can improve its charge/discharge properties. [305] The same observation was made for Na in crystalline Ge (formation energy 0.88 eV), whereas Na in crystalline Sn had a negative (−0.09 eV) formation energy. [175] Crystalline Sn anodes could hence be shown to be more promising NIB anodes than crystalline Ge and Si for future experimental anode studies. AIMD simulations to generate amorphous Si, Ge, and Sn did further show that the Na migration barriers (0.31, 0.27, and 0.26 eV for Si, Ge, and Sn, respectively) are also reduced (as compared to Na migration in the crystalline materials), indicating an alternative materials design route for the alloying anodes. [175] As for the crystalline systems, amorphous Sn shows the most rapid Na diffusion, which from Bader charges have been attributed to the weaker covalent bonding between Sn and Na, as compared to Si and Na. [175] The alloying anodes have been shown to have high storage capacity and facile diffusion, but issues remain regarding their volumetric expansion which could lead to breakdown. To further investigate and prevent this, MD simulations with appropriate force fields or AIMD (which will be limited by the length and time scales these simulations will be possible to conduct with the current technologies) simulations of the thermal and volumetric expansion over anode cycling should be performed. Furthermore, as has been investigated experimentally, composites with graphene [325,326] could both increase electronic conductivity (which is currently a limiting factor of the alloying anodes) and mitigate the volume expansion, and should be pursued also in computational studies.

Conversion-Type Anode Materials
Conversion-type anodes are commonly metal oxides, metal sulfides, or metal selenides. [293,[327][328][329][330][331][332][333] The conversion reaction mechanism can be expressed as A B ( ) B A ynAM yAM x x y n  + + , where AM is the alkali metal ion, A the metal in the conversion type anode, and B is O, N, S, P or Se. [253,[333][334][335][336] This conversion reaction can result in high theoretical capacities, making these anodes highly attractive. [253,334] During charge, AM + moves from the cathode, via the electrolyte, and arrives at the anode surface where a redox reaction leads to the formation of AM n B and A. During discharge, A reacts with AM n B, returning AM + to the cathode. [253,334] Despite having high capacity for LiBs, transition metal oxide anode materials have yet to realize commercial success due to their high voltage, poor cycling stability, and volume expansion, which result in decreased energy density during battery operation. [337]

Strucutral Simulations
The anisotropic structure of the MoS 2 monolayer provides numerous metal adsorption sites. [327,329,330,[338][339][340][341][342] Fast electron transfer and AM ion intercalation/deintercalation is facilitated in the layered bulk. [329] It has been suggested that during charge this anode first undergoes an intercalation type reaction, and then a conversion reaction (forming Mo and AM s S) in the final step. [338,343,344] However, MoS 2 on its own has shown poor cycling behavior and low reversible capacities. Doping with Se or composites with other sulfides and carbon-based materials have been investigated experimentally to improve its stability. [251,329,330] DFT, MD, and AIMD studies of MoS 2 anode materials have focused on the structural stability, and the electronic properties of the MoS 2 structure, forming an understanding of its intrinsic intercalation and surface adsorption properties. Layered MoS 2 exists in five different polytypes; hexagonal 2H (naturally most stable form), tetragonal 1T (high energy, stabilized from liquid exfoliation), rhombohedral 3R, tetragonal 2T (from theoretical studies), and coexistence of 2H and 3R (Figure 10a). [345] DFT studies of Na intercalation into layered MoS 2 with the hexagonal 2H symmetry have shown that Na intercalation leads to a phase transformation to the tetragonal phase. [339,345] At higher Na concentrations, the 3R and 2T polymorphs becomes energetically and geometrically identical to the 2H and 1T polymorphs. [345] A similar observation was made for lithiation of the 2H MoS 2 structure, indicating that the metal concentration directly affects the crystal structure of this anode material and that a phase transition occurs at high metal concentrations due to a softening of the mechanical properties. [346] The incorporation of Li in 2H-MoS 2 leads to a smaller interlayer expansion than for Na intercalation, 11% (to 14 Å) versus 23% (to 15.5 Å). [341,343] K has also shown strong MoS 2 surface adsorption, [347] and stable K intercalation structures. Potassiation of MoS 2 is energetically favorable in the 2H structure up to K 0.875 MoS 2 , and the intercalation energies per K atom are stronger at lower K loading. This high K loading is higher than the experimentally observed K 0.4 MoS 2 . [348] At higher K loading, MoS 2 was shown to reduce to Mo and K x S, in accordance with the previously seen intercalation/conversion reaction for this anode material. [348] As for Na, the surface adsorption energies are stronger (−3.16 to −3.46 eV) [347] than the intercalation (≈−2 eV at low K loading). [349,350] AM intercalation in between the MoS 2 monolayers in bulk MoS 2 incorporates at two distinct sites (octahedral, or tetrahedral), in a similar fashion to Li intercalation in Bi 2 S 3 . [136,345] For the octahedral site, AM is bonded to six S, and for the tetrahedral site it is bonded to four S. [345] The calculated intercalation energies show that Na intercalation at octahedral sites is more energetically favorable than tetragonal sites, whereas Li intercalation is only more energetically favorable at an octahedral site at interlayer distances below 14.7 Å. [341,345] Hence, it was concluded that the weak dispersion interactions in between the MoS 2 layers (leading to the stabilization of the different polymorphs) become negligible in relation to the stronger chemical interaction between the AM and MoS 2 . [345] Na intercalation in Na x MoS 2 at concentrations up to x = 1.0 were shown to be energetically stable in both the 1T and 2H polymorph, with the Na intercalation energies being stronger in the 1T polymorph. [345] The electronic structure of the 1T polymorph from DFT calculations show that the lowest degenerate 4d xy , 4d yz , and 4d xz orbitals are partially occupied, making them readily available for accommodating the AM electron upon bonding. [345,346] Adv. Energy Mater. 2022, 12, 2200662  2200662 (17 of 29) www.advancedsciencenews.com

Migration
Both Li and Na were found to form bonds and interact with the sulfur sites instead of the Mo-sites when intercalated and diffusing in between the MoS 2 layers. Due to these Li-S/Na-S bonds, the interlayer distance between the MoS 2 layers were found to be directly impacting the Li and Na binding and migration energies. [341,343] The Na migration barriers (Figure 10b) in the 1T polymorph was also found to be lower (0.28 eV) than in the 2H polymorph (0.68 eV), suggesting that the 1T polymorph is more suitable as an NIB anode, due to its higher electrochemical performance from voltage simulations, metallic electronic structure, more rapid Na mobility, and higher Na capacity. [345] DFT simulations also showed that the Li and Na diffusion in 2H MoS 2 can be improved by expanding the interlayer distance (originally 12.60 Å) between the MoS 2 layers, showing a way for materials design and engineering to be applied experimentally. [341,343] 1T-MoS 2 was shown experimentally to have faster potassium kinetics in terms of reversible K + intercalation and conversion than 2H MoS 2 . [351] However, to the best of our knowledge, no computational study of the effect of different MoS 2 polymorph on K intercalation has been presented in the literature.
Combining a MoS 2 layer with graphene can further lower the Li and Na migration barrier and improve LIB and NIB performance. [166,352,353] Wu et al. simulated Na migration in between two MoS 2 layers and in between a MoS 2 layer and a graphene layer. In the graphene MoS 2 system, the Na migration energy barrier is 0.21 eV, with the bilayer MoS 2 having a higher Na migration barrier of 0.34 eV. [352] A separate study by Sun et al. showed that the Na migration energy on the MoS 2 layer in the MoS 2 graphene systems has a higher migration energy barrier of 0.35 eV in between two MoS 2 rings, with a smaller migration barrier of 0.07 eV moving Na from on top of a Mo site to the middle of the MoS 2 ring. These two diffusion behaviors (in between layers and on the MoS 2 surface) are likely from these studies to coexist in physical anode materials, with a higher degree of surface diffusion postulated to be due to the stronger Na surface adsorption energy than Na intercalation energy. [353] Graphene can also be combined with a SnS 2 monolayer. Samad et al. modeled the 1T SnS 2 monolayer, both with and without a graphene layer. Na atoms was found to adsorb both at the hollow position (above a sulfur atom) and on top of a Sn atom, with only a small difference in adsorption energy (−1.36, and −1.34 eV, respectively). These energies can be lowered by intercalating Na in between a MoS 2 and a graphene layer, and this was shown to also hold true for the adsorption and intercalation of more than one Na atom. [335]

Novel 2D Anode Materials
In the above sections, computational modeling has mainly been used to investigate anode materials that have already been suggested or tested experimentally to understand the material performance and how atomic scale alterations can improve AMIB performance. Computational modeling is also employed to explore and investigate novel materials that have not been Figure 10. a) shows an illustrated Na migration path in two polymorphs of MoS 2 , where orange is Na, and purple is Mo. b) The corresponding migration energy profiles for the two paths. Reproduced with permission. [345] Copyright 2014, Elsevier. c-g) Combined graphene and MoS 2 /h) SnS 2 heterostructures for NIB anodes. (c) The optimized structure of the MoS 2 /graphene heterointerface from DFT simulations, (d) the Na diffusion path in (c), and (e) the Na diffusion path for Na diffusion between two MoS 2 layers (no graphene). The corresponding NEB migration paths in (f) for (d) and (e) shows that Na migration at the heterointerface is more mobile. (g) The trajectory for Na hopping in the heterointerface in (d) from AIMD simulations at 600 K. Reproduced with permission. [352] Copyright 2018, Royal Society of Chemistry. h) The NEB simulated diffusion path for Na in a SnS 2 /graphene heterostructure (side view left and top view right) following a H1-T1-H3-T3-H1 path as illustrated in the top view figure to the right. Brown spheres are carbon in graphene, blue Na, black Sn, and green S. Reproduced with permission. [335] Copyright 2016, Royal Society of Chemistry.
Adv. Energy Mater. 2022, 12, 2200662  2200662 (18 of 29) www.advancedsciencenews.com tested in real batteries, and to guide the experimental development of anode materials. For this purpose, a plethora of novel 2D materials (Figure 11) have been studied computationally, [354][355][356][357][358][359] with the aim to increase adsorption as an indication of storage and migration of AM. [22,326,[360][361][362][363][364] Recently graphdiyne, which consists of both sp-and sp 2hybridized carbon atoms ( Figure 11) have been investigated computationally. [364] Its structure can be extensively modified by varying the number of sp-and sp 2 -carbons in the simulation cell, making it particularly suitable for computational investigation and screening prior to experimental synthesis. The calculated Na capacity of graphdiyne is NaC 2.57 , which corresponds to a theoretical specific capacity of 497 mAh g −1 . By layering the graphdiyne monolayers, in a manner similar to graphene layering in graphite, a lower capacity (316 mAh g −1 ) is obtained, corresponding to NaC 5.14 . [366] This furthers the argument that the graphdiyne monolayers provide more storage capacity for NIBs per carbon atom. Li adsorption on methyl substituted graphdiyne was further shown to be highly dependent on lattice site, with lattice sites closer to the alkynyl group being more favorable for Li adsorption (adsorption energy −1.34 eV) than away from this site (adsorption energy −1.16 eV). [360] Based on these calculations, the theoretical lithium storage capacity is 1701 mAh g −1 , showing high LIB anode capacity. [360] Upon K insertion in graphdiyne, transformation of sp-hybridized carbon to sp 2 -hybridized carbons, and theoretical capacities of 2870 mAh g −1 in monolayer graphdiyne calculated due to the formation of C 7 K 9 . As for Na, layering of graphdiyne for KIB anodes lowers the theoretical storage capacity of K to 1700 mAh g −1 .
The lower storage capacity of K in the layered structure was assigned to the higher K migration barrier between the graphdiyne layers (0.666 eV) than on the monolayer (0.175 eV). [372] The challenge now lies in synthesizing these monolayer graphdiyne anode materials, and form stable interfaces with the battery electrolytes.
Recently, Zhang et al. explored the possibility of using the carbon nitride 2D materials of C 2 N (Figure 11), C 3 N, and g-C 3 N 4 as LIB anode materials. [356,365] From calculations of Li adsorption energies, C 3 N and g-C 3 N 4 were found to be unsuitable for LIBs. It was however suggested that a nanotube version of g-C 3 N 4 could show improved Li capacity, making this a more suitable anode material for future investigation. [356] The theoretical capacities calculated for the C 2 N material was 671.7 mAh g −1 .
As for Li diffusion through graphene sheets, the Li diffusion through the C 2 N sheet was found to be highly energetically hindered (E m ≈ 7 eV), whereas the metal migration across the material is more rapid. [356] The migration barriers of on-plane Na and K diffusion were reported to be 0.24 eV for Na, significantly higher than on graphene, and 0.10 eV for K (the Li migration barrier for comparison on the same system is 0.25 eV). [373] DFT simulations further showed that C 2 N monolayers could be a promising anode material worthy of experimental investigation with high theoretical specific capacity of 1181.9 mAh g −1 for Li, 599.72 mAh g −1 for Na and 385.81 mAh g −1 for K. [373] A 2D nanoporous graphene was recently synthesized, [374] and shown through DFT simulations to have favorable AMIB properties Figure 11. Computational models of 2D anode materials of C 2 N (Reproduced with permission. [365] Copyright 2017, Wiley-VCH), graphdiyne (Reproduced with permission. [366] Copyright 2017, Elsevier), Zr 2 N (Reproduced with permission. [367] Copyright 2020, Elsevier), Sb nanosheet (Reproduced with permission. [368] Copyright 2019, Elsevier), phosphorene (Reproduced with permission. [301] Copyright 2018, Royal Society of Chemistry), SiGe (Reproduced with permission. [369] Copyright 2019, Elsevier), silicene, germanene, and stanine (Reproduced with permission. [370] Copyright 2016, Royal Society of Chemistry), puckered layer of GeS, GeSe, SiS, and SiSe (Reproduced with permission. [371] Copyright 2016, American Chemical Society), and nanoporous graphene (Reproduced with permission. [358] Copyright 2020 American Chemical Society).
Adv. Energy Mater. 2022, 12, 2200662  2200662 (19 of 29) www.advancedsciencenews.com when doped with heteroatoms. [358,375] This further highlights the importance of heteroatom doping and defect engineering as a tool for anode design and optimization, and the power of computational methods in guiding practical studies.
2D zirconium nitrogen ( Figure 11) was studied for NIB anodes by assessing different stable phases and functionalization (O, OH, and F). [367] Among the considered ZrN, t-Zr 2 N, h-Zr 2 N, t-ZrN 2 , h-ZrN 2 , and p-ZrN 2 , t-Zr 2 N was shown to be the most stable. Functionalization of this 2D material with O improves the Na adsorption energy (−2.042 eV for nonfunctionalized, and −2.505 eV when O-functionalized), whereas OH− and F-functionalization led to weaker Na adsorption (−1.095 and −0.791 eV, respectively). [367] The trend Na adsorption energy is also carried through in the theoretical capacity (C,Zr 2 N = 545.9 mAh g −1 (x AM,max = 4), C,Zr 2 NO 2 = 469.4 mAh g −1 (x AM,max = 4), C,Zr 2 N(OH) = 413.1 mAh g −1 (x AM,max = 3.5)), where C is calculated with Na adsorbing on both sides of the 2D surface. [367] As previously seen for graphene, the stronger adsorption can also lead to higher metal migration barriers, with Na migration at the nonfunctionalized Zr 2 N having E m 0.018 eV, and O-functionalized Zr 2 N 0.13 eV. [367] Hence, this material could be sensitive to electrolyte breakdown and SEI constituents, slowing down NIB anode performance if functional groups are allowed to bind to the surface.
Phosphorene is a single layer of black phosphorus and is widely studied due to its similarity to graphene and graphite in its layered structure. Each phosphorus atom is covalently bonded to three other phosphorus atoms, forming the puckered structure shown in Figure 11e. [355] As for graphene, vacancy defects in the phosphorene plane lead to stronger metal adsorption. [376] When Li atoms are continuously adsorbed on phosphorene with a single vacancy defect, its dangling bonds will be filled first, resulting in a semiconductor to metal transformation, which can improve the electrochemical performance. [376] A similar transformation from semimetallic to metallic conductivity caused by Li adsorption have also been observed in separate studies. [355,[377][378][379][380] It was seen that there are strong noncovalent interactions and extensive charge transfer between the adsorbent and phosphorene, which will improve their corresponding electronic properties. [355,[377][378][379][380] Unfortunately, phosphorene ( Figure 11) is not thermally stable and is highly reactive with air. [141] Hence, studies combining phosphorene with other 2D materials are instigated.
DFT studies of lithium and sodium migration between the hexagonal boron nitride (h-BN) and phosphorene layers and above the surface of the phosphorene layer showed that structural anisotropy has an essential effect on metal migration. [381] This work revealed that the energy barrier along the zigzag direction is far lower than that along the armchair direction, while the other 2D materials (graphene, Ti 3 C 2 ) have higher diffusion barriers due to their isotropy. Moreover, the diffusion barrier of h-BN/phosphorene is lower than that of isolated phosphorene, suggesting the superior charge/discharge capability is improved by the heterostructure. [381] Moreover, the electronic properties of phosphorene/BN and its lithiated/sodiated products have a bandgap of 0.68 V, which would limit electronic conductivity. However, when phosphorene/BN is intercalated with Li or Na, a semiconductor to metal transition occurs. [381] Further systematic DFT studies of the effect of metal loading (which would reflect different stages in the electrochemical cycle) on the electronic conductivity would be required to further optimize these structures for practical anode applications.
2D Si (silicene), Ge (germanene), and Sn (stanine) have been stipulated as possible anode materials based on the high specific capacity of their alloying anode counterparts and the previous successes of graphene. [382][383][384] Making the Si, Ge, and Sn anodes as a 2D material removes the issues relating to battery degradation arising from their high volume expansions, [269,328,370,385,386] although stability issues have been seen under ambient conditions. [359] Silicene, germanene, and stanine have large surface area, high electron mobility, and high mechanical flexibility. [382] As opposed to graphene, the 2D sheets are not flat, but are prone to buckling. [382,386] DFT simulations of silicene, germanene, and stanine predicted that these materials have high charge capacities, strong adsorption of both Li and Na (with Li being more strongly bound at each material and concentration than Na), and Li and Na migration follow a zigzag path, with migration barriers lower than on graphene (Figure 12). [382,384,386] Among these 2D materials, the silicene surface provides the fastest diffusion path for both Li and Na (E m 0.2, and 0.13 eV, respectively). Slightly higher migration energies for Li (0.385 eV), and Na (0.23 eV) were obtained for silicene in a separate study, concluding that Li diffusion on silicene is similar to graphene, whereas Na diffusion on graphene (E m 0.14 eV) would be more rapid. [370] The diffusion path through the sheet (through a hole site) was found to be energetically unfavorable. [382] Building on the favorable anode properties found for the silicene, germanene, and stanine, a later study modeled a 2D material consisting of both Si and Ge for NIBs and KIBs. [369] This material, siligene (Figure 11), was predicted to be thermally stable through DFT formation energy and phonon calculations, and showed similar buckling to silicene and germanene. [369] Surprisingly, adsorption energy calculations showed that both Li and Na have similar adsorption energies at this surface (E ads −0.85, and −0.84 eV for Li, and Na, respectively), whereas K forms a much stronger bond (E ads −1.24 eV). K had the lowest calculated migration barrier, with Na and Li showing similar migration barriers as on silicene and germanene. Voltage simulations further showed that siligene would theoretically be suitable for LIBs, NIBs, and KIBs. [369] Mixing a silicene and graphene layer was investigated by Shi et al., resulting in a material with an interlayer distance of 3.53 Å. [370] This layered silicene graphene heterostructure material was shown to be able to accommodate large Li and Na concentrations (Na intercalation in between bilayer graphene does also become energetically favorable at ≈3.5 Å [85,111,121,142,151,246] ). DFT migration energy simulations also showed that the silicene graphene material retains rapid metal migration with E m 0.36-0.40 eV for Li, and 0.19-0.30 eV for Na. [370] Based on the interlayer distance simulations of K-intercalated graphite, this material would not be expected to be suitable for KIB anodes. DFT stimulations should be employed to further investigate these materials for KIBs, especially as this material has been shown from simulations to have high structural stability during battery cycling. The mechanical stiffness for the combined material (383.45 N m −1 ), and the isolated 2D structures (65.78 N m −1 for silicene and 292.78 N m −1 for graphene) was Adv. Energy Mater. 2022, 12, 2200662  2200662 (20 of 29) www.advancedsciencenews.com calculated. The higher mechanical stiffness of the silicene graphene heterostructure as compared to the isolated materials showed that this anode will preserve its structure during battery cycling and hence have better cycling performance, even though volume expansions of 20% and 40% were observed upon lithiation and sodiation respectively. [370] Building on the research of these rumpled 2D materials, other metals such as Sb have been investigated computationally as rumpled nanosheets for AMIB anodes. Bare (Figure 11d) and hydrogenated Sb nanosheets were shown through B3LYP level DFT simulations to be thermodynamically stable, and to be suitable for NIB anodes. [368] As for the graphdiyne materials, computational materials design has shown to be very efficient in predicting also noncarbon 2D anode materials, and the next challenge is to make these a practical reality.
2D monochalcogenides such as 2D selenides (GeSe, SiSe) and 2D sulfides (GeS, SiS) have also been investigated computationally to explore their possible applications in AMIBs. These 2D monochalcogenides are resistive to oxidation, and are thermodynamically stable under ambient conditions. [141,371] 2D selenides, such as GeSe (Figure 11), have been proposed as anode materials due to their high surface area and superior electronic properties. [141,371] These materials are suggested to be able to be grown through liquid or mechanical exfoliation. DFT simulations of Li, Na, and K have shown that this 2D material has energetically favorable metal adsorption (E ads Li −1.92 eV, Na −1.25 eV, and K −1.36 eV), and rapid metal migration paths in the zigzag direction (E m Li 0.29 eV, Na 0.12 eV, and K 0.11 eV). [141,371] Metal migration could also be conceptually possible in the armchair direction, but leads to higher metal migration energies due to the stronger repulsion from the Ge ions. [371] The higher Li migration energies are due to its stronger adsorption energy. [141,371] 2D GeSe can accommodate a maximum of 8 Na, 6 K, and 2 Li per unit cell, whereas the bulk layered material can only accommodate one AM ion in the unit cell. [141] Hence the theoretical capacities of 177, 707, and 530 mAh g −1 , respectively, from voltage simulations of Li, Na, and K on 2D GeSe are much higher than the theoretical capacity of the layered bulk GeSe anode (89 mAh g −1 for Li, Na, and K). [141] Based on these computational studies, 2D GeSe could be a suitable anode material for NIBs and KIBs, but would have limited success as a LIB anode. Other 2D monochalcogenides such as SiS, SnS, SnS 2 , GeS, SiSe, and SnSe have also been investigated as possible anode materials based on the previous successful modeling of other 2D materials. [328,371] Li has shown stronger adsorption to the silicon-based monochalcogenides, with theoretical specific capacity calculated as 445.7 mAh g −1 for SiS, and 250.44 mAh g −1 for SiSe. [328,371] In summary, DFT simulations allow for exotic 2D materials to be assessed as anode materials for AMIBs. However, the practical applications of these materials in batteries are not straightforward, due to the challenges in synthesizing these materials experimentally. The computational investigation of Zr 2 N anodes also identified practical issues and breakdown mechanisms, and similar treatment and analysis of other 2D materials would be useful.

Summary and Perspectives
High-performance anodes remain one of the most critical optimization and engineering challenges to achieve high-performance alkali metal ion batteries. In this review, we have outlined the current progress on LIB, NIB, and KIB anode materials from computational studies. To date, anode materials have largely been focused on graphite/graphene-based carbon anodes. Based on the success in modeling graphene, a plethora or other 2D anode materials have been suggested and investigated computationally. A common approach to improve anode performance is to couple graphene to other conversion, intercalation, and alloying anodes. The alloying anode materials, whilst successfully modeled remains short of reaching their (c,d) The same for Na. The surface paths follow a zigzag pattern. Reproduced with permission. [382] Copyright 2016, Elsevier.
Adv. Energy Mater. 2022, 12, 2200662  2200662 (21 of 29) www.advancedsciencenews.com potential as an anode material due to large volume expansion. To further optimize these anodes materials for practical applications, additional investigations are needed, and further research into alloying-graphene (or other suitable materials with high electronic conductivity and ability to restrain the volume expansion) composites required.
Performance limiting processes in AMIBs such as irreversible capacity loss, phase stability, and sluggish diffusion, can be probed using computational molecular modeling. The strength of computational methods is to provide insight in the effect of material composition and structure on materials properties, advancing materials design. This review has shown conclusively that atomic scale modeling can predict, verify, and optimize anode materials for alkali-ion batteries. Computational materials design, as highlighted in this review article, is manyfold. One successful approach to increase or tune the intercalation and adsorption ability of anode materials is through defects and dopants. Atomic scale computational studies have been shown to be powerful in directly probing the effect of both bulk and surface defects on capacity, voltage, and diffusion, providing input to experimental studies. Indeed, DFT studies have been further employed to show which defects or heteroatoms should be avoided or limited in an anode material as they would provide strong storage sites, but limit the reversibility and cycling behavior. Similarly, simulations have shown the effect of surface functionalization on anode performance for both carbon-based and MXene anodes, further showing the difference in AM interaction with the same surface functionalization. Moving away from defects, dopants, and functional groups, computational materials design has also been used to directly probe the effect of interlayer distance in intercalation and conversion-type anode materials. By altering the interlayer distance, and indeed also the pore shapes, increased anode performance could be proposed. Similarly, atomic design of the crystal phase, structure, and degree of disorder and amorphousness has been pursued to increase anode performance in both intercalation, alloying, and conversion-type anodes. This has further led to easing of the transference from LIB to KIB anodes, where the ionic radius of the AM ion varies. To move forward, techniques and methodologies to model more complex materials, heterogeneous materials and realistic interfaces are required, in addition to the already very powerful DFT approaches currently in use.
Computational materials design has greatly accelerated the AMIB field. Computational studies can be helpful in screening materials, and present guides to improved materials properties through defect and dopant engineering. Future advances in anode materials need to continue the exploration of novel structures, as well as continue to optimize current materials. More wide integration of experiments with computational methods to gain fundamental understanding of structure, composition, defect, dopant, and electronic and mechanistic features of these materials would be advantageous. Deeper atomic-scale and fundamental insights into local structures, heterogeneities, diffusion mechanisms, electronic conduction, and interfaces with electrolytes will be needed to move to the next frontier. From experimental characterization, there are problems about heat generated during battery cycling, which would require the modeling and understanding of high temperature thermodynamic processes.
The ideal scenario for the computational modeling of any material would be a one-size-fits-all approach, with no empirical parameters and the only input needed being material composition or desired property. For atomic scale modeling, more powerful computing facilities would eradicate many of the challenges currently faced due to the need for approximations and assumptions. As demonstrated in this review, significant advances in computational materials science have been made. Limitations do exist, and improvements and expansion of capabilities are required. One of the areas that would lead to the most impact in DFT and AIMD modeling of anode materials is the development of better performing exchange-correlation functionals. It has been suggested that the new meta-GGA SCAN functionals would improve the modeling of layered cathode materials, [19,387] as it can treat both dispersion and localized states within the functional and as such it would be an interesting topic to pursue for future modeling for the layered anode materials. Extensive testing and benchmarking would be required to evaluate both its accuracy, performance, and computational expense as compared to the current state-of-the-art.
In practice, real systems are often inherently more heterogeneous than that can currently be handled by DFT and MD simulations. MD melt quench methods followed by DFT optimizations have been shown to be excellent to model some glassy materials, but it is heavily reliant on experimental data to validate the generated structures and show that the generated structures are representative of a physical system. Machine learning potentials fitted against experimental pair distribution functions of hard carbon materials have also been reported, but again the length scales currently achievable with DFT simulations, and the absence of MD interatomic potentials that can simulate the electrochemical mechanisms, limits the modeling even of these more realistic models to small motifs (<1000 atoms). Furthermore, the dependence of these machine learning potentials on training sets could limit the predictive powers needed for computational materials design. Nevertheless, such machine learning approach by coupling experimental and computational data to understand the atomic scale mechanisms could be very powerful and should be developed further. In depth reviews of machine learning for materials design can be found in reference [388][389][390][391] The next frontier in atomic scale computational materials design should be scaling methods that allows for larger systems to be modeled with DFT, and for longer AIMD trajectories to be obtained within reasonable computational expense and time. This would allow for both larger and more complex systems to be investigated. Additionally, the possibility to combine different methods for different materials within the same simulation (for example, interfaces between metal and semiconductors/insulators [69] ) would advance the field.
An anode property that is typically not investigated to the same extent as metal ion storage (relating to battery specific capacity and voltage) and metal ion diffusion, but key to a highperforming anode material, is electronic conductivity. Electronic conductivity contributes to high-rate performance and can be calculated from DFT simulations but is dependent on the knowledge of experimental relaxation time which is challenging to obtain. It is however possible to set a relaxation time, and to comparatively investigate the electronic conductivity between different materials. Experimentally, graphene has been Adv. Energy Mater. 2022, 12, 2200662  2200662 (22 of 29) www.advancedsciencenews.com composited (as also evident from the computational studies presented here that have layered different materials with graphene) with both alloying and conversion type anode materials to increase their electronic conductivity. For future anode studies, the electronic properties (which can also be assessed in terms of density of states and band structures), and their thermal variations should be considered, especially in the high specific capacity (here investigated in terms of stable formation of high Li, Na, or K concentrations from intercalation or phase stability simulations) alloying and conversion type anode materials.
Moving toward a working battery setup, it is not enough to design a standalone high-performing anode material. The anode material must form a chemically, thermally, physically, and mechanically stable interphase with electrolytes, whilst at the same time allowing for facile ionic diffusion across this interphase. Additionally, this interphase might vary (forming and reforming) over the electrochemical cycles. One of the main breakdown mechanisms of LIBs, NIBs, and KIBs is related to the anode electrolyte interphase, and computational studies as the next frontier need to expand to not only consider the isolated bulk or surfaces, but also these interphases. This is however a nontrivial task. Homogeneous interphases are currently investigated computationally. More realistic, heterogeneous interphases that consider defects, disorder, amorphous materials (or indeed amorphous domains in the interphase), defects, dopant segregation to surface, electrolyte molecule dissociation, grain boundaries, etc. are yet in their infancy. These simulations would require "large" systems to be modeled with DFT, and for the diffusion processes long AIMD trajectories, bringing the attention to the development of more efficient codes and HPC infrastructure.
Apart from the anode materials reviewed in this paper for AM-ion batteries, it is wort noting that there is increasing interest in using AM metals directly as anodes, [27,[392][393][394] in next-generation batteries such as AM-O 2 batteries, AM-S batteries, and AM solid state batteries, as they could provide the highest theoretical specific capacity. Whilst beyond the scope of the topic of this review, there remain challenges and opportunities for computational materials design to support the development of these next-generation batteries. One challenge with AM metals anodes in AM-O 2 and AM-S batteries with liquid electrolyte is the dendrite formation and growth, and the formation and growth of the SEI. DFT and MD simulations could be used to understand the dendrite formation and growth mechanisms as well as the SEI formation mechanisms at the metal-liquid electrolyte interface. DFT and MD simulations can also be used to support the design of AM metal alloy and interface engineering approaches to stabilize AM metals, [395][396][397] as well as the design of novel electrolyte solutions to stabilize the SEI. [398] One possible solution allowing for the use of AM metal anodes is to exchange the liquid electrolytes to solid ceramic electrolytes. [399,400] Currently, a large research effort, accelerated by DFT and MD simulations, is taking place in designing solid electrolytes with ionic (i.e., Li, Na, and K) conduction comparable to that of the liquid electrolytes. These solid electrolytes were initially expected to mitigate the safety issues related to the use of liquid electrolytes, but recent studies have suggested that lithium dendrite formation could also be an issue at solid electrolyte | lithium metal anode interfaces. [45,[401][402][403][404] To this end, further computational studies are required to pinpoint the nucleation mechanisms at the AM metal|solid electrolyte interface, leading to dendrite formation and interface breakdown, as well as atomic-scale design of solid electrolytes themselves. Moving away from AMIBs, but retaining the goal of using metal anodes, multivalent ion batteries based on Mg, Ca, Zn, and Al ions have been suggested. [9] These battery technologies have the potential to provide higher capacity (due to the multiple electron transfer per metal cation) compared to the single valent AMIBs, and have also been suggested to have improved safety when coupling the metal anodes with liquid electrolytes. [9,280,[405][406][407] Whilst these battery technologies remain less well studied than their AMIB counterparts, computational studies are excellently poised to translate the same materials design strategy used when adopting LIB anodes to NIBs and KIBs to the multivalent metal-ion battery anodes.
Finally, as evident from this review (and not surprising considering the timeline of technology development) KIB anode materials are less studied computationally than NIB and LIB anodes. Hence, there are opportunities for computational studies to screen, validate, optimize, and modify the currently investigate LIB and NIB anode materials for KIBs, and to predict novel ones. The anode materials currently studied computationally for AMIBs should further be evaluated for use in alternative technologies such as divalent batteries, sulfur batteries, air batteries, etc. The current wealth of computational modeling of anode materials have shown to be invaluable for materials design and should continue to be incorporated with experiments based on their success in materials optimization, materials discovery, and materials characterization, and will play an integral part in the future work by also considering nanostructured materials, heterogeneities, composites, and compatibility with electrolytes. Development of a screening methodology with not only anode-specific properties, but also identifying anode properties and materials that will be stable with electrolytes and under different operating conditions will allow the practical advances to accelerate by removing the number of anode-electrolyte couples needed to be explored in a lab.