BDIP-Net: Dual-Interaction Graph Learning for Property Prediction of Bilayer Materials

# BDIP-Net: Dual-Interaction Graph Learning for Property Prediction of Bilayer Materials
Source: [https://arxiv.org/html/2608.14640](https://arxiv.org/html/2608.14640)
,Chen ZhaoBaylor UniversityWacoUSA[chen\_zhao@baylor\.edu](mailto:[email protected]),Jin HuUniversity of ArkansasFayettevilleUSA[jinhu@uark\.edu](mailto:[email protected]),Shui\-Qing YuUniversity of ArkansasFayettevilleUSA[syu@uark\.edu](mailto:[email protected])andXintao WuUniversity of ArkansasFayettevilleUSA[xintaowu@uark\.edu](mailto:[email protected])

###### Abstract\.

Stacked bilayer materials exhibit rich stacking\-dependent properties driven by the interplay between strong intra\-layer bonding and weak inter\-layer van der Waals interactions\. The computational discovery of such materials is challenging because accurate structure generation typically relies on expensive DFT\-based optimization, while existing machine\-learning models often fail to explicitly distinguish different interaction types during property prediction\. To address these challenges, we propose a machine\-learning framework for efficient construction and property prediction of stacked bilayer materials\. The framework employs a MatterSim\-D3\-based structural optimization workflow to generate DFT\-quality bilayer structures from monolayer building blocks and stacking configurations at substantially reduced computational cost\. For property prediction, we introduce BDIP\-Net \(Bilayer Dual\-Interaction Potential Network\), a graph neural network that explicitly models intra\-layer and inter\-layer interactions through interaction\-specific potential representations and adaptive message fusion\. We evaluate the proposed framework on BiDB, HetDB, and SAMBA, encompassing homobilayers, heterobilayers, and twisted bilayer systems\. Results show that the MatterSim\-D3\-based workflow closely reproduces DFT\-PBE\-D3 optimized structures, while BDIP\-Net consistently outperforms existing graph neural network and potential\-based approaches for bilayer property prediction\.

††copyright:none††ccs:Computing methodologies Modeling and simulation††ccs:Computing methodologies Machine learning algorithms## 1\.Introduction

Stacked bilayer materials are formed by vertically assembling two monolayers\. Although the layers are primarily coupled through weak van der Waals interactions, their properties can vary significantly with stacking configuration, interlayer distance, lattice mismatch, and twist angle\. These stacking\-dependent effects can give rise to electronic and optical properties that are absent in the corresponding isolated monolayers\(Liet al\.,[2010](https://arxiv.org/html/2608.14640#bib.bib17)\)\. Density functional theory \(DFT\) with dispersion corrections, such as DFT\-D3\(Grimmeet al\.,[2010](https://arxiv.org/html/2608.14640#bib.bib6)\), is commonly used to optimize bilayer structures and evaluate their properties\. However, full DFT\-D3 structural optimization is computationally expensive because a single monolayer pair can generate numerous stacking configurations, lateral translations, and interlayer distances\. In existing bilayer material resources, the final bilayer structures are typically obtained through multi\-stage optimization workflows involving interlayer\-distance scans, in\-plane registry optimization, and atomic relaxation before DFT property calculations are performed\. The high computational cost of this process limits large\-scale screening of bilayer materials and the construction of datasets for machine learning \(ML\)\.

In this work, we develop a complete workflow for stacked bilayer material property prediction, spanning bilayer structure construction, efficient structural optimization, and ML\-based band gap prediction\. We adopt Machine learning interatomic potentials \(MLIPs\), which learns accurate representations of interatomic interactions while achieving computational efficiencies, and employ MatterSim\(Denget al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib5)\)with DFT\-D3 dispersion corrections\(Grimmeet al\.,[2010](https://arxiv.org/html/2608.14640#bib.bib6)\)as an MLIP\-D3 surrogate to replace costly DFT\-D3 structural optimization\. The resulting MatterSim\-D3 workflow produces optimized bilayer structures at a fraction of the computational cost of conventional DFT\-D3 optimization while maintaining comparable structural quality\. To further exploit the unique interaction characteristics of bilayer materials, we propose BDIP\-Net, an interaction\-aware graph neural network that extends the potential\-based modeling framework of PotNet\(Linet al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib30)\)to stacked bilayer systems\. BDIP\-Net constructs separate intra\-layer and inter\-layer edge sets, encodes them using Coulomb\-based and London\-dispersion\-based potentials, and integrates the resulting messages through an interaction\-message attention mechanism\. By explicitly distinguishing strong intra\-layer interactions from weak inter\-layer van der Waals interactions, BDIP\-Net learns interaction\-specific representations while preserving the physics\-informed advantages of potential\-based message passing\.

We evaluate our framework with three computational stacked 2D materials databases, BiDB\(Pakdelet al\.,[2024](https://arxiv.org/html/2608.14640#bib.bib19)\)which covers homobilayers by optimizing multiple stacking configurations, HetDB which\(Saueret al\.,[2025](https://arxiv.org/html/2608.14640#bib.bib4)\)covers van der Waals heterobilayers, and SAMBA\(Araújoet al\.,[2025](https://arxiv.org/html/2608.14640#bib.bib7)\)that provides over 18,000 twisted homo\- and hetero\-bilayer structures\. All three databases are built upon the Computational 2D Materials Database \(C2DB\)\(Haastrupet al\.,[2018](https://arxiv.org/html/2608.14640#bib.bib32)\)which is a comprehensive repository of first\-principles\-calculated monolayer materials\. Our experiments on these three benchmark bilayer datasets, BiDB, HetDB, and SAMBA, demonstrate that BDIP\-Net consistently outperforms BiMat\-ML, conventional bilayer\-based baselines, PotNet, and SE\-PotNet\. Furthermore, MatterSim\-D3 achieves nearly identical predictive performance to DFT\-PBE\-D3, validating the effectiveness of the proposed efficient optimization workflow\. Together, these results show that combining efficient MLIP\-based structural optimization with interaction\-aware potential\-based graph learning provides an effective and scalable framework for stacked bilayer material property prediction\.

## 2\.Related Work

Graph neural networks have become a widely adopted framework for materials property prediction by representing atoms as nodes and atomic interactions as edges\(Xie and Grossman,[2018](https://arxiv.org/html/2608.14640#bib.bib15); Schüttet al\.,[2017](https://arxiv.org/html/2608.14640#bib.bib3); Chenet al\.,[2019](https://arxiv.org/html/2608.14640#bib.bib10); Qiaoet al\.,[2020](https://arxiv.org/html/2608.14640#bib.bib2)\)\. Most crystal GNNs, e\.g\., SchNet\(Schüttet al\.,[2017](https://arxiv.org/html/2608.14640#bib.bib3)\), CGCNN\(Xie and Grossman,[2018](https://arxiv.org/html/2608.14640#bib.bib15)\), MEGNet\(Chenet al\.,[2019](https://arxiv.org/html/2608.14640#bib.bib10)\), and OrbNet\(Qiaoet al\.,[2020](https://arxiv.org/html/2608.14640#bib.bib2)\), encode interatomic distances using radial basis function \(RBF\) expansions and learn interaction patterns through message passing\. Instead of directly encoding distances, PotNet\(Linet al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib30)\)transforms interatomic distances into potential values and uses the resulting potential\-based features for message passing\. Its experiment results show physically informed interaction modeling can enhance crystal property prediction\. However, PotNet was developed for general crystalline materials and treats all local atomic interactions using a unified representation\. For stacked bilayer materials, the underlying interactions are inherently heterogeneous, consisting of strong intra\-layer bonding and weak inter\-layer van der Waals coupling\. Consequently, directly applying PotNet to bilayer materials may not fully capture the distinct physical characteristics of these interaction types\. Our proposed BDIP\-Net builds upon the potential\-based modeling philosophy of PotNet while explicitly distinguishing intra\-layer and inter\-layer interactions\. By constructing separate interaction graphs and employing interaction\-specific potential encodings together with interaction\-message attention, BDIP\-Net adapts physics\-informed message passing to the unique interaction mechanisms of stacked bilayer materials\.

Recent studies have explored machine learning approaches for predicting the properties of stacked two\-dimensional \(2D\) bilayer materials\. Among them, SE\-PAINN\(Chenet al\.,[2024](https://arxiv.org/html/2608.14640#bib.bib16)\)explicitly distinguishes intra\-layer and inter\-layer interactions by constructing separate edge sets and processing them through dedicated message\-passing branches\. However, SE\-PAINN relies on optimized bilayer structures generated through computationally intensive structural optimization workflows, limiting its scalability for large\-scale screening\. BiMat\-ML\(Vuonget al\.,[2026](https://arxiv.org/html/2608.14640#bib.bib20)\)adopts a different strategy by predicting bilayer properties directly from monolayer structures, stacking configurations, and monolayer properties\. By avoiding explicit bilayer structure construction and structural optimization, BiMat\-ML significantly reduces computational cost\. However, because the model operates on monolayer\-level representations rather than reconstructed bilayer structures, it cannot directly learn atomic interactions arising from the stacked bilayer geometry\. In contrast, our work combines efficient MLIP\-based structural optimization with interaction\-aware graph learning on reconstructed bilayer structures, enabling both scalable structure generation and explicit modeling of intra\-layer and inter\-layer interactions\.

## 3\.Method

![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/workflow3.jpg)Figure 1\.Workflow of the proposed bilayer material construction and property prediction framework\. The workflow consists of two main stages\. In the first stage, bilayer materials are constructed from monolayer structures obtained from C2DB\. Two types of bilayers are considered, homogeneous bilayers formed from two identical monolayers and heterogeneous bilayers formed from two different monolayers\. For each case, the initial bilayer construction uses a stacking configuration consisting of two3×33\\times 3stacking configuration matrices and a translation vector\. The bottom and top stacking configuration matrices are applied to the corresponding layers to produce transformed layers, and the translation vector is applied to the top layer while keeping the bottom layer fixed to form the initial bilayer structure\. The initial bilayer is then optimized using the MatterSim\-D3 structural optimization workflow, includingzz\-scan,x​yxy\-scan, and atomic relaxation, to obtain the structurally optimized bilayer\. In the second stage, the structurally optimized bilayer material is used as input to BDIP\-Net for graph\-level input preparation, message passing convolution layers, and property prediction\.The workflow of the proposed method is illustrated in Figure[1](https://arxiv.org/html/2608.14640#S3.F1)\. It consists of two main stages, including bilayer material construction and bilayer property prediction\.

### 3\.1\.MatterSim\-D3\-based Bilayer Structure Construction

In general, the initial bilayer construction for BiDB, HetDB, and SAMBA starts from DFT\-relaxed monolayer structures available in C2DB\. Before stacking, each layer is transformed using a dataset\-specific3×33\\times 3stacking configuration matrix\. This matrix defines the geometric operation applied to each layer before stacking, including enlarging the monolayer into a supercell, rotating the layer, and optionally flipping the layer along the out\-of\-plane direction\. After these matrix\-based transformations, a translation vector is applied to the top layer while keeping the bottom layer fixed\. This vector includes optional in\-plane translations along thexxandyydirections, and an out\-of\-plane translation along thezzdirection determined by the initial interlayer distance\. We provide the background of stacked 2D materials representation in Appendix[A](https://arxiv.org/html/2608.14640#A1)and the detailed construction processes for BiDB, HetDB, and SAMBA in Appendix[B](https://arxiv.org/html/2608.14640#A2)\.

In the original database workflows, the structural optimization of bilayer structures is performed using first\-principles density functional theory \(DFT\) calculations with explicit treatments of van der Waals \(vdW\) interactions\. All three workflows use DFT\-PBE, based on the Perdew–Burke–Ernzerhof \(PBE\) exchange\-correlation functional\(Perdewet al\.,[1996](https://arxiv.org/html/2608.14640#bib.bib41)\), but they employ different vdW treatments\. BiDB and HetDB use the DFT\-PBE\-D3 framework, where the D3 dispersion correction is added to the PBE total energy to account for long\-range vdW interactions\(Grimmeet al\.,[2010](https://arxiv.org/html/2608.14640#bib.bib6)\)\. The D3 term provides a computationally inexpensive correction for dispersion interactions\. In contrast, SAMBA uses DFT\-PBE with the optB86b vdW functional instead of the D3 correction, where optB86b is used to incorporate vdW interactions calculation\(Klimešet al\.,[2011](https://arxiv.org/html/2608.14640#bib.bib42)\)\. In these DFT\-based workflows, thezz\-scan determines the optimal interlayer distance by minimizing the total energy, thex​yxy\-scan evaluates a9×99\\times 9grid of in\-plane translations to identify the lowest energy stacking registry when applied, and the final atomic relaxation optimizes the structure until the maximum atomic force is below the dataset\-specific thresholdfmaxf\_\{\\max\}\. Although these procedures provide reliable optimized structures, each DFT\-PBE total energy or maximum atomic force calculation is computationally expensive\. Since the scan and relaxation steps require these calculations to be repeated many times for each candidate bilayer, the overall computational cost becomes substantial\.

To reduce this cost, we adopt the MLIP\-D3 structural optimization\(Mishin,[2021](https://arxiv.org/html/2608.14640#bib.bib1)\)\. Instead of repeatedly evaluating the DFT\-PBE total energy and maximum atomic force, we use an MLIP for the PBE energy and force components and add the DFT\-D3 dispersion correction to account for long\-range vdW interactions\. MLIP enables near\-first\-principles accuracy at computational costs that are orders of magnitude lower than DFT\. This choice is computationally efficient because DFT\-D3 has been classified as adding only a small additional computational cost, whereas vdW\-DF methods such as optB86b\-vdW increase the computational time by approximately50%50\\%compared with a regular DFT\-PBE calculation\(Klimešet al\.,[2011](https://arxiv.org/html/2608.14640#bib.bib42); Klimeš and Michaelides,[2012](https://arxiv.org/html/2608.14640#bib.bib43)\)\. Evaluation results in\(Saueret al\.,[2025](https://arxiv.org/html/2608.14640#bib.bib4)\)show that using MatterSim\(Denget al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib5)\)together with the DFT\-D3 correction, namely MatterSim\-D3, gives the best overall performance among the evaluated MLIP\-D3 models in reproducing bilayer structures optimized by DFT\-PBE\-D3, and DFT bandgap calculations based on bilayer structures optimized by MatterSim\-D3 differ from those based on DFT\-PBE\-D3 optimized structures by only about0\.02​eV0\.02~\\mathrm\{eV\}on average\. Therefore, for all bilayer structures, including those from SAMBA, we use MatterSim\-D3 to compute the total energy and maximum atomic force required in thezz\-scan, the optionalx​yxy\-scan, and the relaxation\. The detailed structural optimization process for each dataset is provided in Appendix[C](https://arxiv.org/html/2608.14640#A3)\. The resulting structurally optimized bilayer materialMbiM^\{\\mathrm\{bi\}\}is saved in the Vienna Ab initio Simulation Package \(VASP\) format\(Kresse and Furthmüller,[1996](https://arxiv.org/html/2608.14640#bib.bib40)\)and then used as the input structure for theBDIP\-Netmodel\.

### 3\.2\.BDIP\-Net For Bilayer Property Prediction

![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/bdipnet2.jpg)Figure 2\.Architecture of the BDIP\-Net model for stacked bilayer material property prediction\. The architecture consists of three main stages\. In the first stage, the structurally optimized bilayer material is used to prepare graph\-level inputs\. Atomic features are transformed into initial node embeddings, and adaptive layer neighbor search is applied to construct intra\-layer and inter\-layer edge sets\. Interaction\-specific potential\-based edge embeddings are then computed using a Coulomb\-based potential for intra\-layer edges and a London\-based potential for inter\-layer edges\. In the second stage, BDIP\-Net performs message passing overTTconvolution layers\. In each layer, intra\-layer and inter\-layer messages are constructed and aggregated separately, and interaction\-message attention is used to adaptively fuse the two interaction types when updating node representations\. In the third stage, the final node representations are mean\-pooled into a graph\-level representation, which is passed to an output MLP to predict the target property\. The model is trained end\-to\-end by minimizing the prediction loss\.Figure[2](https://arxiv.org/html/2608.14640#S3.F2)illusrates the overall training and prediction procedure of the proposed BDIP\-Net model for stacked bilayer material property prediction\. BDIP\-Net consists of two main components\. The first component prepares the graph\-level inputs for message passing by finding two interaction\-specific edge sets, namely intra\-layer interaction edges and inter\-layer interaction edges, initializing node embeddings, and computing interaction\-specific potential\-based edge embeddings\. The second component performs node message passing, where intra\-layer and inter\-layer messages are aggregated separately and then adaptively fused to update node representations for property prediction\.

#### 3\.2\.1\.Graph\-Level Input Preparation for Node Message Passing

Given a structurally optimized bilayer materialMbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), whereAbi∈ℝN×daA^\{\\mathrm\{bi\}\}\\in\\mathbb\{R\}^\{N\\times d\_\{a\}\}denotes the one\-hot encoded atomic features ofNNatoms,Pbi∈ℝN×3P^\{\\mathrm\{bi\}\}\\in\\mathbb\{R\}^\{N\\times 3\}denotes the atomic positions, andLbi∈ℝ3×3L^\{\\mathrm\{bi\}\}\\in\\mathbb\{R\}^\{3\\times 3\}denotes the lattice matrix, BDIP\-Net first constructs interaction\-specific edge sets for intra\-layer and inter\-layer interactions\. Adaptive layer neighbor search is then applied separately to each interaction type under periodic boundary conditions\. For each interaction type, the cutoff radius is adaptively increased until every atom has at leastNintraN^\{\\mathrm\{intra\}\}orNinterN^\{\\mathrm\{inter\}\}valid neighbors after the corresponding filtering\. This produces the intra\-layer edge setℰintra\\mathcal\{E\}^\{\\mathrm\{intra\}\}and the inter\-layer edge setℰinter\\mathcal\{E\}^\{\\mathrm\{inter\}\}, preserving the distinction between the two interaction types for subsequent edge embedding and node message passing\. Details of the adaptive layer neighbor search are provided in Appendix[D\.3](https://arxiv.org/html/2608.14640#A4.SS3)\.

After identifying the interaction\-specific edge sets, BDIP\-Net computes potential\-based edge embeddings\. This design is motivated by the PotNet ablation results\(Linet al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib30)\), which show that replacing raw distance\-based features with physically motivated local potential features accounts for most of the performance improvement\. Adding infinite potential\-based features provides only a modest additional improvement over local potential features\. BDIP\-Net therefore adopts only local potential features, retaining the main benefit of the potential\-based representation with a simpler design\. Specifically, BDIP\-Net follows PotNet by transforming the interatomic distancedi​jd\_\{ij\}of each intra\-layer edge\(i,j\)\(i,j\)using a Coulomb\-based potential\. Since inter\-layer interactions in stacked bilayer materials are typically dominated by weak van der Waals forces, they are commonly modeled using the Lennard–Jones potential\(Jones,[1924](https://arxiv.org/html/2608.14640#bib.bib39)\)\. BDIP\-Net instead adopts a London\-dispersion\-based potential\(Wagner and Schreiner,[2015](https://arxiv.org/html/2608.14640#bib.bib38)\)as a simpler formulation for capturing the dominant attractived−6d^\{\-6\}component of these interactions\. The scalar potentials areVCoulomb​\(i,j\)=−ϵCoulomb/di​jV\_\{\\mathrm\{Coulomb\}\}\(i,j\)=\-\\epsilon^\{\\mathrm\{Coulomb\}\}/d\_\{ij\}andVLondon​\(i,j\)=−ϵLondon/\(di​j\)6V\_\{\\mathrm\{London\}\}\(i,j\)=\-\\epsilon^\{\\mathrm\{London\}\}/\(d\_\{ij\}\)^\{6\}, whereϵCoulomb\\epsilon^\{\\mathrm\{Coulomb\}\}andϵLondon\\epsilon^\{\\mathrm\{London\}\}are scaling coefficients\.

\(1\)ei​jintra\\displaystyle e\_\{ij\}^\{\\mathrm\{intra\}\}=MLPeintra​\(RBFCoulomb​\(−ϵCoulomb/di​j\)\)∈ℝdh,\\displaystyle=\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{intra\}\}\\\!\\left\(\\mathrm\{RBF\}^\{\\mathrm\{Coulomb\}\}\\\!\\left\(\-\\epsilon^\{\\mathrm\{Coulomb\}\}/d\_\{ij\}\\right\)\\right\)\\in\\mathbb\{R\}^\{d\_\{h\}\},ei​jinter\\displaystyle e\_\{ij\}^\{\\mathrm\{inter\}\}=MLPeinter​\(RBFLondon​\(−ϵLondon/\(di​j\)6\)\)∈ℝdh,\\displaystyle=\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{inter\}\}\\\!\\left\(\\mathrm\{RBF\}^\{\\mathrm\{London\}\}\\\!\\left\(\-\\epsilon^\{\\mathrm\{London\}\}/\(d\_\{ij\}\)^\{6\}\\right\)\\right\)\\in\\mathbb\{R\}^\{d\_\{h\}\},
whereRBFCoulomb\\mathrm\{RBF\}^\{\\mathrm\{Coulomb\}\}andRBFLondon\\mathrm\{RBF\}^\{\\mathrm\{London\}\}denote separate RBF expansions, whileMLPeintra\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{intra\}\}andMLPeinter\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{inter\}\}are interaction\-specific edge networks\. Details of the RBF expansions are provided in Appendix[D\.2](https://arxiv.org/html/2608.14640#A4.SS2)\. Next, the atomic node embeddings are initialized as𝒱\(0\)=Linear𝒱​\(Abi\)∈ℝN×dh\\mathcal\{V\}^\{\(0\)\}=\\mathrm\{Linear\}\_\{\\mathcal\{V\}\}\(A^\{\\mathrm\{bi\}\}\)\\in\\mathbb\{R\}^\{N\\times d\_\{h\}\}, whereLinear𝒱\\mathrm\{Linear\}\_\{\\mathcal\{V\}\}denotes a linear projection\.

#### 3\.2\.2\.Message Passing and Property Prediction

This subsection describes the second component of BDIP\-Net, which updates node representations through message passing and then predicts the target bilayer material property\. Specifically, BDIP\-Net appliesTTmessage\-passing convolution layers to update the initial node embeddings obtained from Section[3\.2\.1](https://arxiv.org/html/2608.14640#S3.SS2.SSS1)\. At each layer, the representation of each node is updated by aggregating information from two interaction\-specific branches, namely an intra\-layer branch and an inter\-layer branch\. The intra\-layer branch aggregates edge messages from neighboring nodes within the same layer, while the inter\-layer branch aggregates edge messages from neighboring nodes across different layers\. Within each branch, each node receives edge messages from its neighboring nodes, and BDIP\-Net computes a scalar attention weight for each edge message to measure its importance during interaction\-message aggregation\. The weighted edge messages in each branch are then aggregated to produce one intra\-layer interaction\-message and one inter\-layer interaction message for each node\. These two interaction\-messages are further adaptively fused by interaction\-message attention, which computes scalar attention weights to learn the relative contribution of the intra\-layer and inter\-layer interaction messages when updating each node representation\. The use of scalar attention weights for both edge messages and interaction\-specific messages is conceptually inspired by the slot\-based attention idea in SlotGAT\(Zhouet al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib21)\), where separate representation slots are maintained and their relative importance is learned by attention\. However, unlike SlotGAT, where slots correspond to node types in heterogeneous graphs, BDIP\-Net applies scalar attention weights to edge messages within each interaction branch and to interaction\-specific messages across the intra\-layer and inter\-layer branches\. For each interaction typeη∈\{intra,inter\}\\eta\\in\\\{\\mathrm\{intra\},\\mathrm\{inter\}\\\}and each edge\(i,j\)∈ℰη\(i,j\)\\in\\mathcal\{E\}^\{\\eta\}, BDIP\-Net constructs the pairwise interaction representation as

\(2\)wi​j\(t\),η=vi\(t\)⊕vj\(t\)⊕ei​jη∈ℝ3​dh,w\_\{ij\}^\{\(t\),\\eta\}=v\_\{i\}^\{\(t\)\}\\oplus v\_\{j\}^\{\(t\)\}\\oplus e\_\{ij\}^\{\\eta\}\\in\\mathbb\{R\}^\{3d\_\{h\}\},where⊕\\oplusdenotes feature concatenation andei​jηe\_\{ij\}^\{\\eta\}denotes the edge feature associated with interaction typeη\\eta\. The pairwise interaction representation is then transformed into an interaction\-specific edge message:

\(3\)zi​j\(t\),η=MLPz\(t\),η​\(wi​j\(t\),η\)∈ℝdh,z\_\{ij\}^\{\(t\),\\eta\}=\\mathrm\{MLP\}\_\{z\}^\{\(t\),\\eta\}\\left\(w\_\{ij\}^\{\(t\),\\eta\}\\right\)\\in\\mathbb\{R\}^\{d\_\{h\}\},whereMLPz\(t\),η\\mathrm\{MLP\}\_\{z\}^\{\(t\),\\eta\}maps the3​dh3d\_\{h\}\-dimensional pairwise interaction representation to adhd\_\{h\}\-dimensional edge message for interaction typeη\\eta\. The attention weight of each edge is computed from the transformed edge message:

\(4\)αi​j\(t\),η=softmaxj∈𝒩\(i\)η⁡\(MLPα\(t\),η​\(zi​j\(t\),η\)\)∈ℝ,\\alpha\_\{ij\}^\{\(t\),\\eta\}=\\operatorname\{softmax\}\_\{j\\in\\mathcal\{N\}^\{\\eta\}\_\{\(i\)\}\}\\left\(\\mathrm\{MLP\}\_\{\\alpha\}^\{\(t\),\\eta\}\\left\(z\_\{ij\}^\{\(t\),\\eta\}\\right\)\\right\)\\in\\mathbb\{R\},whereMLPα\(t\),η\\mathrm\{MLP\}\_\{\\alpha\}^\{\(t\),\\eta\}maps eachdhd\_\{h\}\-dimensional edge message to an attention score, and the softmax is normalized over the neighbors of nodeiiassociated with interaction typeη\\eta\. The interaction\-specific message of nodeiiis then aggregated as

\(5\)mi\(t\),η=∑j∈𝒩\(i\)ηαi​j\(t\),η​zi​j\(t\),η∈ℝdh\.m\_\{i\}^\{\(t\),\\eta\}=\\sum\_\{j\\in\\mathcal\{N\}^\{\\eta\}\_\{\(i\)\}\}\\alpha\_\{ij\}^\{\(t\),\\eta\}z\_\{ij\}^\{\(t\),\\eta\}\\in\\mathbb\{R\}^\{d\_\{h\}\}\.
After the intra\-layer and inter\-layer messages are obtained, BDIP\-Net applies interaction\-message attention fusion to adaptively combine them for each atom\. The interaction\-message attention weights are computed as

\(6\)\[βi\(t\),intraβi\(t\),inter\]=softmax⁡\(MLPβ\(t\)​\(\[mi\(t\),intrami\(t\),inter\]\)\)∈ℝ2,\\begin\{bmatrix\}\\beta\_\{i\}^\{\(t\),\\mathrm\{intra\}\}\\\\ \\beta\_\{i\}^\{\(t\),\\mathrm\{inter\}\}\\end\{bmatrix\}=\\operatorname\{softmax\}\\left\(\\mathrm\{MLP\}\_\{\\beta\}^\{\(t\)\}\\left\(\\begin\{bmatrix\}m\_\{i\}^\{\(t\),\\mathrm\{intra\}\}\\\\ m\_\{i\}^\{\(t\),\\mathrm\{inter\}\}\\end\{bmatrix\}\\right\)\\right\)\\in\\mathbb\{R\}^\{2\},whereMLPβ\(t\)\\mathrm\{MLP\}\_\{\\beta\}^\{\(t\)\}is shared across the intra\-layer and inter\-layer interaction\-messages at layertt\. The MLP maps the twodhd\_\{h\}\-dimensional interaction\-messages to two attention scores, and the softmax normalizes the scores over the two message slots\. The node representation is then updated by adding the weighted interaction messages to the previous\-layer representation:

\(7\)vi\(t\+1\)=vi\(t\)\+βi\(t\),intra​mi\(t\),intra\+βi\(t\),inter​mi\(t\),inter∈ℝdh\.v\_\{i\}^\{\(t\+1\)\}=v\_\{i\}^\{\(t\)\}\+\\beta\_\{i\}^\{\(t\),\\mathrm\{intra\}\}m\_\{i\}^\{\(t\),\\mathrm\{intra\}\}\+\\beta\_\{i\}^\{\(t\),\\mathrm\{inter\}\}m\_\{i\}^\{\(t\),\\mathrm\{inter\}\}\\in\\mathbb\{R\}^\{d\_\{h\}\}\.
AfterTTmessage\-passing layers, the final node representations\{vi\(T\)\}i∈𝒱\\\{v\_\{i\}^\{\(T\)\}\\\}\_\{i\\in\\mathcal\{V\}\}are averaged to obtain the graph\-level representation:

\(8\)Z𝒢=1\|𝒱\|​∑i∈𝒱vi\(T\)∈ℝdh\.Z\_\{\\mathcal\{G\}\}=\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{i\\in\\mathcal\{V\}\}v\_\{i\}^\{\(T\)\}\\in\\mathbb\{R\}^\{d\_\{h\}\}\.
The graph\-level representationZ𝒢Z\_\{\\mathcal\{G\}\}is passed through an output MLP to predict the target bilayer property asY^=MLPY^​\(Z𝒢\)\\hat\{Y\}=\\mathrm\{MLP\}\_\{\\hat\{Y\}\}\(Z\_\{\\mathcal\{G\}\}\)\. The entire framework is trained end\-to\-end by minimizing the mean absolute error lossℒ=\|Y^−Y\|\\mathcal\{L\}=\|\\hat\{Y\}\-Y\|\. The outputMLPY^\\mathrm\{MLP\}\_\{\\hat\{Y\}\}uses a Linear–ShiftedSoftplus–Linear network architecture, while all other MLPs in BDIP\-Net use a Linear–SiLU–Linear network architecture, where ShiftedSoftplus and SiLU denote the activation functions between the two linear layers\. MLPs with different subscripts or superscripts, including those for different interaction types and different convolution layers, are independently parameterized\. The pseudocode of our BDIP\-Net is shown in Algorithm[8](https://arxiv.org/html/2608.14640#alg8)of Appendix[D\.1](https://arxiv.org/html/2608.14640#A4.SS1)\.

## 4\.Experiments

### 4\.1\.Experiment Settings

![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/band_gap_box_plot.png)Figure 3\.Bandgap distribution in BiDB, HetDB, and SAMBA\.Datasets\.We follow Bimat\-ML\(Vuonget al\.,[2026](https://arxiv.org/html/2608.14640#bib.bib20)\)and use the same processed BiDB and HetDB datasets\. The SAMBA dataset contains more than 18,000 twisted bilayer structures generated from 63 monolayers in the C2DB database\. We use 980 bilayers \(144 homobilayers and 836 heterobilayers\) with available DFT\-optimized structures and bandgap values\. For all three datasets, we use bandgap as the target property\. Fig\.[3](https://arxiv.org/html/2608.14640#S4.F3)shows the composition and bandgap distributions of the three datasets\. BiDB contains 6,683 homobilayers with a mean bandgap of0\.930\.93eV, while HetDB contains 336 heterobilayers with a mean bandgap of0\.760\.76eV\. SAMBA contains 146 homobilayers and 834 heterobilayers, with mean bandgaps of0\.660\.66and0\.280\.28eV, respectively\. The three datasets provide complementary coverage of aligned homobilayers, heterobilayers, and twisted bilayer systems\. The differences in dataset size, bilayer type, bandgap distribution, and bilayer\-monolayer bandgap relationships form an ideal benchmark for evaluation\.

Cross\-validation strategy\.We perform 4\-fold cross\-validation by splitting each dataset at the monolayer\-pair level rather than the bilayer level\. All bilayers constructed from the same monolayer pair are assigned to the same fold\. For homogeneous bilayers, the pair is represented as\(M,M\)\(M,M\), whereas for heterogeneous bilayers it is represented as\(MA,MB\)\(M\_\{A\},M\_\{B\}\)\. This strategy prevents the same monolayer pair from appearing in both the training and test sets, thereby evaluating the models’ ability to generalize to bilayers constructed from unseen monolayer pairs\.

BDIP\-Net\.Our BDIP\-Net uses the MatterSim\-D3 structural optimization procedure described in Section[3\.1](https://arxiv.org/html/2608.14640#S3.SS1)\. We also report the results from DFT\-PBE\-D3 where the optimized bilayer structures provided by BiDB, HetDB, and SAMBA are directly converted into the VASP format and used as model inputs\. For BiDB, HetDB, and SAMBA, we use a consistent BDIP\-Net hyperparameter configuration, with only the batch size adjusted for HetDB\. The model consists of three convolutional layers \(T=3T=3\) and uses adaptive layer neighbor search withNintra=16N^\{\\mathrm\{intra\}\}=16andNinter=16N^\{\\mathrm\{inter\}\}=16, starting from initial search radii ofRintra=Rinter=4R^\{\\mathrm\{intra\}\}=R^\{\\mathrm\{inter\}\}=4\. Both potential scaling factors are set toϵCoulomb=ϵLondon=1\\epsilon^\{\\mathrm\{Coulomb\}\}=\\epsilon^\{\\mathrm\{London\}\}=1, and the hidden dimension is set todh=128d\_\{h\}=128for both node and edge embeddings\. The model is trained for 500 epochs using the Adam optimizer with a learning rate of2×10−42\\times 10^\{\-4\}\. The batch size is set to 128 for BiDB and SAMBA and 16 for HetDB\.

Baselines\.We compare with Bimat\-ML\(Vuonget al\.,[2026](https://arxiv.org/html/2608.14640#bib.bib20)\), SE\-CGCNN\(Chenet al\.,[2024](https://arxiv.org/html/2608.14640#bib.bib16)\), SE\-MEGNET\(Chenet al\.,[2024](https://arxiv.org/html/2608.14640#bib.bib16)\), SE\-PAINN\(Chenet al\.,[2024](https://arxiv.org/html/2608.14640#bib.bib16)\)and PotNet\(Linet al\.,[2023](https://arxiv.org/html/2608.14640#bib.bib30)\)\. BDIP\-Net and the bilayer\-based baselines are evaluated under the two structural optimization settings described above, as indicated in the Optimization column of the result tables\. In our experiments, Bimat\-ML uses only monolayer structures and stacking configuration information as inputs, without using monolayer properties, ensuring a fair comparison with the other models\. Since Bimat\-ML does not take an optimized bilayer structure as input, its entry in the Optimization column is marked as N/A\. All experiments in this work are conducted on an NVIDIA H200 GPU with 140 GB of memory\.

### 4\.2\.Experiment Results

Tables[1](https://arxiv.org/html/2608.14640#S4.T1)summarizes the single\-domain bandgap prediction results on BiDB, HetDB, and SAMBA, respectively\. Each table reports the MAE, MSE, RMSE, andR2R^\{2\}of BDIP\-Net and the baseline models using bilayer structures optimized by DFT\-PBE\-D3 and MatterSim\-D3\. Bimat\-ML, which does not use optimized bilayer structures or monolayer property inputs, is included as a structure\-free baseline\.

Table 1\.Regression performance on the BiDB, HetDB, and SAMBA datasets for bandgap prediction\.DatasetModelOptimizationMAE↓\\downarrowMSE↓\\downarrowRMSE↓\\downarrowR2↑\\uparrowBiDBBimat\-MLN/A0\.35±0\.050\.35\\pm 0\.050\.35±0\.080\.35\\pm 0\.080\.59±0\.060\.59\\pm 0\.060\.68±0\.060\.68\\pm 0\.06SE\-CGCNNDFT\-PBE\-D30\.38±0\.050\.38\\pm 0\.050\.49±0\.190\.49\\pm 0\.190\.69±0\.090\.69\\pm 0\.090\.68±0\.060\.68\\pm 0\.06MatterSim\-D30\.38±0\.050\.38\\pm 0\.050\.50±0\.180\.50\\pm 0\.180\.70±0\.130\.70\\pm 0\.130\.66±0\.080\.66\\pm 0\.08SE\-MEGNETDFT\-PBE\-D30\.36±0\.050\.36\\pm 0\.050\.37±0\.130\.37\\pm 0\.130\.60±0\.110\.60\\pm 0\.110\.67±0\.060\.67\\pm 0\.06MatterSim\-D30\.36±0\.060\.36\\pm 0\.060\.37±0\.120\.37\\pm 0\.120\.60±0\.100\.60\\pm 0\.100\.67±0\.080\.67\\pm 0\.08SE\-PAINNDFT\-PBE\-D30\.36±0\.040\.36\\pm 0\.040\.40±0\.080\.40\\pm 0\.080\.63±0\.070\.63\\pm 0\.070\.64±0\.060\.64\\pm 0\.06MatterSim\-D30\.36±0\.060\.36\\pm 0\.060\.40±0\.070\.40\\pm 0\.070\.63±0\.060\.63\\pm 0\.060\.63±0\.070\.63\\pm 0\.07PotNetDFT\-PBE\-D30\.33±0\.020\.33\\pm 0\.020\.36±0\.060\.36\\pm 0\.060\.60±0\.060\.60\\pm 0\.060\.66±0\.030\.66\\pm 0\.03MatterSim\-D30\.33±0\.030\.33\\pm 0\.030\.36±0\.070\.36\\pm 0\.070\.60±0\.060\.60\\pm 0\.060\.65±0\.040\.65\\pm 0\.04BDIP\-NetDFT\-PBE\-D30\.30±0\.04\{0\.30\\pm 0\.04\}0\.32±0\.10\{0\.32\\pm 0\.10\}0\.57±0\.100\.57\\pm 0\.100\.72±0\.050\.72\\pm 0\.05MatterSim\-D30\.30±0\.03\\mathbf\{0\.30\\pm 0\.03\}0\.32±0\.10\\mathbf\{0\.32\\pm 0\.10\}0\.56±0\.09\\mathbf\{0\.56\\pm 0\.09\}0\.73±0\.05\\mathbf\{0\.73\\pm 0\.05\}HetDBBimat\-MLN/A0\.16±0\.010\.16\\pm 0\.010\.09±0\.040\.09\\pm 0\.040\.29±0\.060\.29\\pm 0\.060\.76±0\.120\.76\\pm 0\.12SE\-CGCNNDFT\-PBE\-D30\.12±0\.020\.12\\pm 0\.020\.03±0\.020\.03\\pm 0\.020\.18±0\.020\.18\\pm 0\.020\.91±0\.020\.91\\pm 0\.02MatterSim\-D30\.12±0\.010\.12\\pm 0\.010\.03±0\.020\.03\\pm 0\.020\.17±0\.040\.17\\pm 0\.040\.92±0\.050\.92\\pm 0\.05SE\-MEGNETDFT\-PBE\-D30\.19±0\.020\.19\\pm 0\.020\.09±0\.030\.09\\pm 0\.030\.29±0\.060\.29\\pm 0\.060\.76±0\.080\.76\\pm 0\.08MatterSim\-D30\.19±0\.020\.19\\pm 0\.020\.09±0\.040\.09\\pm 0\.040\.30±0\.060\.30\\pm 0\.060\.75±0\.070\.75\\pm 0\.07SE\-PAINNDFT\-PBE\-D30\.11±0\.020\.11\\pm 0\.020\.03±0\.020\.03\\pm 0\.020\.17±0\.020\.17\\pm 0\.020\.91±0\.020\.91\\pm 0\.02MatterSim\-D30\.11±0\.010\.11\\pm 0\.010\.03±0\.020\.03\\pm 0\.020\.17±0\.050\.17\\pm 0\.050\.90±0\.030\.90\\pm 0\.03PotNetDFT\-PBE\-D30\.11±0\.020\.11\\pm 0\.020\.04±0\.030\.04\\pm 0\.030\.21±0\.060\.21\\pm 0\.060\.89±0\.070\.89\\pm 0\.07MatterSim\-D30\.11±0\.020\.11\\pm 0\.020\.05±0\.020\.05\\pm 0\.020\.21±0\.050\.21\\pm 0\.050\.89±0\.050\.89\\pm 0\.05BDIP\-NetDFT\-PBE\-D30\.09±0\.02\{0\.09\}\\pm 0\.020\.02±0\.02\{0\.02\}\\pm 0\.020\.14±0\.05\{0\.14\}\\pm 0\.050\.94±0\.05\{0\.94\}\\pm 0\.05MatterSim\-D30\.09±0\.01\\mathbf\{0\.09\\pm 0\.01\}0\.02±0\.02\\mathbf\{0\.02\\pm 0\.02\}0\.15±0\.06\\mathbf\{0\.15\\pm 0\.06\}0\.93±0\.05\\mathbf\{0\.93\\pm 0\.05\}SAMBABimat\-MLN/A0\.18±0\.030\.18\\pm 0\.030\.13±0\.090\.13\\pm 0\.090\.35±0\.120\.35\\pm 0\.120\.66±0\.120\.66\\pm 0\.12SE\-CGCNNDFT\-PBE\-D30\.16±0\.020\.16\\pm 0\.020\.09±0\.080\.09\\pm 0\.080\.29±0\.110\.29\\pm 0\.110\.71±0\.110\.71\\pm 0\.11MatterSim\-D30\.16±0\.020\.16\\pm 0\.020\.09±0\.060\.09\\pm 0\.060\.29±0\.100\.29\\pm 0\.100\.70±0\.090\.70\\pm 0\.09SE\-MEGNETDFT\-PBE\-D30\.15±0\.020\.15\\pm 0\.020\.08±0\.050\.08\\pm 0\.050\.28±0\.080\.28\\pm 0\.080\.73±0\.040\.73\\pm 0\.04MatterSim\-D30\.15±0\.020\.15\\pm 0\.020\.08±0\.050\.08\\pm 0\.050\.28±0\.080\.28\\pm 0\.080\.73±0\.060\.73\\pm 0\.06SE\-PAINNDFT\-PBE\-D30\.15±0\.020\.15\\pm 0\.020\.09±0\.070\.09\\pm 0\.070\.27±0\.100\.27\\pm 0\.100\.73±0\.090\.73\\pm 0\.09MatterSim\-D30\.15±0\.020\.15\\pm 0\.020\.09±0\.060\.09\\pm 0\.060\.28±0\.110\.28\\pm 0\.110\.73±0\.090\.73\\pm 0\.09PotNetDFT\-PBE\-D30\.13±0\.020\.13\\pm 0\.020\.07±0\.090\.07\\pm 0\.090\.26±0\.100\.26\\pm 0\.100\.74±0\.090\.74\\pm 0\.09MatterSim\-D30\.13±0\.020\.13\\pm 0\.020\.09±0\.060\.09\\pm 0\.060\.26±0\.100\.26\\pm 0\.100\.73±0\.080\.73\\pm 0\.08BDIP\-NetDFT\-PBE\-D30\.10±0\.02\{0\.10\}\\pm 0\.020\.05±0\.040\.05\\pm 0\.040\.21±0\.070\.21\\pm 0\.070\.83±0\.03\{0\.83\}\\pm 0\.03MatterSim\-D30\.10±0\.01\\mathbf\{0\.10\\pm 0\.01\}0\.04±0\.02\\mathbf\{0\.04\\pm 0\.02\}0\.20±0\.05\\mathbf\{0\.20\\pm 0\.05\}0\.83±0\.03\\mathbf\{0\.83\\pm 0\.03\}

Table 2\.P\-values of pairwise comparisons against BDIP\-Net over combined 3 datasets\.ModelMAEMSERMSER2Bimat\-ML0\.0000\.0000\.0000\.000SE\-CGCNN0\.0000\.0060\.0000\.001SE\-MEGNET0\.0000\.0000\.0000\.000SE\-PAINN0\.0000\.0010\.0010\.000PotNet0\.0000\.0010\.0000\.000BDIP\-Net\-0\.0020\.0020\.0020\.001#### 4\.2\.1\.Prediction Performance

As shown in Table[1](https://arxiv.org/html/2608.14640#S4.T1), our BDIP\-Net achieves the best performance across all four metrics \(MAE, MSE, RMSE, andR2R^\{2\}\) and and all three datasets \(BiDB, HetDB, and SAMBA\)\. Specifically, compared with the best\-performed baseline PotNet, BDIP\-Net increases theR2R^\{2\}from0\.650\.65to0\.730\.73on BiDB, from0\.890\.89to0\.930\.93on HetDB, from0\.730\.73to0\.830\.83on SAMBA\.

Statistical significance\.We perform one\-sided pairedtt\-tests to compare BDIP\-Net with each baseline model for MAE, MSE, RMSE, andR2R^\{2\}, with the alternative hypothesis that BDIP\-Net achieves lower MAE, MSE, and RMSE and a higherR2R^\{2\}\. For each comparison, the four\-fold results from BiDB, HetDB, and SAMBA are combined to form 12 paired observations\. All comparisons between BDIP\-Net and the baseline models yieldp<0\.01p<0\.01, as shown in Table[2](https://arxiv.org/html/2608.14640#S4.T2), indicating that the improvements of BDIP\-Net are statistically significant\.

#### 4\.2\.2\.Structural Optimization and Its Effect on Prediction

Table[1](https://arxiv.org/html/2608.14640#S4.T1)also reports the bandgap prediction results obtained using DFT\-PBE\-D3\-optimized structure\. Both BDIP\-Net and baseline models exhibit only minor variations when DFT\-PBE\-D3\-optimized structures are replaced with MatterSim\-D3\-optimized structures\. Overall, the small differences between the two optimization settings indicate that MatterSim\-D3 preserves the structural information required for bandgap prediction and achieves performance comparable to that obtained using DFT\-PBE\-D3\-optimized structures while significantly reducing the computational cost\.

#### 4\.2\.3\.Ablation Study

We conduct an ablation study on the message\-passing operation and interaction\-message attention module of BDIP\-Net by comparing the original model with an ablated variant, denoted as BDIP\-Net\-\. In the intra\-layer and inter\-layer branches of BDIP\-Net\-, the original message\-passing operation, which assigns a single scalar softmax attention weight to each neighboring edge message, is replaced with the PotNet\-style message\-passing operation based on a vector of feature\-wise sigmoid gates\. In addition, the scalar softmax interaction\-message attention module, which assigns a scalar attention weight to each aggregated intra\-layer and inter\-layer message, is removed, and the two messages are combined by direct summation\. This comparison evaluates the overall effect of replacing the proposed message\-passing operation and scalar softmax interaction\-message attention module with PotNet\-style message passing and direct message summation\.

Table 3\.Ablation study of BDIP\-Net\-and BDIP\-Net using MatterSim\-D3\-optimized structures\.DatasetModelMAE↓\\downarrowMSE↓\\downarrowRMSE↓\\downarrowR2↑\\uparrowBiDBBDIP\-Net\-0\.32±0\.030\.32\\pm 0\.030\.34±0\.120\.34\\pm 0\.120\.59±0\.110\.59\\pm 0\.110\.69±0\.040\.69\\pm 0\.04BDIP\-Net0\.30±0\.03\\mathbf\{0\.30\\pm 0\.03\}0\.32±0\.11\\mathbf\{0\.32\\pm 0\.11\}0\.56±0\.10\\mathbf\{0\.56\\pm 0\.10\}0\.73±0\.05\\mathbf\{0\.73\\pm 0\.05\}HetDBBDIP\-Net\-0\.10±0\.010\.10\\pm 0\.010\.03±0\.020\.03\\pm 0\.020\.18±0\.060\.18\\pm 0\.060\.91±0\.070\.91\\pm 0\.07BDIP\-Net0\.09±0\.01\\mathbf\{0\.09\\pm 0\.01\}0\.02±0\.02\\mathbf\{0\.02\\pm 0\.02\}0\.15±0\.06\\mathbf\{0\.15\\pm 0\.06\}0\.93±0\.05\\mathbf\{0\.93\\pm 0\.05\}SAMBABDIP\-Net\-0\.13±0\.040\.13\\pm 0\.040\.08±0\.060\.08\\pm 0\.060\.26±0\.110\.26\\pm 0\.110\.75±0\.110\.75\\pm 0\.11BDIP\-Net0\.10±0\.01\\mathbf\{0\.10\\pm 0\.01\}0\.04±0\.02\\mathbf\{0\.04\\pm 0\.02\}0\.20±0\.05\\mathbf\{0\.20\\pm 0\.05\}0\.83±0\.03\\mathbf\{0\.83\\pm 0\.03\}

Table[3](https://arxiv.org/html/2608.14640#S4.T3)shows that BDIP\-Net consistently outperforms BDIP\-Net\-across all three datasets\. Specifically, theR2R^\{2\}values increase from0\.690\.69to0\.730\.73on BiDB, from0\.910\.91to0\.930\.93on HetDB, and from0\.750\.75to0\.830\.83on SAMBA, respectively\. Similar to the comparisons with the other baseline models, one\-sided pairedtt\-tests on the combined four\-fold results from BiDB, HetDB, and SAMBA show that the improvements of BDIP\-Net over BDIP\-Net\-are statistically significant across all evaluation metrics \(p<0\.01p<0\.01\)\.

#### 4\.2\.4\.Runtime Analysis and Computational Efficiency

Structural Optimization Runtime\.BDIP\-Net adopts the efficient MatterSim\-D3 to performs explicit structural optimization rather than using expensive DFT\-PBE\-D3\-based energy and force calculations\. The complete structural optimization takes an average of 4\.5 seconds per material on BiDB, 3\.5 seconds on HetDB, and 13\.5 seconds on SAMBA\. The runtime differences among the three datasets mainly arise from their different optimization workflows and relaxation criteria\. The structural optimization time is longest for SAMBA because its workflow additionally includes fullzz\-scan,x​yxy\-scan, and relaxation procedures for twisted bilayers\. In contrast, obtaining aDFT\-PBE\-D3\-optimized structure typically requires hours of DFT\-based structural optimization\.

Training and Inference Runtime\.On BiDB,BDIP\-Netrequires575​s575~\\mathrm\{s\}for training, comparable to PotNet \(550​s550~\\mathrm\{s\}\) and SE\-MEGNET \(475​s475~\\mathrm\{s\}\), but longer than SE\-CGCNN \(400​s400~\\mathrm\{s\}\), SE\-PAINN \(300​s300~\\mathrm\{s\}\), and BiMat\-ML \(285​s285~\\mathrm\{s\}\)\. On HetDB, the training time ofBDIP\-Netis40​s40~\\mathrm\{s\}, compared with35​s35~\\mathrm\{s\}for PotNet,25​s25~\\mathrm\{s\}for SE\-MEGNET,20​s20~\\mathrm\{s\}for SE\-CGCNN, and12​s12~\\mathrm\{s\}for both SE\-PAINN and BiMat\-ML\. On SAMBA,BDIP\-Netrequires80​s80~\\mathrm\{s\}while PotNet, SE\-MEGNET, SE\-CGCNN, SE\-PAINN, and BiMat\-ML requires70​s70~\\mathrm\{s\},65​s65~\\mathrm\{s\},60​s60~\\mathrm\{s\},50​s50~\\mathrm\{s\}, and48​s48~\\mathrm\{s\}, respectively\. Despite its modestly higher training cost, BDIP\-Net achieves the best prediction accuracy across all three datasets\. All methods exhibit similar inference efficiency, with an average per\-sample inference time of approximately 0\.05 s\. Detailed training runtime comparisons are provided in Appendix[E\.2](https://arxiv.org/html/2608.14640#A5.SS2)\(Figure[6](https://arxiv.org/html/2608.14640#A5.F6)\)\.

#### 4\.2\.5\.Cross\-dataset Domain Generalization

We use HetDB and a subset of BiDB \(denoted as BiDB∗\) as the source\-domain training datasets, and the full SAMBA as the target\-domain test set\. BiDB contains 6,683 homobilayers, while HetDB contains 336 heterobilayers\. The SAMBA test set contains 980 bilayers, including 146 homobilayers and 834 heterobilayers\. Directly combining the full BiDB dataset with HetDB would make the training data heavily dominated by homobilayers\. BiDB∗contains 249 bilayers formed from 39 monolayers overlapping with SAMBA\. BiDB and HetDB have 39 and 13 monolayers overlapping with SAMBA, respectively, and the complete lists are provided in the Appendix[E\.3\.1](https://arxiv.org/html/2608.14640#A5.SS3.SSS1)\. Since PotNet achieves the best overall single\-domain performance among the baseline models , it is selected as the representative baseline for comparison\. To evaluate how prior exposure to constituent monolayers affects cross\-dataset generalization, we partition the SAMBA test set based on monolayer overlap with the training data, irrespective of layer ordering\. The test set comprises 147 bilayers with no monolayer overlap, 386 with one\-monolayer overlap, and 447 with two\-monolayer overlap\. To further distinguish previously observed bilayer compositions from novel combinations of known monolayers, the two\-monolayer\-overlap subset is divided into 121 exact\-pair\-overlap bilayers, in which the same monolayer pair appears in the training set \(regardless of layer ordering\), and 326 bilayers containing known monolayers but unseen pairings\.

Table 4\.BiDB∗\+HetDB\-to\-SAMBA domain generalization performance\.SettingNumber ofBilayersModelMAE↓\\downarrowMSE↓\\downarrowRMSE↓\\downarrowR2↑\\uparrowFull SAMBA980BDIP\-Net0\.270\.160\.400\.40PotNet0\.340\.220\.470\.15No\-monolayer\-overlap147BDIP\-Net0\.430\.380\.62\-0\.62PotNet0\.480\.340\.59\-0\.47One\-monolayer\-overlap386BDIP\-Net0\.310\.180\.43\-0\.35PotNet0\.410\.290\.54\-1\.13Two\-monolayer\-overlap447BDIP\-Net0\.170\.060\.250\.82PotNet0\.230\.120\.350\.67Exact\-pair\-overlap121BDIP\-Net0\.090\.020\.130\.98PotNet0\.090\.020\.140\.97Two\-monolayer\-overlap excludingexact\-pair\-overlap326BDIP\-Net0\.200\.080\.290\.54PotNet0\.290\.160\.400\.12

![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/box_plot_subset_bidb_hetdb_samba.png)Figure 4\.BiDB∗\+HetDB\-to\-SAMBA Domain Generalization Performance Breakdown![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/scatter_plot_examples.jpg)Figure 5\.BiDB∗\+HetDB\-to\-SAMBA Domain Generalization AE DistributionPrediction Performance\.As shown in Table[4](https://arxiv.org/html/2608.14640#S4.T4), BDIP\-Net achieves better prediction performance than PotNet on the full SAMBA test set of 980 bilayers, with an MAE of0\.270\.27, an MSE of0\.160\.16, an RMSE of0\.400\.40, and anR2R^\{2\}of0\.400\.40, compared with0\.340\.34,0\.220\.22,0\.470\.47, and0\.150\.15, respectively, for PotNet\. Figs\.[4](https://arxiv.org/html/2608.14640#S4.F4)and[5](https://arxiv.org/html/2608.14640#S4.F5)show that the prediction performance of both models improves as the degree of monolayer overlap increases\. On the 147 no\-monolayer\-overlap bilayers, BDIP\-Net and PotNet achieve MAEs of0\.430\.43and0\.480\.48, respectively\. On the 386 one\-monolayer\-overlap bilayers, the corresponding MAEs are0\.310\.31and0\.410\.41\. On the 447 two\-monolayer\-overlap bilayers, BDIP\-Net achieves an MAE of0\.170\.17, compared with0\.230\.23for PotNet\. The best results are obtained on the 121 exact\-pair\-overlap bilayers, where both models achieve an MAE of0\.090\.09\. After removing these exact\-pair\-overlap bilayers, BDIP\-Net still achieves an MAE of0\.200\.20, an RMSE of0\.290\.29, and anR2R^\{2\}of0\.540\.54on the remaining 326 two\-monolayer\-overlap bilayers, compared with0\.290\.29,0\.400\.40, and0\.120\.12, respectively, for PotNet\. These results demonstrate that BDIP\-Net provides stronger cross\-dataset generalization to unseen bilayer combinations of known monolayers\.

Representative Prediction Cases\.Table[5](https://arxiv.org/html/2608.14640#S4.T5)presents representative bilayers to illustrate how prediction difficulty varies with monolayer overlap and the relationship between bilayer and monolayer bandgaps\. In the two\-monolayer\-overlap cases, BDIP\-Net is more accurate than PotNet whether the bilayer bandgap is close to one monolayer bandgap or differs from both\. For one\-monolayer\-overlap cases, BDIP\-Net performs well when the bilayer bandgap is close to the overlapping monolayer bandgap, while both models show larger errors when it is closer to the non\-overlapping monolayer bandgap or differs from both\. In the no\-monolayer\-overlap case, both models have large prediction errors\.

Table 5\.Representative bilayers with different degrees of monolayer overlap in the BiDB∗\+HetDB\-to\-SAMBA cross\-domain generalization experiment\.BilayerBottomLayerTopLayerTwistAngleBottomBandgapTopBandgapBilayerBandgapBDIP\-NetPrediction \(AE\)PotNetPrediction \(AE\)Ga2Se/2ZrSe2Ga2Se2ZrSe20\.00∘0\.00^\{\\circ\}2\.020\.200\.210\.24 \(0\.03\)0\.55 \(0\.34\)PdS2/WTe2PdS2WTe260\.00∘60\.00^\{\\circ\}1\.170\.810\.130\.12 \(0\.01\)0\.32 \(0\.19\)As/14\{\}\_\{14\}/Hg4Te4As2HgTe19\.11∘19\.11^\{\\circ\}1\.390\.050\.080\.07 \(0\.01\)0\.45 \(0\.37\)B12N/12\{\}\_\{12\}/B7P7BNBP49\.11∘49\.11^\{\\circ\}4\.680\.910\.912\.08 \(1\.17\)2\.59 \(1\.68\)P6/Sn4S4P2SnS30\.00∘30\.00^\{\\circ\}1\.982\.321\.700\.61 \(1\.09\)0\.81 \(0\.89\)BP/BPBPBP60\.00∘60\.00^\{\\circ\}0\.910\.910\.402\.63 \(2\.23\)2\.64 \(2\.24\)

Note:Underlined bottom\- and top\-layer monolayers overlap with monolayers in the training dataset\.

HetDB\-to\-SAMBA\.We further evaluate the more challenging HetDB\-to\-SAMBA setting, where only the 336 HetDB heterobilayers are used for training\. Since the number of monolayers overlapping with SAMBA decreases from 39 to 13, the overall performance of both BDIP\-Net and PotNet degrades substantially\. Nevertheless, BDIP\-Net still generalizes better than PotNet to unseen bilayer combinations of known monolayers, particularly on the two\-monolayer\-overlap subset after excluding exact\-pair overlap\. Detailed results are provided in Appendix[E\.3\.2](https://arxiv.org/html/2608.14640#A5.SS3.SSS2)\.

## 5\.Conclusion

In this work, we developed a pipeline for stacked bilayer material property prediction, from bilayer structure construction and efficient MatterSim\-D3\-based structural optimization to ML\-based property prediction\. Compared with Bimat\-ML, which directly predicts bilayer properties from monolayer structures and stacking information, our workflow includes an additional bilayer structure construction and structural optimization step\. However, this step is performed efficiently using MatterSim\-D3, making the overall bilayer construction and optimization process significantly faster than workflows based on full DFT\-based structural optimization\. Moreover, the reconstructed bilayer graph enables GNN\-based models to explicitly learn intra\-layer and inter\-layer atomic interactions\. We further proposed BDIP\-Net, an interaction\-aware GNN that separates intra\-layer and inter\-layer interactions and fuses their messages through the interaction\-message attention module\. Experiments on BiDB, HetDB, and SAMBA show that BDIP\-Net consistently outperforms the baseline models under both structural optimization settings\. In addition, when bilayer structures optimized by MatterSim\-D3 are used as model inputs, the prediction performance is nearly identical to that obtained using the corresponding DFT\-optimized structures\. The cross\-dataset generalization experiments further show that BDIP\-Net provides stronger performance than PotNet on unseen bilayer combinations formed from monolayers observed during training, while prediction becomes more challenging when one or both target monolayers are absent from the training data\. These results demonstrate that the proposed workflow substantially reduces the structural optimization cost while maintaining accurate property prediction and supporting cross\-dataset generalization for stacked bilayer materials\.

## Code Availability

## Limitations and Ethical Considerations

This work is evaluated on publicly available computational materials databases and therefore inherits any biases, inaccuracies, or limited material diversity present in these datasets, which may affect model generalization\. Although the proposed framework substantially reduces the computational cost of bilayer structure generation and property prediction, its predictions should complement rather than replace high\-fidelity first\-principles calculations and experimental validation\. The study does not involve human participants or personal data; consequently, issues of data privacy and informed consent are not applicable\.

## Generative AI Usage

During the preparation of this manuscript, the authors used ChatGPT solely to assist with language editing\. The generative AI tool was not used for the conception of the research, experimental design, data analysis, interpretation of results, or generation of scientific content\. All scientific content, analyses, and conclusions were developed and verified by the authors, who take full responsibility for the accuracy and integrity of the work\.

## Acknowledgements

This work was supported in part by the National Institute of General Medical Sciences of National Institutes of Health under award P20GM139768, and the Arkansas Integrative Metabolic Research Center at the University of Arkansas\.

## References

- A\. L\. Araújo, P\. H\. Sophia, F\. Crasto de Lima, and A\. Fazzio \(2025\)A high\-throughput framework and database for twisted 2d van der waals bilayers\.npj Computational Materials\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p3.1)\.
- C\. Chen, W\. Ye, Y\. Zuo, C\. Zheng, and S\. P\. Ong \(2019\)Graph networks as a universal machine learning framework for molecules and crystals\.Chemistry of Materials31\(9\),pp\. 3564–3572\.Cited by:[§2](https://arxiv.org/html/2608.14640#S2.p1.1)\.
- X\. Chen, G\. Ru, and W\. Qi \(2024\)Structural embedding methods for machine learning models accelerate research on stacked 2d materials\.The Journal of Physical Chemistry C128\(37\),pp\. 15512–15521\.Cited by:[§2](https://arxiv.org/html/2608.14640#S2.p2.1),[§4\.1](https://arxiv.org/html/2608.14640#S4.SS1.p4.1)\.
- B\. Deng, Y\. Zhong, J\. Li,et al\.\(2023\)MatterSim: a deep learning atomistic model for materials science\.arXiv preprint arXiv:2306\.XXXX\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p2.1),[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p3.5)\.
- S\. Grimme, J\. Antony, S\. Ehrlich, and H\. Krieg \(2010\)A consistent and accurate ab initio parametrization of density functional dispersion correction \(dft\-d\) for the 94 elements h–pu\.The Journal of Chemical Physics132\(15\),pp\. 154104\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p1.1),[§1](https://arxiv.org/html/2608.14640#S1.p2.1),[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p2.4)\.
- S\. Haastrup, M\. Strange, M\. Pandey, T\. Deilmann, P\. S\. Schmidt, N\. F\. Hinsche, M\. N\. Gjerding, D\. Torelli, P\. M\. Larsen, A\. C\. Riis\-Jensen, J\. Gath, K\. W\. Jacobsen, J\. J\. Mortensen, T\. Olsen, and K\. S\. Thygesen \(2018\)The computational 2d materials database: high\-throughput modeling and discovery of atomically thin crystals\.2D Materials5\(4\),pp\. 042002\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p3.1)\.
- J\. E\. Jones \(1924\)On the determination of molecular fields\. ii\. from the equation of state of a gas\.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences106\(738\),pp\. 463–477\.External Links:[Document](https://dx.doi.org/10.1098/rspa.1924.0082)Cited by:[§3\.2\.1](https://arxiv.org/html/2608.14640#S3.SS2.SSS1.p2.7)\.
- J\. Klimeš, D\. R\. Bowler, and A\. Michaelides \(2011\)Van der waals density functionals applied to solids\.Physical Review B83\(19\),pp\. 195131\.Cited by:[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p2.4),[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p3.5)\.
- J\. Klimeš and A\. Michaelides \(2012\)Perspective: advances and challenges in treating van der waals dispersion forces in density functional theory\.The Journal of Chemical Physics137\(12\),pp\. 120901\.Cited by:[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p3.5)\.
- G\. Kresse and J\. Furthmüller \(1996\)Efficient iterative schemes for ab initio total\-energy calculations using a plane\-wave basis set\.Physical review B54\(16\),pp\. 11169\.Cited by:[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p3.5)\.
- G\. Li, A\. Luican, J\. Lopes dos Santos, A\. Castro Neto, A\. Reina, J\. Kong, and E\. Andrei \(2010\)Observation of van hove singularities in twisted graphene layers\.Nature physics6\(2\),pp\. 109–113\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p1.1)\.
- Y\. Lin, K\. Yan, Y\. Luo, Y\. Liu, X\. Qian, and S\. Ji \(2023\)Efficient approximations of complete interatomic potentials for crystal property prediction\.InProceedings of the 40th International Conference on Machine Learning,Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p2.1),[§2](https://arxiv.org/html/2608.14640#S2.p1.1),[§3\.2\.1](https://arxiv.org/html/2608.14640#S3.SS2.SSS1.p2.7),[§4\.1](https://arxiv.org/html/2608.14640#S4.SS1.p4.1)\.
- Y\. Mishin \(2021\)Machine\-learning interatomic potentials for materials science\.Acta Materialia214,pp\. 116980\.Cited by:[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p3.5)\.
- S\. Pakdel, A\. Rasmussen, A\. Taghizadeh, M\. Kruse, T\. Olsen, and K\. S\. Thygesen \(2024\)High\-throughput computational stacking reveals emergent properties in natural van der waals bilayers\.Nature Communications15\(1\),pp\. 932\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p3.1)\.
- J\. P\. Perdew, K\. Burke, and M\. Ernzerhof \(1996\)Generalized gradient approximation made simple\.Physical Review Letters77\(18\),pp\. 3865–3868\.Cited by:[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p2.4)\.
- Z\. Qiao, M\. Welborn, A\. Anandkumar, F\. R\. Manby, and T\. F\. Miller \(2020\)OrbNet: deep learning for quantum chemistry using symmetry\-adapted atomic\-orbital features\.The Journal of chemical physics153\(12\)\.Cited by:[§2](https://arxiv.org/html/2608.14640#S2.p1.1)\.
- M\. O\. Sauer, P\. M\. Lyngby, and K\. S\. Thygesen \(2025\)Dispersion\-corrected machine learning potentials for 2d van der waals materials\.Physical Review Materials9\(7\),pp\. 074007\.Cited by:[§1](https://arxiv.org/html/2608.14640#S1.p3.1),[§3\.1](https://arxiv.org/html/2608.14640#S3.SS1.p3.5)\.
- K\. Schütt, P\. Kindermans, H\. E\. Sauceda Felix, S\. Chmiela, A\. Tkatchenko, and K\. Müller \(2017\)Schnet: a continuous\-filter convolutional neural network for modeling quantum interactions\.Advances in neural information processing systems30\.Cited by:[§2](https://arxiv.org/html/2608.14640#S2.p1.1)\.
- A\. Vuong, M\. Van, C\. Zhao, and X\. Wu \(2026\)Property prediction of stacked bilayer materials: a multimodal learning approach\.InIJCAI,Cited by:[§E\.1](https://arxiv.org/html/2608.14640#A5.SS1.p2.1),[§E\.1](https://arxiv.org/html/2608.14640#A5.SS1.p3.1),[§2](https://arxiv.org/html/2608.14640#S2.p2.1),[§4\.1](https://arxiv.org/html/2608.14640#S4.SS1.p1.4),[§4\.1](https://arxiv.org/html/2608.14640#S4.SS1.p4.1)\.
- J\. P\. Wagner and P\. R\. Schreiner \(2015\)London dispersion in molecular chemistry—reconsidering steric effects\.Angewandte Chemie International Edition54\(42\),pp\. 12274–12296\.Cited by:[§3\.2\.1](https://arxiv.org/html/2608.14640#S3.SS2.SSS1.p2.7)\.
- T\. Xie and J\. C\. Grossman \(2018\)Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties\.Physical review letters120\(14\),pp\. 145301\.Cited by:[§2](https://arxiv.org/html/2608.14640#S2.p1.1)\.
- Z\. Zhou, J\. Shi, R\. Yang, Y\. Zou, and Q\. Li \(2023\)SlotGAT: slot\-based message passing for heterogeneous graphs\.InProceedings of the 40th International Conference on Machine Learning,Cited by:[§3\.2\.2](https://arxiv.org/html/2608.14640#S3.SS2.SSS2.p1.3)\.

## Appendix ABackground on Stacked 2D Materials Representation

A crystal structure is represented as a periodic repetition of a unit cell in three\-dimensional \(3D\) space\. Specifically, a crystal is described byM=\(A,P,L\)M=\(A,P,L\), whereA=\[𝐚1,…,𝐚N\]T∈ℝN×daA=\[\\mathbf\{a\}\_\{1\},\\dots,\\mathbf\{a\}\_\{N\}\]^\{T\}\\in\\mathbb\{R\}^\{N\\times d\_\{a\}\}is the atom feature matrix,P=\[𝐩1,…,𝐩N\]T∈ℝN×3P=\[\\mathbf\{p\}\_\{1\},\\dots,\\mathbf\{p\}\_\{N\}\]^\{T\}\\in\\mathbb\{R\}^\{N\\times 3\}is the matrix of Cartesian coordinates of atoms in the unit cell, andL=\[𝐥1,𝐥2,𝐥3\]T∈ℝ3×3L=\[\\mathbf\{l\}\_\{1\},\\mathbf\{l\}\_\{2\},\\mathbf\{l\}\_\{3\}\]^\{T\}\\in\\mathbb\{R\}^\{3\\times 3\}is the lattice matrix describing how the unit cell repeats in three directions through the lattice vectors𝐥1,𝐥2,𝐥3\\mathbf\{l\}\_\{1\},\\mathbf\{l\}\_\{2\},\\mathbf\{l\}\_\{3\}\. The lattice vectors define the lattice parametersa=‖𝐥1‖a=\\\|\\mathbf\{l\}\_\{1\}\\\|,b=‖𝐥2‖b=\\\|\\mathbf\{l\}\_\{2\}\\\|,c=‖𝐥3‖c=\\\|\\mathbf\{l\}\_\{3\}\\\|, and anglesγ=∠​\(𝐥1,𝐥2\)\\gamma=\\angle\(\\mathbf\{l\}\_\{1\},\\mathbf\{l\}\_\{2\}\),α=∠​\(𝐥2,𝐥3\)\\alpha=\\angle\(\\mathbf\{l\}\_\{2\},\\mathbf\{l\}\_\{3\}\),β=∠​\(𝐥1,𝐥3\)\\beta=\\angle\(\\mathbf\{l\}\_\{1\},\\mathbf\{l\}\_\{3\}\)\. The infinite crystal structure can then be represented as

P^=\{𝐩^i\|𝐩^i=𝐩i\+k1​𝐥1\+k2​𝐥2\+k3​𝐥3,k1,k2,k3∈ℤ,1≤i≤N\},\\hat\{P\}=\\left\\\{\\hat\{\\mathbf\{p\}\}\_\{i\}\\,\\middle\|\\,\\hat\{\\mathbf\{p\}\}\_\{i\}=\\mathbf\{p\}\_\{i\}\+k\_\{1\}\\mathbf\{l\}\_\{1\}\+k\_\{2\}\\mathbf\{l\}\_\{2\}\+k\_\{3\}\\mathbf\{l\}\_\{3\},\\;k\_\{1\},k\_\{2\},k\_\{3\}\\in\\mathbb\{Z\},\\;1\\leq i\\leq N\\right\\\},A^=\{𝐚^i\|𝐚^i=𝐚i,1≤i≤N\}\.\\hat\{A\}=\\left\\\{\\hat\{\\mathbf\{a\}\}\_\{i\}\\,\\middle\|\\,\\hat\{\\mathbf\{a\}\}\_\{i\}=\\mathbf\{a\}\_\{i\},\\;1\\leq i\\leq N\\right\\\}\.
Here,P^\\hat\{P\}contains all periodic images of atomic positions generated by lattice translations, whileA^\\hat\{A\}assigns the corresponding atomic features to each periodic image\.

A two\-dimensional \(2D\) material in the Computational 2D Materials Database \(C2DB\) is represented within this framework as a crystal structure with periodicity restricted to the in\-plane directions and is described byM=\(A,P,L\)M=\(A,P,L\)\. The lattice is defined by the vectors𝐥1,𝐥2,𝐥3\\mathbf\{l\}\_\{1\},\\mathbf\{l\}\_\{2\},\\mathbf\{l\}\_\{3\}, where𝐥1\\mathbf\{l\}\_\{1\}and𝐥2\\mathbf\{l\}\_\{2\}span the in\-plane Bravais lattice, while the out\-of\-plane vector satisfiesα=β=90∘\\alpha=\\beta=90^\{\\circ\}andc=t\+vc=t\+v, such that𝐥3=\(0,0,c\)\\mathbf\{l\}\_\{3\}=\(0,0,c\), wherettdenotes the thickness of the material andvvis a vacuum region introduced to prevent spurious interactions between periodic images along the out\-of\-plane direction\. The atomic positions are given by the Cartesian coordinate matrixP=\[𝐩1,…,𝐩N\]TP=\[\\mathbf\{p\}\_\{1\},\\dots,\\mathbf\{p\}\_\{N\}\]^\{T\}\. The infinite 2D material structure can then be represented as

P^=\{𝐩^i\|𝐩^i=𝐩i\+k1​𝐥1\+k2​𝐥2,k1,k2∈ℤ,1≤i≤N\},\\hat\{P\}=\\left\\\{\\hat\{\\mathbf\{p\}\}\_\{i\}\\,\\middle\|\\,\\hat\{\\mathbf\{p\}\}\_\{i\}=\\mathbf\{p\}\_\{i\}\+k\_\{1\}\\mathbf\{l\}\_\{1\}\+k\_\{2\}\\mathbf\{l\}\_\{2\},\\;k\_\{1\},k\_\{2\}\\in\\mathbb\{Z\},\\;1\\leq i\\leq N\\right\\\},A^=\{𝐚^i\|𝐚^i=𝐚i,1≤i≤N\}\.\\hat\{A\}=\\left\\\{\\hat\{\\mathbf\{a\}\}\_\{i\}\\,\\middle\|\\,\\hat\{\\mathbf\{a\}\}\_\{i\}=\\mathbf\{a\}\_\{i\},\\;1\\leq i\\leq N\\right\\\}\.
Here,P^\\hat\{P\}contains all periodic images of atomic positions generated by in\-plane lattice translations, whileA^\\hat\{A\}assigns the corresponding atomic features to each periodic image\.

A stacked two\-dimensional material in BiDB is constructed by combining two identical monolayer structures to form a homogeneous bilayer\. LetMmono=\(Amono,Pmono,Lmono\)M^\{\\mathrm\{mono\}\}=\(A^\{\\mathrm\{mono\}\},P^\{\\mathrm\{mono\}\},L^\{\\mathrm\{mono\}\}\)denote the monolayer structure, and define the bottom and top layers asMbotM^\{\\mathrm\{bot\}\}andMtopM^\{\\mathrm\{top\}\}, respectively\. Prior to stacking, a set of transformations may be applied toMtopM^\{\\mathrm\{top\}\}, including in\-plane rotation, translation, and an optional flip operation\. These transformations are parameterized by a stacking configuration matrixSSand applied to the atomic positions of the top layer asPtop←Ptop​SP^\{\\mathrm\{top\}\}\\leftarrow P^\{\\mathrm\{top\}\}S\. The two layers are then stacked along the out\-of\-plane direction\. The optimal interlayer distanced∗d\_\{\*\}is selected by minimizing the total energy of the stacked structure along the out\-of\-plane direction\. The bilayer lattice is constructed by preserving the in\-plane lattice vectors and redefining the out\-of\-plane lattice vector asLbi=\[𝐥1mono,𝐥2mono,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\[\\mathbf\{l\}\_\{1\}^\{\\mathrm\{mono\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{mono\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\]^\{T\}, where𝐥3bi=\(0,0,cbi\)\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,c^\{\\mathrm\{bi\}\}\)andcbi=2​tmono\+d\+vc^\{\\mathrm\{bi\}\}=2t^\{\\mathrm\{mono\}\}\+d\+v\. After stacking with the optimized interlayer distance, the bilayer structure is relaxed to optimize the atomic positions in all three spatial directions \(xx,yy, andzz\), allowing further reduction of the total energy and yielding a stable bilayer material\. The resulting bilayer is represented asMbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), whereAbiA^\{\\mathrm\{bi\}\}andPbiP^\{\\mathrm\{bi\}\}are obtained by combining the atomic features and positions ofMbotM^\{\\mathrm\{bot\}\}andMtopM^\{\\mathrm\{top\}\}\.

For HetDB and SAMBA, the bilayer structure is constructed using a commensurate\-supercell procedure\. LetMtop=\(Atop,Ptop,Ltop\)M^\{\\mathrm\{top\}\}=\(A^\{\\mathrm\{top\}\},P^\{\\mathrm\{top\}\},L^\{\\mathrm\{top\}\}\)andMbot=\(Abot,Pbot,Lbot\)M^\{\\mathrm\{bot\}\}=\(A^\{\\mathrm\{bot\}\},P^\{\\mathrm\{bot\}\},L^\{\\mathrm\{bot\}\}\)denote the two input monolayers, which may be identical for homobilayers or different for heterobilayers\. Given a pair of commensurate supercell matricesStopS^\{\\mathrm\{top\}\}andSbotS^\{\\mathrm\{bot\}\}, selected according to the dataset\-specific lattice matching criteria, the corresponding supercellsMtop,scM^\{\\mathrm\{top,sc\}\}andMbot,scM^\{\\mathrm\{bot,sc\}\}are first generated\. The in\-plane lattice vectors of the two supercells are then aligned by averaging the in\-plane components ofLtop,scL^\{\\mathrm\{top,sc\}\}andLbot,scL^\{\\mathrm\{bot,sc\}\}to define a common bilayer lattice, with𝐥1bi=12​\(𝐥1top,sc\+𝐥1bot,sc\)\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\}=\\frac\{1\}\{2\}\(\\mathbf\{l\}\_\{1\}^\{\\mathrm\{top,sc\}\}\+\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bot,sc\}\}\)and𝐥2bi=12​\(𝐥2top,sc\+𝐥2bot,sc\)\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\}=\\frac\{1\}\{2\}\(\\mathbf\{l\}\_\{2\}^\{\\mathrm\{top,sc\}\}\+\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bot,sc\}\}\)\. The bilayer lattice is defined asLbi=\[𝐥1bi,𝐥2bi,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\[\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\]^\{T\}, where𝐥3bi=\(0,0,v\+d\+tmono,top\+tmono,bot\)\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,v\+d\+t^\{\\mathrm\{mono,top\}\}\+t^\{\\mathrm\{mono,bot\}\}\)\. The two supercells are then stacked along the out\-of\-plane direction\. The optimal interlayer distanced∗d\_\{\*\}is first selected by minimizing the total energy of the stacked structure along the out\-of\-plane direction\. Depending on the dataset\-specific construction workflow, an in\-planex​yxy\-scan may also be performed after determiningd∗d\_\{\*\}\. The in\-plane stacking configuration is then selected as the lateral shift that gives the lowest total energy of the bilayer structure\. After that the bilayer is structurally relaxed to optimize the atomic positions in all three spatial directions \(xx,yy, andzz\), yielding a stable bilayer material\. The final bilayer is represented asMbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), whereAbiA^\{\\mathrm\{bi\}\}andPbiP^\{\\mathrm\{bi\}\}are obtained by combining the atomic features and positions ofMbot,scM^\{\\mathrm\{bot,sc\}\}andMtop,scM^\{\\mathrm\{top,sc\}\}\.

## Appendix BInitial Bilayer Structure Construction

This appendix provides the detailed process used to construct the initial stacked bilayer structures for BiDB, HetDB, and SAMBA before structural optimization\.

### B\.1\.BiDB

Algorithm 1Initial Bilayer Construction in BiDB\(ConstructBIDB\)\(\\mathrm\{ConstructBIDB\}\)1:Monolayer crystal

Mmono=\(Amono,Pmono,Lmono\)M^\{\\mathrm\{mono\}\}=\(A^\{\\mathrm\{mono\}\},P^\{\\mathrm\{mono\}\},L^\{\\mathrm\{mono\}\}\), bottom and top stacking configuration matrices

SbotS^\{\\mathrm\{bot\}\}and

StopS^\{\\mathrm\{top\}\}, initial interlayer distance

dd, in\-plane translations

\(δx,δy\)\(\\delta\_\{x\},\\delta\_\{y\}\), monolayer thickness

tmonot^\{\\mathrm\{mono\}\}, vacuum

vv
2:Initial bilayer structure

MinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}
3:Define the bottom and top stacking configuration matrices:

Sbot=\[c1botc3bot0c2botc4bot000c5bot\],Stop=\[c1topc3top0c2topc4top000c5top\]S^\{\\mathrm\{bot\}\}=\\begin\{bmatrix\}c\_\{1\}^\{\\mathrm\{bot\}\}&c\_\{3\}^\{\\mathrm\{bot\}\}&0\\\\ c\_\{2\}^\{\\mathrm\{bot\}\}&c\_\{4\}^\{\\mathrm\{bot\}\}&0\\\\ 0&0&c\_\{5\}^\{\\mathrm\{bot\}\}\\end\{bmatrix\},\\qquad S^\{\\mathrm\{top\}\}=\\begin\{bmatrix\}c\_\{1\}^\{\\mathrm\{top\}\}&c\_\{3\}^\{\\mathrm\{top\}\}&0\\\\ c\_\{2\}^\{\\mathrm\{top\}\}&c\_\{4\}^\{\\mathrm\{top\}\}&0\\\\ 0&0&c\_\{5\}^\{\\mathrm\{top\}\}\\end\{bmatrix\}
4:Construct the bottom and top monolayer copies:

Mbot=Mmono,Mtop=MmonoM^\{\\mathrm\{bot\}\}=M^\{\\mathrm\{mono\}\},\\qquad M^\{\\mathrm\{top\}\}=M^\{\\mathrm\{mono\}\}
5:Compute the fractional coordinates of the bottom and top layers:

Fbot=Pmono​\(Lmono\)−1,Ftop=Pmono​\(Lmono\)−1F^\{\\mathrm\{bot\}\}=P^\{\\mathrm\{mono\}\}\(L^\{\\mathrm\{mono\}\}\)^\{\-1\},\\qquad F^\{\\mathrm\{top\}\}=P^\{\\mathrm\{mono\}\}\(L^\{\\mathrm\{mono\}\}\)^\{\-1\}
6:Apply the bottom and top stacking configuration matrices:

Fbot=\(Fbot​Sbot\)mod1,Ftop=\(Ftop​Stop\)mod1F^\{\\mathrm\{bot\}\}=\(F^\{\\mathrm\{bot\}\}S^\{\\mathrm\{bot\}\}\)\\bmod 1,\\qquad F^\{\\mathrm\{top\}\}=\(F^\{\\mathrm\{top\}\}S^\{\\mathrm\{top\}\}\)\\bmod 1
7:Define the bilayer thickness and lattice matrix:

tbi=2​tmono\+d,𝐥3bi=\(0,0,v\+tbi\),Lbi=\[𝐥1mono,𝐥2mono,𝐥3bi\]T∈ℝ3×3t^\{\\mathrm\{bi\}\}=2t^\{\\mathrm\{mono\}\}\+d,\\quad\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,\\,v\+t^\{\\mathrm\{bi\}\}\),\\quad L^\{\\mathrm\{bi\}\}=\\begin\{bmatrix\}\\mathbf\{l\}\_\{1\}^\{\\mathrm\{mono\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{mono\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\end\{bmatrix\}^\{T\}\\in\\mathbb\{R\}^\{3\\times 3\}
8:Calculate the fractional out\-of\-plane layer translation:

δ​z=zmaxbot\+d‖𝐥3bi‖−zmintop\\delta z=z\_\{\\max\}^\{\\mathrm\{bot\}\}\+\\dfrac\{d\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\-z\_\{\\min\}^\{\\mathrm\{top\}\}
9:Define the layer translation vector

𝝉\\boldsymbol\{\\tau\}:

𝝉=\[δxδyδ​z\]\\boldsymbol\{\\tau\}=\\begin\{bmatrix\}\\delta\_\{x\}&\\delta\_\{y\}&\\delta z\\end\{bmatrix\}
10:Apply the layer translation to the top layer:

Ftop=\(Ftop\+𝟏​𝝉\)mod1,𝟏∈ℝN×1F^\{\\mathrm\{top\}\}=\\left\(F^\{\\mathrm\{top\}\}\+\\mathbf\{1\}\\boldsymbol\{\\tau\}\\right\)\\bmod 1,\\quad\\mathbf\{1\}\\in\\mathbb\{R\}^\{N\\times 1\}
11:Construct the bilayer fractional coordinates:

Fbi=Fbot∪FtopF^\{\\mathrm\{bi\}\}=F^\{\\mathrm\{bot\}\}\\cup F^\{\\mathrm\{top\}\}
12:Convert the bilayer fractional coordinates back to Cartesian coordinates:

Pbi=Fbi​LbiP^\{\\mathrm\{bi\}\}=F^\{\\mathrm\{bi\}\}L^\{\\mathrm\{bi\}\}
13:Construct the bilayer atomic features:

Abi=Amono∪AmonoA^\{\\mathrm\{bi\}\}=A^\{\\mathrm\{mono\}\}\\cup A^\{\\mathrm\{mono\}\}
14:return

Minitbi=\(Abi,Pbi,Lbi\)M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)

The initial bilayer construction workflow in the BiDB dataset follows Algorithm[1](https://arxiv.org/html/2608.14640#alg1)\. Given a monolayer crystalMmonoM^\{\\mathrm\{mono\}\}, bottom and top stacking configuration matricesSbotS^\{\\mathrm\{bot\}\}andStopS^\{\\mathrm\{top\}\}, an initial interlayer distancedd, an in\-plane translation\(δx,δy\)\(\\delta\_\{x\},\\delta\_\{y\}\), and the monolayer thicknesstmonot^\{\\mathrm\{mono\}\}, the goal is to construct the corresponding initial bilayer structureMinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}, whereMmono=\(Amono,Pmono,Lmono\)M^\{\\mathrm\{mono\}\}=\(A^\{\\mathrm\{mono\}\},P^\{\\mathrm\{mono\}\},L^\{\\mathrm\{mono\}\}\),Minitbi=\(Abi,Pbi,Lbi\)\.M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)\.

The stacking configuration matrices for the bottom and top layers are written as

\(9\)Sbot=\[c1botc3bot0c2botc4bot000c5bot\],Stop=\[c1topc3top0c2topc4top000c5top\],S^\{\\mathrm\{bot\}\}=\\begin\{bmatrix\}c\_\{1\}^\{\\mathrm\{bot\}\}&c\_\{3\}^\{\\mathrm\{bot\}\}&0\\\\ c\_\{2\}^\{\\mathrm\{bot\}\}&c\_\{4\}^\{\\mathrm\{bot\}\}&0\\\\ 0&0&c\_\{5\}^\{\\mathrm\{bot\}\}\\end\{bmatrix\},\\qquad S^\{\\mathrm\{top\}\}=\\begin\{bmatrix\}c\_\{1\}^\{\\mathrm\{top\}\}&c\_\{3\}^\{\\mathrm\{top\}\}&0\\\\ c\_\{2\}^\{\\mathrm\{top\}\}&c\_\{4\}^\{\\mathrm\{top\}\}&0\\\\ 0&0&c\_\{5\}^\{\\mathrm\{top\}\}\\end\{bmatrix\},wherec1,c2,c3,c4c\_\{1\},c\_\{2\},c\_\{3\},c\_\{4\}define the in\-plane transformation andc5c\_\{5\}defines the out\-of\-plane flip operation\. In the current version of BiDB, the bottom layer is kept unchanged, while the stacking information is applied to the top layer\. Therefore, the bottom\-layer stacking configuration matrix is set to

Sbot=I3=\[100010001\]\.S^\{\\mathrm\{bot\}\}=I\_\{3\}=\\begin\{bmatrix\}1&0&0\\\\ 0&1&0\\\\ 0&0&1\\end\{bmatrix\}\.
The user\-provided stacking information consists of the top\-layer stacking configuration matrixStopS^\{\\mathrm\{top\}\}and the in\-plane translations\(δx,δy\)\(\\delta\_\{x\},\\delta\_\{y\}\)\. From the input monolayer, two monolayer copies are first constructed asMbot=MmonoM^\{\\mathrm\{bot\}\}=M^\{\\mathrm\{mono\}\}andMtop=MmonoM^\{\\mathrm\{top\}\}=M^\{\\mathrm\{mono\}\}\. The atomic positions of both copies are then converted into fractional coordinates with respect to the monolayer lattice, whereFbot=Pmono​\(Lmono\)−1F^\{\\mathrm\{bot\}\}=P^\{\\mathrm\{mono\}\}\(L^\{\\mathrm\{mono\}\}\)^\{\-1\}andFtop=Pmono​\(Lmono\)−1F^\{\\mathrm\{top\}\}=P^\{\\mathrm\{mono\}\}\(L^\{\\mathrm\{mono\}\}\)^\{\-1\}\. The bottom and top stacking configuration matrices are applied separately, followed by modulo 1:

\(10\)Fbot=\(Fbot​Sbot\)mod1,Ftop=\(Ftop​Stop\)mod1\.F^\{\\mathrm\{bot\}\}=\\left\(F^\{\\mathrm\{bot\}\}S^\{\\mathrm\{bot\}\}\\right\)\\bmod 1,\\qquad F^\{\\mathrm\{top\}\}=\\left\(F^\{\\mathrm\{top\}\}S^\{\\mathrm\{top\}\}\\right\)\\bmod 1\.Thus, the bottom layer remains unchanged, while the top layer is transformed according to the user\-providedStopS^\{\\mathrm\{top\}\}\. The bilayer thickness is first defined astbi=2​tmono\+dt^\{\\mathrm\{bi\}\}=2t^\{\\mathrm\{mono\}\}\+d\. The bilayer out\-of\-plane lattice vector is then given by𝐥3bi=\(0,0,v\+tbi\)\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,v\+t^\{\\mathrm\{bi\}\}\), and the bilayer lattice matrix is defined asLbi=\[𝐥1mono,𝐥2mono,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\\begin\{bmatrix\}\\mathbf\{l\}\_\{1\}^\{\\mathrm\{mono\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{mono\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\end\{bmatrix\}^\{T\}\. The fractional out\-of\-plane layer translationδ​z\\delta zis computed to impose the initial interlayer distancedd\. First, the top surface of the bottom layer and the bottom surface of the top layer are computed as

\(11\)zmaxbot=max⁡\(f13,bot,…,fN3,bot\),zmintop=min⁡\(f13,top,…,fN3,top\)\.z\_\{\\max\}^\{\\mathrm\{bot\}\}=\\max\\left\(f\_\{1\}^\{3,\\mathrm\{bot\}\},\\dots,f\_\{N\}^\{3,\\mathrm\{bot\}\}\\right\),\\qquad z\_\{\\min\}^\{\\mathrm\{top\}\}=\\min\\left\(f\_\{1\}^\{3,\\mathrm\{top\}\},\\dots,f\_\{N\}^\{3,\\mathrm\{top\}\}\\right\)\.Then, the fractional out\-of\-plane translation is given by

\(12\)δ​z=zmaxbot\+d‖𝐥3bi‖−zmintop,\\delta z=z\_\{\\max\}^\{\\mathrm\{bot\}\}\+\\dfrac\{d\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\-z\_\{\\min\}^\{\\mathrm\{top\}\},wherefibot=\(fi1,bot,fi2,bot,fi3,bot\)f\_\{i\}^\{\\mathrm\{bot\}\}=\(f\_\{i\}^\{1,\\mathrm\{bot\}\},f\_\{i\}^\{2,\\mathrm\{bot\}\},f\_\{i\}^\{3,\\mathrm\{bot\}\}\)andfitop=\(fi1,top,fi2,top,fi3,top\)f\_\{i\}^\{\\mathrm\{top\}\}=\(f\_\{i\}^\{1,\\mathrm\{top\}\},f\_\{i\}^\{2,\\mathrm\{top\}\},f\_\{i\}^\{3,\\mathrm\{top\}\}\)denote the fractional coordinates of atomiiin the bottom and top layers, respectively\. The layer translation vector is then defined using the user\-provided in\-plane translations and the computed out\-of\-plane translation:

\(13\)𝝉=\[δxδyδ​z\]\.\\boldsymbol\{\\tau\}=\\begin\{bmatrix\}\\delta\_\{x\}&\\delta\_\{y\}&\\delta z\\end\{bmatrix\}\.The layer translation is applied only to the top layer:

\(14\)Ftop=\(Ftop\+𝟏​𝝉\)mod1,𝟏∈ℝN×1\.F^\{\\mathrm\{top\}\}=\\left\(F^\{\\mathrm\{top\}\}\+\\mathbf\{1\}\\boldsymbol\{\\tau\}\\right\)\\bmod 1,\\qquad\\mathbf\{1\}\\in\\mathbb\{R\}^\{N\\times 1\}\.
The bilayer fractional coordinates are constructed by combining the bottom layer and the translated top layer:Fbi=Fbot∪FtopF^\{\\mathrm\{bi\}\}=F^\{\\mathrm\{bot\}\}\\cup F^\{\\mathrm\{top\}\}\. These fractional coordinates are converted back to Cartesian coordinates using the bilayer lattice:Pbi=Fbi​LbiP^\{\\mathrm\{bi\}\}=F^\{\\mathrm\{bi\}\}L^\{\\mathrm\{bi\}\}\. Finally, the bilayer atomic features are constructed asAbi=Amono∪AmonoA^\{\\mathrm\{bi\}\}=A^\{\\mathrm\{mono\}\}\\cup A^\{\\mathrm\{mono\}\}, and the initial bilayer structure is given byMinitbi=\(Abi,Pbi,Lbi\)M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)\.

### B\.2\.HetDB

The HetDB initial bilayer construction consists of two main steps\. First, an optimal pair of compatible two\-dimensional supercell transformation matrices,CtopC^\{\\mathrm\{top\}\}andCbotC^\{\\mathrm\{bot\}\}, is identified for the two input monolayers using Algorithm[2](https://arxiv.org/html/2608.14640#alg2)\. Second, the selected two\-dimensional matrix pair is used as input to Algorithm[3](https://arxiv.org/html/2608.14640#alg3), where it is extended to the corresponding three\-dimensional supercell transformation matrices,StopS^\{\\mathrm\{top\}\}andSbotS^\{\\mathrm\{bot\}\}, by adding the out\-of\-plane flip operation\. The resulting three\-dimensional transformation matrix pair is then used to construct the initial heterogeneous bilayer structureMinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\.

Algorithm 2Supercell Pair Search in HetDB Dataset\(SupercellPairSearch\)\(\\mathrm\{SupercellPairSearch\}\)1:Monolayers

Mtop,MbotM^\{\\mathrm\{top\}\},M^\{\\mathrm\{bot\}\}; search bound

cmaxc\_\{\\max\}
2:Optimal pair

\(Copttop,Coptbot\)\(C^\{\\mathrm\{top\}\}\_\{\\mathrm\{opt\}\},C^\{\\mathrm\{bot\}\}\_\{\\mathrm\{opt\}\}\)
3:Generate

𝒞\(n\)\\mathcal\{C\}^\{\(n\)\}using Eq\. \([15](https://arxiv.org/html/2608.14640#A2.E15)\), where the entries

ci\(n\)c\_\{i\}^\{\(n\)\}of each

C\(n\)C^\{\(n\)\}are constrained by

cmaxc\_\{\\max\},

n∈\{top,bot\}n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}
4:Sort

𝒞\(n\)\\mathcal\{C\}^\{\(n\)\}by increasing

\|det\(C\(n\)\)\|\|\\det\(C^\{\(n\)\}\)\|,

n∈\{top,bot\}n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}
5:Initialize

\(Copttop,Coptbot\)←∅\(C^\{\\mathrm\{top\}\}\_\{\\mathrm\{opt\}\},C^\{\\mathrm\{bot\}\}\_\{\\mathrm\{opt\}\}\)\\leftarrow\\emptyset
6:Initialize

\(Nopt,ρopt,ϕopt\)←\(∞,∞,∞\)\(N\_\{\\mathrm\{opt\}\},\\rho\_\{\\mathrm\{opt\}\},\\phi\_\{\\mathrm\{opt\}\}\)\\leftarrow\(\\infty,\\infty,\\infty\)
7:for

Cbot∈𝒞botC^\{\\mathrm\{bot\}\}\\in\\mathcal\{C\}^\{\\mathrm\{bot\}\}do

8:for

Ctop∈𝒞topC^\{\\mathrm\{top\}\}\\in\\mathcal\{C\}^\{\\mathrm\{top\}\}do

9:Compute total number of atoms

NNusing Eq\. \([16](https://arxiv.org/html/2608.14640#A2.E16)\)

10:if

N\>NmaxN\>N\_\{\\max\}then

11:break

12:endif

13:Construct supercell lattices

L\(n\),scL^\{\(n\),\\mathrm\{sc\}\}using Eq\. \([17](https://arxiv.org/html/2608.14640#A2.E17)\)

14:Compute twist angle

ϕ\\phiusing Eq\. \([18](https://arxiv.org/html/2608.14640#A2.E18)\)

15:Compute

ρtop,ρbot\\rho^\{\\mathrm\{top\}\},\\rho^\{\\mathrm\{bot\}\}and

ρ\\rhousing Eq\. \([19](https://arxiv.org/html/2608.14640#A2.E19)\)

16:Compute

χtop,χbot\\chi^\{\\mathrm\{top\}\},\\chi^\{\\mathrm\{bot\}\}using Eq\. \([20](https://arxiv.org/html/2608.14640#A2.E20)\)

17:if

ϕ∉\[0∘,90∘\]\\phi\\notin\[0^\{\\circ\},90^\{\\circ\}\]or

ρ\>ρmax\\rho\>\\rho\_\{\\max\}or

χtop∉\[15∘,165∘\]\\chi^\{\\mathrm\{top\}\}\\notin\[15^\{\\circ\},165^\{\\circ\}\]or

χbot∉\[15∘,165∘\]\\chi^\{\\mathrm\{bot\}\}\\notin\[15^\{\\circ\},165^\{\\circ\}\]then

18:continue

19:endif

20:Construct common lattice

Lsc,commonL^\{\\mathrm\{sc\},\\mathrm\{common\}\}using Eq\. \([21](https://arxiv.org/html/2608.14640#A2.E21)\)

21:Compute deformation tensors

s\(n\)s^\{\(n\)\}and strain tensors

ε\(n\)\\varepsilon^\{\(n\)\}using Eq\. \([22](https://arxiv.org/html/2608.14640#A2.E22)\)–\([23](https://arxiv.org/html/2608.14640#A2.E23)\)

22:Compute maximum strain

εmax\\varepsilon\_\{\\text\{max\}\}using Eq\. \([24](https://arxiv.org/html/2608.14640#A2.E24)\)

23:if

εmax\>εthreshold\\varepsilon\_\{\\text\{max\}\}\>\\varepsilon\_\{\\text\{threshold\}\}then

24:continue

25:endif

26:if

\(N,ρ,ϕ\)<\(Nopt,ρopt,ϕopt\)\(N,\\rho,\\phi\)<\(N\_\{\\mathrm\{opt\}\},\\rho\_\{\\mathrm\{opt\}\},\\phi\_\{\\mathrm\{opt\}\}\)then

27:

\(Copttop,Coptbot\)←\(Ctop,Cbot\)\(C^\{\\mathrm\{top\}\}\_\{\\mathrm\{opt\}\},C^\{\\mathrm\{bot\}\}\_\{\\mathrm\{opt\}\}\)\\leftarrow\(C^\{\\mathrm\{top\}\},C^\{\\mathrm\{bot\}\}\)
28:

\(Nopt,ρopt,ϕopt\)←\(N,ρ,ϕ\)\(N\_\{\\mathrm\{opt\}\},\\rho\_\{\\mathrm\{opt\}\},\\phi\_\{\\mathrm\{opt\}\}\)\\leftarrow\(N,\\rho,\\phi\)
29:endif

30:endfor

31:endfor

32:return

\(Copttop,Coptbot\)\(C^\{\\mathrm\{top\}\}\_\{\\mathrm\{opt\}\},C^\{\\mathrm\{bot\}\}\_\{\\mathrm\{opt\}\}\)

The supercell pair search in HetDB dataset follows Algorithm[2](https://arxiv.org/html/2608.14640#alg2)\. The procedure begins by extracting the primitive in\-plane lattice matricesL\(n\)∈ℝ2×2L^\{\(n\)\}\\in\\mathbb\{R\}^\{2\\times 2\}forn∈\{top,bot\}n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}from the input monolayers, whereL\(n\)=\[𝐥1\(n\),𝐥2\(n\)\]TL^\{\(n\)\}=\[\\mathbf\{l\}\_\{1\}^\{\(n\)\},\\mathbf\{l\}\_\{2\}^\{\(n\)\}\]^\{T\}contains only the in\-plane lattice vectors, while the out\-of\-plane lattice vector𝐥3\(n\)\\mathbf\{l\}\_\{3\}^\{\(n\)\}is retained but not used in the supercell pair search procedure\. For each layern∈\{top,bot\}n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}, all possible integer two\-dimensional supercell transformation matrices are generated as

\(15\)C\(n\)=\(c1\(n\)c2\(n\)c3\(n\)c4\(n\)\),ci\(n\)∈ℤ,\|ci\(n\)\|≤cmax\.C^\{\(n\)\}=\\begin\{pmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}\\end\{pmatrix\},\\quad c\_\{i\}^\{\(n\)\}\\in\\mathbb\{Z\},\\quad\|c\_\{i\}^\{\(n\)\}\|\\leq c\_\{\\max\}\.whereci\(n\)c\_\{i\}^\{\(n\)\}are the integer coefficients defining the in\-plane supercell transformation andcmaxc\_\{\\max\}bounds their magnitude\. This defines the candidate sets𝒞\(n\)\\mathcal\{C\}^\{\(n\)\},n∈\{top,bot\}n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}\. The determinant\|det\(C\(n\)\)\|\|\\det\(C^\{\(n\)\}\)\|gives the scaling factor for the number of atoms, such that the supercell contains\|det\(C\(n\)\)\|\|\\det\(C^\{\(n\)\}\)\|times the number of atomsNmono\(n\)N^\{\\mathrm\{\(n\)\}\}\_\{\\mathrm\{mono\}\}in the monolayer unit cell\. The candidate sets𝒞\(n\)\\mathcal\{C\}^\{\(n\)\}are then sorted in increasing order of\|det\(C\(n\)\)\|\|\\det\(C^\{\(n\)\}\)\|to prioritize smaller supercells\. For any candidate pair\(Ctop,Cbot\)\(C^\{\\mathrm\{top\}\},C^\{\\mathrm\{bot\}\}\), the total number of atoms is computed as

\(16\)N=\|det\(Ctop\)\|​Nmonotop\+\|det\(Cbot\)\|​Nmonobot\.N=\|\\det\(C^\{\\mathrm\{top\}\}\)\|\\,N^\{\\mathrm\{top\}\}\_\{\\mathrm\{mono\}\}\+\|\\det\(C^\{\\mathrm\{bot\}\}\)\|\\,N^\{\\mathrm\{bot\}\}\_\{\\mathrm\{mono\}\}\.Due to this ordering, for a fixedCbotC^\{\\mathrm\{bot\}\}in the outer loop, the total number of atomsNNincreases monotonically asCtopC^\{\\mathrm\{top\}\}in the inner loop progresses\. Therefore, ifN\>NmaxN\>N\_\{\\max\}, the current candidate pair is discarded and the inner loop overCtopC^\{\\mathrm\{top\}\}is terminated early\.

For the remaining candidates, the corresponding in\-plane supercell lattices are constructed as

\(17\)L\(n\),sc=C\(n\)​L\(n\),L\(n\),sc=\[𝐥1\(n\),sc,𝐥2\(n\),sc\]T,n∈\{top,bot\}\.L^\{\(n\),\\mathrm\{sc\}\}=C^\{\(n\)\}L^\{\(n\)\},\\quad L^\{\(n\),\\mathrm\{sc\}\}=\[\\mathbf\{l\}\_\{1\}^\{\(n\),\\mathrm\{sc\}\},\\mathbf\{l\}\_\{2\}^\{\(n\),\\mathrm\{sc\}\}\]^\{T\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}\.whereL\(n\),sc∈ℝ2×2L^\{\(n\),\\mathrm\{sc\}\}\\in\\mathbb\{R\}^\{2\\times 2\}denotes the supercell lattice of layernn\. The twist angle is computed from the first lattice vectors of the two layers:

\(18\)ϕ=tan−1⁡\(𝐥1top,sc×𝐥1bot,sc𝐥1top,sc⋅𝐥1bot,sc\)\.\\phi=\\tan^\{\-1\}\\left\(\\frac\{\\mathbf\{l\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\_\{1\}\\times\\mathbf\{l\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\_\{1\}\}\{\\mathbf\{l\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\_\{1\}\\cdot\\mathbf\{l\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\_\{1\}\}\\right\)\.
The internal angleχ\(n\)\\chi^\{\(n\)\}is computed as

\(19\)χ\(n\)=cos−1⁡\(𝐥1\(n\),sc⋅𝐥2\(n\),sc‖𝐥1\(n\),sc‖​‖𝐥2\(n\),sc‖\)\.\\chi^\{\(n\)\}=\\cos^\{\-1\}\\left\(\\frac\{\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{1\}\\cdot\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{2\}\}\{\\\|\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{1\}\\\|\\,\\\|\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{2\}\\\|\}\\right\)\.
The lattice vector norm ratioρ\(n\)\\rho^\{\(n\)\}is computed as

\(20\)ρ\(n\)=max⁡\(‖𝐥1\(n\),sc‖‖𝐥2\(n\),sc‖,‖𝐥2\(n\),sc‖‖𝐥1\(n\),sc‖\),ρ=max⁡\(ρtop,ρbot\)\.\\rho^\{\(n\)\}=\\max\\left\(\\frac\{\\\|\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{1\}\\\|\}\{\\\|\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{2\}\\\|\},\\frac\{\\\|\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{2\}\\\|\}\{\\\|\\mathbf\{l\}^\{\(n\),\\mathrm\{sc\}\}\_\{1\}\\\|\}\\right\),\\quad\\rho=\\max\(\\rho^\{\\mathrm\{top\}\},\\rho^\{\\mathrm\{bot\}\}\)\.
Candidate pairs violating the geometric constraintsϕ∈\[0∘,90∘\]\\phi\\in\[0^\{\\circ\},90^\{\\circ\}\],χtop,χbot∈\[15∘,165∘\]\\chi^\{\\mathrm\{top\}\},\\chi^\{\\mathrm\{bot\}\}\\in\[15^\{\\circ\},165^\{\\circ\}\], andρ≤ρmax\\rho\\leq\\rho\_\{\\max\}are discarded\. For valid candidates, a common in\-plane lattice is constructed as

\(21\)Lsc,common=12​\(Ltop,sc\+Lbot,sc\),L^\{\\mathrm\{sc\},\\mathrm\{common\}\}=\\frac\{1\}\{2\}\\left\(L^\{\\mathrm\{top\},\\mathrm\{sc\}\}\+L^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\\right\),which defines the shared in\-plane lattice of the bilayer prior to introducing the out\-of\-plane lattice vector\.

The deformation tensor of each layer is then computed by mapping its unstrained supercell lattice onto the common lattice:

\(22\)𝒔\(n\)=\(L\(n\),sc\)−1​Lsc,common−I,n∈\{top,bot\},\\boldsymbol\{s\}^\{\(n\)\}=\\left\(L^\{\(n\),\\mathrm\{sc\}\}\\right\)^\{\-1\}L^\{\\mathrm\{sc\},\\mathrm\{common\}\}\-I,\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\},whereIIis the2×22\\times 2identity matrix\. The corresponding finite Lagrangian strain tensor is defined as

\(23\)𝜺\(n\)=12​\[𝒔\(n\)\+\(𝒔\(n\)\)T\+𝒔\(n\)​\(𝒔\(n\)\)T\]\.\\boldsymbol\{\\varepsilon\}^\{\(n\)\}=\\frac\{1\}\{2\}\\left\[\\boldsymbol\{s\}^\{\(n\)\}\+\\left\(\\boldsymbol\{s\}^\{\(n\)\}\\right\)^\{T\}\+\\boldsymbol\{s\}^\{\(n\)\}\\left\(\\boldsymbol\{s\}^\{\(n\)\}\\right\)^\{T\}\\right\]\.The maximum strain of the candidate pair is obtained from the eigenvalues of the strain tensors:

\(24\)εmax=maxn∈\{top,bot\}⁡maxi⁡\|λi​\(𝜺\(n\)\)\|\.\\varepsilon\_\{\\max\}=\\max\_\{n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}\}\\max\_\{i\}\\left\|\\lambda\_\{i\}\\left\(\\boldsymbol\{\\varepsilon\}^\{\(n\)\}\\right\)\\right\|\.Ifεmax\>εthreshold\\varepsilon\_\{\\max\}\>\\varepsilon\_\{\\mathrm\{threshold\}\}, the candidate pair is discarded\. Otherwise, it is compared with the current best solution using the lexicographic order of\(N,ρ,ϕ\)\(N,\\rho,\\phi\), where smaller values are preferred\. Finally, the best candidate pair is selected as\(Copttop,Coptbot\)\(C^\{\\mathrm\{top\}\}\_\{\\mathrm\{opt\}\},C^\{\\mathrm\{bot\}\}\_\{\\mathrm\{opt\}\}\)\.

Algorithm 3Initial Bilayer Construction in HetDB\(ConstructHETDB\)\(\\mathrm\{ConstructHETDB\}\)1:Monolayers

Mtop,MbotM^\{\\mathrm\{top\}\},M^\{\\mathrm\{bot\}\}; two\-dimensional supercell transformation matrices

Ctop,CbotC^\{\\mathrm\{top\}\},C^\{\\mathrm\{bot\}\}; flip operations

c5top,c5botc\_\{5\}^\{\\mathrm\{top\}\},c\_\{5\}^\{\\mathrm\{bot\}\}; monolayer thicknesses

tmono,top,tmono,bott^\{\\mathrm\{mono,top\}\},t^\{\\mathrm\{mono,bot\}\}; vacuum region

vv, interlayer distance

dd
2:Initial bilayer structure

MinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}
3:Extend the two\-dimensional supercell transformation matrices to three\-dimensional matrices:

C\(n\)=\[c1\(n\)c2\(n\)c3\(n\)c4\(n\)\],S\(n\)=\[c1\(n\)c2\(n\)0c3\(n\)c4\(n\)000c5\(n\)\],n∈\{top,bot\}\.C^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}\\end\{bmatrix\},\\quad S^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}&0\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}&0\\\\ 0&0&c\_\{5\}^\{\(n\)\}\\end\{bmatrix\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}\.
4:Construct supercells:

M\(n\),sc=\(A\(n\),sc,P\(n\),sc,L\(n\),sc\)=MakeSupercell​\(M\(n\),S\(n\)\)\.M^\{\(n\),\\mathrm\{sc\}\}=\(A^\{\(n\),\\mathrm\{sc\}\},P^\{\(n\),\\mathrm\{sc\}\},L^\{\(n\),\\mathrm\{sc\}\}\)=\\mathrm\{MakeSupercell\}\(M^\{\(n\)\},S^\{\(n\)\}\)\.
5:Compute common in\-plane lattice

Lsc,commonL^\{\\mathrm\{sc\},\\mathrm\{common\}\}using Eq\. \([21](https://arxiv.org/html/2608.14640#A2.E21)\)

6:Define bilayer lattice vectors:

𝐥1bi=𝐥1sc,common,𝐥2bi=𝐥2sc,common,𝐥3bi=\(0,0,v\+d\+tmono,top\+tmono,bot\)\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{1\}^\{\\mathrm\{sc\},\\mathrm\{common\}\},\\quad\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{2\}^\{\\mathrm\{sc\},\\mathrm\{common\}\},\\quad\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,v\+d\+t^\{\\mathrm\{mono,top\}\}\+t^\{\\mathrm\{mono,bot\}\}\)
7:Construct bilayer lattice:

Lbi=\[𝐥1bi,𝐥2bi,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\\begin\{bmatrix\}\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\end\{bmatrix\}^\{T\}
8:Convert to fractional coordinates:

F\(n\),sc=P\(n\),sc​\(Lbi\)−1,n∈\{top,bot\}F^\{\(n\),\\mathrm\{sc\}\}=P^\{\(n\),\\mathrm\{sc\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}
9:Compute current interlayer distance:

dcurrent=\(zmeantop−zmeanbot\)​∥𝐥3bi∥d\_\{\\mathrm\{current\}\}=\\left\(z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{top\}\}\-z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{bot\}\}\\right\)\\lVert\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\rVert
10:Compute fractional out\-of\-plane translation:

δ​z=d−dcurrent∥𝐥3bi∥\\delta z=\\frac\{d\-d\_\{\\mathrm\{current\}\}\}\{\\lVert\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\rVert\}
11:Define the layer translation vector without in\-plane translations:

𝝉=\[00δ​z\]\\boldsymbol\{\\tau\}=\\begin\{bmatrix\}0&0&\\delta z\\end\{bmatrix\}
12:Apply the layer translation to the top layer:

Ftop,sc←\(Ftop,sc\+𝟏​𝝉\)mod1,𝟏∈ℝNtop,sc×1F^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\leftarrow\\left\(F^\{\\mathrm\{top\},\\mathrm\{sc\}\}\+\\mathbf\{1\}\\boldsymbol\{\\tau\}\\right\)\\bmod 1,\\qquad\\mathbf\{1\}\\in\\mathbb\{R\}^\{N^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\times 1\}
13:Convert back to Cartesian coordinates:

P\(n\),sc=F\(n\),sc​Lbi,n∈\{top,bot\}P^\{\(n\),\\mathrm\{sc\}\}=F^\{\(n\),\\mathrm\{sc\}\}L^\{\\mathrm\{bi\}\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}
14:Combine layers:

Pbi=Ptop,sc∪Pbot,sc,Abi=Atop,sc∪Abot,scP^\{\\mathrm\{bi\}\}=P^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup P^\{\\mathrm\{bot\},\\mathrm\{sc\}\},\\quad A^\{\\mathrm\{bi\}\}=A^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup A^\{\\mathrm\{bot\},\\mathrm\{sc\}\}
15:return

Mbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)

The process of constructing an initial bilayer material in the HetDB dataset is summarized in Algorithm[3](https://arxiv.org/html/2608.14640#alg3)\. Given the top and bottom monolayersMtopM^\{\\mathrm\{top\}\}andMbotM^\{\\mathrm\{bot\}\}, the two\-dimensional supercell transformation matricesCtopC^\{\\mathrm\{top\}\}andCbotC^\{\\mathrm\{bot\}\}, the flip operationsc5topc\_\{5\}^\{\\mathrm\{top\}\}andc5botc\_\{5\}^\{\\mathrm\{bot\}\}, the monolayer thicknesses, the vacuum spacingvv, and the target interlayer distancedd, the goal is to construct the heterogeneous bilayer materialMinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\. The in\-plane translations are fixed as\(δ​x,δ​y\)=\(0,0\)\(\\delta x,\\delta y\)=\(0,0\)\. We adopt the same definition of the flip operation as in BiDB, where the in\-plane directions are preserved and the out\-of\-plane direction is either kept or reversed\. In HetDB, two configurations are considered: keeping both layers unflipped or flipping both layers simultaneously\.

The two\-dimensional supercell transformation matrices are first extended to three\-dimensional supercell transformation matrices by including the flip operation along the out\-of\-plane direction:

\(25\)S\(n\)=\[c1\(n\)c2\(n\)0c3\(n\)c4\(n\)000c5\(n\)\],n∈\{top,bot\},S^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}&0\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}&0\\\\ 0&0&c\_\{5\}^\{\(n\)\}\\end\{bmatrix\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\},wherec1\(n\),c2\(n\),c3\(n\),c4\(n\)c\_\{1\}^\{\(n\)\},c\_\{2\}^\{\(n\)\},c\_\{3\}^\{\(n\)\},c\_\{4\}^\{\(n\)\}are from the two\-dimensional supercell transformation matrixC\(n\)C^\{\(n\)\}andc5\(n\)=±1c\_\{5\}^\{\(n\)\}=\\pm 1controls the out\-of\-plane flip\. The functionMakeSupercell\\mathrm\{MakeSupercell\}is then applied to each monolayer structureM\(n\)=\(A\(n\),P\(n\),L\(n\)\)M^\{\(n\)\}=\(A^\{\(n\)\},P^\{\(n\)\},L^\{\(n\)\}\), whereA\(n\)A^\{\(n\)\}is the atomic features,P\(n\)P^\{\(n\)\}contains the Cartesian atomic coordinates, andL\(n\)L^\{\(n\)\}is the lattice matrix\. The corresponding supercell is constructed as

\(26\)M\(n\),sc=\(A\(n\),sc,P\(n\),sc,L\(n\),sc\)=MakeSupercell​\(M\(n\),S\(n\)\)\.M^\{\(n\),\\mathrm\{sc\}\}=\(A^\{\(n\),\\mathrm\{sc\}\},P^\{\(n\),\\mathrm\{sc\}\},L^\{\(n\),\\mathrm\{sc\}\}\)=\\mathrm\{MakeSupercell\}\(M^\{\(n\)\},S^\{\(n\)\}\)\.
Inside this function, the fractional coordinates are first computed from the Cartesian coordinates as

\(27\)F\(n\)=P\(n\)​\(L\(n\)\)−1,F^\{\(n\)\}=P^\{\(n\)\}\\left\(L^\{\(n\)\}\\right\)^\{\-1\},and the supercell lattice is obtained as

\(28\)L\(n\),sc=S\(n\)​L\(n\)\.L^\{\(n\),\\mathrm\{sc\}\}=S^\{\(n\)\}L^\{\(n\)\}\.The determinant ofS\(n\)S^\{\(n\)\}gives the number of primitive cells contained in the supercell, so the number of atoms becomes

\(29\)N\(n\),sc=\|detS\(n\)\|​N\(n\)=\|c1\(n\)​c4\(n\)−c2\(n\)​c3\(n\)\|​N\(n\),N^\{\(n\),\\mathrm\{sc\}\}=\\left\|\\det S^\{\(n\)\}\\right\|N^\{\(n\)\}=\\left\|c\_\{1\}^\{\(n\)\}c\_\{4\}^\{\(n\)\}\-c\_\{2\}^\{\(n\)\}c\_\{3\}^\{\(n\)\}\\right\|N^\{\(n\)\},sincec5\(n\)=±1c\_\{5\}^\{\(n\)\}=\\pm 1only changes the out\-of\-plane orientation\. The supercell fractional coordinates are generated by translating each𝐟i\(n\)∈F\(n\)\\mathbf\{f\}\_\{i\}^\{\(n\)\}\\in F^\{\(n\)\}with integer vectors𝐭=\(t1,t2,t3\)\\mathbf\{t\}=\(t\_\{1\},t\_\{2\},t\_\{3\}\), wheret1,t2,t3∈ℤt\_\{1\},t\_\{2\},t\_\{3\}\\in\\mathbb\{Z\}\. For a two\-dimensional material, translations are applied only in the in\-plane directions, sot3=0t\_\{3\}=0and𝐭=\(t1,t2,0\)\\mathbf\{t\}=\(t\_\{1\},t\_\{2\},0\)\. Only translated images inside the supercell are retained untilN\(n\),scN^\{\(n\),\\mathrm\{sc\}\}atoms are generated:

\(30\)F\(n\),sc=\{𝐟i,𝐭\(n\),sc=\(𝐟i\(n\)\+𝐭\)​\(S\(n\)\)−1\|𝐟i,𝐭\(n\),sc∈\[0,1\)3\}\.F^\{\(n\),\\mathrm\{sc\}\}=\\left\\\{\\mathbf\{f\}\_\{i,\\mathbf\{t\}\}^\{\(n\),\\mathrm\{sc\}\}=\\left\(\\mathbf\{f\}\_\{i\}^\{\(n\)\}\+\\mathbf\{t\}\\right\)\\left\(S^\{\(n\)\}\\right\)^\{\-1\}\\,\\middle\|\\,\\mathbf\{f\}\_\{i,\\mathbf\{t\}\}^\{\(n\),\\mathrm\{sc\}\}\\in\[0,1\)^\{3\}\\right\\\}\.
The retained fractional coordinates are converted back to Cartesian coordinates as

\(31\)P\(n\),sc=F\(n\),sc​L\(n\),sc\.P^\{\(n\),\\mathrm\{sc\}\}=F^\{\(n\),\\mathrm\{sc\}\}L^\{\(n\),\\mathrm\{sc\}\}\.The corresponding atomic features are copied fromA\(n\)A^\{\(n\)\}, formingA\(n\),scA^\{\(n\),\\mathrm\{sc\}\}\. Next, the common in\-plane latticeLsc,commonL^\{\\mathrm\{sc\},\\mathrm\{common\}\}is computed from the top and bottom supercell lattices using Eq\. \([21](https://arxiv.org/html/2608.14640#A2.E21)\)\. The bilayer out\-of\-plane length is first computed astbi=v\+d\+tmono,top\+tmono,bott^\{\\mathrm\{bi\}\}=v\+d\+t^\{\\mathrm\{mono,top\}\}\+t^\{\\mathrm\{mono,bot\}\}, which includes the top\-layer thickness, bottom\-layer thickness, interlayer distance, and vacuum spacing\. The bilayer lattice is then constructed by using the common lattice for the in\-plane directions, with𝐥1bi=𝐥1sc,common\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{1\}^\{\\mathrm\{sc\},\\mathrm\{common\}\},𝐥2bi=𝐥2sc,common\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{2\}^\{\\mathrm\{sc\},\\mathrm\{common\}\}, and𝐥3bi=\(0,0,tbi\)\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,t^\{\\mathrm\{bi\}\}\)\. The bilayer lattice matrix is then defined asLbi=\[𝐥1bi,𝐥2bi,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\\begin\{bmatrix\}\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\end\{bmatrix\}^\{T\}\. The atomic positions of both supercells are converted into fractional coordinates with respect to the bilayer lattice:

\(32\)F\(n\),sc=P\(n\),sc​\(Lbi\)−1,n∈\{top,bot\}\.F^\{\(n\),\\mathrm\{sc\}\}=P^\{\(n\),\\mathrm\{sc\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}\.The current interlayer distance is computed from the mean fractionalzz\-coordinates of the top and bottom layers:

\(33\)dcurrent=\(zmeantop−zmeanbot\)​∥𝐥3bi∥\.d\_\{\\mathrm\{current\}\}=\\left\(z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{top\}\}\-z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{bot\}\}\\right\)\\lVert\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\rVert\.To match the target interlayer distancedd, the required fractional out\-of\-plane translation is

\(34\)δ​z=d−dcurrent∥𝐥3bi∥\.\\delta z=\\frac\{d\-d\_\{\\mathrm\{current\}\}\}\{\\lVert\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\rVert\}\.The layer translation vector without in\-plane translations is then defined as

\(35\)𝝉=\[00δ​z\]\.\\boldsymbol\{\\tau\}=\\begin\{bmatrix\}0&0&\\delta z\\end\{bmatrix\}\.The layer translation is applied only to the top layerFtop,scF^\{\\mathrm\{top\},\\mathrm\{sc\}\}, following the same operation as in Eq\. \([14](https://arxiv.org/html/2608.14640#A2.E14)\)\. Finally, the translated fractional coordinates are converted back to Cartesian coordinates asPbot,sc=Fbot,sc​LbiP^\{\\mathrm\{bot\},\\mathrm\{sc\}\}=F^\{\\mathrm\{bot\},\\mathrm\{sc\}\}L^\{\\mathrm\{bi\}\}andPtop,sc=Ftop,sc​LbiP^\{\\mathrm\{top\},\\mathrm\{sc\}\}=F^\{\\mathrm\{top\},\\mathrm\{sc\}\}L^\{\\mathrm\{bi\}\}\. The two layers are then combined asPbi=Ptop,sc∪Pbot,scP^\{\\mathrm\{bi\}\}=P^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup P^\{\\mathrm\{bot\},\\mathrm\{sc\}\}andAbi=Atop,sc∪Abot,scA^\{\\mathrm\{bi\}\}=A^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup A^\{\\mathrm\{bot\},\\mathrm\{sc\}\}, yielding the initial heterogeneous bilayer structureMinitbi=\(Abi,Pbi,Lbi\)M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)\.

### B\.3\.SAMBA

For a given pair of monolayersMtopM^\{\\mathrm\{top\}\}andMbotM^\{\\mathrm\{bot\}\}, Algorithm[4](https://arxiv.org/html/2608.14640#alg4)summarizes the SAMBA initial bilayer construction\. Candidate two\-dimensional supercell transformation matrix pairs are first generated and filtered according to geometric compatibility constraints, including lattice\-area mismatch, internal\-angle mismatch, and lattice\-vector mismatch\. Each retained two\-dimensional matrix pair is then used as the input to Algorithm[5](https://arxiv.org/html/2608.14640#alg5), where it is extended to the corresponding three\-dimensional supercell transformation matrix pair by adding the out\-of\-plane flip operation\. The resulting three\-dimensional transformation matrix pair is then used to construct a candidate initial bilayer structure\. After all valid matrix pairs have been processed, duplicate bilayer structures are removed using RMSD\-based filtering, yielding the final set of non\-duplicated initial SAMBA bilayer structuresℳinitbi\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\.

Algorithm 4Initial Bilayer Structures Generation in SAMBA1:Monolayers

Mtop,MbotM^\{\\mathrm\{top\}\},M^\{\\mathrm\{bot\}\}; search bound

cmaxc\_\{\\max\}; thresholds

ϵarea\\epsilon\_\{\\mathrm\{area\}\},

ϵχ%\\epsilon\_\{\\chi\}^\{\\%\},

ϵχ∘\\epsilon\_\{\\chi\}^\{\\circ\},

ϵℓ\\epsilon\_\{\\ell\}, initial interlayer distance

dd
2:Set of initial bilayer structures

ℳinitbi\\mathcal\{M\}^\{\\mathrm\{bi\}\}\_\{\\mathrm\{init\}\}
3:Initialize

ℳinitbi←∅\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\\leftarrow\\emptyset
4:Generate candidate two\-dimensional supercell transformation matrices

𝒞\(n\)\\mathcal\{C\}^\{\(n\)\}using Eq\. \([15](https://arxiv.org/html/2608.14640#A2.E15)\), with

\|ci\(n\)\|≤cmax\|c\_\{i\}^\{\(n\)\}\|\\leq c\_\{\\max\}
5:for

Cbot∈𝒞botC^\{\\mathrm\{bot\}\}\\in\\mathcal\{C\}^\{\\mathrm\{bot\}\}do

6:for

Ctop∈𝒞topC^\{\\mathrm\{top\}\}\\in\\mathcal\{C\}^\{\\mathrm\{top\}\}do

7:Construct supercell lattices

L\(n\),scL^\{\(n\),\\mathrm\{sc\}\}using Eq\. \([17](https://arxiv.org/html/2608.14640#A2.E17)\)

8:Compute area mismatch

Δarea\\Delta\_\{\\mathrm\{area\}\}using Eq\. \([36](https://arxiv.org/html/2608.14640#A2.E36)\)

9:Compute internal angles

χtop,χbot\\chi^\{\\mathrm\{top\}\},\\chi^\{\\mathrm\{bot\}\}using Eq\. \([37](https://arxiv.org/html/2608.14640#A2.E37)\)

10:Compute angle mismatches

Δχ%,Δχ∘\\Delta\_\{\\chi\}^\{\\%\},\\Delta\_\{\\chi\}^\{\\circ\}using Eq\. \([38](https://arxiv.org/html/2608.14640#A2.E38)\)

11:Compute lattice\-vector mismatches

Δ𝐥,1,Δ𝐥,2\\Delta\_\{\\mathbf\{l\},1\},\\Delta\_\{\\mathbf\{l\},2\}using Eq\. \([39](https://arxiv.org/html/2608.14640#A2.E39)\)

12:if

Δarea\>ϵarea\\Delta\_\{\\mathrm\{area\}\}\>\\epsilon\_\{\\mathrm\{area\}\}or

Δχ%\>ϵχ%\\Delta\_\{\\chi\}^\{\\%\}\>\\epsilon\_\{\\chi\}^\{\\%\}or

Δχ∘\>ϵχ∘\\Delta\_\{\\chi\}^\{\\circ\}\>\\epsilon\_\{\\chi\}^\{\\circ\}or

Δ𝐥,1\>ϵℓ\\Delta\_\{\\mathbf\{l\},1\}\>\\epsilon\_\{\\ell\}or

Δ𝐥,2\>ϵℓ\\Delta\_\{\\mathbf\{l\},2\}\>\\epsilon\_\{\\ell\}then

13:continue

14:endif

15:Compute twist angle

ϕ\\phiusing Eq\. \([41](https://arxiv.org/html/2608.14640#A2.E41)\)

16:

Minitbi←ConstructSAMBA​\(Mtop,Mbot,Ctop,Cbot,d\)M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\\leftarrow\\mathrm\{ConstructSAMBA\}\(M^\{\\mathrm\{top\}\},M^\{\\mathrm\{bot\}\},C^\{\\mathrm\{top\}\},C^\{\\mathrm\{bot\}\},\\;d\)
17:Add

MinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}to

ℳinitbi\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}
18:endfor

19:endfor

20:Remove duplicate structures in

ℳbi\\mathcal\{M\}^\{\\mathrm\{bi\}\}via RMSD\-based filtering

21:return

ℳinitbi\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}

LetLtop=\[𝐥1top,𝐥2top\]T∈ℝ2×2L^\{\\mathrm\{top\}\}=\[\\mathbf\{l\}\_\{1\}^\{\\mathrm\{top\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{top\}\}\]^\{T\}\\in\\mathbb\{R\}^\{2\\times 2\}andLbot=\[𝐥1bot,𝐥2bot\]T∈ℝ2×2L^\{\\mathrm\{bot\}\}=\[\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bot\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bot\}\}\]^\{T\}\\in\\mathbb\{R\}^\{2\\times 2\}denote the primitive in\-plane lattice matrices of the top and bottom monolayers, respectively\. For each monolayer pair, the initial bilayer candidate set is first initialized asℳinitbi←∅\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\\leftarrow\\emptyset\. The candidate sets of two\-dimensional supercell transformation matrices,𝒞top\\mathcal\{C\}^\{\\mathrm\{top\}\}and𝒞bot\\mathcal\{C\}^\{\\mathrm\{bot\}\}, are then generated using Eq\. \([15](https://arxiv.org/html/2608.14640#A2.E15)\) with the search boundcmaxc\_\{\\max\}\. For each matrix pairCtop∈𝒞topC^\{\\mathrm\{top\}\}\\in\\mathcal\{C\}^\{\\mathrm\{top\}\}andCbot∈𝒞botC^\{\\mathrm\{bot\}\}\\in\\mathcal\{C\}^\{\\mathrm\{bot\}\}, the corresponding supercell lattice candidatesLtop,scL^\{\\mathrm\{top\},\\mathrm\{sc\}\}andLbot,scL^\{\\mathrm\{bot\},\\mathrm\{sc\}\}are computed using Eq\. \([17](https://arxiv.org/html/2608.14640#A2.E17)\)\. For each candidate pair, the area mismatch is defined as

\(36\)Δarea=\|‖𝐥1top,sc×𝐥2top,sc‖−‖𝐥1bot,sc×𝐥2bot,sc‖\|‖𝐥1top,sc×𝐥2top,sc‖\.\\Delta\_\{\\mathrm\{area\}\}=\\frac\{\\left\|\\\|\\mathbf\{l\}\_\{1\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\times\\mathbf\{l\}\_\{2\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\\|\-\\\|\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\\times\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\\\|\\right\|\}\{\\\|\\mathbf\{l\}\_\{1\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\times\\mathbf\{l\}\_\{2\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\\|\}\.
The internal angle of each supercell is defined as

\(37\)χ\(n\)=cos−1⁡\(𝐥1\(n\),sc⋅𝐥2\(n\),sc‖𝐥1\(n\),sc‖​‖𝐥2\(n\),sc‖\)\.\\chi^\{\(n\)\}=\\cos^\{\-1\}\\left\(\\frac\{\\mathbf\{l\}\_\{1\}^\{\(n\),\\mathrm\{sc\}\}\\cdot\\mathbf\{l\}\_\{2\}^\{\(n\),\\mathrm\{sc\}\}\}\{\\\|\\mathbf\{l\}\_\{1\}^\{\(n\),\\mathrm\{sc\}\}\\\|\\,\\\|\\mathbf\{l\}\_\{2\}^\{\(n\),\\mathrm\{sc\}\}\\\|\}\\right\)\.
The relative and absolute internal angle mismatches are defined as

\(38\)Δχ%=\|χtop−χbot\|χtop,Δχ∘=\|χtop−χbot\|\.\\Delta\_\{\\chi\}^\{\\%\}=\\frac\{\|\\chi^\{\\mathrm\{top\}\}\-\\chi^\{\\mathrm\{bot\}\}\|\}\{\\chi^\{\\mathrm\{top\}\}\},\\qquad\\Delta\_\{\\chi\}^\{\\circ\}=\|\\chi^\{\\mathrm\{top\}\}\-\\chi^\{\\mathrm\{bot\}\}\|\.
The lattice\-vector mismatch is defined as

\(39\)Δ𝐥,i=\|‖𝐥itop,sc‖−‖𝐥ibot,sc‖\|‖𝐥itop,sc‖,i∈\{1,2\}\.\\Delta\_\{\\mathbf\{l\},i\}=\\frac\{\\left\|\\\|\\mathbf\{l\}\_\{i\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\\|\-\\\|\\mathbf\{l\}\_\{i\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\\\|\\right\|\}\{\\\|\\mathbf\{l\}\_\{i\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\\|\},\\quad i\\in\\\{1,2\\\}\.
Letϵarea\\epsilon\_\{\\mathrm\{area\}\},ϵχ%\\epsilon\_\{\\chi\}^\{\\%\},ϵχ∘\\epsilon\_\{\\chi\}^\{\\circ\}, andϵℓ\\epsilon\_\{\\ell\}denote the corresponding thresholds\. A candidate pair is discarded if

\(40\)Δarea\>ϵarea​or​Δχ%\>ϵχ%​or​Δχ∘\>ϵχ∘​or​Δ𝐥,1\>ϵℓ​or​Δ𝐥,2\>ϵℓ\.\\Delta\_\{\\mathrm\{area\}\}\>\\epsilon\_\{\\mathrm\{area\}\}\\;\\;\\text\{or\}\\;\\;\\Delta\_\{\\chi\}^\{\\%\}\>\\epsilon\_\{\\chi\}^\{\\%\}\\;\\;\\text\{or\}\\;\\;\\Delta\_\{\\chi\}^\{\\circ\}\>\\epsilon\_\{\\chi\}^\{\\circ\}\\;\\;\\text\{or\}\\;\\;\\Delta\_\{\\mathbf\{l\},1\}\>\\epsilon\_\{\\ell\}\\;\\;\\text\{or\}\\;\\;\\Delta\_\{\\mathbf\{l\},2\}\>\\epsilon\_\{\\ell\}\.
Only pairs satisfying all constraints are retained\. For these valid pairs, the twist angle is computed as

\(41\)ϕ=tan−1⁡\(𝐥1top,sc×𝐥1bot,sc𝐥1top,sc⋅𝐥1bot,sc\)\.\\phi=\\tan^\{\-1\}\\left\(\\frac\{\\mathbf\{l\}\_\{1\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\times\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\}\{\\mathbf\{l\}\_\{1\}^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cdot\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\}\\right\)\.
The retained two\-dimensional matrix pairs\(Ctop,Cbot\)\(C^\{\\mathrm\{top\}\},C^\{\\mathrm\{bot\}\}\)are then used as inputs to Algorithm[5](https://arxiv.org/html/2608.14640#alg5), where they are extended to the corresponding three\-dimensional supercell transformation matrices\(Stop,Sbot\)\(S^\{\\mathrm\{top\}\},S^\{\\mathrm\{bot\}\}\)\. In the current SAMBA construction, the out\-of\-plane direction is kept unchanged for all materials\. Therefore, the out\-of\-plane transformation is fixed asc5top=c5bot=1c\_\{5\}^\{\\mathrm\{top\}\}=c\_\{5\}^\{\\mathrm\{bot\}\}=1, and no layer is flipped along the out\-of\-plane direction\. For each retained matrix pair, a candidate bilayer structure is constructed using Algorithm[5](https://arxiv.org/html/2608.14640#alg5)with fixed parametersddand\(δx,δy\)=\(0,0\)\(\\delta\_\{x\},\\delta\_\{y\}\)=\(0,0\), and the resulting structureMinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}is added to the candidate setℳinitbi\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\. After all retained matrix pairs have been processed, duplicate structures inℳinitbi\\mathcal\{M\_\{\\mathrm\{init\}\}\}^\{\\mathrm\{bi\}\}are removed using RMSD\-based geometric matching, with optimal alignment performed via the Kabsch algorithm and atomic correspondence determined using the Hungarian algorithm, yielding the final non\-duplicated setℳinitbi\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}for the givenMtopM^\{\\mathrm\{top\}\}andMbotM^\{\\mathrm\{bot\}\}monolayers\.

Algorithm 5Initial Bilayer Construction\(ConstructSAMBA\)\(\\mathrm\{ConstructSAMBA\}\)1:Monolayers

Mtop,MbotM^\{\\mathrm\{top\}\},M^\{\\mathrm\{bot\}\}; two\-dimensional supercell transformation matrices

Ctop,CbotC^\{\\mathrm\{top\}\},C^\{\\mathrm\{bot\}\}; monolayer thicknesses

tmono,top,tmono,bott^\{\\mathrm\{mono,top\}\},t^\{\\mathrm\{mono,bot\}\}; vacuum region

vv; interlayer distance

dd
2:Initial bilayer structure

MinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}
3:Extend the two\-dimensional supercell transformation matrices to three\-dimensional matrices:

C\(n\)=\[c1\(n\)c2\(n\)c3\(n\)c4\(n\)\],S\(n\)=\[c1\(n\)c2\(n\)0c3\(n\)c4\(n\)0001\],n∈\{top,bot\}\.C^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}\\end\{bmatrix\},\\quad S^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}&0\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}&0\\\\ 0&0&1\\end\{bmatrix\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}\.
4:Construct supercells:

M\(n\),sc=\(A\(n\),sc,P\(n\),sc,L\(n\),sc\)=MakeSupercell​\(M\(n\),S\(n\)\)\.M^\{\(n\),\\mathrm\{sc\}\}=\(A^\{\(n\),\\mathrm\{sc\}\},P^\{\(n\),\\mathrm\{sc\}\},L^\{\(n\),\\mathrm\{sc\}\}\)=\\mathrm\{MakeSupercell\}\(M^\{\(n\)\},S^\{\(n\)\}\)\.
5:Compute common in\-plane lattice

Lsc,commonL^\{\\mathrm\{sc\},\\mathrm\{common\}\}using Eq\. \([21](https://arxiv.org/html/2608.14640#A2.E21)\)

6:Define bilayer lattice vectors:

𝐥1bi=𝐥1sc,common,𝐥2bi=𝐥2sc,common,𝐥3bi=\(0,0,v\+d\+tmono,top\+tmono,bot\)\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{1\}^\{\\mathrm\{sc\},\\mathrm\{common\}\},\\quad\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{2\}^\{\\mathrm\{sc\},\\mathrm\{common\}\},\\quad\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,\\,v\+d\+t^\{\\mathrm\{mono,top\}\}\+t^\{\\mathrm\{mono,bot\}\}\)
7:Construct bilayer lattice:

Lbi=\[𝐥1bi,𝐥2bi,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\\begin\{bmatrix\}\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\end\{bmatrix\}^\{T\}
8:Convert atomic coordinates to fractional coordinates:

F\(n\),sc=P\(n\),sc​\(Lbi\)−1,n∈\{top,bot\}F^\{\(n\),\\mathrm\{sc\}\}=P^\{\(n\),\\mathrm\{sc\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}
9:Compute the fractional out\-of\-plane translation:

δ​z=zmaxbot\+d‖𝐥3bi‖−zmintop\\delta z=z\_\{\\max\}^\{\\mathrm\{bot\}\}\+\\frac\{d\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\-z\_\{\\min\}^\{\\mathrm\{top\}\}
10:Define the layer translation vector without in\-plane translations:

𝝉=\[00δ​z\]\\boldsymbol\{\\tau\}=\\begin\{bmatrix\}0&0&\\delta z\\end\{bmatrix\}
11:Apply the layer translation to the top layer:

Ftop,sc←\(Ftop,sc\+𝟏​𝝉\)mod1,𝟏∈ℝNtop,sc×1F^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\leftarrow\\left\(F^\{\\mathrm\{top\},\\mathrm\{sc\}\}\+\\mathbf\{1\}\\boldsymbol\{\\tau\}\\right\)\\bmod 1,\\qquad\\mathbf\{1\}\\in\\mathbb\{R\}^\{N^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\times 1\}
12:Convert back to Cartesian coordinates:

P\(n\),sc=F\(n\),sc​Lbi,n∈\{top,bot\}P^\{\(n\),\\mathrm\{sc\}\}=F^\{\(n\),\\mathrm\{sc\}\}L^\{\\mathrm\{bi\}\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\}
13:Combine layers:

Pbi=Ptop,sc∪Pbot,sc,Abi=Atop,sc∪Abot,scP^\{\\mathrm\{bi\}\}=P^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup P^\{\\mathrm\{bot\},\\mathrm\{sc\}\},\\quad A^\{\\mathrm\{bi\}\}=A^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup A^\{\\mathrm\{bot\},\\mathrm\{sc\}\}
14:return

Minitbi=\(Abi,Pbi,Lbi\)M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)

The process of constructing an initial bilayer structure in the SAMBA dataset is summarized in Algorithm[5](https://arxiv.org/html/2608.14640#alg5)\. Given two selected two\-dimensional supercell transformation matrices,CtopC^\{\\mathrm\{top\}\}andCbotC^\{\\mathrm\{bot\}\}, for the top and bottom monolayers, the corresponding top and bottom supercells are first constructed\. In the current SAMBA construction, the out\-of\-plane direction is kept unchanged for all materials\. Therefore, each integer2×22\\times 2supercell transformation matrixC\(n\)C^\{\(n\)\}is extended to a3×33\\times 3supercell transformation matrix by preserving the out\-of\-plane direction:

\(42\)C\(n\)=\[c1\(n\)c2\(n\)c3\(n\)c4\(n\)\],S\(n\)=\[c1\(n\)c2\(n\)0c3\(n\)c4\(n\)0001\],n∈\{top,bot\},C^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}\\end\{bmatrix\},\\quad S^\{\(n\)\}=\\begin\{bmatrix\}c\_\{1\}^\{\(n\)\}&c\_\{2\}^\{\(n\)\}&0\\\\ c\_\{3\}^\{\(n\)\}&c\_\{4\}^\{\(n\)\}&0\\\\ 0&0&1\\end\{bmatrix\},\\quad n\\in\\\{\\mathrm\{top\},\\mathrm\{bot\}\\\},whereci\(n\)∈ℤc\_\{i\}^\{\(n\)\}\\in\\mathbb\{Z\}define the in\-plane supercell transformation fromC\(n\)C^\{\(n\)\}, and the out\-of\-plane transformation is fixed asc5\(n\)=1c\_\{5\}^\{\(n\)\}=1\. The top and bottom supercells are then constructed using the sameMakeSupercell\\mathrm\{MakeSupercell\}function adopted in the HetDB construction procedure\. Specifically,Mtop,scM^\{\\mathrm\{top\},\\mathrm\{sc\}\}is computed fromMmono,topM^\{\\mathrm\{mono,top\}\}andStopS^\{\\mathrm\{top\}\}, whileMbot,scM^\{\\mathrm\{bot\},\\mathrm\{sc\}\}is computed fromMmono,botM^\{\\mathrm\{mono,bot\}\}andSbotS^\{\\mathrm\{bot\}\}, following Eq\. \([26](https://arxiv.org/html/2608.14640#A2.E26)\)\.

The common in\-plane latticeLsc,commonL^\{\\mathrm\{sc\},\\mathrm\{common\}\}is computed fromLtop,scL^\{\\mathrm\{top\},\\mathrm\{sc\}\}andLbot,scL^\{\\mathrm\{bot\},\\mathrm\{sc\}\}using Eq\. \([21](https://arxiv.org/html/2608.14640#A2.E21)\)\. The bilayer out\-of\-plane length is first computed astbi=v\+d\+tmono,top\+tmono,bott^\{\\mathrm\{bi\}\}=v\+d\+t^\{\\mathrm\{mono,top\}\}\+t^\{\\mathrm\{mono,bot\}\}, which includes the top\-layer thickness, bottom\-layer thickness, interlayer distance, and vacuum spacing\. The bilayer lattice is then constructed by using the common lattice for the in\-plane directions, with𝐥1bi=𝐥1sc,common\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{1\}^\{\\mathrm\{sc\},\\mathrm\{common\}\},𝐥2bi=𝐥2sc,common\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\}=\\mathbf\{l\}\_\{2\}^\{\\mathrm\{sc\},\\mathrm\{common\}\}, and𝐥3bi=\(0,0,tbi\)\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}=\(0,0,\\,t^\{\\mathrm\{bi\}\}\)\. The bilayer lattice matrix is then defined asLbi=\[𝐥1bi,𝐥2bi,𝐥3bi\]TL^\{\\mathrm\{bi\}\}=\\begin\{bmatrix\}\\mathbf\{l\}\_\{1\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{2\}^\{\\mathrm\{bi\}\},\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\end\{bmatrix\}^\{T\}\. The atomic positions are converted to fractional coordinates with respect to the bilayer lattice, givingFtop,sc=Ptop,sc​\(Lbi\)−1F^\{\\mathrm\{top\},\\mathrm\{sc\}\}=P^\{\\mathrm\{top\},\\mathrm\{sc\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\}andFbot,sc=Pbot,sc​\(Lbi\)−1F^\{\\mathrm\{bot\},\\mathrm\{sc\}\}=P^\{\\mathrm\{bot\},\\mathrm\{sc\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\}\. The fractional out\-of\-plane translation is then computed from the given initial interlayer distancedd, following Eq\. \([12](https://arxiv.org/html/2608.14640#A2.E12)\) used in BiDB\. Without in\-plane translations, the layer translation vector𝝉\\boldsymbol\{\\tau\}is defined using Eq\. \([35](https://arxiv.org/html/2608.14640#A2.E35)\)\. The fractional coordinates of the top layerFtop,scF^\{\\mathrm\{top\},\\mathrm\{sc\}\}are then updated using the same layer\-translation operation as in Eq\. \([14](https://arxiv.org/html/2608.14640#A2.E14)\)\. Finally, the coordinates are converted back to Cartesian form asPtop,sc=Ftop,sc​LbiP^\{\\mathrm\{top\},\\mathrm\{sc\}\}=F^\{\\mathrm\{top\},\\mathrm\{sc\}\}L^\{\\mathrm\{bi\}\}andPbot,sc=Fbot,sc​LbiP^\{\\mathrm\{bot\},\\mathrm\{sc\}\}=F^\{\\mathrm\{bot\},\\mathrm\{sc\}\}L^\{\\mathrm\{bi\}\}\. The two layers are then combined asPbi=Ptop,sc∪Pbot,scP^\{\\mathrm\{bi\}\}=P^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup P^\{\\mathrm\{bot\},\\mathrm\{sc\}\}andAbi=Atop,sc∪Abot,scA^\{\\mathrm\{bi\}\}=A^\{\\mathrm\{top\},\\mathrm\{sc\}\}\\cup A^\{\\mathrm\{bot\},\\mathrm\{sc\}\}, yielding the initial bilayer structureMinitbi=\(Abi,Pbi,Lbi\)M\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\)\.

## Appendix CStructural Optimization

After the dataset\-specific initial bilayer structuresMinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}are constructed, all initial structures are optimized using the MatterSim\-D3 framework\. Starting from an initial constructed bilayer structureMinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}, we first perform azz\-scan to determine the optimal interlayer distancedoptd\_\{\\mathrm\{opt\}\}, as summarized in Algorithm[7](https://arxiv.org/html/2608.14640#alg7)\. During thezz\-scan, the in\-plane translations are fixed asδ​x=0\\delta x=0andδ​y=0\\delta y=0, and the interlayer distanceddis varied using the layer\-translation procedure in Algorithm[6](https://arxiv.org/html/2608.14640#alg6)\. For SAMBA, afterdoptd\_\{\\mathrm\{opt\}\}is obtained, an additionalx​yxy\-scan is performed to select the lowest\-energy lateral stacking position\. Finally, the selected bilayer structure is relaxed using MatterSim\-D3 until the maximum atomic force is below the dataset\-specific maximum force threshold, yielding the optimized bilayer structureMbiM^\{\\mathrm\{bi\}\}\.

For BiDB, we apply a unifiedzz\-scan procedure to both magnetic and non\-magnetic systems\. This differs from the original BiDB generation workflow, where non\-magnetic systems are optimized using a SciPy\-based energy minimization and magnetic systems are optimized using a stepwisezz\-scan with energy normalized by the total number of electrons\. In our workflow, thezz\-scan uses the total MatterSim\-D3 energy without normalization\. The scan starts fromdstart=5\.0​Åd\_\{\\text\{start\}\}=5\.0~\\text\{\\AA \}and decreases todend=0​Åd\_\{\\text\{end\}\}=0~\\text\{\\AA \}\. The coarse and fine step sizes are set toscoarse=0\.25​Ås\_\{\\text\{coarse\}\}=0\.25~\\text\{\\AA \}andsfine=0\.05​Ås\_\{\\text\{fine\}\}=0\.05~\\text\{\\AA \}, respectively\. The coarse scan first identifies an approximate minimum\-energy distance\. Starting from this distance, the fine scan further decreases the interlayer distance and is terminated once the energy increases relative to the previous step\. The previous distance is selected asdoptd\_\{\\mathrm\{opt\}\}\. The final BFGS relaxation is performed until the maximum atomic force is below0\.01​eV/Å0\.01~\\mathrm\{eV/\\text\{\\AA \}\}, producingMbiM^\{\\mathrm\{bi\}\}\.

For HetDB, each initial heterogeneous bilayer structureMbiM^\{\\mathrm\{bi\}\}is optimized using the same MatterSim\-D3\-basedzz\-scan and BFGS relaxation workflow\. Thezz\-scan starts fromdstart=10\.0​Åd\_\{\\text\{start\}\}=10\.0~\\text\{\\AA \}and decreases todend=0​Åd\_\{\\text\{end\}\}=0~\\text\{\\AA \}, with coarse and fine step sizes ofscoarse=0\.25​Ås\_\{\\text\{coarse\}\}=0\.25~\\text\{\\AA \}andsfine=0\.05​Ås\_\{\\text\{fine\}\}=0\.05~\\text\{\\AA \}, respectively\. The coarse scan is used to locate the approximate minimum\-energy distance, and the fine scan is then performed from this distance with step sizesfines\_\{\\text\{fine\}\}\. The fine scan stops when the energy increases relative to the previous step, and the previous distance is selected as the optimal interlayer distancedoptd\_\{\\mathrm\{opt\}\}\. The final BFGS relaxation is performed until the maximum atomic force is below0\.05​eV/Å0\.05~\\mathrm\{eV/\\text\{\\AA \}\}, producingMbiM^\{\\mathrm\{bi\}\}\.

For SAMBA, each non\-duplicated initial bilayer structureMinitbi∈ℳinitbiM\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}\\in\\mathcal\{M\}\_\{\\mathrm\{init\}\}^\{\\mathrm\{bi\}\}is optimized separately\. The interlayer distance is first optimized using thezz\-scan procedure in Algorithm[7](https://arxiv.org/html/2608.14640#alg7), where the layer\-translation procedure in Algorithm[6](https://arxiv.org/html/2608.14640#alg6)is applied withδ​x=0\\delta x=0andδ​y=0\\delta y=0\. Thezz\-scan starts fromdstart=5\.0​Åd\_\{\\text\{start\}\}=5\.0~\\text\{\\AA \}and decreases todend=0​Åd\_\{\\text\{end\}\}=0~\\text\{\\AA \}, with coarse and fine step sizesscoarse=0\.5​Ås\_\{\\text\{coarse\}\}=0\.5~\\text\{\\AA \}andsfine=0\.05​Ås\_\{\\text\{fine\}\}=0\.05~\\text\{\\AA \}, respectively\. The coarse scan first locates an approximate minimum\-energy distance\. The fine scan then starts from this distance and decreasesdduntil the MatterSim\-D3 total energy increases relative to the previous step, after which the previous distance is selected asdoptd\_\{\\mathrm\{opt\}\}\. After obtainingdoptd\_\{\\mathrm\{opt\}\}, an additionalx​yxy\-scan is performed to optimize the lateral stacking registry\. This scan also uses the layer\-translation procedure in Algorithm[6](https://arxiv.org/html/2608.14640#alg6), but with the interlayer distance fixed atdoptd\_\{\\mathrm\{opt\}\}\. Specifically, a9×99\\times 9grid of in\-plane translations\(δ​x,δ​y\)\(\\delta x,\\delta y\), uniformly sampled from0to5/65/6along both in\-plane directions, is evaluated\. The translation pair with the lowest MatterSim\-D3 total energy is selected as\(δ​x,δ​y\)opt\(\\delta x,\\delta y\)\_\{\\mathrm\{opt\}\}\. The final BFGS relaxation is performed until the maximum atomic force is below0\.01​eV/Å0\.01~\\mathrm\{eV/\\text\{\\AA \}\}, producingMbiM^\{\\mathrm\{bi\}\}\.

Algorithm 6Layer Translation\(LayerTranslation\)\(\\mathrm\{LayerTranslation\}\)1:Bilayer material

Mbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), dataset type

𝒟\\mathcal\{D\}, in\-plane translations

\(δ​x,δ​y\)\(\\delta x,\\delta y\), interlayer distance

dd
2:Translated bilayer material

MtranslatedbiM^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}
3:Compute the layer\-separation cutoff

zcutz\_\{\\mathrm\{cut\}\}using Eq\. \([43](https://arxiv.org/html/2608.14640#A3.E43)\)

4:Split the bilayer positions using

zcutz\_\{\\mathrm\{cut\}\}:

Pbot=\{pibi∈Pbi\|pi3,bi<zcut\},Ptop=\{pibi∈Pbi\|pi3,bi≥zcut\}P^\{\\mathrm\{bot\}\}=\\left\\\{p\_\{i\}^\{\\mathrm\{bi\}\}\\in P^\{\\mathrm\{bi\}\}\\,\\middle\|\\,p\_\{i\}^\{3,\\mathrm\{bi\}\}<z\_\{\\mathrm\{cut\}\}\\right\\\},\\qquad P^\{\\mathrm\{top\}\}=\\left\\\{p\_\{i\}^\{\\mathrm\{bi\}\}\\in P^\{\\mathrm\{bi\}\}\\,\\middle\|\\,p\_\{i\}^\{3,\\mathrm\{bi\}\}\\geq z\_\{\\mathrm\{cut\}\}\\right\\\}
5:Convert the positions to fractional coordinates:

Fbot=Pbot​\(Lbi\)−1,Ftop=Ptop​\(Lbi\)−1F^\{\\mathrm\{bot\}\}=P^\{\\mathrm\{bot\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\},\\qquad F^\{\\mathrm\{top\}\}=P^\{\\mathrm\{top\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\}
6:Compute the required fractional out\-of\-plane layer translation

δ​z\\delta z:

δ​z=\{zmaxbot\+d‖𝐥3bi‖−zmintop,𝒟∈\{BiDB,SAMBA\},d−\(zmeantop−zmeanbot\)​‖𝐥3bi‖‖𝐥3bi‖,𝒟=HetDB\.\\delta z=\\begin\{cases\}z\_\{\\max\}^\{\\mathrm\{bot\}\}\+\\dfrac\{d\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\-z\_\{\\min\}^\{\\mathrm\{top\}\},&\\mathcal\{D\}\\in\\\{\\mathrm\{BiDB\},\\mathrm\{SAMBA\}\\\},\\\\\[8\.0pt\] \\dfrac\{d\-\\left\(z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{top\}\}\-z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{bot\}\}\\right\)\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\},&\\mathcal\{D\}=\\mathrm\{HetDB\}\.\\end\{cases\}
7:Define the layer translation vector:

𝝉=\[δ​xδ​yδ​z\]\\boldsymbol\{\\tau\}=\\begin\{bmatrix\}\\delta x&\\delta y&\\delta z\\end\{bmatrix\}
8:Apply the layer translation to the top layer:

Ftranslatedtop←\(Ftop\+𝟏​𝝉\)mod1,𝟏∈ℝNtop×1F^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}\\leftarrow\\left\(F^\{\\mathrm\{top\}\}\+\\mathbf\{1\}\\boldsymbol\{\\tau\}\\right\)\\bmod 1,\\qquad\\mathbf\{1\}\\in\\mathbb\{R\}^\{N\_\{\\mathrm\{top\}\}\\times 1\}
9:Convert the translated top layer back to Cartesian coordinates:

Ptranslatedtop=Ftranslatedtop​LbiP^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}=F^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}L^\{\\mathrm\{bi\}\}
10:Construct the translated bilayer:

Ptranslatedbi=Pbot∪PtranslatedtopP^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}=P^\{\\mathrm\{bot\}\}\\cup P^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}Mtranslatedbi=\(Abi,Ptranslatedbi,Lbi\)M^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\},L^\{\\mathrm\{bi\}\}\)
11:return

MtranslatedbiM^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}

The layer translation workflow follows Algorithm[6](https://arxiv.org/html/2608.14640#alg6)\. This procedure can be used in both thezz\-scan and thex​yxy\-scan to update the atomic coordinates of a bilayer structure\. Given a bilayer materialMbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), in\-plane translations\(δ​x,δ​y\)\(\\delta x,\\delta y\), and an interlayer distancedd, the goal is to translateMbiM^\{\\mathrm\{bi\}\}intoMtranslatedbiM^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}\. In thezz\-scan, the in\-plane translations are fixed asδ​x=0\\delta x=0andδ​y=0\\delta y=0, whileddis varied to find the optimal interlayer distancedoptd\_\{\\mathrm\{opt\}\}\. Afterdoptd\_\{\\mathrm\{opt\}\}is obtained,ddis fixed todoptd\_\{\\mathrm\{opt\}\}, and the same layer\-translation procedure can be used for thex​yxy\-scan by varying\(δ​x,δ​y\)\(\\delta x,\\delta y\)\. In the layer\-translation algorithm, the Cartesianzz\-coordinate thresholdzcutz\_\{\\mathrm\{cut\}\}used to separate the bottom and top layers is first defined as the midpoint between the lowest and highest atomic Cartesianzz\-coordinates of the bilayer:

\(43\)zcut=mini⁡pi3,bi\+maxi⁡pi3,bi2,z\_\{\\mathrm\{cut\}\}=\\frac\{\\min\_\{i\}p\_\{i\}^\{3,\\mathrm\{bi\}\}\+\\max\_\{i\}p\_\{i\}^\{3,\\mathrm\{bi\}\}\}\{2\},wherepi3,bip\_\{i\}^\{3,\\mathrm\{bi\}\}denotes the Cartesianzz\-coordinate of atomiiinPbiP^\{\\mathrm\{bi\}\}\. The bilayer atomic positions are then split into the bottom and top layers usingzcutz\_\{\\mathrm\{cut\}\}:

\(44\)Pbot=\{pibi∈Pbi\|pi3,bi<zcut\},Ptop=\{pibi∈Pbi\|pi3,bi≥zcut\}\.P^\{\\mathrm\{bot\}\}=\\left\\\{p\_\{i\}^\{\\mathrm\{bi\}\}\\in P^\{\\mathrm\{bi\}\}\\,\\middle\|\\,p\_\{i\}^\{3,\\mathrm\{bi\}\}<z\_\{\\mathrm\{cut\}\}\\right\\\},\\qquad P^\{\\mathrm\{top\}\}=\\left\\\{p\_\{i\}^\{\\mathrm\{bi\}\}\\in P^\{\\mathrm\{bi\}\}\\,\\middle\|\\,p\_\{i\}^\{3,\\mathrm\{bi\}\}\\geq z\_\{\\mathrm\{cut\}\}\\right\\\}\.Both layers are then converted into fractional coordinates using the bilayer latticeLbiL^\{\\mathrm\{bi\}\}, givingFbot=Pbot​\(Lbi\)−1F^\{\\mathrm\{bot\}\}=P^\{\\mathrm\{bot\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\}andFtop=Ptop​\(Lbi\)−1F^\{\\mathrm\{top\}\}=P^\{\\mathrm\{top\}\}\(L^\{\\mathrm\{bi\}\}\)^\{\-1\}\. To impose the target interlayer distancedd, the fractional out\-of\-plane translationδ​z\\delta zis computed according to the dataset construction protocol:

\(45\)δ​z=\{zmaxbot\+d‖𝐥3bi‖−zmintop,𝒟∈\{BiDB,SAMBA\},d−\(zmeantop−zmeanbot\)​‖𝐥3bi‖‖𝐥3bi‖,𝒟=HetDB\.\\delta z=\\begin\{cases\}z\_\{\\max\}^\{\\mathrm\{bot\}\}\+\\dfrac\{d\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\-z\_\{\\min\}^\{\\mathrm\{top\}\},&\\mathcal\{D\}\\in\\\{\\mathrm\{BiDB\},\\mathrm\{SAMBA\}\\\},\\\\\[8\.0pt\] \\dfrac\{d\-\\left\(z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{top\}\}\-z\_\{\\mathrm\{mean\}\}^\{\\mathrm\{bot\}\}\\right\)\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\}\{\\\|\\mathbf\{l\}\_\{3\}^\{\\mathrm\{bi\}\}\\\|\},&\\mathcal\{D\}=\\mathrm\{HetDB\}\.\\end\{cases\}wherezmintopz\_\{\\min\}^\{\\mathrm\{top\}\}denotes the minimum fractionalzz\-coordinate of the top layer,zmaxbotz\_\{\\max\}^\{\\mathrm\{bot\}\}denotes the maximum fractionalzz\-coordinate of the bottom layer, andzmeantopz\_\{\\mathrm\{mean\}\}^\{\\mathrm\{top\}\}andzmeanbotz\_\{\\mathrm\{mean\}\}^\{\\mathrm\{bot\}\}denote the mean fractionalzz\-coordinates of the top and bottom layers, respectively\. The layer translation vector𝝉\\boldsymbol\{\\tau\}is defined using Eq\. \([13](https://arxiv.org/html/2608.14640#A2.E13)\)\. The bottom layer is kept fixed, while the layer translation is applied only to the top layer following Eq\. \([14](https://arxiv.org/html/2608.14640#A2.E14)\)\. The translated top\-layer fractional coordinates are converted back to Cartesian coordinates using the same bilayer lattice asPtranslatedtop=Ftranslatedtop​LbiP^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}=F^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}L^\{\\mathrm\{bi\}\}\. Finally, the translated bilayer positions are constructed asPtranslatedbi=Pbot∪PtranslatedtopP^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}=P^\{\\mathrm\{bot\}\}\\cup P^\{\\mathrm\{top\}\}\_\{\\mathrm\{translated\}\}, and the translated bilayer material is given byMtranslatedbi=\(Abi,Ptranslatedbi,Lbi\)M^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\},L^\{\\mathrm\{bi\}\}\)\.

Algorithm 7z\-scan\(z−scan\)\(\\mathrm\{z\-scan\}\)1:Input:Bilayer material

Mbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), dataset type

𝒟\\mathcal\{D\}, search range

dstartd\_\{\\text\{start\}\},

dendd\_\{\\text\{end\}\}, coarse step size

scoarses\_\{\\text\{coarse\}\}, fine step size

sfines\_\{\\text\{fine\}\}, MLIP model

MatterSim​\-​D3\\mathrm\{MatterSim\\mbox\{\-\}D3\}
2:Output:Optimal interlayer distance

doptd\_\{\\mathrm\{opt\}\}
3:Stage 1: Coarse scan

4:

dcurrent←dstartd\_\{\\text\{current\}\}\\leftarrow d\_\{\\text\{start\}\}
5:

dopt←dcurrentd\_\{\\mathrm\{opt\}\}\\leftarrow d\_\{\\text\{current\}\}
6:

E​\(dopt\)←\+∞E\(d\_\{\\mathrm\{opt\}\}\)\\leftarrow\+\\infty
7:while

dcurrent\>dendd\_\{\\text\{current\}\}\>d\_\{\\text\{end\}\}do

8:

Mtranslatedbi←LayerTranslation​\(Mbi,𝒟,0,0,dcurrent\)M^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}\\leftarrow\\mathrm\{LayerTranslation\}\(M^\{\\mathrm\{bi\}\},\\mathcal\{D\},0,0,d\_\{\\text\{current\}\}\)
9:

E​\(dcurrent\)←MatterSim​\-​D3​\(Mtranslatedbi\)E\(d\_\{\\text\{current\}\}\)\\leftarrow\\mathrm\{MatterSim\\mbox\{\-\}D3\}\(M^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}\)
10:if

E​\(dcurrent\)<E​\(dopt\)E\(d\_\{\\text\{current\}\}\)<E\(d\_\{\\mathrm\{opt\}\}\)then

11:

E​\(dopt\)←E​\(dcurrent\)E\(d\_\{\\mathrm\{opt\}\}\)\\leftarrow E\(d\_\{\\text\{current\}\}\)
12:

dopt←dcurrentd\_\{\\mathrm\{opt\}\}\\leftarrow d\_\{\\text\{current\}\}
13:endif

14:

dcurrent←dcurrent−scoarsed\_\{\\text\{current\}\}\\leftarrow d\_\{\\text\{current\}\}\-s\_\{\\text\{coarse\}\}
15:endwhile

16:Stage 2: Fine scan

17:

dcurrent←doptd\_\{\\text\{current\}\}\\leftarrow d\_\{\\mathrm\{opt\}\}
18:

dprevious←dcurrentd\_\{\\text\{previous\}\}\\leftarrow d\_\{\\text\{current\}\}
19:

E​\(dprevious\)←E​\(dopt\)E\(d\_\{\\text\{previous\}\}\)\\leftarrow E\(d\_\{\\mathrm\{opt\}\}\)
20:whileTruedo

21:

dcurrent←dcurrent−sfined\_\{\\text\{current\}\}\\leftarrow d\_\{\\text\{current\}\}\-s\_\{\\text\{fine\}\}
22:

Mtranslatedbi←LayerTranslation​\(Mbi,𝒟,0,0,dcurrent\)M^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}\\leftarrow\\mathrm\{LayerTranslation\}\(M^\{\\mathrm\{bi\}\},\\mathcal\{D\},0,0,d\_\{\\text\{current\}\}\)
23:

E​\(dcurrent\)←MatterSim​\-​D3​\(Mtranslatedbi\)E\(d\_\{\\text\{current\}\}\)\\leftarrow\\mathrm\{MatterSim\\mbox\{\-\}D3\}\(M^\{\\mathrm\{bi\}\}\_\{\\mathrm\{translated\}\}\)
24:if

E​\(dcurrent\)\>E​\(dprevious\)E\(d\_\{\\text\{current\}\}\)\>E\(d\_\{\\text\{previous\}\}\)then

25:break

26:endif

27:

dprevious←dcurrentd\_\{\\text\{previous\}\}\\leftarrow d\_\{\\text\{current\}\}
28:

E​\(dprevious\)←E​\(dcurrent\)E\(d\_\{\\text\{previous\}\}\)\\leftarrow E\(d\_\{\\text\{current\}\}\)
29:endwhile

30:

dopt←dpreviousd\_\{\\mathrm\{opt\}\}\\leftarrow d\_\{\\text\{previous\}\}
31:return

doptd\_\{\\mathrm\{opt\}\}

Thezz\-scan procedure implemented in this work is summarized in Algorithm[7](https://arxiv.org/html/2608.14640#alg7)\. Because the originalzz\-scan implementations are dataset\-specific and may include details that are not fully specified in the published descriptions, our implementation may not exactly match the original code used to generate each dataset\. Starting from an initial bilayer structureMbiM^\{\\mathrm\{bi\}\}, the algorithm searches for the optimal interlayer distancedoptd\_\{\\mathrm\{opt\}\}by evaluating theMatterSim​\-​D3\\mathrm\{MatterSim\\mbox\{\-\}D3\}total energy over a sequence of candidate distances\. For each candidate distancedcurrentd\_\{\\text\{current\}\}, the layer\-translation procedure in Algorithm[6](https://arxiv.org/html/2608.14640#alg6)is applied with fixed in\-plane translationsδ​x=0\\delta x=0andδ​y=0\\delta y=0to update the separation between the two layers\. The scan is performed in two stages\. First, a coarse scan fromdstartd\_\{\\text\{start\}\}todendd\_\{\\text\{end\}\}with step sizescoarses\_\{\\text\{coarse\}\}is used to locate the approximate minimum\-energy distance\. Then, a fine scan starts from this distance and decreasesddwith step sizesfines\_\{\\text\{fine\}\}\. The fine scan continues as long as the energy does not increase relative to the previous step\. Once the energy increases, the scan is terminated, and the previous distance, corresponding to the lowest energy encountered in the fine scan, is selected as the optimal interlayer distancedoptd\_\{\\mathrm\{opt\}\}\.

## Appendix DBDIP\-Net For Bilayer Property Prediction

### D\.1\.BDIP\-Net Algorithm

Algorithm 8BDIP\-Net Algorithm for Bilayer Property Prediction1:Structurally optimized bilayer material

Mbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\); initial cutoffs

Rintra,RinterR^\{\\mathrm\{intra\}\},R^\{\\mathrm\{inter\}\}; required neighbor counts

Nintra,NinterN^\{\\mathrm\{intra\}\},N^\{\\mathrm\{inter\}\}; scaling coefficients

ϵCoulomb,ϵLondon\\epsilon^\{\\mathrm\{Coulomb\}\},\\epsilon^\{\\mathrm\{London\}\}; ground\-truth property

YY
2:Predicted bilayer property

Y^\\hat\{Y\}
3:

ℰintra←AdaptiveLayerNeighborSearch​\(Mbi,Rintra,Nintra,intra\)\\mathcal\{E\}^\{\\mathrm\{intra\}\}\\leftarrow\\mathrm\{AdaptiveLayerNeighborSearch\}\\left\(M^\{\\mathrm\{bi\}\},R^\{\\mathrm\{intra\}\},N^\{\\mathrm\{intra\}\},\\mathrm\{intra\}\\right\)
4:

ℰinter←AdaptiveLayerNeighborSearch​\(Mbi,Rinter,Ninter,inter\)\\mathcal\{E\}^\{\\mathrm\{inter\}\}\\leftarrow\\mathrm\{AdaptiveLayerNeighborSearch\}\\left\(M^\{\\mathrm\{bi\}\},R^\{\\mathrm\{inter\}\},N^\{\\mathrm\{inter\}\},\\mathrm\{inter\}\\right\)
5:foreach edge

\(i,j\)∈ℰintra\(i,j\)\\in\\mathcal\{E\}^\{\\mathrm\{intra\}\}do

6:

ei​jintra←MLPeintra​\(RBFCoulomb​\(−ϵCoulomb/di​j\)\)e\_\{ij\}^\{\\mathrm\{intra\}\}\\leftarrow\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{intra\}\}\\left\(\\mathrm\{RBF\}^\{\\mathrm\{Coulomb\}\}\\left\(\-\\epsilon^\{\\mathrm\{Coulomb\}\}/d\_\{ij\}\\right\)\\right\)
7:endfor

8:foreach edge

\(i,j\)∈ℰinter\(i,j\)\\in\\mathcal\{E\}^\{\\mathrm\{inter\}\}do

9:

ei​jinter←MLPeinter​\(RBFLondon​\(−ϵLondon/\(di​j\)6\)\)e\_\{ij\}^\{\\mathrm\{inter\}\}\\leftarrow\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{inter\}\}\\left\(\\mathrm\{RBF\}^\{\\mathrm\{London\}\}\\left\(\-\\epsilon^\{\\mathrm\{London\}\}/\(d\_\{ij\}\)^\{6\}\\right\)\\right\)
10:endfor

11:

𝒱\(0\)←Linear𝒱​\(Abi\)\\mathcal\{V\}^\{\(0\)\}\\leftarrow\\mathrm\{Linear\}\_\{\\mathcal\{V\}\}\\left\(A^\{\\mathrm\{bi\}\}\\right\)
12:for

t=0t=0to

T−1T\-1do

13:foreach atom

i∈𝒱i\\in\\mathcal\{V\}do

14:foreach interaction type

η∈\{intra,inter\}\\eta\\in\\\{\\mathrm\{intra\},\\mathrm\{inter\}\\\}do

15:

mi\(t\),η←𝟎m\_\{i\}^\{\(t\),\\eta\}\\leftarrow\\mathbf\{0\}
16:foreach edge

\(i,j\)∈ℰiη\(i,j\)\\in\\mathcal\{E\}^\{\\eta\}\_\{i\}do

17:

wi​j\(t\),η←vi\(t\)⊕vj\(t\)⊕ei​jηw\_\{ij\}^\{\(t\),\\eta\}\\leftarrow v\_\{i\}^\{\(t\)\}\\oplus v\_\{j\}^\{\(t\)\}\\oplus e\_\{ij\}^\{\\eta\}
18:

zi​j\(t\),η←MLPz\(t\),η​\(wi​j\(t\),η\)z\_\{ij\}^\{\(t\),\\eta\}\\leftarrow\\mathrm\{MLP\}\_\{z\}^\{\(t\),\\eta\}\\left\(w\_\{ij\}^\{\(t\),\\eta\}\\right\)
19:

αi​j\(t\),η←softmaxj∈𝒩\(i\)η⁡\(MLPα\(t\),η​\(zi​j\(t\),η\)\)\\alpha\_\{ij\}^\{\(t\),\\eta\}\\leftarrow\\operatorname\{softmax\}\_\{j\\in\\mathcal\{N\}^\{\\eta\}\_\{\(i\)\}\}\\left\(\\mathrm\{MLP\}\_\{\\alpha\}^\{\(t\),\\eta\}\\left\(z\_\{ij\}^\{\(t\),\\eta\}\\right\)\\right\)
20:

mi\(t\),η←mi\(t\),η\+αi​j\(t\),η​zi​j\(t\),ηm\_\{i\}^\{\(t\),\\eta\}\\leftarrow m\_\{i\}^\{\(t\),\\eta\}\+\\alpha\_\{ij\}^\{\(t\),\\eta\}z\_\{ij\}^\{\(t\),\\eta\}
21:endfor

22:endfor

23:

\[βi\(t\),intraβi\(t\),inter\]←softmax⁡\(MLPβ\(t\)​\(\[mi\(t\),intrami\(t\),inter\]\)\)\\begin\{bmatrix\}\\beta\_\{i\}^\{\(t\),\\mathrm\{intra\}\}\\\\ \\beta\_\{i\}^\{\(t\),\\mathrm\{inter\}\}\\end\{bmatrix\}\\leftarrow\\operatorname\{softmax\}\\left\(\\mathrm\{MLP\}\_\{\\beta\}^\{\(t\)\}\\left\(\\begin\{bmatrix\}m\_\{i\}^\{\(t\),\\mathrm\{intra\}\}\\\\ m\_\{i\}^\{\(t\),\\mathrm\{inter\}\}\\end\{bmatrix\}\\right\)\\right\)
24:

vi\(t\+1\)←vi\(t\)\+βi\(t\),intra​mi\(t\),intra\+βi\(t\),inter​mi\(t\),interv\_\{i\}^\{\(t\+1\)\}\\leftarrow v\_\{i\}^\{\(t\)\}\+\\beta\_\{i\}^\{\(t\),\\mathrm\{intra\}\}m\_\{i\}^\{\(t\),\\mathrm\{intra\}\}\+\\beta\_\{i\}^\{\(t\),\\mathrm\{inter\}\}m\_\{i\}^\{\(t\),\\mathrm\{inter\}\}
25:endfor

26:endfor

27:

Z𝒢←1\|𝒱\|​∑i∈𝒱vi\(T\)Z\_\{\\mathcal\{G\}\}\\leftarrow\\frac\{1\}\{\|\\mathcal\{V\}\|\}\\sum\_\{i\\in\\mathcal\{V\}\}v\_\{i\}^\{\(T\)\}
28:

Y^←MLPY^​\(Z𝒢\)\\hat\{Y\}\\leftarrow\\mathrm\{MLP\}\_\{\\hat\{Y\}\}\\left\(Z\_\{\\mathcal\{G\}\}\\right\)
29:

ℒ←\|Y^−Y\|\\mathcal\{L\}\\leftarrow\|\\hat\{Y\}\-Y\|
30:Update all trainable parameters by minimizing

ℒ\\mathcal\{L\}
31:return

Y^\\hat\{Y\}

Algorithm[8](https://arxiv.org/html/2608.14640#alg8)summarizes the complete training and prediction procedure of BDIP\-Net for bilayer material property prediction\. The algorithm takes as input a structurally optimized bilayer materialMbi=\(Abi,Pbi,Lbi\)M^\{\\mathrm\{bi\}\}=\(A^\{\\mathrm\{bi\}\},P^\{\\mathrm\{bi\}\},L^\{\\mathrm\{bi\}\}\), the initial intra\-layer and inter\-layer cutoff radiiRintraR^\{\\mathrm\{intra\}\}andRinterR^\{\\mathrm\{inter\}\}, the required neighbor countsNintraN^\{\\mathrm\{intra\}\}andNinterN^\{\\mathrm\{inter\}\}, the potential scaling coefficientsϵCoulomb\\epsilon^\{\\mathrm\{Coulomb\}\}andϵLondon\\epsilon^\{\\mathrm\{London\}\}, and the ground\-truth propertyYYduring training\. The output is the predicted bilayer propertyY^\\hat\{Y\}\. Lines 1–2 construct the interaction\-specific edge setsℰintra\\mathcal\{E\}^\{\\mathrm\{intra\}\}andℰinter\\mathcal\{E\}^\{\\mathrm\{inter\}\}using the adaptive layer neighbor search procedure\. Lines 3–8 compute the intra\-layer and inter\-layer edge embeddingsei​jintrae\_\{ij\}^\{\\mathrm\{intra\}\}andei​jintere\_\{ij\}^\{\\mathrm\{inter\}\}from the corresponding Coulomb\-based and London\-based potentials\. Line 9 projects the atomic features into the initial node embeddings𝒱\(0\)=\{vi\(0\)\}i=1N\\mathcal\{V\}^\{\(0\)\}=\\\{v\_\{i\}^\{\(0\)\}\\\}\_\{i=1\}^\{N\}\. Lines 10–24 performTTmessage\-passing convolution layers\. Within each convolution layer, the message\-passing operations are performed separately for each interaction typeη∈\{intra,inter\}\\eta\\in\\\{\\mathrm\{intra\},\\mathrm\{inter\}\\\}\. Lines 15–16 construct the interaction\-specific representationwi​j\(t\),ηw\_\{ij\}^\{\(t\),\\eta\}fromvi\(t\)v\_\{i\}^\{\(t\)\},vj\(t\)v\_\{j\}^\{\(t\)\}, andei​jηe\_\{ij\}^\{\\eta\}, and transformwi​j\(t\),ηw\_\{ij\}^\{\(t\),\\eta\}into the corresponding edge messagezi​j\(t\),ηz\_\{ij\}^\{\(t\),\\eta\}\. Lines 17–18 compute the interaction\-specific edge attention coefficientαi​j\(t\),η\\alpha\_\{ij\}^\{\(t\),\\eta\}fromzi​j\(t\),ηz\_\{ij\}^\{\(t\),\\eta\}, useαi​j\(t\),η\\alpha\_\{ij\}^\{\(t\),\\eta\}to weight the corresponding edge messagezi​j\(t\),ηz\_\{ij\}^\{\(t\),\\eta\}, and aggregate the weighted edge messages over neighbors of interaction typeη\\etato obtain the interaction\-specific node messagemi\(t\),ηm\_\{i\}^\{\(t\),\\eta\}\. Consequently, each node obtains an intra\-layer messagemi\(t\),intram\_\{i\}^\{\(t\),\\mathrm\{intra\}\}and an inter\-layer messagemi\(t\),interm\_\{i\}^\{\(t\),\\mathrm\{inter\}\}\. Lines 21–22 compute the interaction attention weightsβi\(t\),intra\\beta\_\{i\}^\{\(t\),\\mathrm\{intra\}\}andβi\(t\),inter\\beta\_\{i\}^\{\(t\),\\mathrm\{inter\}\}frommi\(t\),intram\_\{i\}^\{\(t\),\\mathrm\{intra\}\}andmi\(t\),interm\_\{i\}^\{\(t\),\\mathrm\{inter\}\}, respectively, use them to weight and fuse the two interaction\-specific node messages, and add the fused message to the residual node representationvi\(t\)v\_\{i\}^\{\(t\)\}to obtain the updated node embeddingvi\(t\+1\)v\_\{i\}^\{\(t\+1\)\}\. After theTTmessage\-passing convolution layers, Line 25 performs mean pooling over the final node embeddings\{vi\(T\)\}i∈𝒱\\\{v\_\{i\}^\{\(T\)\}\\\}\_\{i\\in\\mathcal\{V\}\}to obtain the graph representationZ𝒢Z\_\{\\mathcal\{G\}\}\. Line 26 transformsZ𝒢Z\_\{\\mathcal\{G\}\}into the predicted bilayer propertyY^\\hat\{Y\}\. During training, Lines 27–28 compute the MAE lossℒ\\mathcal\{L\}fromY^\\hat\{Y\}andYYand update all trainable parameters through backpropagation\. During inference, Line 29 returns the predicted bilayer propertyY^\\hat\{Y\}after the forward pass\.

### D\.2\.Interaction\-Specific Radial Basis Function Expansions

BDIP\-Net uses separate radial basis function \(RBF\) expansions for the Coulomb and London dispersion potential values\. Each expansion maps a scalar potential value to aded\_\{e\}\-dimensional feature vector\. The RBF centers are defined over the potential ranges induced by the distance intervald∈\[dmin,dmax\]d\\in\[d\_\{\\min\},d\_\{\\max\}\], wheredmin=1d\_\{\\min\}=1anddmax=8d\_\{\\max\}=8in our experiments\. For the Coulomb potential, the RBF centers are uniformly placed over

\[−ϵCoulombdmax,−ϵCoulombdmin\]\.\\left\[\-\\frac\{\\epsilon^\{\\mathrm\{Coulomb\}\}\}\{d\_\{\\max\}\},\-\\frac\{\\epsilon^\{\\mathrm\{Coulomb\}\}\}\{d\_\{\\min\}\}\\right\]\.For the London dispersion potential, the RBF centers are uniformly placed over

\[−ϵLondon\(dmax\)6,−ϵLondon\(dmin\)6\]\.\\left\[\-\\frac\{\\epsilon^\{\\mathrm\{London\}\}\}\{\(d\_\{\\max\}\)^\{6\}\},\-\\frac\{\\epsilon^\{\\mathrm\{London\}\}\}\{\(d\_\{\\min\}\)^\{6\}\}\\right\]\.
The resulting RBF\-expanded features are mapped todhd\_\{h\}\-dimensional edge embeddings by the independently parameterized networksMLPeintra\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{intra\}\}andMLPeinter\\mathrm\{MLP\}\_\{e\}^\{\\mathrm\{inter\}\}\.

### D\.3\.Adaptive Layer Neighbor Search

Algorithm 9Adaptive Layer Neighbor Search1:Structurally optimized bilayer material

Mbi=\(A,P,L\)M^\{\\mathrm\{bi\}\}=\(A,P,L\); initial cutoff radius

RR; required number of neighbors

NN; interaction type

η∈\{intra,inter\}\\eta\\in\\\{\\mathrm\{intra\},\\mathrm\{inter\}\\\}
2:Interaction\-specific edge set

ℰη\\mathcal\{E\}^\{\\eta\}
3:Compute the layer\-separation cutoff

zcutz^\{\\mathrm\{cut\}\}using Eq\. \([43](https://arxiv.org/html/2608.14640#A3.E43)\)

4:foreach atom

i∈Ai\\in Ado

5:if

zi\>zcutz\_\{i\}\>z^\{\\mathrm\{cut\}\}then

6:

ℓi←1\\ell\_\{i\}\\leftarrow 1
7:else

8:

ℓi←0\\ell\_\{i\}\\leftarrow 0
9:endif

10:endfor

11:

𝑒𝑛𝑜𝑢𝑔ℎ←False\\mathit\{enough\}\\leftarrow\\texttt\{False\}
12:whilenot

𝑒𝑛𝑜𝑢𝑔ℎ\\mathit\{enough\}do

13:Compute

𝒩i←FindNeighbors​\(Mbi,R\)\\mathcal\{N\}\_\{i\}\\leftarrow\\mathrm\{FindNeighbors\}\(M^\{\\mathrm\{bi\}\},R\)for each atom

i∈Ai\\in A
14:if

η=intra\\eta=\\mathrm\{intra\}then

15:foreach atom

i∈Ai\\in Ado

16:

𝒩i←\{\(j,𝐤,di​j\)∈𝒩i∣ℓi=ℓj\}\\mathcal\{N\}\_\{i\}\\leftarrow\\\{\(j,\\mathbf\{k\},d\_\{ij\}\)\\in\\mathcal\{N\}\_\{i\}\\mid\\ell\_\{i\}=\\ell\_\{j\}\\\}
17:endfor

18:elseif

η=inter\\eta=\\mathrm\{inter\}then

19:foreach atom

i∈Ai\\in Ado

20:

𝒩i←\{\(j,𝐤,di​j\)∈𝒩i∣ℓi≠ℓj\}\\mathcal\{N\}\_\{i\}\\leftarrow\\\{\(j,\\mathbf\{k\},d\_\{ij\}\)\\in\\mathcal\{N\}\_\{i\}\\mid\\ell\_\{i\}\\neq\\ell\_\{j\}\\\}
21:endfor

22:endif

23:

𝑒𝑛𝑜𝑢𝑔ℎ←True\\mathit\{enough\}\\leftarrow\\texttt\{True\}
24:foreach atom

i∈Ai\\in Ado

25:Sort

𝒩i\\mathcal\{N\}\_\{i\}by ascending

di​jd\_\{ij\}
26:if

\|𝒩i\|≥N\|\\mathcal\{N\}\_\{i\}\|\\geq Nthen

27:

di\(N\)←distance to theN\-th neighbor in​𝒩id\_\{i\}^\{\(N\)\}\\leftarrow\\text\{distance to the $N$\-th neighbor in \}\\mathcal\{N\}\_\{i\}
28:

𝒩ifiltered←\{\(j,𝐤,di​j\)∈𝒩i∣di​j≤di\(N\)\}\\mathcal\{N\}\_\{i\}^\{\\mathrm\{filtered\}\}\\leftarrow\\\{\(j,\\mathbf\{k\},d\_\{ij\}\)\\in\\mathcal\{N\}\_\{i\}\\mid d\_\{ij\}\\leq d\_\{i\}^\{\(N\)\}\\\}
29:else

30:

𝒩ifiltered←𝒩i\\mathcal\{N\}\_\{i\}^\{\\mathrm\{filtered\}\}\\leftarrow\\mathcal\{N\}\_\{i\}
31:

𝑒𝑛𝑜𝑢𝑔ℎ←False\\mathit\{enough\}\\leftarrow\\texttt\{False\}
32:endif

33:endfor

34:ifnot

𝑒𝑛𝑜𝑢𝑔ℎ\\mathit\{enough\}then

35:

R←2​RR\\leftarrow 2R
36:endif

37:endwhile

38:

ℰη←∅\\mathcal\{E\}^\{\\eta\}\\leftarrow\\emptyset
39:foreach atom

i∈Ai\\in Ado

40:foreach

\(j,𝐤,di​j\)∈𝒩ifiltered\(j,\\mathbf\{k\},d\_\{ij\}\)\\in\\mathcal\{N\}\_\{i\}^\{\\mathrm\{filtered\}\}do

41:Add edge

\(i,j\)𝐤\(i,j\)\_\{\\mathbf\{k\}\}to

ℰη\\mathcal\{E\}^\{\\eta\}
42:endfor

43:endfor

44:return

ℰη\\mathcal\{E\}^\{\\eta\}

Algorithm[9](https://arxiv.org/html/2608.14640#alg9)describes an adaptive neighbor search for a specified interaction type applied to a structurally optimized bilayer materialMbi=\(A,P,L\)M^\{\\mathrm\{bi\}\}=\(A,P,L\)\. The algorithm first computes the layer\-separation cutoffzcutz^\{\\mathrm\{cut\}\}using Eq\. \([43](https://arxiv.org/html/2608.14640#A3.E43)\) and assigns each atom to either the top or bottom layer\. Atoms withzi\>zcutz\_\{i\}\>z^\{\\mathrm\{cut\}\}are assigned to the top layer, whereas atoms withzi≤zcutz\_\{i\}\\leq z^\{\\mathrm\{cut\}\}are assigned to the bottom layer\. It then performs a neighbor search with an initial cutoff radiusRRand filters the resulting neighbors according to the selected interaction typeη\\eta\. Forη=intra\\eta=\\mathrm\{intra\}, only neighbors satisfyingℓi=ℓj\\ell\_\{i\}=\\ell\_\{j\}are retained, whereas forη=inter\\eta=\\mathrm\{inter\}, only neighbors satisfyingℓi≠ℓj\\ell\_\{i\}\\neq\\ell\_\{j\}are retained\. The filtered neighbor list of each atom is sorted by distance\. If an atom has at leastNNvalid neighbors, the distance to itsNN\-th nearest valid neighbor is used as an adaptive cutoff, and all valid neighbors within this distance are retained\. Consequently, the complete neighbor shell at theNN\-th distance is preserved, including any additional neighbors located at the same distance as theNN\-th neighbor\. Otherwise, the cutoff radius is doubled and the neighbor search is repeated until every atom has at leastNNvalid neighbors\. The resulting filtered neighbor lists are then used to construct the interaction\-specific edge setℰη\\mathcal\{E\}^\{\\eta\}\.

## Appendix EAdditional Experimental Results

### E\.1\.Datasets

We conduct our evaluation on BiDB, HetDB, and SAMBA because they provide complementary coverage of aligned homobilayers, heterobilayers, and twisted bilayer systems, together forming a diverse benchmark for evaluating our MatterSim\-D3\-based structural optimization workflow and the interaction\-aware graph learning capability of BDIP\-Net\.

Van der Waals Bilayer Database\.We follow Bimat\-ML\(Vuonget al\.,[2026](https://arxiv.org/html/2608.14640#bib.bib20)\)and use the same processed BiDB dataset\. The dataset contains homogeneous bilayers derived from monolayers in the C2DB database\. After removing 250 bilayers associated with 10 monolayers lacking CIF files, the processed dataset contains 10,899 valid bilayer structures\. We use the bandgap as the target prediction property and exclude bilayers without valid bandgap values, resulting in a final dataset of 6,683 bilayer materials formed from 940 unique monolayers\.

Van der Waals 2D Heterostructure Database\.We follow Bimat\-ML\(Vuonget al\.,[2026](https://arxiv.org/html/2608.14640#bib.bib20)\)and use the same HetDB dataset\. The dataset comprises 336 heterogeneous bilayer materials constructed from 38 distinct monolayers in the C2DB database\. Each bilayer is formed by stacking two different monolayers, and no duplicate bilayers are generated from the same monolayer pair under different stacking configurations\. Bandgap values are available for all 336 bilayers and are used as the target property\.

Simulation and Automated Methods for Bilayer Analysis\.We use the SAMBA dataset, which contains more than 18,000 twisted bilayer structures generated from 63 monolayers in the C2DB database\. Due to the high computational cost of DFT calculations, only 980 bilayer structures have undergone DFT structural optimization and property calculations, including 144 homogeneous bilayer materials and 836 heterogeneous bilayer materials\. Bandgap values are available for all 980 DFT\-calculated bilayers and are used as the target property\.

Table 6\.Bilayer types and bandgap relationships in the BiDB, HetDB, and SAMBA datasets\.DatasetBilayerTypeNumber ofBilayersMean / Min / MaxBandgap \(eV\)Monolayer Bandgap RangeClose to Monolayer BandgapsWithinOutsideCloseNeitherBiDBHomo66830\.93 / 0\.00 / 5\.76N/AN/A47231960HetDBHetero3360\.76 / 0\.00 / 3\.3576260149187SAMBAHomo1460\.66 / 0\.00 / 4\.45N/AN/A7571Hetero8340\.28 / 0\.00 / 2\.01435399607227

Dataset statistics\.Table[6](https://arxiv.org/html/2608.14640#A5.T6)summarizes the composition and bandgap distributions of the three datasets\. BiDB contains 6,683 homobilayers with a mean bandgap of0\.930\.93eV, while HetDB contains 336 heterobilayers with a mean bandgap of0\.760\.76eV\. SAMBA contains 146 homobilayers and 834 heterobilayers, with mean bandgaps of0\.660\.66and0\.280\.28eV, respectively\. The SAMBA heterobilayers are therefore concentrated at lower bandgap values, whereas the homobilayers cover a wider bandgap range\. In this statistical analysis, we use a threshold of0\.20\.2eV and consider a bilayer bandgap to be close to one of the two monolayer bandgaps forming the bilayer when their absolute difference does not exceed this threshold\. Under this criterion, 4,723 BiDB homobilayers are close to their corresponding monolayer bandgaps, while 1,960 are not\. For heterobilayers, the monolayer bandgap range is defined as the interval between the smaller and larger monolayer bandgaps\. In HetDB, 76 bilayers lie within this range and 260 lie outside it\. Among all HetDB samples, 149 are close to the smaller monolayer bandgap, while none are close to the larger one\. In SAMBA, 435 heterobilayers lie within the monolayer bandgap range and 399 lie outside it\. Among the 605 samples close to a monolayer bandgap, 603 are close to the smaller bandgap and only 2 are close to the larger one\. These results show that, after stacking, when a heterobilayer bandgap remains close to one of its monolayer bandgaps, it tends to be closer to the smaller monolayer bandgap\. For the SAMBA homobilayers, 75 samples are close to the corresponding monolayer bandgap, while 71 are not\. Overall, the differences in dataset size, bilayer type, bandgap distribution, and bilayer\-monolayer bandgap relationships introduce substantial distributional variation and make cross\-dataset prediction more challenging\.

### E\.2\.Running Time

Model Training and Inference Runtime\.

![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/training_time.jpg)Figure 6\.Training Runtime Comparison Across 3 DatasetsWe compare the training and inference runtimes ofBDIP\-Net, PotNet and the baseline models on BiDB, HetDB, and SAMBA\. Since the training and inference runtimes are nearly identical for structures optimized byDFT\-PBE\-D3andMatterSim\-D3, we report a single runtime for each dataset and model\. The training runtimes are summarized in Figure[6](https://arxiv.org/html/2608.14640#A5.F6)\. On BiDB,BDIP\-Netrequires575​s575~\\mathrm\{s\}for training, which is comparable to PotNet \(550​s550~\\mathrm\{s\}\) and SE\-MEGNET \(475​s475~\\mathrm\{s\}\), but slower than SE\-CGCNN \(400​s400~\\mathrm\{s\}\), SE\-PAINN \(300​s300~\\mathrm\{s\}\), and BiMat\-ML \(285​s285~\\mathrm\{s\}\)\. On HetDB,BDIP\-Netrequires40​s40~\\mathrm\{s\}, compared with35​s35~\\mathrm\{s\}for PotNet,25​s25~\\mathrm\{s\}for SE\-MEGNET,20​s20~\\mathrm\{s\}for SE\-CGCNN, and12​s12~\\mathrm\{s\}for both SE\-PAINN and BiMat\-ML\. On SAMBA,BDIP\-Netrequires80​s80~\\mathrm\{s\}, compared with70​s70~\\mathrm\{s\}for PotNet,65​s65~\\mathrm\{s\}for SE\-MEGNET,60​s60~\\mathrm\{s\}for SE\-CGCNN,50​s50~\\mathrm\{s\}for SE\-PAINN, and48​s48~\\mathrm\{s\}for BiMat\-ML\. Across all datasets and models, the average per\-sample inference time is approximately0\.05​s0\.05~\\mathrm\{s\}\. Although the training and inference runtimes are nearly identical when usingMatterSim\-D3\-optimized andDFT\-PBE\-D3\-optimized structures as model inputs, the time required to obtain these two types of input structures is significantly different\. Starting from the same initial bilayer structure,MatterSim\-D3can generate a structurally optimized bilayer structure in only a few seconds, whereas obtaining aDFT\-PBE\-D3\-optimized structure typically requires hours of DFT\-based structural optimization\. Thus,BDIP\-Netcombined withMatterSim\-D3\-based structure generation provides a computationally practical pipeline for bilayer property prediction\.

### E\.3\.Cross\-Domain Generalization

#### E\.3\.1\.Monolayer Overlap

Table[7](https://arxiv.org/html/2608.14640#A5.T7)lists the 63 unique monolayers that appear in the SAMBA dataset\. Monolayers overlapping with the BiDB dataset are marked with the subscript\(B\)\(B\), while those overlapping with both the BiDB and HetDB datasets are marked with\(B,H\)\(B,H\)\. This annotation is used to identify monolayer overlap between the training and test datasets in the cross\-domain generalization analysis presented in the main paper\.

Table 7\.List of monolayers in the SAMBA dataset\. Subscripts\(B\)\(B\)and\(B,H\)\(B,H\)indicate monolayers overlapping with BiDB only and with both BiDB and HetDB, respectively\.As2\\mathrm\{As\_\{2\}\}As2​S3\\mathrm\{As\_\{2\}S\_\{3\}\}As\\mathrm\{As\}Bi2​Te3\(B\)\\mathrm\{Bi\_\{2\}Te\_\{3\}\}\_\{\\scriptscriptstyle\(B\)\}As4​S6Orpiment\\mathrm\{As\_\{4\}S\_\{6\}^\{Orpiment\}\}Bi2\\mathrm\{Bi\_\{2\}\}Bi2​Se3\\mathrm\{Bi\_\{2\}Se\_\{3\}\}Bi2​SeTe2\(B\)\\mathrm\{Bi\_\{2\}SeTe\_\{2\}\}\_\{\\scriptscriptstyle\(B\)\}BN\(B,H\)\\mathrm\{BN\}\_\{\\scriptscriptstyle\(B,H\)\}As4​S6Anorpiment\\mathrm\{As\_\{4\}S\_\{6\}^\{Anorpiment\}\}BP\\mathrm\{BP\}C2\(B,H\)\\mathrm\{C\_\{2\}\}\_\{\\scriptscriptstyle\(B,H\)\}CdTe\\mathrm\{CdTe\}Ga2​S2\(B,H\)\\mathrm\{Ga\_\{2\}S\_\{2\}\}\_\{\\scriptscriptstyle\(B,H\)\}Ga2​Se2\(B,H\)\\mathrm\{Ga\_\{2\}Se\_\{2\}\}\_\{\\scriptscriptstyle\(B,H\)\}Ga2​Te2\(B\)\\mathrm\{Ga\_\{2\}Te\_\{2\}\}\_\{\\scriptscriptstyle\(B\)\}Ge2\\mathrm\{Ge\_\{2\}\}GeS\(B\)\\mathrm\{GeS\}\_\{\\scriptscriptstyle\(B\)\}Hf2​Te6\\mathrm\{Hf\_\{2\}Te\_\{6\}\}HgTe\(B\)\\mathrm\{HgTe\}\_\{\\scriptscriptstyle\(B\)\}In2​Se2\(B,H\)\\mathrm\{In\_\{2\}Se\_\{2\}\}\_\{\\scriptscriptstyle\(B,H\)\}In2​Se3\\mathrm\{In\_\{2\}Se\_\{3\}\}IrTe21​T\\mathrm\{IrTe\_\{2\}^\{1T\}\}MoS22​H\(B,H\)\\mathrm\{MoS\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B,H\)\}MoSe22​H\(B,H\)\\mathrm\{MoSe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B,H\)\}MoTe22​H\(B,H\)\\mathrm\{MoTe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B,H\)\}NbS22​H\(B\)\\mathrm\{NbS\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B\)\}NbSe22​H\(B\)\\mathrm\{NbSe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B\)\}NbTe22​H\(B\)\\mathrm\{NbTe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B\)\}NiS21​T\(B\)\\mathrm\{NiS\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}NiSe21​T\(B\)\\mathrm\{NiSe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}NiTe21​T\(B\)\\mathrm\{NiTe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}P2\\mathrm\{P\_\{2\}\}P4\\mathrm\{P\_\{4\}\}Pd2​Se4\(B\)\\mathrm\{Pd\_\{2\}Se\_\{4\}\}\_\{\\scriptscriptstyle\(B\)\}ReSe22​H\\mathrm\{ReSe\_\{2\}^\{2H\}\}PdS21​T\(B\)\\mathrm\{PdS\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}PdSe21​T\(B\)\\mathrm\{PdSe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}PdTe21​T\(B\)\\mathrm\{PdTe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}Pd4​Se6​Hg2\\mathrm\{Pd\_\{4\}Se\_\{6\}Hg\_\{2\}\}PtS21​T\(B\)\\mathrm\{PtS\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}PtSe21​T\(B,H\)\\mathrm\{PtSe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B,H\)\}PtTe21​T\(B\)\\mathrm\{PtTe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}Pt4​Se6​Hg2\\mathrm\{Pt\_\{4\}Se\_\{6\}Hg\_\{2\}\}ReS22​H\\mathrm\{ReS\_\{2\}^\{2H\}\}Sb2\\mathrm\{Sb\_\{2\}\}Sb2​Te3\(B\)\\mathrm\{Sb\_\{2\}Te\_\{3\}\}\_\{\\scriptscriptstyle\(B\)\}Si2\\mathrm\{Si\_\{2\}\}SnS21​T\(B\)\\mathrm\{SnS\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}SnS\(B\)\\mathrm\{SnS\}\_\{\\scriptscriptstyle\(B\)\}SnSe21​T\(B\)\\mathrm\{SnSe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}TaS22​H\(B\)\\mathrm\{TaS\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B\)\}TaSe22​H\(B\)\\mathrm\{TaSe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B\)\}Ti2​S6\\mathrm\{Ti\_\{2\}S\_\{6\}\}TiS21​T\\mathrm\{TiS\_\{2\}^\{1T\}\}TiSe21​T\(B\)\\mathrm\{TiSe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}WS22​H\(B,H\)\\mathrm\{WS\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B,H\)\}WSe22​H\(B,H\)\\mathrm\{WSe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B,H\)\}WTe22​H\(B,H\)\\mathrm\{WTe\_\{2\}^\{2H\}\}\_\{\\scriptscriptstyle\(B,H\)\}Zn2​Se2\(B\)\\mathrm\{Zn\_\{2\}Se\_\{2\}\}\_\{\\scriptscriptstyle\(B\)\}Zr2​Te6\\mathrm\{Zr\_\{2\}Te\_\{6\}\}ZrS21​T\(B\)\\mathrm\{ZrS\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B\)\}ZrSe21​T\(B,H\)\\mathrm\{ZrSe\_\{2\}^\{1T\}\}\_\{\\scriptscriptstyle\(B,H\)\}

#### E\.3\.2\.HetDB\-to\-SAMBA

In this setting, the models are trained on the full HetDB dataset consisting of 336 heterobilayers and evaluated on the full SAMBA test set of 980 bilayers\. Among the monolayers forming the HetDB heterobilayers, 13 overlap with monolayers in SAMBA\. Both BDIP\-Net and PotNet are trained with a learning rate of2×10−42\\times 10^\{\-4\}and a batch size of 128, while their remaining model configurations follow those described in Section[4\.1](https://arxiv.org/html/2608.14640#S4.SS1)\. Similar to the previous setting, the SAMBA test set is divided according to the degree of monolayer overlap with the training data\. The no\-monolayer\-overlap, one\-monolayer\-overlap, and two\-monolayer\-overlap subsets contain 626, 284, and 70 bilayers, respectively\. The two\-monolayer\-overlap subset is further divided into 36 exact\-pair\-overlap bilayers and 34 bilayers with two\-monolayer overlap excluding exact\-pair overlap\. Table[8](https://arxiv.org/html/2608.14640#A5.T8)summarizes the prediction performance of BDIP\-Net and PotNet under this setting\.

Table 8\.HetDB\-to\-SAMBA domain generalization prediction performance of BDIP\-Net and PotNet\.SettingNumber ofBilayersModelMAE↓\\downarrowMSE↓\\downarrowRMSE↓\\downarrowR2↑\\uparrowFull SAMBA980BDIP\-Net0\.560\.500\.71\-0\.93PotNet0\.500\.390\.63\-0\.50No\-monolayer\-overlap626BDIP\-Net0\.650\.610\.78\-2\.71PotNet0\.570\.460\.68\-1\.77One\-monolayer\-overlap284BDIP\-Net0\.460\.380\.61\-0\.85PotNet0\.430\.320\.56\-0\.56Two\-monolayer\-overlap70BDIP\-Net0\.120\.030\.160\.97PotNet0\.150\.100\.310\.89Exact\-pair\-overlap36BDIP\-Net0\.100\.020\.140\.95PotNet0\.070\.010\.110\.97Two\-monolayer\-overlap excludingexact\-pair\-overlap34BDIP\-Net0\.140\.030\.180\.98PotNet0\.240\.190\.430\.86

![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/box_plot_samba_hetdb.png)Figure 7\.HetDB\-to\-SAMBA Domain Generalization Breakdown Analysis![Refer to caption](https://arxiv.org/html/2608.14640v1/figures/scatter_plot_monolayer_overlap_hetdb_samba_sorted.png)Figure 8\.HetDB\-to\-SAMBA Domain Generalization AE DistributionPrediction Performance\.As shown in Table[8](https://arxiv.org/html/2608.14640#A5.T8), PotNet achieves better prediction performance than BDIP\-Net on the full SAMBA test set of 980 bilayers, with an MAE of0\.500\.50, an MSE of0\.390\.39, an RMSE of0\.630\.63, and anR2R^\{2\}of−0\.50\-0\.50, compared with0\.560\.56,0\.500\.50,0\.710\.71, and−0\.93\-0\.93, respectively, for BDIP\-Net\. Figs\.[7](https://arxiv.org/html/2608.14640#A5.F7)and[8](https://arxiv.org/html/2608.14640#A5.F8)show that the prediction performance of both models improves as the degree of monolayer overlap increases\. On the 626 no\-monolayer\-overlap bilayers, PotNet achieves a lower MAE than BDIP\-Net \(0\.570\.57vs\.0\.650\.65\)\. On the 284 one\-monolayer\-overlap bilayers, the corresponding MAEs are0\.430\.43and0\.460\.46\. In contrast, BDIP\-Net performs better on the 70 two\-monolayer\-overlap bilayers, achieving an MAE of0\.120\.12, compared with0\.150\.15for PotNet\. The best results are obtained on the 36 exact\-pair\-overlap bilayers, where PotNet achieves a slightly lower MAE \(0\.070\.07vs\.0\.100\.10\)\. However, after removing these exact\-pair\-overlap bilayers, BDIP\-Net achieves an MAE of0\.140\.14, an RMSE of0\.180\.18, and anR2R^\{2\}of0\.980\.98on the remaining 34 two\-monolayer\-overlap bilayers, compared with0\.240\.24,0\.430\.43, and0\.860\.86, respectively, for PotNet\. These results indicate that although PotNet performs better overall, BDIP\-Net provides stronger cross\-dataset generalization to unseen bilayer combinations of known monolayers\. The overall HetDB\-to\-SAMBA performance is nevertheless limited because most SAMBA bilayers belong to the no\-monolayer\-overlap and one\-monolayer\-overlap subsets\.