Ab initio study of lithium diffusion at a layered carbon/amorphous silicon interface

Abstract

Silicon–carbon composites are widely investigated as anode materials for lithium-ion batteries because they offer improved electrochemical performance compared to the respective pristine elements. In this study, we comprehensively examine Li diffusion in model interface systems composed of γ-graphdiyne and amorphous silicon using ab initio molecular dynamics simulations and nudged elastic band calculations. We analyze the structure–property relationships between the Li diffusion coefficient – which is found to be one to two orders of magnitude higher in γ-graphdiyne than in amorphous Si – and the local environment and charge state of Li in both materials. Migration of lithium across the silicon–carbon interface exhibits only moderate energy barriers and therefore does not limit the overall diffusion kinetics of the composite; instead, bulk amorphous Si constitutes the primary kinetic bottleneck. We conclude that nanostructure design should prioritize controlling the size and shape of Si particles rather than targeting the nature of the interface.

---

1 Introduction

Recent advances in electric mobility and renewable-energy production have provided considerable impetus for battery research. These developments demand high-capacity energy storage, highlighting the need to enhance the performance of lithium-ion batteries (LIBs), the most widely adopted technology. While novel materials are being explored for every battery component, particular emphasis has been placed on discovering new electrode materials to increase energy density.^1,2

Currently, most commercial LIBs employ graphite anodes, which form Li intercalation compounds via the reversible insertion of Li ions into their layered structure. During lithiation, the material's structure and volume undergo relatively minor changes, expanding by only 13% to yield LiC₆. This modest expansion enables very good cycling performance, with minimal capacity loss during repeated (de) lithiation.^3,4 In addition to graphite, other carbon allotropes – such as graphene and amorphous carbon – have been investigated as negative electrode materials.^5–7 Although some carbons exceed graphite's gravimetric energy density of 372 mAh g⁻¹,³ certain non-carbonaceous materials offer substantially higher theoretical values; silicon, for example, reaches 3579 mAh g⁻¹.^8–10

Fully lithiated silicon, Li₁₅Si₄, can store roughly ten times more Li than graphite, making it a promising candidate for high-performance anodes in next-generation LIBs.^8–10 During electrochemical lithiation, both crystalline and amorphous Si first form amorphous Li_xSi,¹¹ which subsequently crystallizes to Li₁₅Si₄ at x = 3.75, after which no further Li insertion occurs.^12–14 Delithiation proceeds through amorphous Li_xSi and finally recovers amorphous Si.^12,13 However, upon full lithiation, the volume of the initially crystalline Si expands by a factor of 3.8, inducing severe mechanical stress and particle cracking. Fractured particles eventually lose contact, leading to Li trapping and rapid capacity fade, i.e. poor cyclability. In addition, Si exhibits low intrinsic electrical and ionic conductivity.^8–10

Ongoing research addresses these cyclability issues by employing Si nanostructures and exploring composite materials. In particular, Si nanocomposites with different carbon allotropes have been intensively studied, demonstrating promising electrochemical performance. Si/C composites are synthesized via methods such as pyrolysis, chemical or physical vapor deposition (CVD, PVD), mechanical milling, electrospinning, or combinations thereof.^8–10 Such approaches produce a wide variety of nanoscale architectures, which are typically characterized structurally by SEM, TEM, XRD, XPS, and Raman spectroscopy.^15–23 Many Si/C composite anodes exhibit significantly improved cycling stability: the soft carbon matrix buffers the volume expansion of Si particles during lithiation, preventing mechanical stress, particle cracking, and contact loss of the material. Additionally, graphitic carbon domains enhance both electrical and ionic conductivity throughout the electrode.^8–10

Li diffusion in crystalline (c-Si) and amorphous silicon (a-Si), as well as in fully lithiated (c-Li₁₅Si₄) and partially lithiated Li–Si compounds (a-Li_xSi), has been investigated both experimentally^14,24–26 and computationally.^26–41 Likewise, Li diffusivity in graphite has been studied through both experimental^42–44 and theoretical^45–47 works. While atomistic Li migration pathways and energy barriers in bulk Si and graphite are now well characterized, the impact of Si/C interfaces on Li diffusion in composite anodes remains an open question. Interfaces can formally be seen as two-dimensional crystal defects, which considerably affect physical material properties such as binding affinities and diffusion constants. Experimental techniques often have difficulties to discriminate unambiguously between bulk and interface regions as the origin of particular spectroscopic observations; in contrast to this, atomistic simulations are well-suited to isolate these effects.

Computational works report enhanced Li mobility at Si/C interfaces compared to bulk silicon.^48–56 However, most studies methodologically focus on either molecular dynamics (MD) simulations^48–53,56 or nudged elastic band (NEB) calculations,⁵⁴ with only one study combining both techniques.⁵⁵ Moreover, although two different carbon allotropes, graphene^{48,50,54–56} and amorphous C,^49,51,53 have been modeled, nearly all interface models employ crystalline Si,^{49–51,53–56} despite amorphous Si being the relevant phase in anodes.

In this work, we construct model systems from γ-graphdiyne and amorphous silicon to investigate Li diffusion in Si/C composites. To our knowledge, this is the first study to explore Li mobility at a Si/C interface incorporating a graphyne allotrope. Moreover, recognizing that Si anodes transition from crystalline to amorphous during initial lithiation/delithiation, we use amorphous rather than crystalline Si, contrary to most previous work. We employ a combination of ab initio MD simulations and NEB calculations to identify meaningful Li migration paths and characterize the underlying diffusion mechanisms. Our analysis delivers a comprehensive picture of Li mobility, local structural environment, and partial charges in bulk C and Si, as well as at the interface. Finally, we single out the parameters governing Li diffusivity and discuss their implications for the design of Si/C composite materials.

2 Computational methods

2.1 Density functional theory parameters

The CP2K^57,58 program package was employed for all density functional theory (DFT) calculations.^59,60 Within CP2K's Quickstep (QS) module,^58,61 the Gaussian plane wave (GPW) method⁶² solves the Kohn–Sham (KS) equations⁶⁰ by combining an atom-centered Gaussian basis for the orbitals with plane waves for the electron density. A DZVP-MOLOPT-SR-GTH⁶³ basis set was used in conjunction with GTH-PBE pseudopotentials.^64–66 The plane-wave energy cutoff was set to 350 Ry, with a relative cutoff of 40 Ry. Since all systems studied here are sufficiently large, only the Γ point was sampled. Exchange–correlation (XC) effects were described using the Perdew–Burke–Ernzerhof (PBE) functional,^67,68 and dispersion interactions were treated via the DFT-D3 scheme⁶⁹ with Becke–Johnson damping⁷⁰ and revised damping parameters.⁷¹ To achieve self-consistent field (SCF) convergence to 10⁻⁶, orbital transformation (OT)^72,73 was employed with a DIIS minimizer and the FULL_KINETIC preconditioner. Periodic boundary conditions (PBCs) were imposed in all three Cartesian directions. Because our models are slabs, they are truly periodic only in x and y direction, with x and y cell parameters determined as described in section 2.2. In the z direction, a 20 Å vacuum region was added to the slab thickness to prevent interactions between periodic images.

2.2 γ-Graphdiyne slab generation

Model systems for Si/C interfaces were constructed using γ-graphdiyne (γ-GDY) as the carbon component. Graphynes, which are a relatively new class of carbon allotropes and remain challenging to synthesize,^74–76 have been demonstrated experimentally to constitute promising anode materials for LIBs.^76–78 In this study, γ-GDY was employed as a layered carbon model exhibiting three-dimensional Li diffusivity. Although graphite is often incorporated into experimentally studied Si/C composites, Li migration through a defect-free graphite sheet is hindered by an energy barrier of approximately 8 eV,⁴⁶ rendering the process effectively unobservable within 100 ps MD simulations at moderate temperatures and confining Li mobility to two dimensions between layers. To facilitate three-dimensional diffusion, defects could be introduced into the graphene sheets of graphite; here, however, a layered carbon with intrinsically large pores was chosen. Comparable migration barriers for Li diffusion both parallel to and through graphyne⁷⁹ and graphdiyne⁸⁰ layers have been reported by previous DFT-NEB investigations.

ABC-type stacking of γ-GDY was determined by Matsuoka et al. via transmission electron microscopy (TEM), selected-area electron diffraction (SAED), and 2D grazing-incidence wide-angle X-ray scattering (2D GIWAXS) measurements, with lattice parameters of a = b = 9.6 Å and c = 10.2 Å (interlayer distance 3.4 Å), and angles α = β = 90° and γ = 120°.⁸¹ Based on these experimental values, a single γ-GDY layer was initially constructed as a 2 × 2 supercell containing 72 atoms. To determine the optimal lattice constant for subsequent simulations, geometry optimizations of this single layer were performed using a Broyden–Fletcher–Goldfarb–Shanno (BFGS) minimizer.^82–85 Atomic positions were fully relaxed under fixed cell dimensions x = y, beginning from the experimentally reported value of 19.2 Å.⁸¹ Relaxations were carried out over the range 18.2 Å to 20.2 Å in 0.1 Å increments (energy minimum at 18.9 Å), then refined between 18.85 Å and 18.95 Å in 0.01 Å increments (min. at 18.94 Å), and finally between 18.935 Å and 18.945 Å in 0.001 Å increments. From these calculations, the minimum-energy lattice parameter for the single 2 × 2 γ-GDY layer was found to be x = y = 18.937 Å, which is in reasonable agreement with the experimental value of 19.2 Å.

A slab was then constructed from three γ-GDY layers, arranged in ABC stacking with an initial interlayer distance of 3.4 Å as determined experimentally.⁸¹ Atomic positions were first relaxed within x = y = 18.937 Å under the constraint that the z coordinates were held fixed. In a second relaxation step, all atomic coordinates were fully relaxed without any constraints, resulting in a final interlayer spacing of 3.29 Å. The resulting 216-atom γ-GDY slab was used to build the Si/C interfaces described in section 2.4.

2.3 Amorphous silicon slab generation

For the silicon component of our interface models, amorphous Si (a-Si) was generated independently from the graphdiyne slab by melting and then rapidly quenching a box with initially crystalline Si, using force field MD (FFMD) simulations with the LAMMPS program package.⁸⁶ Eleven force fields were rigorously assessed for their ability to reproduce experimentally determined properties of amorphous Si. This benchmarking procedure is described in detail in Section S1 of the SI. The second parametrization of the empirical interatomic potential proposed by Tersoff⁸⁷ was selected for the simulations.

Starting with the conventional unit cell of crystalline Si (c-Si, a = 5.4304 Å, α = 90°),⁸⁸ a 3 × 3 × 2 supercell containing 144 atoms was constructed. After initial geometry optimization with the Tersoff potential, the structure was equilibrated in a microcanonical (NVE) ensemble for 500 ps with the temperature maintained at 300 K by velocity rescaling. The temperature was then increased to 3400 K in a canonical (NVT) ensemble to melt c-Si over 2 ns (+1.55 K ps⁻¹). Temperature control was provided by a Nosé–Hoover chain thermostat^89–91 with a damping parameter of 100 fs. Next, the density of the Si melt was determined at 1 atm in a constant NPT simulation of 144 ns. Pressure control was achieved using a Nosé–Hoover chain barostat with a damping parameter of 2 ps. The system volume was averaged over the final 100 ps of the NPT simulation, yielding a density of 2.31 g cm⁻³ from a volume of 2909 ų. To match the γ-GDY slab described in section 2.2, the cubic system was deformed to dimensions of x = y = 18.937 Å and z = 9.365 Å (computed according to the determined density), with α = β = 90° and γ = 120°. This deformation was applied during 5 ns of NVT simulation. During the following 5 ns, the z dimension was expanded to 30 Å to create a slab with approximately 20 Å of vacuum space between periodic images. The molten Si slab was equilibrated for 225 ns and then rapidly quenched to 300 K over 310 ps (−10 K ps⁻¹) to produce an amorphous silicon slab. The a-Si slab was then equilibrated for an additional 225 ns, and after geometry optimization at the DFT level (see section 2.1), the fully relaxed atomic coordinates were used for Si/C interface construction as described in section 2.4. All FFMD simulations were performed with a time step of 1 fs under PBCs in all directions. A second a-Si slab was generated using the same protocol, with the melt equilibration extended by 100 ns.

2.4 Interface construction, AIMD and NEB methodology

In sections 2.2 and 2.3, the preparation of a γ-GDY slab and two a-Si slabs, each having dimensions x = y = 18.937 Å, α = β = 90° and γ = 120°, was described. The slab thicknesses in the z direction were 6.58 Å for γ-GDY, and 9.46 Å and 8.74 Å for the first and second a-Si slabs, respectively. These two subsystems were then combined to form the Si/C interface models. By varying the initial separation, defined as the difference in z coordinates between the closest Si atom and the C layer, and performing subsequent geometry optimizations, several interface configurations exhibiting different numbers of Si–C bonds were obtained. Following the same energy minimization protocol with decreasing step size as detailed in section 2.2 for finding x and y, three low-energy initial Si–C slab separations were identified for each Si slab: 0.254 Å, 0.544 Å, and 2.251 Å for the first slab; and 0.636 Å, 0.842 Å, and 3.427 Å for the second. Note that these are separations prior to relaxation. At each of these six separations, one γ-GDY slab was joined with one of the two a-Si slabs, resulting in six distinct 360-atom Si/C interface models. This ensemble of models improves the statistical robustness of our simulations and permits correlation of physical properties with interfacial separation.

After full relaxation of the atomic coordinates, four Li atoms were placed in the interfacial γ-GDY layer at the maximum possible Li–Li separation of 9.47 Å to minimize electrostatic interactions. The geometry was then reoptimized. All DFT calculations involving four Li atoms were performed with a total charge of +4 to represent an ionic state. The relaxed atomic configurations of the single γ-GDY layer, the γ-GDY slab, the two a-Si slabs, and the six Si/C interface model systems with and without Li atoms are given as SI⁹² (.xyz files). Additionally, the preparation procedure described in sections 2.2 to 2.4 and the relaxed Li/Si/C models are illustrated in Fig. S1 and S2 of the SI, respectively.

Ab initio MD (AIMD) simulations of the relaxed Li/Si/C systems were carried out in an NVT ensemble at 400 K. Temperature was controlled with a Nosé–Hoover chain thermostat, and a time step of 1 fs was used. The equilibration protocol comprised three stages: (1) massive thermostatting for 1 ps with a time constant of 10 fs; (2) massive thermostatting for 1 ps with a time constant of 100 fs; and (3) global thermostatting for 1 ps with a time constant of 100 fs.

After equilibration, the four Li atoms were pulled gradually from γ-GDY into a-Si one after another by applying a harmonic restraint k·(z − z_F)² on the z coordinate of Li. Initially, z_F was set to z₀–0.5 Å, where z₀ is the Li z coordinate after equilibration. Every 500 fs, z_F was reduced by an additional 0.5 Å, that is, z_F = z₀ − j·0.5 Å, with the step number j. This caused the spring force F_z = 2k·(z − z_F) to increase progressively. The force constant k was chosen as 0.001E_h/a₀² to ensure that Li velocities did not exceed the maximum equilibrium values of approximately 0.0015a₀E_h/ℏ. Thermostatting during the pulling simulations followed the same global scheme as in equilibration stage (3). With this approach, migration of a single Li atom from C into Si took between 6 ps and 25 ps.

Climbing image⁹³ nudged elastic band^94–97 (CI-NEB) calculations were carried out for all Li diffusion pathways identified in the AIMD trajectories. The initial and final atomic configurations of each path were fully relaxed prior to mapping the path with eight replicas. During the NEB calculation, the atomic coordinates of the first and last images were held fixed.

2.5 Analysis and visualization of ab initio data

For the analysis, mean square displacements (MSDs) and diffusion coefficients D were calculated using TRAVIS.^98,99 D was determined from the MSD(τ) according to the following relation:

The MSD of the particle type of interest, in this case Li, in an MD simulation is a function of the correlation time τ. Denoting the position vector of particle i at time t by r⃑_i(t) and at time t + τ by r⃑_i(t + τ), the MSD is obtained by averaging |r⃑_i(t + τ) − r⃑_i(t)|² over all particles i and over all time origins t. D is then given by the slope of MSD(τ) in the limit of large τ. In practice, D is extracted by performing a linear regression of MSD(τ) versus τ within the regime where the MSD grows linearly.

When a harmonic restraint is applied to Li, the motion of the Li atom becomes anisotropic and is characterized by the mobility μ instead of the diffusion coefficient D:

In this equation, v is the particle velocity, E is the hypothetical electric field strength, q is the particle charge (assumed to be +1 for Li), F is the spring force defined above, z_i and z_f are the initial and final z coordinates of the particle, and Δt is the observation time interval. For each 500 fs interval, z is approximated by z_i, which keeps F constant and allows μ to be calculated. To obtain a single mobility value for the simulation, μ is averaged over all such intervals. Finally, diffusivity and mobility are connected by the Einstein–Smoluchowski equation:

Li coordination in MD trajectories is assessed using radial distribution functions (RDFs), which are computed with TRAVIS.^98,99 An RDF g_ab expresses the probability of finding a pair of particles – a reference particle (Li) and an observed particle (C or Si) – separated by a distance r, relative to the probability expected for a uniform distribution of the observed particles. RDFs are defined by the following equation:

In this equation, i and j label the reference and observed particles, respectively, and N_a and N_b are their total numbers. The vectors r⃑_i(t) and r⃑_j(t) give the positions of particles i and j at simulation time t, and V is the volume of the simulation cell. The discretized δ function is equal to 1/2w when the separation |r⃑_i(t) − r⃑_j(t)| falls within the interval [r − w, r + w], where w is the bin width; otherwise, it is zero.

Mulliken partial atomic charges¹⁰⁰ were calculated using the implementation in CP2K. Structural models were visualized with VMD.¹⁰¹

3 Results and discussion

3.1 Li mobility, diffusion mechanism and migration barriers

To investigate Li diffusion at the interfaces between γ-GDY and a-Si, we performed AIMD simulations following the protocol outlined in section 2.4 and depicted in Fig. 1(a). Throughout the MD trajectories, the Li/Si/C model systems remained structurally stable: the atomic configurations of both the carbon and silicon slabs showed no significant changes, and the Si–C interfacial distance remained essentially constant. The Li diffusion coefficients D_Li at 400 K were calculated from MSDs (eqn (1)) for free diffusion within the bulk materials, and from eqn (2) and (3) for Li pulled across the interface by a harmonic restraint. The corresponding results are presented in Fig. 1(b).

Fig. 1 (a) Schematic illustration of the AIMD simulation methodology described in section 2.4. One Li atom is pulled from γ-GDY into a-Si by applying a spring force F_z, while the other three atoms move freely. (b) Diffusion coefficients D_Li (in Ų ps⁻¹) at 400 K for Li atoms moving freely in bulk C (black) and bulk Si (blue), determined from various AIMD simulations via the MSDs (eqn (1)), and for Li atoms pulled across the interface (green), obtained using eqn (2) and (3). Median values are indicated by solid lines and labeled, and mean values are shown by dashed lines.

Li diffusivities in bulk γ-GDY range from 0 to 0.768 Ų ps⁻¹, with a median of 0.173 Ų ps⁻¹. In bulk a-Si the median D_Li is 0.004 Ų ps⁻¹, with individual values between 0 and 0.429 Ų ps⁻¹. The diffusion coefficients from Li mobilities across the interface range from 0.013 to 0.381 Ų ps⁻¹, with a median of 0.065 Ų ps⁻¹. The average Li diffusivity values are 0.186, 0.027, and 0.083 Ų ps⁻¹, respectively, and are therefore slightly higher than the corresponding medians, indicating right-skewed distributions. In bulk γ-GDY the standard deviation is 0.160 Ų ps⁻¹, which is more than twice the values of 0.071 and 0.076 Ų ps⁻¹ in bulk a-Si and across the interface, respectively.

Although Li diffusion coefficients in all three regions range from very small values – even zero in the bulk materials – to comparatively high outliers, the distribution of values shows greater variation in bulk γ-GDY, indicating more heterogeneous Li mobility in this material. Overall, the median Li diffusivity in bulk γ-GDY is one to two orders of magnitude higher than in bulk a-Si, which shows that there is a significant drop in D_Li when a Li atom diffuses from carbon into silicon in a Si/C composite anode material. The Li mobility across the interface is intermediate between the two bulk regions, demonstrating that the interface does not impede Li migration and that the bulk Si component ultimately governs the overall Li diffusion kinetics.

To validate our results, we compare the computed Li diffusivities, ≈10⁻⁵ cm² s⁻¹ in bulk γ-GDY and ≈10⁻⁷ cm² s⁻¹ in bulk a-Si, with values reported in the literature. For graphite, a common component of Si/C composite anode materials, the experimentally measured in-plane diffusivity is 4.4 × 10⁻⁶ cm² s⁻¹ (Persson et al.¹⁰²), while a value of ≈10⁻⁷ cm² s⁻¹ is reported from simulations (Persson et al.¹⁰²). Our D_Li in γ-GDY is therefore roughly one order of magnitude larger, which is consistent with Li migration barriers in γ-GDY being ≈0.1 eV lower than in graphite (see below). In a-Si, the Li diffusion coefficient has been measured as 5.1 × 10⁻¹² cm² s⁻¹ (Ding et al.²⁴) and estimated from first principles as ≈10⁻¹² cm² s⁻¹ (Tritsaris et al.³²). Note that Yan et al. further demonstrate that D_Li in computational studies of a-Si can span a very wide range from 10⁻¹⁶ to 10⁻⁸ cm² s⁻¹,³⁶ with our values again at least one order of magnitude larger. In summary, although the AIMD-derived diffusivities reported here may somewhat overestimate more realistic experimental and KMC values, most likely because of the limited timescales accessible to AIMD and different measurement/simulation temperatures, our simulations reproduce the substantial (several orders of magnitude) decrease in D_Li from carbon to silicon that is observed in both experiment and modelling.

To further elucidate Li diffusion across γ-GDY/a-Si interfaces at the atomistic level, we analyzed the diffusion mechanism in our AIMD simulations. We found that Li migration from C into Si proceeds by hopping between relatively stable, low-energy positions. These hops can be grouped into six distinct pathways constituting the following migration mechanism, illustrated for one example Li atom in Fig. 2(a) and in Fig. S3 in the SI: (1) between the upper and middle graphdiyne layers, (2) through the middle graphdiyne layer, (3) between the middle and lower graphdiyne layers, (4) across the Si/C interface, (5) from the interface into a-Si, and (6) within a-Si. Hereafter, we refer to these as migration paths (1) to (6), where paths (1) to (3) occur within bulk C, path (4) spans the interface, and paths (5) and (6) occur within bulk Si. Although each Li atom could in principle traverse all six paths, some simulations omit certain paths, particularly paths (1) and (6), depending on the initial and final Li positions. In a few cases, there is an additional stable Li site within a pathway, leading to two hops along the same path during a simulation. Conversely, when no stable site exists between two paths, they occur consecutively at the same time. For all subsequent analyses, any merged paths are treated as a single data point for each constituent pathway.

Fig. 2 (a) Exemplary illustration of the Li diffusion mechanism across the interface, showing six distinct migration pathways indicated by arrows and color-coded. (b) Energy barriers for Li migration ΔE^M (in eV) along diffusion paths (1) to (6) from NEB calculations. Data points for C → Si migration are shown as ▽, and those for Si → C migration as △. Median values for C → Si are indicated by ↓ and labeled, and median values for Si → C are shown by ↑. (c) Li migration distances a (in Å) for paths (1) to (6) from NEB calculations. Median values are indicated by solid lines and labeled, and mean values are shown by dashed lines.

The Li hops observed in our AIMD simulations were further examined using NEB calculations as described in section 2.4. Each NEB calculation directly provides two energy barriers for Li migration ΔE^M (in C → Si and Si → C directions), as well as a series of fully relaxed atomic configurations along the path, which are given as SI⁹² (.xyz files) together with the respective energy profiles (.txt files). Additionally, the energy profiles are illustrated for one example Li atom in Fig. S4 in the SI. From the initial and final positions of the hopping Li atom, we extracted the migration distance a. The energy barriers and migration distances for all Li hops identified across the various AIMD simulations are presented in Fig. 2(b) and (c).

In the C → Si direction, the minimum Li migration barriers for paths (1) to (6) are all near zero, while the maximum barriers increase significantly from 0.31 to 0.55, 0.85, 1.51, 1.91 and 1.88 eV, respectively. The median barriers increase accordingly: 0.09, 0.17, 0.23, 0.36, 0.48 and 0.69 eV. In the Si → C direction, the minimum barriers are zero for paths (1) to (5) and 0.10 eV for path (6), with maximum values of 0.09, 0.38, 0.38, 1.69, 1.69 and 1.68 eV, respectively. The median values again increase sequentially from 0.03 to 0.04, 0.06, 0.08, 0.15 and 0.54 eV. The corresponding Li hopping distances span 0.90–3.67 Å for path (1) (median ã = 1.95 Å), 1.57–6.29 Å for path (2) (ã = 2.72 Å), 2.08–6.92 Å for path (3) (ã = 2.76 Å), 1.60–6.92 Å for path (4) (ã = 3.18 Å), 1.49–6.45 Å for path (5) (ã = 3.09 Å) and 1.27–4.94 Å for path (6) (ã = 3.26 Å).

For comparison, Tritsaris et al. reported energy barriers for Li migration in bulk a-Si ranging from 0.1 to 2.4 eV, with an average of 0.58 eV.³² For graphdiyne, Sun et al. calculated minimum barriers of 0.18 eV for in-plane migration and 0.17 eV for out-of-plane diffusion.⁸⁰ Compared with the corresponding values for path (6) in a-Si and paths (1) and (2) in graphdiyne, these results indicate that our slab models are sufficiently large to approximate bulk diffusion. Furthermore, given reported migration barriers of 0.28 to 0.30 eV for in-plane diffusion in graphite,¹⁰² graphdiyne can be regarded as a reasonable model system for Li migration in layered carbon materials.

To evaluate Li mobility based on the observed migration barriers, the maximum values can safely be neglected, since high-energy Li hops rarely occur according to the Boltzmann factor exp(−ΔE^M/kT). Although very-low-energy barriers were observed for each path in our simulations, their actual availability depends strongly on the local structural environment of each Li atom. To assess the overall diffusivity of an ensemble of Li atoms, it is therefore most useful to discuss the median values of ΔE^M. These median barriers increase continuously in the C → Si direction, in agreement with the results presented above for D_Li, and reflect a progressive immobilization of Li atoms as they migrate from γ-GDY across the interface into a-Si. Conversely, the barriers decrease continuously in the Si → C direction. It is noteworthy that for each path the ΔE^M values are asymmetric between the two directions, indicating that Li atoms on average have higher potential energies in the final configurations than in the initial ones. This means that migration from C into Si is thermodynamically endothermic. From the median barriers, the overall energy gain of a Li atom between paths (1) and (6) is roughly 1.12 eV. This asymmetry is largest for paths (4) and (5), which together account for 0.61 eV, indicating that the endothermic step mainly occurs at the interface when Li atoms leave one material and enter the other.

Regarding Li hopping distances, these increase steadily from paths (1) to (4) and then remain relatively constant for paths (5) and (6). Therefore, the increase in migration barriers can be attributed primarily to the longer hop distances for the first four paths (see section 3.4), which describe bulk graphdiyne migration and interfacial migration. In contrast, bulk a-Si migration involves additional structural effects that further increase the energy barriers. These effects are discussed in the following sections. Overall, Li migration in bulk γ-GDY (paths (1) to (3)) consists of hops with median barriers of approximately 0.1–0.2 eV and distances of 1.9–2.8 Å, followed by migration across the interface (path (4)) with barriers around 0.4 eV and distances around 3.2 Å, and finally by migration in bulk a-Si (paths (5) and (6)) with barriers of 0.5–0.7 eV and distances of 3.1–3.3 Å. The barriers in the reverse migration direction are 0.1–0.3 eV lower.

In summary, in bulk graphdiyne, Li diffusivity is high because migration barriers are low and hopping distances are relatively short. In contrast, in bulk amorphous silicon, migration barriers are high and hopping distances are longer, yielding Li diffusion coefficients that are one to two orders of magnitude lower than in the carbon component. Li migration from C into Si is endothermic and proceeds via a six-step diffusion mechanism in our systems. Li mobilities are intermediate between those of the two bulk materials, reflecting a continuous increase in migration barriers in the C → Si diffusion direction and indicating progressive Li immobilization. Migration across the Si/C interface itself involves relatively long hopping distances but only moderate energy barriers, so it is not the rate-limiting step for overall diffusion. From the macroscopic perspective of a Si/C composite anode material, and assuming that results for γ-GDY can be transferred to other layered carbons, Li atoms migrate relatively freely within the carbon component, while the silicon component governs the overall diffusion kinetics. Because the interface does not constitute a diffusion bottleneck, its nature does not need to be specifically targeted in material design.

3.2 Local structural environment of Li

To further relate Li diffusion at γ-GDY/a-Si interfaces to structural properties, we analyzed the local bonding environment of Li by calculating Li–C and Li–Si RDFs from our MD simulations using eqn (4). The first maximum of each RDF r^[1]_Li–X corresponds to the most frequent Li–X separation observed in the trajectory and can therefore be interpreted as the bond distance. The values of these first maxima for Li atoms diffusing freely in bulk C and in bulk Si are presented in Fig. 3(b) for all of our AIMD simulations.

Fig. 3 (a) Schematic illustration of the AIMD simulation methodology described in section 2.4. One Li atom is pulled from γ-GDY into a-Si by applying a spring force F_z, while the other three atoms move freely. (b) Li coordination, given by the position of the first peaks of the Li–C and Li–Si RDFs r^[1]_Li–X (in Å) at 400 K, for Li atoms moving freely in bulk C (black) and bulk Si (blue). RDFs were determined from various AIMD simulations via eqn (4). Median values are indicated by solid lines and labeled, and mean values are shown by dashed lines. (c) Median Mulliken charges of Li q̃_Li (in units of e) at 400 K for Li atoms moving freely in bulk C (black) and bulk Si (blue), and for Li atoms pulled across the interface (green). Values were obtained from various AIMD simulations. Median values are indicated by solid lines and labeled, and mean values are shown by dashed lines.

Li–C distances range from 2.35 to 2.58 Å, with both the median and the mean at 2.44 Å and a standard deviation of 0.05 Å. Li–Si distances span 2.48 to 3.32 Å, with a median of 2.68 Å. The mean value of 2.76 Å reflects a right-skewed distribution caused by a few outliers at larger separations, consistent with a larger standard deviation of 0.21 Å. These outliers arise from Li atoms that migrated to the lower vacuum surface of the Si slab and are only weakly bonded there, which also leads to the higher D_Li values observed in Fig. 1(b). Such atoms are not truly bonded in bulk Si and should be excluded from discussion, as they are artifacts of our model system setup.

Considering the minimum values of 2.35 Å for Li–C and 2.48 Å for Li–Si, alongside the medians of 2.44 Å and 2.68 Å, respectively, and noting that the covalent radius of sp²/sp-hybridized carbon (0.73/0.69 Å (ref. 103)) is roughly 0.4 Å smaller than that of silicon (1.11 Å (ref. 103)), it is evident that Li–Si bonds are effectively about 0.2 Å shorter than Li–C bonds in our systems. This difference could imply either a stronger, more covalent Li–Si interaction or that Li is more “squeezed” into a-Si than into γ-GDY. The results presented in section 3.1, showing that Li positions in Si are energetically less favorable than in C, support the latter interpretation. Thus, Li is likely accommodated without strain between the graphdiyne layers, whereas the a-Si lattice must be deformed more to host Li atoms, resulting in energetically less favorable, shorter Li–Si distances.

Thus far we have discussed the local environment of Li based on RDFs – that is, statistical results from our complete AIMD trajectories. Because Li hops are rare events in these trajectories, changes in the local structure during those hops are not captured by the RDFs. To analyze these events in more detail, we used the fully relaxed atomic configurations from our NEB calculations of Li migration pathways to determine the closest Li–C and Li–Si distances r^min_Li–X in the initial, transition state (TS), and final structures of each hop, assuming that these distances influence other relevant properties such as migration barriers. The extracted minimum Li–X distances are shown in Fig. 4(b).

Fig. 4 (a) Exemplary illustration of the Li diffusion mechanism across the interface, showing six distinct migration pathways indicated by arrows and color-coded. (b) Li coordination, expressed as the minimum Li–C and Li–Si distances r^min_Li–X (in Å) along diffusion paths (1) to (6) from NEB calculations. Data points for the initial and final atomic configurations are shown as ▷ and ◁, respectively, and those for the transition state configurations are shown as ◇. Median values for the transition state are indicated by and labeled, and median values for the initial and final configurations are shown by , respectively. (c) Mulliken charges of Li q_Li (in units of e) along diffusion paths (1) to (6) from NEB calculations. The symbols and median-value indicators are the same as in panel (b).

The median values of the closest Li–C and Li–Si distances along each migration path, using the notation initial structure → TS → final structure, are as follows: path (1), 2.32 → 2.31 → 2.31 Å; path (2), 2.33 → 2.30 → 2.33 Å; path (3), 2.32 → 2.30 → 2.32 Å; path (4), 2.32 → 2.27 → 2.28 Å; path (5), 2.32 → 2.44 → 2.50 Å; and path (6), 2.56 → 2.32 → 2.42 Å. Fig. S5 in the SI presents separate analyses for Li–C and Li–Si distances, showing that the minimum Li–X distance corresponds to a Li–C distance from the initial structure of path (1) through that of path (5), and to a Li–Si distance from the TS structure of path (5) onward.

For paths (1) to (4), corresponding to Li migration in bulk C and across the interface, all minimum distances are Li–C distances of about 2.3 Å and vary very little along each path. This indicates that Li hopping is not associated with unusually small interatomic distances at the TS, but instead proceeds readily between the graphdiyne layers without straining them. For path (5), the minimum Li–X distance is Li–C in the initial structure but becomes Li–Si in both the TS and the final structure. Given the 0.4 Å difference in covalent radii between C and Si, two points emerge: (i) the Li–X distances in the TS and final structure are effectively 0.2 to 0.3 Å shorter than in the initial structure, consistent with Li being “squeezed” into the a-Si lattice as seen in the RDF results; and (ii) the Li–Si distance is even slightly shorter in the TS than in the final structure, indicating that the lattice is strained more during migration between stable sites. The same effect is observed for path (6), where the Li–Si distance in the TS is 0.1 to 0.2 Å smaller than at the stable Li sites. When comparing the TS interatomic distances of paths (1) to (4) with those of paths (5) and (6), and correcting for the difference in covalent radii of Si and C, it becomes clear that in Si-dominated transition states the atoms are effectively about 0.3 to 0.4 Å closer than in C-dominated ones. This demonstrates that Li atoms must squeeze through the a-Si lattice not only relative to the stable sites, but also in comparison with migration through graphdiyne.

In summary, our results show that the carbon layers of γ-GDY accommodate Li atoms without significant strain, characterizing it as an intercalation compound comparable to graphite. Li–C interactions are generally weaker, as reflected by effectively longer Li–C distances that remain essentially unchanged during facile Li migration between the layers, indicating that the carbon lattice has little effect on the transition state energy. In bulk a-Si, Li insertion causes substantial deformation of the silicon lattice, leading to stronger Li–Si interactions and overall shorter Li–Si distances, considering the larger atomic size of Si. At the transition state during Li migration, Li–Si distances become even shorter than in the stable Li sites, indicating increased lattice strain and a pronounced influence of the Si environment on the transition state energy, which is discussed further in section 3.4. Li migration across the Si/C interface is governed by the same weak Li–C interactions with the local structural environment as in bulk graphdiyne, consistent with the moderate interfacial migration barriers.

3.3 Li charge dynamics

As a final physical quantity to characterize Li diffusion at γ-GDY/a-Si interfaces, we focused on the Li partial atomic charges. We calculated Mulliken charges along our AIMD trajectories using CP2K. While Mulliken charges are not real physical observables and should not be interpreted as rigorous charge states of atoms, they do give a complementary perspective on the chemistry of Li–C and Li–Si bonds. Note that the physical trends discussed hereafter are also observed qualitatively for other calculation methods for partial atomic charges that we have tested. In Fig. 3(c), the median Li charge from each trajectory q̃_Li is shown, distinguishing Li atoms diffusing freely in bulk C, those diffusing in bulk Si, and those moving across the interface. The median Li charges in bulk C range from +0.44 e to +0.54 e, with a median of +0.53 e, an average of +0.52 e, and a standard deviation of 0.02 e. In bulk Si, they range from +0.06 e to +0.73 e, with a median of +0.31 e, a mean of +0.35 e, and a standard deviation of 0.14 e. For Li atoms pulled across the interface, the charge values span +0.22 e to +0.51 e, with a median of +0.40 e, an average of +0.39 e, and a standard deviation of 0.07 e. As noted in section 3.2 for the Li–Si distances in Fig. 3(b), the very positive Li charges are outliers caused by loosely bound Li atoms at the vacuum surface of a-Si; these can be excluded from further discussion.

Observing the trends in q̃_Li, Li Mulliken charges decrease by about 0.2 e when diffusing from C into Si. This aligns with the results from section 3.2, and can be interpreted as stronger Li–Si interaction compared to Li–C. In addition to Li being more “squeezed” into the a-Si lattice than between graphdiyne layers, this can also be understood through electronegativity: Si is less electronegative than C, so the Li–Si bond is less polarized and Li less positively charged. Moreover, Li charges in bulk C are much more narrowly distributed than in bulk Si. This agrees with the Li–C and Li–Si distance data in Fig. 3(b), indicating that γ-GDY provides a structurally more homogeneous environment for Li than a-Si does. In contrast, Fig. 1(b) shows that Li diffusivities in γ-GDY are more heterogeneous, reflecting the combination of generally significantly lower migration barriers in that material and the limited time scale of the AIMD simulations.

To atomistically resolve Li charge dynamics during migration from C into Si, we computed Mulliken charges for the initial, transition state (TS), and final structures of all Li hops studied with NEB calculations. Fig. 4(c) presents the results for q_Li along diffusion paths (1) to (6). Using the notation initial structure → TS → final structure, the median Li charges (in e) are as follows: path (1), +0.53 → +0.53 → +0.50; path (2), +0.50 → +0.51 → +0.49; path (3), +0.48 → +0.50 → +0.48; path (4), +0.48 → +0.48 → +0.42; path (5), +0.41 → +0.29 → +0.30; and path (6), +0.33 → +0.25 → +0.29.

The trends in these values closely resemble those observed for the shortest Li–X distances, see Fig. 4(b). This underscores the clear link between these two quantities, which we examine further in section 3.4. For migration in bulk C, i.e. along paths (1) to (3), Li charges remain near +0.5 e, indicating little change in the weak Li–C interactions between stable sites and the TS. Along path (4), a first modest charge drop of about 0.05 e occurs between the TS and the final structure, even though Li–C interactions still dominate rather than Li–Si as discussed in section 3.2. Migration across the interface thus accounts for roughly one quarter of the total ≈0.2 e charge decrease seen between C and Si in both AIMD and NEB results. Along path (5), the second and largest charge drop of about 0.1 e takes place between the initial structure and the TS, reflecting the shift from primarily Li–C to primarily Li–Si interactions. Along path (6), the TS Li charge is about 0.05 e lower than at the stable sites, consistent with the shorter Li–Si distances in the TS. Overall, these observations confirm the two effects described in section 3.2: (i) Li charges are generally less positive in a-Si than in γ-GDY, corresponding to stronger Li–Si interactions, and (ii) only in bulk Si does the TS exhibit a smaller Li charge than the stable sites, a consequence of greater lattice deformation during migration.

In summary, Li charges in bulk graphdiyne are more homogeneous and generally more positive, showing minimal variation during migration. In contrast, Li charges in bulk a-Si are more heterogeneous and less positive, which reflects both the electronegativity concept and greater deformation of the Si lattice, and these charges decrease even further in the highly strained transition state during migration. The decrease in Li charge observed when migrating from C into Si proceeds in two steps. The first one involves a modest reduction as Li hops across the interface, and the second one involves a larger decrease as it diffuses from the interface into the a-Si lattice.

3.4 Structure–property relationships

So far, we have separately discussed several physical quantities that characterize Li diffusion at the γ-GDY/a-Si interfaces in sections 3.1 to 3.3, but we have only suggested the links between them. We analyzed all these properties for correlations. In this section, we present the correlations we found and establish several structure–property relationships for Li mobility in these systems. First, in Fig. 5, the energy barriers for Li migration ΔE^M, averaged over both directions for each hop, are plotted against the migration distances a. It is evident that for paths (1) through (6) – that is, for diffusion in the bulk materials as well as across the interface – the migration barrier is positively correlated with the hopping distance. Thus, the longer a Li atom must migrate between two stable sites, the higher the energy barrier tends to be, regardless of the component of the Si/C composite anode material.

Fig. 5 Li migration barriers ΔE^M (in eV), averaged over both C → Si and Si → C directions, plotted against Li hopping distances a (in Å) for (a) path (1), (b) path (2), (c) path (3), (d) path (4), (e) path (5), and (f) path (6).

Next, the Li migration barriers ΔE^M are plotted against the minimum Li–Si distances in the transition state r^min_Li–Si,TS for paths (3) to (6) in Fig. 6(a), and against the Li charges in the transition state q_Li,TS for paths (5) and (6) in Fig. 6(b). For migration in bulk a-Si, corresponding to paths (5) and (6), the energy barrier is negatively correlated with both the minimum Li–Si distance and the Li charge in the TS. In other words, when the TS is dominated by Li–Si interactions, as discussed in section 3.2, shorter interatomic distances and lower positive charge on Li both make the hop less energetically favourable. This trend could be explained by lattice deformation: the more the Li atom must “squeeze” through the Si network to pass the TS between two stable sites, the greater the strain imposed on the lattice.

Fig. 6 Li migration barriers ΔE^M (in eV), averaged over both C → Si and Si → C directions, plotted in (a) against the minimum Li–Si distances in the transition state r^min_Li–Si,TS (in Å) for paths (3) to (6), and in (b) against the Li Mulliken charges in the transition state q_Li,TS (in units of e) for paths (5) and (6).

To examine whether Li–Si distance and Li charge are correlated, q̃_Li in bulk a-Si from the AIMD simulations is plotted against the corresponding r^[1]_Li–Si in Fig. 7(a), and q_Li in the transition state is plotted against the respective r^min_Li–Si for paths (3) to (6) in Fig. 7(b). Both plots show a positive correlation between the two quantities for stable Li sites (AIMD) and TS (NEB): the closer Li and Si atoms are, the less positive becomes the Li charge. Although correlation does not imply causation, the consistency of these results suggests that there is an actual causal relationship between these two properties, and that their influence on the migration barriers discussed above reflects the same underlying effect. We therefore conclude that shorter Li–Si distances lead to less positive Li charges, and in the TS this implies higher migration barriers because of increased lattice strain.

Fig. 7 (a) Median Li Mulliken charges q̃_Li (in units of e) for Li atoms moving freely in bulk Si at 400 K, plotted against the position of the first peaks of the Li–Si RDFs r^[1]_Li–Si (in Å). (b) Li Mulliken charges in the transition state q_Li,TS (in units of e), plotted against the minimum Li–Si distances in the transition state r^min_Li–Si,TS (in Å) for paths (3) to (6).

Finally, the minimum Li–Si distances in the TS r^min_Li–Si,TS are plotted against the Li hopping distances a for paths (4) to (6) in Fig. 8. These plots reveal an additional correlation for Li migration both across the interface and within bulk a-Si: as the hopping distance increases, the minimum Li–Si distance in the TS tends to decrease. Because both of these effects were shown above to lead to higher migration barriers, it is reasonable to assume a causal link between them. Although this relationship is not immediately obvious, we suggest a probabilistic interpretation: the farther a Li atom must migrate between two stable sites in the a-Si network, the more likely it is to pass very close to a Si atom, which then defines the TS.

Fig. 8 Minimum Li–Si distances in the transition state r^min_Li–Si,TS (in Å), plotted against the Li hopping distances a (in Å) for (a) path (4), (b) path (5), and (c) path (6).

To summarize, we have identified several structure–property relationships governing Li mobility in our graphdiyne/a-Si interface model systems. First, migration barriers increase with hopping distance in every component of our systems: longer hops cause higher energy barriers. Second, in bulk a-Si, and to a lesser extent across the interface, longer hopping distances lead to smaller Li–Si separations in the transition state accompanied by a lower positive charge on Li, which further raises the migration barrier. We attribute these trends to variations in lattice strain imposed by the migrating Li atom on the a-Si network.

Conclusions

We have thoroughly investigated Li diffusion in Si/C composite anodes, represented by a suitable model system based on amorphous silicon and the carbon allotrope γ-graphdiyne. By combining ab initio molecular dynamics simulations with nudged elastic band calculations, we have characterized Li mobility and examined structure–property relationships with the local chemical environment and charge state of Li atoms. In the γ-GDY component, Li diffusivity is high due to generally low migration barriers. Li atoms intercalate between carbon layers without inducing significant strain, as reflected by more positive Li charges. During facile migration, the weak Li–C interactions change only minimally, so the carbon layers have little effect on the transition state energy. In contrast, the penetration of Li into the amorphous Si phase induces substantial lattice strain, accompanied by less positive Li charges. This stronger Li–Si interaction becomes even more pronounced in the highly strained transition state during an elementary Li migration step, where it determines the transition state energy. This results in high migration barriers and Li diffusion coefficients that are one to two orders of magnitude lower than those in the carbon component.

For Li migration from carbon into silicon, the Li mobility is reduced compared to the carbon phase, reflecting progressively higher energy barriers and increasing Li immobilization. During this process, the Li charge decreases in two steps: (i) while hopping across the interface, and (ii) during migration from the interface into the a-Si. Li hopping across the interface itself is governed by weak Li–C interactions in the transition state, and thus, exhibits only moderate energy barriers and does not limit the overall diffusion rate. Assuming that γ-GDY is a suitable model system for other layered carbons such as graphite, Li diffusion in Si/C composite anode materials remains relatively unhindered within the carbon material. The carbon–silicon interface does not constitute a diffusion bottleneck, so composite particle design does not need to focus specifically on its nature. Instead, because the silicon material governs the diffusion kinetics of the composite, nanostructural design should aim to minimize diffusion distances within bulk a-Si – for example, by reducing Si particle size in at least one dimension. This suggests the synthesis of Si/C anode materials with Si nanoparticles, nanowires or nanofilms.

We plan to extend the time and length scales accessible to materials modelling of Si/C composite anodes by employing machine-learned force fields, kinetic Monte Carlo simulations, and, ultimately, continuum modelling. Our goal is to develop multiscale models that provide a realistic description of Li diffusion kinetics in these materials under battery operating conditions. A key step in this direction is to incorporate the effects of further Li insertion, together with the associated volume expansion and mechanical stress, on Li diffusion, since the present study focused on the very initial phase of lithiation with only isolated Li atoms in a-Si. Moreover, to move beyond the limited AIMD time and length scales will require accounting for spatial and temporal variability in kinetic properties, for example, by introducing dispersion factors for their distributions in a-Si.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Data for this article, including fully relaxed atomic configurations of the studied systems and migration paths (.xyz files), as well as energy profiles for the latter (.txt files) are available at Zenodo at https://doi.org/10.5281/zenodo.16759836.

Supplementary information: additional figures of the studied systems, migration paths, energy profiles, and analyses of Li–C and Li–Si distances. See DOI: https://doi.org/10.1039/d5lf00225g.

Acknowledgements

This project was funded by the Deutsche Forschungsgemeinschaft, project ID 435886714.

References

1. J. Tarascon and M. Armand, Nature, 2001, 414, 359–367.
2. M. Armand and J. Tarascon, Nature, 2008, 451, 652–657.
3. R. A. Huggins, Advanced Batteries - Materials Science Aspects, Springer, Boston, 1st edn, 2009.
4. S. Schweidler, L. de Biasi, A. Schiele, P. Hartmann, J. Brezesinski and T. Janek, J. Phys. Chem. C, 2018, 122, 8829–8835.
5. N. A. Kaskhedikar and J. Maier, Adv. Mater., 2009, 21, 2664–2680.
6. J. Hassoun, F. Bonaccorso, M. Agostini, M. Angelucci, M. G. Betti, R. Cingolani, M. Gemmi, C. Mariani, S. Panero, V. Pellegrini and B. Scrosati, Nano Lett., 2014, 14, 4901–4906.
7. Y. Wang, W. Tian, L. Wang, H. Zhang, J. Liu, T. Peng, L. Pan, X. Wang and M. Wu, ACS Appl. Mater. Interfaces, 2018, 10, 5577–5585.
8. U. Kasavajjula, C. Wang and A. J. Appleby, J. Power Sources, 2007, 163, 1003–1039.
9. X. Liu, X. Zhu and D. Pan, R. Soc. Open Sci., 2018, 5, 172370.
10. X. Li, M. Zhang, S. Yuan and C. Lu, ChemElectroChem, 2020, 7, 4289–4302.
11. P. Limthongkul, Y. Jang, N. J. Dudney and Y. Chiang, Acta Mater., 2003, 51, 1103–1113.
12. M. N. Obrovac and L. Christensen, Electrochem. Solid-State Lett., 2004, 7, A93–A96.
13. J. Li and J. R. Dahn, J. Electrochem. Soc., 2007, 154, A156–A161.
14. M. Gu, Z. Wang, J. G. Connell, D. E. Perea, L. J. Lauhon, F. Gao and C. Wang, ACS Nano, 2013, 7, 6303–6309.
15. Y. Xu, G. Yin, Y. Ma, P. Zuo and X. Cheng, J. Mater. Chem., 2010, 20, 3216–3220.
16. Y.-X. Yin, S. Xin, L.-J. Wan, C.-J. Li and Y.-G. Guo, J. Phys. Chem. C, 2011, 115, 14148–14154.
17. B. Wang, X. Li, X. Zhang, B. Luo, Y. Zhang and L. Zhi, Adv. Mater., 2013, 25, 3560–3565.
18. X. Y. Zhou, J. J. Tang, J. Yang, J. Xie and L. L. Ma, Electrochim. Acta, 2013, 87, 663–668.
19. H. Jo, J. Kim, D.-T. Nguyen, K. K. Kang, D.-M. Jeon, A.-R. Yang and S.-W. Song, J. Phys. Chem. C, 2016, 120, 22466–22475.
20. S. Guo, X. Hu, Y. Hou and Z. Wen, ACS Appl. Mater. Interfaces, 2017, 9, 42084–42092.
21. S. Zhu, J. Zhou, Y. Guan, W. Cai, Y. Zhao, Y. Zhu, L. Zhu, Y. Zhu and Y. Qian, Small, 2018, 14, 1802457.
22. W. M. Dose, M. J. Piernas-Muñoz, V. A. Maroni, S. E. Trask, I. Bloom and C. S. Johnson, Chem. Commun., 2018, 54, 3586–3589.
23. G. Hou, B. Cheng, Y. Yang, Y. Du, Y. Zhang, B. Li, J. He, Y. Zhou, D. Yi, N. Zhao, Y. Bando, D. Golberg, J. Yao, X. Wang and F. Yuan, ACS Nano, 2019, 13, 10179–10190.
24. N. Ding, J. Xu, Y. Yao, G. Wegner, X. Fang, C. Chen and I. Lieberwirth, Solid State Ionics, 2009, 180, 222–225.
25. S. Dupke, T. Langer, R. Pöttgen, M. Winter, S. Passerini and H. Eckert, Phys. Chem. Chem. Phys., 2012, 14, 6496–6508.
26. Z. Wang, M. Gu, Y. Zhou, X. Zu, J. G. Connell, J. Xiao, D. Perea, L. J. Lauhon, J. Bang, S. Zhang, C. Wang and F. Gao, Nano Lett., 2013, 13, 4511–4516.
27. W. Wan, Q. Zhang, Y. Cui and E. Wang, J. Phys.: Condens. Matter, 2010, 22, 415501.
28. H. Kim, K. E. Kweon, C.-Y. Chou, J. G. Ekerdt and G. S. Hwang, J. Phys. Chem. C, 2010, 114, 17942–17946.
29. C.-Y. Chou, H. Kim and G. S. Hwang, J. Phys. Chem. C, 2011, 115, 20018–20026.
30. P. Johari, Y. Qi and V. B. Shenoy, Nano Lett., 2011, 11, 5494–5500.
31. H. Kim, C.-Y. Chou, J. G. Ekerdt and G. S. Hwang, J. Phys. Chem. C, 2011, 115, 2514–2521.
32. G. A. Tritsaris, K. Zhao, O. U. Okeke and E. Kaxiras, J. Phys. Chem. C, 2012, 116, 22212–22216.
33. O. I. Malyi, T. L. Tan and S. Manzhos, Appl. Phys. Express, 2013, 6, 027301.
34. P. Haldar and A. Chatterjee, Energy Procedia, 2014, 54, 310–319.
35. Z. Wang, Q. Su, H. Deng, W. He, J. Lin and Y. Q. Fu, J. Mater. Chem. A, 2014, 2, 13976–13982.
36. X. Yan, A. Gouissem and P. Sharma, Mech. Mater., 2015, 91, 306–312.
37. Z. Wang, Q. Su, H. Deng and Y. Fu, ChemElectroChem, 2015, 2, 1292–1297.
38. J. Moon, B. Lee, M. Cho and K. Cho, J. Power Sources, 2016, 328, 558–566.
39. G. Wang, B. Xu, J. Shi, M. Wu, H. Su and C. Ouyang, Nanoscale, 2019, 11, 14042–14049.
40. C. Kirsch, C. Dreßler and D. Sebastiani, J. Phys. Chem. C, 2022, 126, 12136–12149.
41. C. Kirsch, C. Dreßler and D. Sebastiani, JPhys Energy, 2025, 7, 025003.
42. M. D. Levi and D. Aurbach, J. Phys. Chem. B, 1997, 101, 4641–4647.
43. T. L. Kulova, A. M. Skundin, E. A. Nizhnikovskii and A. V. Fesenko, Russ. J. Electrochem., 2006, 42, 259–262.
44. H.-J. Guo, X.-H. Li, X.-M. Zhang, H.-Q. Wang, Z.-X. Wang and W.-J. Peng, Xinxing Tan Cailiao, 2007, 22, 7–10.
45. K. Toyoura, Y. Koyama, A. Kuwabara, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 214303.
46. S. Thinius, M. M. Islam, P. Heitjans and T. Bredow, J. Phys. Chem. C, 2014, 118, 2273–2280.
47. E. G. Leggesse, C.-L. Chen and J.-C. Jiang, Carbon, 2016, 103, 209–216.
48. C.-Y. Chou and G. S. Hwang, J. Phys. Chem. C, 2013, 117, 9598–9604.
49. Z. Zhang, N. Liao, H. Zhou and W. Xue, Acta Mater., 2019, 178, 173–178.
50. F. Zhou, N. Liao, M. Zhang and W. Xue, Appl. Surf. Sci., 2019, 463, 610–615.
51. X. Gao, W. Lu and J. Xu, ACS Appl. Mater. Interfaces, 2021, 13, 21362–21370.
52. S. B. Olou'ou Guifo, J. E. Mueller, D. van Duin, M. K. Talkhoncheh, A. C. T. van Duin, D. Henriques and T. Markus, J. Phys. Chem. C, 2023, 127, 2818–2834.
53. A. Lashani Zand, A. Niksirat, Z. Sanaee and M. Pourfath, ACS Omega, 2023, 8, 44698–44707.
54. K. Odbadrakh, N. W. McNutt, D. M. Nicholson, O. Rios and D. J. Keffer, Appl. Phys. Lett., 2014, 105, 053906.
55. M. Wang, C. X. Yang, X. Y. Leng, Y. Chen, S. B. Yang, W. Li, W. Hong and Y. Xu, Int. J. Quantum Chem., 2024, 124, e27343.
56. L. Si, R. Huang, F. Liu, Q. Chen, P. Ouyang and Z. Dou, ACS Omega, 2025, 10, 19502–19509.
57. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
58. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
59. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
60. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
61. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
62. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–488.
63. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
64. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
65. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
66. M. Krack, Theor. Chem. Acc., 2005, 114, 145–152.
67. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
68. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
69. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
70. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
71. D. G. A. Smith, L. A. Burns, K. Patkowski and C. D. Sherrill, J. Phys. Chem. Lett., 2016, 7, 2197–2203.
72. J. VandeVondele and J. Hutter, J. Chem. Phys., 2003, 118, 4365–4369.
73. V. Weber, J. VandeVondele, J. Hutter and A. M. N. Niklasson, J. Chem. Phys., 2008, 128, 084113.
74. G. Li, Y. Li, H. Liu, Y. Guo, Y. Li and D. Zhu, Chem. Commun., 2010, 46, 3256–3258.
75. Q. Li, Y. Li, Y. Chen, L. Wu, C. Yang and X. Cui, Carbon, 2018, 136, 248–254.
76. C. Yang, Y. Li, Y. Chen, Q. Li, L. Wu and X. Cui, Small, 2019, 15, 1804710.
77. S. Zhang, H. Liu, C. Huang, G. Cui and Y. Li, Chem. Commun., 2015, 51, 1834–1837.
78. C. Huang, S. Zhang, H. Liu, Y. Li, G. Cui and Y. Li, Nano Energy, 2015, 11, 481–489.
79. H. Zhang, M. Zhao, X. He, Z. Wang, X. Zhang and X. Liu, J. Phys. Chem. C, 2011, 115, 8845–8850.
80. C. Sun and D. J. Searles, J. Phys. Chem. C, 2012, 116, 26222–26226.
81. R. Matsuoka, R. Sakamoto, K. Hoshiko, S. Sasaki, H. Masunaga, K. Nagashio and H. Nishihara, J. Am. Chem. Soc., 2017, 139, 3145–3152.
82. C. G. Broyden, IMA J. Appl. Math., 1970, 6, 76–90.
83. R. Fletcher, Comput. J., 1970, 13, 317–322.
84. D. Goldfarb, Math. Comput., 1970, 24, 23–26.
85. D. F. Shanno, Math. Comput., 1970, 24, 647–656.
86. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
87. J. Tersoff, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 9902–9905.
88. B. N. Dutta, Phys. Status Solidi B, 1962, 2, 984–987.
89. S. Nosé, Mol. Phys., 1984, 52, 255–268.
90. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
91. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
92. C. Kirsch and D. Sebastiani, Zenodo, 2025.
93. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
94. G. Mills and H. Jónsson, Phys. Rev. Lett., 1994, 72, 1124–1127.
95. G. Mills, H. Jónsson and G. K. Schenter, Surf. Sci., 1995, 324, 305–337.
96. B. J. Berne, G. Ciccotti and D. F. Coker, Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific, Singapore, 1st edn, 1998.
97. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
98. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.
99. M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, J. Chem. Phys., 2020, 152, 164105.
100. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
101. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
102. K. Persson, V. A. Sethuraman, L. J. Hardwick, Y. Hinuma, Y. S. Meng, A. van der Ven, V. Srinivasan, R. Kostecki and G. Ceder, J. Phys. Chem. Lett., 2010, 1, 1176–1180.
103. B. Cordero, V. Gómez, A. E. Platero-Prats, M. Revés, J. Echeverría, E. Cremades, F. Barragán and S. Alvarez, Dalton Trans., 2008, 2832–2838.

---