Developing Potential Energy Surfaces for Graphene-Based 2D–3D Interfaces From Modified High-Dimensional Neural Networks for Applications in Energy Storage

Intricately engineered device designs are put forth very often in recent years under the broad spectrum of advanced technology to target sustainability and enhanced efficiency. Complexity revolving around such device designs calls for the discovery of new multicomponent systems or new materials altogether. The concept of engineering materials has been applied in the vast sense in multicomponent systems to enhance functionality based on structure–function relationship [1]. For instance, in energy storage alone, the development of electrode design as core-shell, layered, or coating is continuously practiced for enhancing the electrode capacity and cycle life [2–6]. Such arrangements of different materials together in engineered devices create unique interfaces whose equilibrium and properties remain of intense research interest.

Two-dimensional (2D) material-based heterostructures (2D + nD, n = 0,1,2,3) have attracted enormous interest in various fields [7–9], including batteries [10,11]. The discovery of graphene in 2004 [12] opened a new space for mixed-dimensional materials where heterostructures based on graphene-like 2D materials layered with three-dimensional (3D) bulk materials have exhibited intriguing properties [7,13]. The mechanical, electrical, optical, magnetic, and thermal properties exhibited by these heterostructures have applications in batteries, optics, solar technologies, and electronics [14–22].

Arrangement of 2D–3D materials can be achieved either during epitaxial growth of one material over the other or with post-growth modification techniques [13]. Despite the development of multiple synthesis modes, a complete characterization of these complex nanoscale interfaces is yet a challenge. Scanning tunneling microscopy (STM) and atomic force microscopy (AFM) have enabled studying the interface between ordered 2D material heterostructures [23]. Yet, their scope is not expandable to polymorphing interfaces where 3D bulk undergoes lattice distortions and phase transitions to stabilize itself over a 2D substrate.

The fundamental understanding of interfacial properties in these systems is critical for sustaining the desired electro-chemo-mechanical behavior. Quantitative computational modeling efforts with ab-initio methods such as density functional theory (DFT) indicate the influence of structural polymorphism in 3D bulk on interfacial properties such as interfacial strength, electron transfer, and conductivity [24,25]. However, studies in this domain remain less explored due to the high computational cost involved in ab-initio methods. Moreover, molecular simulations of interfacial properties need reliable inter-atomic potential, limited to only a few elements. Due to these restrictions, there has been no significant advancement in interface studies at the atomic/molecular level. These limitations are expected to be overcome with the emergence of artificial intelligence (AI)-based models in materials.

Machine learning (ML)-based potential energy surface (PES) can describe complex systems at low computational cost and with close to first principles accuracy. These methods rely on a large amount of DFT data of materials (structures, energies, and forces) to efficiently explore the chemical space with respect to the target properties during training. Atomic structures of materials are represented by appropriate descriptors and fed to the neural network (NN) algorithm to generate PES that is invariant to translational, rotational, and permutation of homonuclear atoms [26,27]. Such PES are independent of any physical parameters and approximations, unlike empirical force fields, and hence are more accurate if descriptors describe the atomic local environments efficiently. Several ML techniques have been employed for PES development: linear regression [28,29], Gaussian approximation [30,31], high-dimensional neural networks (HDNN) [32], and graph neural networks (GNN) [33,34]. Of particular interest are HDNN using atom-centered symmetry functions (ACSF or SF) as input descriptors [26] that have been successful for a wide range of materials due to their generality [35–38]. Consequently, a lot of effort has been spent on refining the original HDNN given by Behler and Parrinello[32] to develop more efficient approaches [39–41].

The development of ML-based PES is still in its nascent stage and has barely explored materials’ domain beyond small organic molecules [42–44], single component condensed systems [32,45,46], bicomponent systems [47–49], and crystalline ordered motifs [50]. ML research in heterostructure territory is mostly limited to searching and predicting new ordered heterostructures with targeted properties [51–53]. The recent attempts to model 3D-3D interface grain boundaries demonstrate the scope of ML models in capturing mechanics at the nanoscale when trained on high accuracy data [54]. However, the primary bottleneck ML needs to scale in modeling complex materials is the tradeoff between quantity and quality of training data. As the complexity of the problem increases, acquiring ample data either computationally or empirically becomes a challenge.

The present work explores a mixed-dimensional heterostructure (2D–3D) interface formed by graphene (Gr) and tin (Sn) for the lattice distortions by utilizing DFT. Geometry optimization of Sn bulk over Gr results in the lattice rearrangements in interfacial Sn atomic layers. The extensive computation undertaken to study these interfaces with DFT has been further used to develop an ML-based PES for Gr and Sn interfaces. Here, the development of ML models could reinforce the atomic realization of these complex interfaces. This is one of the first attempts at modeling complex polymorphing 2D–3D interfaces with ML. We use a modified approach to the HDNN algorithm to develop PES for Gr|Sn interfaces which exhibit good transferability to new Gr|Sn interface systems despite limited training data. The design and utilization of ML-based potentials can be extended to expedite interface simulations in the future. Their scope of utilization in overcoming existing battery electrode design issues has been presented as a comprehensive discussion in the last section.

A large part of this assumption relies on the accuracy in describing atomic interactions in the local chemical neighborhoods. Appropriate descriptors that can convert cartesian coordinates to numeric vectors describing atomic interactions with translational, rotational, and permutational invariance are critical for the PES development. The descriptor defines the chemical environment of each atom within a specific cutoff radius (R_c). A large R_c value ensures that all energetically relevant interactions such as covalent, electrostatic, dispersion, and van der Waals interactions are included. Once the descriptor obtains fingerprints of atomic local space, separate feedforward neural networks are used for each atom in the system to express atomic energy contributions E_i. These E_i are then added to yield the total energy of the system E_total.

High-dimensional neural networks model in the present work (Fig. 1(b)) is a modified version of second-generation HDNN for bicomponent systems (also called BPNN for Behler and Parrinello Neural Networks, Fig. 1(a)). Figure 1 differentiates both HDNN models for a system of six atoms, three Sn and three C. In the BPNN model, chemical species are separated by a distinct set of weights. Atomic input features defined by descriptors (G_i) are fed to atomic neural networks (ann), and two sets of ann weights are fitted corresponding to each chemical species in the system. Weights corresponding to Sn atoms (set-a, red-ann) are identical, as are the weights corresponding to C atoms (set-b, yellow-ann). This ensures the invariance of total energy against the interchanging of two identical atoms within set-a and set-b. In addition, this permits easy size extrapolation of the model. If another atom is added to the system, additional ann corresponding to the species can be appended to the architecture and added to the total energy expression (represented in Eq. (1)). Complete details of second-generation HDNN (BPNN) for multicomponent systems can be found in Ref. [55]. In contrast, modified HDNN in Fig. 1(b) identifies chemical species with an added input feature rather than relying on separate trained weights. Thus, weights of all ann (set-k) are trained to be identical for the Sn–C system. All ann architectures and weights are suited for a bicomponent system of Sn–C rather than an individual species. This approach permits easy system size extrapolation.

Here, G_i² is a sum of Gaussians multiplied by the cutoff function. The width of the Gaussian and the center of the Gaussian can be defined by parameters η and R_s. A non-zero R_s value can shift the center of Gaussian away from the reference atom. Therefore, it is preferably set to 0. The parameter η is a Gaussian exponent responsible for indicating reduced interaction strength with increasing distance between the two atoms. Parameters ζ and λ in the function G_i⁴ control the angular resolution and cosine function, respectively. λ usually takes value either +1 or −1 for inverting the cosine function maxima from θ_ijk = 0 deg to θ_ijk = 180 deg [26]. The most preferred value for ζ is 1 as it provides sufficient coverage centered at 0 deg (when λ = 1). Higher values can increase angular resolution close to the center at a reduced coverage cost [57]. Multiple parameter set for each symmetry function (SF) type are needed to cover different portions of the chemical environment. Values of these parameters define the high-dimensional input vector representing the local environment of each atom in the material system. It is advisable that 100–150 G_i be used for the bicomponent system, such as the interface systems, to capture the full dimensionality of the system. Redundancy of the information that can arise due to a large number of G_i has been recognized to be not a problem for HDNN [26]. SF are uniquely beneficial as input vectors for HDNN because input vector size is independent of the actual number of neighboring atoms within a set cutoff radius R_c [36]. We used DScribe package in python to calculate SF for interface systems with the parameters defined in Table 1 [58].

We use a range of η values to capture the full dimensionality of the structures. The presented parameter set is chosen based on the benchmarking studies on descriptors for bicomponent bulk systems having elemental makeup similar to our structures [35,59,60]. λ has both values +1 and −1 corresponding to both centers in the SF set. To attain high angular resolution as well as complete coverage for intended interface systems, we use higher values for ζ in addition to 1 (ζ = 1, 2, 4). These SF set yielded the best results in our comparative evaluation and served as the foundation for further assessment of our model. As described in Sec. 2.1, we preferred to use large R_c for sufficient atomic interaction coverage. Here, R_c = 8.9 Å is found to be sufficient for optimum coverage of atom’s local environment (validation presented in supplemental information section I available in the Supplemental Materials on the ASME Digital Collection). By using the SF parameter set in Table 1, the chemical environment of each atom was represented by 162 input features. It is important to note here that the DScribe package used for the conversion of cartesian coordinates to SF considers the atomic number (Z) of chemical species in the environment of a central atom by appending the value to the SF (G² and G⁴). However, it does not consider an atomic number of central species in any way. To overcome this drawback, an additional input feature was added to the SF input feature for each atom which was its own atomic number. Thus, each atom was represented by 163 input features in our study.

Coherent interface models were created between ordered single layer Gr and Sn allotropes with an optimum interfacial gap of 3.5 Å. Gr contains 60 sp² hybridized carbon atoms in all the interface structures, whereas the size of Sn bulk varies from atom size of 16 to 64. Sn is known for having many energetically similar allotropes, with alpha (α-Sn) and beta (β-Sn) being the prominent ones. The interfaces are set in the periodic x–y plane across Sn (100) miller indices. These interface structures were modeled with a vacuum of 15 Å in z dimensions to circumvent the periodic influences, followed by DFT optimization to obtain relaxed strain-free interface configurations. All crystalline Sn bulks were derived from the materials project database [61]. Amorphous Sn was created with computational quenching of α-Sn₆₄ following the methodology discussed in our earlier works [24,62,63]. Both Gr and Sn structures were DFT optimized individually before interfacing. All DFT calculations are done using Vienna Ab initio Simulation Package (VASP) [64]. Projector-augmented-wave (PAW) potentials are used to mimic inert core electrons, while the valence electrons are represented by plane-wave basis set [65,66]. Plane-wave energy cutoff and convergence tolerance for all relaxations are 550 eV and 10⁻⁶ eV, respectively. The generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) exchange-correlation function is taken into account [67] with inclusivity of vdW correction to incorporate the effect of weak long-range van der Waals (vdW) forces [68]. The energy minimizations are done by conjugate gradient method with Hellmann-Feynman forces less than 0.02 eV/Å. Considering the vacuum slab structure of all interfaces, gamma-centered k-meshes 4 × 4 × 1 are used for good convergence.

Data for training and testing are systematically selected to meet the primary objective of the study, which is developing PES for Gr|Sn interfaces with the least possible computation necessary and best possible transferability. While material databases are a reliable source for most material’s structural data for learning-based workflows, they are significantly lacking in the domain of interfacial configurations. Tracing the equilibration of interfaces where exists the possibility of prominent lattice distortions, proved to be computationally expensive. It took approximately 700 h with 72 CPU cores to finish the complete relaxation of a single Gr|α-Sn interface system. To initiate 2D–3D interface-based structural analysis in the future, there is a need to develop machine learning-centered cheaper and faster workflows.

We optimized multiple unequilibrated Sn and Gr interface structures and divided them into training and testing datasets. The training dataset consists of five Gr|Sn interfaces: crystalline α-Sn_{(32 and 64)}, β-Sn_{(16 and 18)}, and amorphous Sn₆₄ interfaced with Gr (shown in supplemental information section II available in the Supplemental Materials on the ASME Digital Collection). The training dataset is built from convergence iterations of DFT simulations which covers the trajectory of minima search for five interfaces starting from the initial non-equilibrated structure. This scheme ensured that non-equilibrated and intermittent interface structures were as much part of the learning process as the relaxed structures. The intermediate DFT iterations of five Gr|Sn interface structures accounted for 5789 structures with their reference energies for training.

To analyze the performance of our model and test the transferability of PES in a sequential order of unfamiliarity, we use four carefully contemplated test interface structures (Fig. 2). The first test structure (T1) is the very Gr|β-Sn₁₆ interface used in the training dataset, except that the orientation of Sn bulk is slightly shifted over the Gr surface in T1. It is an example of a known interface structure and unknown orientation. The second test structure (T2) is again a familiar interface with increased Sn bulk size (Gr|β-Sn₁₆ → Gr|β-Sn₃₂). The objective of T2 is to test the system size extrapolation capabilities of the PES (see Fig. S2 available in the Supplemental Materials on the ASME Digital Collection).

In contrast to T1 and T2, the third test structure (T3) is the completely unfamiliar interface. A new Sn bulk (mp-949028 from Materials Project Database) is interfaced with Gr. The fourth and final test interface (T4) is derived from T3 by creating divacancy defects in the interfacing Gr. This change adds complexities of defects and surface adsorption in the interface structure, which are not noted in the earlier test interfaces. Differences in test structures from the training dataset are summarized in Table 2. The initial test interface structures were created like training interfaces and subjected to DFT optimizations. Since the current machine learning scheme does not include automated equilibration, the ability of PES to predict energies of intermediate configurations as structures search for global minima is assessed by predicting energies of test interfaces between initial and final configurations. The variations in the test interfaces during DFT optimization can be noted in the iteration snapshots presented in supplemental information section III available in the Supplemental Materials on the ASME Digital Collection. While minimal Sn alignment changes are observed in initial—final T2 and T3 structures, major structural transformative and surface defects are seen in T1 and T4, respectively.

This section discusses the atomic specifications of the Gr|Sn interface in optimized structures, which predominantly set apart these interfaces from Gr-based interfaces reported previously [7,13]. Discussion of interface phenomenon is important for designing interfaces and controlling heterostructure properties in applied technologies. The Gr|Sn interface systems are optimized by DFT to obtain relaxed strain-free interface configurations. Presented interfaces structurally resemble interfaces assembled in devices post-growth rather than interfaces originated during direct epitaxial growth of Gr on a substrate. Final interfacial configurations of Gr|Sn are depicted in Fig. 3 for two prominent Sn allotropes, α-Sn, and β-Sn, respectively. α-Sn is a diamond cubic crystal and β-Sn is a body-centered tetragonal crystal (Fig. 3(b)), which are two solid allotropes of Sn commonly in use. At temperatures below room temperature (286 K), α-Sn is the stable phase, which transitions to its β configurations rather quickly as temperature rises [69]. The sensitivity of temperature conditions symbolizes the significance of solid Sn α«b transitions for practical applications. This sensitivity elevates in the interfacial conditions with large lattice mismatch.

Differences in materials and lattice constants imply strained conditions in the buffer layer of Sn at the Gr|Sn interface, which conditions the Sn bulk towards a possible phase change. Consequently, the observed lattice constant of interfacial Sn (c = 4.5 Å) is different from the rest of the α-Sn bulk (c = 4.7 Å) in Fig. 3(a). This indicates a possible phase transformation from α-Sn to β-Sn in the buffer layer at Gr|α-Sn interface. However, these structural transformations of Sn are at a few layers limit at the surfaces and do not proliferate to central regions of the bulk where α conformations are retained (Fig. 4). Surface relaxation of independent α-Sn slab did not show any distortions, which eradicates the possibility of this structural reconstruction at Gr|α-Sn interface as mere surface defects. We repeated the simulation with increased vacuum in z dimensions to ensure structural distortions in the top layer are not due to periodic influences (shown in supplemental information section IV available in the Supplemental Materials on the ASME Digital Collection). Our observations indicate that this surface hardening of α-Sn is nucleated due to the presence of the Gr interface. In contrast to Gr|α-Sn interface, no structural distortions are noted in the relaxed Gr|β-Sn in Fig. 3(c), indicating the preference of the Gr interface toward β-Sn.

Structural changes in Sn at the Gr interface are also significantly impacted by Sn bulk size. Figures 3(d) and 3(e) exhibit Gr|Sn interfaces with smaller α-Sn and β-Sn bulks. Complete distortion of α-Sn₃₂ bulk is noted in Fig. 3(d) with an increased density (see supporting information section V available in the Supplemental Materials on the ASME Digital Collection). Likewise, β-Sn₁₆ rearranged over Gr surface as a single atomic layer in Fig. 3(e) with near-atomic distances of 3.15 Å. A drop in the Sn bulk size causes a reduction in the dimensions of Sn bulk and brings all the Sn atoms to the surface, much like 2D materials. These observations in small Sn crystals fall in line with prior experimental evidence of Sn nanocrystals becoming denser over graphene surfaces [70]. The relative differences in the material surfaces and charge analysis of Gr|Sn interfaces strongly indicate the weak van der Waals forces as the foundation of the formed interfaces. Charge analysis was performed in the said interfaces using Bader charge scripts by the Henkelman group [71]. The net electron exchange across the Gr|Sn interface is less than 1 e⁻¹ for all interfaces denoting negligible covalent interaction.

Sn is a well-known high-capacity anode for lithium ion batteries (LIB) and sodium ion batteries (NIB) with prominent shortcomings from phase transitions (β ↔ α) and volumetric strains. The presence of Gr substrate for Sn anode provenly scales down the volume expansion associated with mechanical failures [62,70]. Our analysis suggests it can possibly minimize frequented phase transitions from β ↔ α due to its preference for β-Sn. With our computational study, we attempt to closely understand the swift structural transformations in Gr|Sn interfaces. However, there is scope for further experimental X-ray diffraction analysis of Sn crystals interfaced with Gr, which can validate the presented observations.

HDNN is traditionally trained on reference energy per atom. However, electronic simulations do not provide actual energy per atom for reference, and energy per atom is deduced based on the total energy of the system [26]. Atomic energies are realized by dividing the total energy of the system by the number of atoms. This scheme gives equivalent atomic energies for all the atoms in a system. Such an approximation is suitable for a single component condensed system where atoms are present in a single phase. However, in Sec. 3.1, it can be observed that Gr|Sn interface structures have multiple phases with strained interfacial Sn atoms. Assuming interfacial Sn atoms will have higher atomic energies than sub-interfacial Sn atoms, training distinct atomic chemical environments in interface systems on uniform atomic energies is not considered a suitable approach.

The batch size for the training is kept at 10 (Size), and the accuracy of the energy prediction after each epoch is measured in terms of root-mean-square error (RMSE). Models are trained until accuracy metrics RMSE of at least 0.002 eV/Atom has been achieved. This amounted to 5000 epochs. The performance of the trained model is validated (validation split = 10%) and then tested on test structures T1 and T2. The performance of the trained model on a 10% validation split results in RMSE 0.0042 eV/Atom. This result is comparable to earlier reported deep learning studies on condensed phase systems [47,73], thereby concluding that this approach effectively develops PES for interfaces if target interfaces are similar to the training data. Next, trained PES are further used to predict atomic energies of the test structures T1 and T2. The results are summarized in Table 3. While the model performs well on validation split (test data randomly separated out of training data), its performance for new interfaces has been poor, indicating an overfitting case. A potential reason for this performance could be the training process where different atomic environments in a single system are assigned the same atomic energies. This renders the model inefficient in differentiating atomic energies of the atoms with different chemical neighborhoods.

The second training approach considers the total energy of the system as the final output of modified HDNN. In the last layer of HDNN, atomic energies E_i from all ann are added to yield E_total. This required fixing the atomic size (number of atoms in each system) in the training data to be 300. Input features for any system having less than 300 atoms have been extrapolated by zero padding. Inputs (G_i) to anns for each interface system are shuffled during each epoch. G_i is a one-dimensional array that encodes information about the central atom species and its local chemical environment. Hence, permutations at this stage do not change either the atomic energy or combined total energy. Loss is calculated from the total energies of the system. This allowed model to assign atomic energies based on the chemical space of each atom. Nested ann in modified HDNN has three hidden layers with 100–50–10 nodes. Hyperparameters of ann (nodes, activation function) remain the same as described before. Weights and biases are optimized through a supervised learning process using Adam optimizer and a learning rate of 0.00001. The loss function after each epoch was determined by the mean squared error of predicted total energies from reference DFT energies as per Eq. (5). The batch size for the training was kept one. The accuracy of the energy prediction after each epoch is measured in terms of RMSE, and the model is trained on 5789 Gr|Sn interface structures until the total energy RMSE of at least 0.2 eV is achieved.

Once trained, new PES is used to predict energies of test structures obtained from T1 (familiar interface, different orientation). Between unequilibrated and equilibrated T1 structures, there are approximately 260 structural configurations. All of which were used as test data. Figure 5 compares system energies of T1 structures obtained from DFT (E_DFT) with energies predicted by new PES (E_predict). Both system energies match closely with the RMSE value 0.0901 eV. Slight error is noted for non-equilibrated structures (below 50 DFT iterative structures), which further reduces as the structure stabilizes. Because HDNN is fitted to the total energies of the structure, we note that E_predict is bordering on E_DFT values but is not equivalent to the exact values even though the T1 interface is very close to trained structural configuration.

The primary objective of the targeted potentials (PES) is the ability to predict energies of new Gr|Sn interfaces and avail the least computation necessary during the development. The HDNN model trained on structures from 5 Gr|Sn interfaces has been tested on new interfaces T1–T4 described in Sec. 2.4. The results are summarized in Table 4. Between unequilibrated and equilibrated system configurations, there are approximately 260–400 structural configurations for each test structure used for testing. Predicted energies of new interfaces by modified HDNN weights have smaller RMSE values (eV/Atom). RMSE values for T1 and T2 in Table 4 are significantly lower than RMSE values noted in Table 3. This clearly indicates that the deep neural network model designed for such multiphasic interfaces should train across total energies rather than atomic energies to gather complete system information. The difference noted in energies of completely unfamiliar interfaces T3 and T4 is also relatively small. In the absence of accurate empirical potentials in literature, the developed PES demonstrates acceptable performance for new Gr|Sn interfaces constituting structural distortions.

Mathematical models like machine learning or deep learning can have powerful predictive accuracy when trained on large datasets. Data requirements are major barriers that delay the adoption of ML/DL techniques in artificial intelligence-assisted material development and discovery. The lack of a sizable dataset can be compensated with advanced sampling techniques such as heterogeneous dataset, random sampling, and stochastic surface walking [35,74]. Each of these methods concentrates on a different aspect of PES development. In this work, we use the geometry optimization arc of interfaces as a dataset, which includes unequilibrated, metastable, and equilibrated structures. The advisable course of using ab-initio molecular dynamics (AIMD) simulations to explore new equilibria for interfaces prior to DFTs is skipped due to the computational demand of the undertaking. The time to run a single AIMD simulation with 500 steps and 2 fs timestep on Gr|Sn interface ranged between 400 and –700 h on a CPU with 72 cores during our initial tests. Despite the limited dataset, our model predicted energies of new interface structures with lower RMSE when compared with the traditional approach. This can be said especially for the interface structures such as T1 and T2, where the Sn phase is familiar with reference training data. By avoiding incorrect decomposition of total energy into atomic energies, the model can obtain a high degree of accuracy and transferability even with a limited dataset.

To overcome the limitations of data, a heterogenous dataset was created consisting of individual Gr and Sn structures and distribution of Sn clusters adsorbed on the surface of Gr (2D–1D interface). This new heterogeneous dataset consisted of 9646 structures (see supporting information Fig. S6 available in the Supplemental Materials on the ASME Digital Collection). However, the modified HDNN model trained on heterogeneous datasets performed poorly for test interface structures (T1–T4). Since RMSE values from the heterogenous dataset approach were higher than the values in Table 4, these results have not been shown in the present work. The failure of the heterogenous dataset approach to capture 2D–3D interface structures accurately is primarily on account of the unique microstructural characteristics of 2D–3D interface from its individual 2D–1D interface counterparts. While the present study emphasizes a correct description of atomic energies in neural networks, there is further scope to develop successful data sampling approaches for 2D–3D interfaces.

The current state of existing ML models is still miles away from completely replacing DFT for complex interface systems. However, for the purpose of molecular dynamics simulations, modified HDNN models trained with appropriate atomic energy decomposition could be sufficient. Applications of ML-based PES depend on the dataset sampling approach adopted. While AIMD generated dataset trained model could be successfully used for simulations where more than one equilibrium is searched for, DFT-based dataset could be sufficient to explore one global minimum for the structures with a certain tradeoff between accuracy and computation times. Ongoing advancements in training and sampling techniques can possibly overcome the challenges associated with ML-assisted interface studies.

The conventional bulk materials-based batteries (Fig. 6(a)) have practical issues to meet the ever-increasing energy demand [75]. The two most critical chemo-mechanical failure modes in batteries are—(i) interface failure leading to electrical isolation of active electrode particles [76] and (ii) mechanical failure of active materials [77–79]. The active electrode particles (e.g., Si) have to be in contact with the metal current collector (Fig. 6(b1)), such as Ni, to ensure a uniform electron exchange [80]. During Li intercalation, active particles undergo substantial volume expansion [81]. For example, Si undergoes 300% volume expansion upon lithiation [82]. Since the metal current collectors (e.g., Cu and Ni) act as non-slippery surfaces, the volume expansion/contraction of active particles generates excessive interfacial stress (Fig. 6(b1)), leading to fracture of active particles, causing battery failures [62].

Another critical interface is between the polymer binder and active particles (Fig. 6(b2)) [83]. The polymer binder network keeps active particles adhered together and ensures continued electrical contact throughout the electrode [84]. However, volume expansion/contraction of active particles causes excessive interfacial stress at polymer-active particle interface, leading to detachment and electrical isolation of active particles (Fig. 6(b2)) [83]. Besides interface failure, another prominent reason for battery failure is the active particles’ fracture (Fig. 6(c)) [85]. Again volume expansion/contraction of Si anode particle upon lithiation/delithiation causes its formation of cracks [86], leading to battery failure [87]. These practical problems in battery electrodes need to be addressed by strategizing the electrode design.

Two-dimensional material-based heterostructures are promising candidates to solve these burning issues [88–90]. Electrical isolation of active particles can be avoided by replacing polymer binders with conductive and flexible 2D material such as MXenes, which can form an omnipresent electron-conducting network within the electrode [91] (Fig. 7(a1)). Next, the issue of the current collector and active particle interface failure (Fig. 7(b1)) can be addressed with two options: (i) Fig. 7(a2): addition of 2DM such as graphene as “coating” on the current collector to make it a “slippery” interface [24,62], (ii) Fig. 7(a3): completely replace the current collector with 2DM such as MXenes [92,93]. On the other hand, the problem of active materials failure (Fig. 6(c)) can be solved by using 2D materials-based anode (Fig. 7(b1)) [10] or 3D active materials integrated with 2D materials (Fig. 7(b2)) [7,94].

In all the 2D materials-based cases discussed, interface plays a critical role in electrochemical performance, i.e., energy and power density, volumetric capacity and so on.. The mechanical integrity of these interfaces dictates the long-term performance of energy storage systems [95–97]. Computational modeling methods such as DFT and MD simulations have been good alternatives to expensive experimental characterization for developing a deeper understanding of complex interfacial characteristics in batteries. Work by Basu et al.[62] on recognized benefits of slippery graphene surface at current collector end in combating stresses in Si anode upon lithiation is an example of the comprehensive scope of simulation studies. However, MD-based methodology to study interfacial stress cannot be extended to new and heavy metal electrodes such as Sn, Se, and more, due to the lack of appropriate forcefields. Presented work lays the foundation for developing the futuristic AI-based potentials to study a wide variety of 2D materials-based interfaces. Although the present study considers only graphene-based interfaces, the presented approach can be implemented in other 2D materials-based systems.

In summary, we performed optimization of different graphene and Sn-based 2D–3D interfaces resulting in a unique dataset, a kind lacking in existing databases. Our DFT simulation results show lattice distortions in Sn interfacing with graphene in great atomic details and highlight the preferred stability of β-Sn over graphene as opposed to α-Sn. One of the best application cases of Gr|Sn interface systems is in Sodium-ion batteries, where the presence of graphene interface can alleviate mechanical stresses upon Na intercalation in otherwise high capacity-low stability Sn anode. Usage of graphene-based heterostructures is undisputedly vast, and the need to model structural-functional aspects of such interfaces is an emergent need. We present the development of ML-based PES that can predict the energies of complex graphene-Sn 2D|3D interface systems with good accuracies that could be used to replace expensive ab-initio methods in future modeling efforts. Applicability of high-dimensional neural networks (HDNN) developed by Behler and Parrinello that utilize atom-centered symmetry functions as structural descriptors has been shown. The widely used approach to calculate loss function on atomic energies showed good performance on validation split but failed to predict energies of the new interface systems. To overcome this, we modify the HDNN model to enable training on the total energies of the system rather than atomic energies. This latter approach significantly improves the performance of PES in predicting the total energies of new Gr|Sn interface systems that constitute the test interfaces. The primary reason for this enhanced performance is the freedom model gained to assign atomic energies based on the atomic chemical environment. This opens the possibility for the more accurate evaluation of atomic energies and forces from ML models, allowing the scope for automated equilibration.

The work is supported by National Science Foundation (NSF) Civil, Mechanical and Manufacturing Innovation (CMMI) Program, Award Number # 1911900. The authors acknowledge the Extreme Science and Engineering Discovery Environment (XSEDE) for the computational facilities (Award Number—DMR180013). V.S. acknowledges Mr. Joy Datta for fruitful discussion.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.