Physics-informed Hamiltonian learning for large-scale optoelectronic property prediction

Domain: biology
Published: 2026-03-20
Editorial outcome: minor_revision

Abstract
Predicting optoelectronic properties of large-scale atomistic systems under realistic conditions is crucial for rational materials design, yet computationally prohibitive with first-principles simulations. Recent neural network models have shown promise in overcoming these challenges, but typically require large datasets and lack physical interpretability. Physics-inspired approximate models offer greater data efficiency and intuitive understanding, but often sacrifice accuracy and transferability. Here we present HAMSTER, a physics-informed machine learning framework for predicting the quantum-mechanical Hamiltonian of complex chemical systems. Starting from an approximate model encoding essential physical effects, HAMSTER captures the critical influence of dynamic environments on Hamiltonians using only few explicit first-principles calculations. We demonstrate our approach on halide perovskites, achieving accurate prediction of optoelectronic properties across temperature and compositional variations, and scalability to systems containing tens of thousands of atoms. This work highlights the power of physics-informed Hamiltonian learning for accurate and interpretable optoelectronic property prediction in large, complex systems.

Introduction
Quantum-mechanical predictions of large-scale chemical systems are of key importance to accelerating progress in diverse fields ranging from quantum biology1 and materials science2 to drug design3. In optoelectronic materials, which are key components in photovoltaic and lightning applications, computations of large-scale systems are essential to capture real-world factors such as the presence of defects and dynamic disorder that strongly influence device behavior through charge transport and recombination processes. Ideally, one would predict the electronic structure of such large and disordered systems at finite temperature directly from first-principles. Achieving this would move quantum-mechanical modeling beyond the prevailing paradigm, which assumes an ideal lattice at 0 K, Bloch-wave electronic states, and perturbative treatments of critical effects such as electron-phonon or electron-defect interactions. Computational capabilities that go beyond such frameworks and accurately capture thermal effects on optoelectronic properties represent a crucial step towards realistic materials modeling under operating conditions, ultimately advancing the design of next-generation optoelectronic technologies.

Progress in machine learning force fields (MLFFs) suggests promising routes to overcome steep limitations in the modeling of large-scale molecules and materials. In particular, MLFFs enable first-principles accuracy of molecular dynamics (MD) simulations at dramatically reduced computational cost4–14. While MLFFs can accurately predict atomic trajectories at specified thermodynamic conditions, they typically lack the ability to provide insights into electronic structure, an increasingly important demand given the recent availability of highly accurate MD trajectories of large-scale quantum systems. Relying on explicit first-principles electronic-structure theory, such as density functional theory (DFT), remains impractical because of its unfavorable scaling of computational cost with increasing system size.

Recent work15,16 suggests that equivariant machine learning (ML) provides a powerful framework to sidestep explicit first-principles calculations of large-scale chemical systems by learning and predicting the quantum-mechanical Hamiltonian. Such Hamiltonian learning frameworks require training an ML model with first-principles data, from which one can then predict the Hamiltonian matrix using atomic coordinates as input to the model. Omitting explicit first-principles calculations, one can then, in principle, calculate many quantum-mechanical observables of large-scale chemical systems.

Relevant in this context, learning the Hamiltonian is challenging because it depends not only on a system’s atomic structure but also on the chosen representation, so there is no single unique mapping from structure to matrix elements. However, equivariant ML models respecting fundamental physical transformations, such as rotations and translations, are a natural choice for Hamiltonian learning. Equivariant neural networks (ENNs) were recently developed as particularly promising frameworks for Hamiltonian learning17–19. Although ENNs improve data efficiency by encoding symmetry and therefore require fewer samples than comparable non-equivariant models, they may still demand substantial amounts of training data to capture complex dynamical behavior, such as that exhibited in multi-component materials at finite temperature. Related to this limitation, there exists a lack of consistency in how Hamiltonian learning models are represented in a consistent basis, which restricts their applicability to specific DFT frameworks requiring careful benchmarking20.

To alleviate some of these challenges, Gu et al. recently combined a tight-binding (TB) model with deep learning to predict Hamiltonians, demonstrating accurate electronic-property predictions for large-scale simulations using DFT eigenvalues as training data21. Although this suggests that incorporating approximate physical models into the ML workflow can substantially improve prediction accuracy, there is still little understanding of how to best integrate physics-aware approaches with emerging ML architectures or hybrid multiscale simulations for more complex material classes. Thus, there is a need for new approaches seeking to combine the efficiency of physics-based approximations with the flexibility of ML while maintaining transferability across material classes.

In this work, we address these critical limitations by developing a physics-informed ML framework for predicting the quantum-mechanical Hamiltonian of large-scale chemical systems. Our approach is physically motivated and requires little training data, yet it achieves high predictive accuracy and transferability. We augment the physical model with an ML scheme that captures the key effect of the dynamic environment on the electronic Hamiltonian. We find that our integrated model, HAMSTER (Hamiltonian-learning Approach for Multiscale Simulations using a Transferable and Efficient Representation), achieves first-principles accuracy with only a small number of explicit first-principles calculations, offering substantially greater data efficiency than existing frameworks. We assess the performance of HAMSTER on halide perovskite semiconductors, a challenging case due to their soft lattices, dynamic disorder, and chemical complexity, well beyond the systems typically considered in prior work. We show that the model yields accurate electronic-property predictions, consistent with first-principles and experimental results, across diverse compositions and temperatures for large-scale systems containing tens of thousands of atoms.

Results
The objective of Hamiltonian learning is to predict the Hamiltonian matrix using the atomic coordinates as input. For chemical systems containing many atoms, the Hamiltonian is a large matrix that encodes rich quantum-mechanical interactions. Predicting it from scratch with ML is a challenging task that requires fitting many parameters, which is computationally expensive. Our goal is to enhance the data efficiency of Hamiltonian learning by incorporating underlying physical principles into the model prior to the ML procedure. To achieve it, we suggest a \Delta ML approach, where the ML model learns differences between results from a physical model and the ground truth instead of initializing the fitting procedure from scratch (see Fig. 1a).

In our HAMSTER approach, we therefore start from the well-known TB model that requires a small amount of first-principles data as input and provides an approximate yet physically motivated description of the Hamiltonian. Specifically, assuming periodic boundary conditions, we write the matrix elements of the TB Hamiltonian in k-space as

H_{ij} = \sum_{R} e^{ik \cdot R} \langle \chi_i(r) | H_{eff} | \chi_j(r - R) \rangle := H_{ij}^R, (1)

where \chi_i is the i-th electron orbital, R a real-space lattice vector, and H_{eff} an effective independent-particle Hamiltonian. The indices i and j are composite indices that specify both the orbital and the atomic site. To identify different types of interactions, it is useful to write the effective Hamiltonian as

H_{eff} = -\frac{1}{2} \nabla^2 + \sum_{k=1}^N V_k(r - r_k) + H_{soc}, (2)

where H_{soc} encapsulates spin-orbit coupling (SOC) effects (see “Methods” section), and the potential is expressed in terms of contributions by atoms located at sites k22.

Using Eq. (2) in Eq. (1), we can identify different types of matrix elements: (i) i = j corresponds to on-site matrix elements; (ii) i = j \neq k to environment effects for on-site elements; (iii) i \neq j = k or i = k \neq j corresponds to off-site matrix elements; (iv) i \neq j \neq k to environment effects for off-site elements; and (v) SOC matrix elements. Contributions due to (ii) and (iv) are often neglected in common TB models, which is known as the two-center approximation23–25. Using this approximation, in our TB model26, we compute the two-center contributions as

H_{ij}^R = \sum_{v} C_{ij}^v V_v(\Delta r) V_v^0, (3)

where C_{ij}^v takes into account the relative orientation in the v-th orbital, V_v(\Delta r) accounts for the distance between the orbitals, and V_v^0 are parameters.

At finite temperatures, atoms in molecules and materials are constantly in motion, which leads to fluctuations in their local environments. Starting from our TB model, we therefore augment the Hamiltonian in a \Delta ML step27 (see Fig. 1a) to capture the effects of each atom’s dynamic environment via interaction terms (ii) and (iv). To this end, each non-zero matrix element in the Hamiltonian of Eq. (3), determined by a cutoff radius r_{cut}, is associated with a descriptor vector, x_{ij}^R, whose specific construction is detailed below. From the full set of descriptor vectors, we select a subset of size N_p using the k-means clustering algorithm. These serve as kernel support points and are used to predict the corrections for remaining or unseen matrix elements (see “Methods” for details on the sampling procedure).

We employ a radial basis function kernel model, where each kernel support point, x_n, is associated with a learnable parameter, h_n. Given a trial environment descriptor vector, x_{ij}^R, the ML correction for the respective matrix element in Eq. (3) is calculated as

\delta H_{ij}^R = \sum_{n=1}^{N_p} h_n \exp \left( - \frac{||x_n - x_{ij}^R||^2}{2\sigma^2} \right), (4)

where the kernel width, \sigma, acts as a hyperparameter that has to be specified before optimization, as is the case with N_p. With H_{eff} defined, one may proceed with kernel ridge regression, using first-principles Hamiltonian matrix elements as ground-truth values. Instead, we propose using energy eigenvalues obtained from first-principles calculations as training targets and optimize the model parameters via gradient descent (see “Methods” section for more details).

A key remaining part of the model is selecting descriptor vectors for Hamiltonian matrix elements. Various descriptors were proposed in prior studies16,28–30. A significant advantage is the choice of a physical model for the Hamiltonian, which already accounts for changes in the distance and relative orientation between orbitals i and j (see Eq. (3)). Then, the ML correction only needs to predict finer dynamic modifications of matrix elements via Eq. (4). To be physically meaningful, the descriptors should vary smoothly with changes in atomic positions and respect the symmetries of the Hamiltonian. Accordingly, we adopt a relatively simple descriptor that encodes both the interaction type and surrounding atomic environment (see Fig. 1b): we approximate the environment of H_{ij} as a combination of local environments around atoms i and j, treated independently. The local environment of each atom i is computed as

h_{env,i} = \sum_{k \neq i} f_{cut}(|r_i - r_k|) \int d^3r \chi_i(r - r_i) \chi_k(r - r_k), (5)

where the sum runs over all atoms and orbitals, indexed by k, within a cutoff radius, r_{cut}. The smooth cutoff function, f_{cut}(r) = \frac{1}{2}[\cos(\pi r) + 1], ensures a gradual decay of contributions from distant atoms. The local environment is identical for all lattice translation vectors, R, as a consequence of translational symmetry. Note that the integrals in Eq. (5) reduce to those in Eq. (3) when V_v = 1, and therefore need to be evaluated only once for both the TB and ML contributions.

To provide the model with information about the interaction type, we augment the descriptor vector with several structural features: the equilibrium distance between atoms i and j, the angular orientation of their respective orbitals relative to the connecting vector, r_{ij}, and the atomic numbers Z_i and Z_j. These quantities can either be kept fixed during atomic motion or updated dynamically, thereby providing the model with the flexibility to accommodate additional physical effects such as thermal lattice expansion. This design ensures that the ML model focuses on learning changes in the local environment, thereby improving data efficiency by reducing the complexity of the learning task. The descriptor vector then reads

x_{ij}^R = [Z_i, Z_j, \Delta r_{ij}^R, \theta_i^R, \theta_j^R, \phi_{ij}^R, h_{env,i}, h_{env,j}]^T. (6)

Here, \Delta r_{ij}^R is the distance between atoms i and j with the translation vector R applied to atom j, \theta_x^R (x = i, j) is the angle between the orientation of the orbital on atom x and the bonding axis, and \phi_{ij}^R is the angle between the orbital orientations on atoms i and j.

We apply the HAMSTER approach to bulk GaAs as a test case, starting from the above-described TB model (see ref. 26 for full details). Figure 1c shows that training the \Delta ML method requires only a few explicit first-principles calculations of MD-generated structures at a temperature of 400 K for the validation error, calculated as mean absolute error across the eigenvalues (MAE, see “Methods” section), to converge. Figure 1d shows that systematic deviations in the TB eigenvalue spectrum at 400 K are markedly reduced when adding the \Delta ML correction, yielding sub-50 meV accuracy for many of the eigenvalues in comparison to DFT. These results suggest that the primary reason for the inaccuracies of the TB model stems from its omission of dynamic environment effects, which are effectively addressed through the addition of the ML correction in a data-efficient way.

Having established the effectiveness of the HAMSTER approach for a relatively simple test case, we now turn to the significantly more complex halide perovskite semiconductors. These materials have attracted significant attention due to their promising optoelectronic properties and applications in photovoltaics and light-emitting devices32,33. At the same time, they exhibit complex vibrational properties, characterized by anharmonic, overdamped phonons and unusual electron-phonon coupling34–37, which are generally challenging to capture in atomistic models38,39. Furthermore, the halide p-orbitals in the ABX3 perovskite crystal structure are embedded in two distinctly different environments, as these orbitals are oriented either along a bonding axis or perpendicular to it. These aspects render the computational modeling of perovskite optoelectronic properties at finite temperature particularly challenging. We consider cesium lead bromide (CsPbBr3) in its cubic phase as a prototypical inorganic halide perovskite material and investigate the performance of HAMSTER.

We first construct a TB model for CsPbBr3 as a physically meaningful approximation of the Hamiltonian. The model includes bromine p-orbitals, leads s- and p-orbitals, and cesiums s-orbitals in the basis, as these dominantly contribute to the electronic structure close to the valence band maximum (VBM) and conduction band minimum (CBM). Furthermore, we explicitly account for SOC in the TB Hamiltonian, given that it plays a crucial role in the electronic structure of halide perovskites. Further details regarding the model and its optimization are described in the “Methods” section.

To determine an appropriate training set size for \Delta ML, we train our model on an increasing number of structures drawn from an MD trajectory at a temperature of 425 K. Sampling an equal number of kernel support points from each structure, we compare the use of r_{cut} = 5.5 \AA and r_{cut} = 6.2 \AA (see Eq. (5)) in two separate optimizations to assess the impact of the interaction range on model performance.

The learning curves (see Fig. 2a) show that the validation error plateaus after using around 10 training structures, reaching an accuracy for the eigenvalues that drops below 50 meV. While this level of accuracy is comparable to that reported in a recent study21 using a neural-network scheme for elemental and binary semiconductor compounds, we emphasize that our approach needs far fewer training structures for the considerably more complex halide perovskite CsPbBr3, well below the few hundred typically required by neural-network-based Hamiltonian learning models17,21,31. Remarkably, for this material, HAMSTER still generalizes reasonably well with as few as two training structures, remaining within 0.03 eV of the training error, highlighting the markedly superior data efficiency of our approach. The mild overfitting at very low training set sizes, notable through the small gap between training and validation error, disappears as the training set increases.

Furthermore, we observe that the prediction errors decrease with increasing cutoff radius (see Fig. 2a). This is expected, as a larger cutoff captures more interactions, increasing the expressiveness of the model by introducing additional non-zero Hamiltonian elements. However, this improvement in accuracy comes at the cost of increased computational effort, as the Hamiltonian becomes denser as r_{cut} increases. Thus, selecting the cutoff radius requires balancing accuracy against computational efficiency.

We analyze the structure of the descriptor space to gain insight into the model’s increased expressiveness when a larger r_{cut} is used. To this end, we perform a principal component analysis (PCA, see “Methods”) on the sampled descriptors for both cutoff radii, which projects the eight-dimensional descriptor vectors onto a two-dimensional subspace while preserving as much of their variance as possible. This enables a clear visualization of their distribution (see Fig. 2b). In this representation, descriptor vectors that encode types of interactions lie close together. Since each descriptor component corresponds to a specific orbital interaction in the Hamiltonian (see Eq. (6)), the PCA projection groups and separates regions of descriptor space according to the orbital couplings they represent.

For each interaction type, we observe multiple distinct clusters, which can be associated with specific, physically-interpretable orbital coupling types, such as sp\sigma or pp\pi. When increasing the cutoff radius, new clusters emerge that correspond to additional interactions between atoms in neighboring unit cells, including Pb-Pb, Cs-Cs, and Br-Br couplings. This not only demonstrates the ability of the ML model to distinguish between interaction types but also its capability to learn physically-interpretable corrections according to additional interactions, which further enhances accuracy.

Having trained and benchmarked HAMSTER at its reference temperature, we now examine its transferability across temperatures to assess robustness at varying thermal conditions. Thermal variations of optoelectronic properties in halide perovskites are notoriously difficult to capture computationally. A prime example is the band gap, whose anomalously weak temperature dependence defies the trends seen in conventional semiconductors and poses a significant challenge for predictive modeling40–42.

Before turning to large-scale simulations, we apply the model to MD trajectories of a 2 \times 2 \times 2 supercell of CsPbBr3 for benchmarking against DFT data, using r_{cut} = 5.5 \AA in our model. Figure 3a compares the temperature dependence of the band gap in CsPbBr3 as obtained with HAMSTER and DFT. Although the model was trained solely on structures from an MD trajectory at 425 K, its predictions remain in close agreement with DFT also at higher temperatures. Across all three considered temperatures, the difference between the DFT-computed band gap and the value predicted by HAMSTER remains below 50 meV, and the MAEs of energy eigenvalues are 0.047 eV (425 K), 0.049 eV (525 K), and 0.057 eV (625 K). Furthermore, we find that the fluctuations of instantaneous band gaps, gauged by the standard deviation at each temperature, are similar in DFT and HAMSTER. These results underscore the transferability of our model over a considerable temperature range within the cubic phase of CsPbBr3.

Figure 3b presents the average residual for each band obtained in the pristine TB model (left panel) in comparison to HAMSTER (right panel) at a temperature of 425 K. Without accounting for the effect of the dynamic environment, the energy of certain bands appears systematically over- or underestimated, indicating inherent limitations of the TB model in capturing band-specific dynamic fluctuations in the electronic structure. Incorporating the environment-dependent interactions via the \Delta ML step, these systematic deviations are significantly reduced, with residuals now appearing much less amplified and more randomly distributed around zero. A notable exception remains for the uppermost conduction bands, where a consistent mismatch persists. This can be attributed to their reduced weighting during optimization (see “Methods” section). Overall, these results demonstrate the strength of the physics-informed Hamiltonian learning framework, which combines the ability of TB models to capture the essential features of the electronic structure with the flexibility of ML to incorporate subtle but crucial environmental and multi-center effects.

A key advantage of our approach is that it enables simulations of system sizes far beyond the practical reach of explicit first-principles methods such as DFT. Modeling of larger systems can reveal distinct dynamical behavior and, by reducing finite-size effects, provide a more realistic description of the material. To explore this regime, we train an MLFF for CsPbBr3 (see “Methods”) and perform MD simulations for a 16 \times 16 \times 16 supercell of CsPbBr3 consisting of 20,480 atoms. Using the HAMSTER model, we subsequently predict Hamiltonians and compute band gaps by extracting eigenvalues at the VBM and CBM via the Lanczos algorithm as implemented in ARPACK43. Figure 3a shows that the mean value of the band gap at 425 K closely matches the one obtained for the smaller supercell. At the same time, the sizable band-gap fluctuations seen for the smaller supercells are almost entirely averaged out, yielding a much narrower distribution and a more realistic description of the material at finite temperature. This improved statistical convergence allows us to assess the temperature dependence of the band gap. While the absolute value of the band gap is affected by the known inaccuracies of PBE-level DFT, the resulting temperature dependence can still be meaningfully compared to experiment. For the 16 \times 16 \times 16 supercell, we obtain a slope of 3.95 \times 10^{-4} eV K^{-1}, in good agreement with the experimental value of 3.41 \times 10^{-4} eV K^{-1} reported in ref. 44. We note, however, that thermal expansion is not included in the simulation and that the temperature range considered here extends beyond that of the experimental study, both of which are expected to affect the band gap values.

To demonstrate the transferability of our workflow to other perovskites, we apply the HAMSTER model to MAPbBr3. This choice is motivated by a recent study showing that in the related perovskite FAPbI3, the supercell size was particularly important for reaching accurate band gap predictions, primarily due to molecular asymmetry, requiring simulation cells with at least 6000 atoms45.

We optimize the model parameters following the same procedure as for CsPbBr3 (see “Methods” section), again using r_{cut} = 5.5 \AA. However, for this material, we do not include orbitals stemming from the cation and perform classical force-field MD calculations (see “Methods” section). When applied to an independent test set consisting of 100 MD snapshots of a 2 \times 2 \times 2 supercell at 300 K, the model reaches an MAE of 0.060 eV relative to DFT. This demonstrates that the model maintains a level of accuracy comparable to that achieved for CsPbBr3. Figure 4a compares the band gaps of MAPbBr3 obtained from DFT and HAMSTER across different supercell sizes at 300 K. Leveraging the computational efficiency of HAMSTER, we extend our simulations to 16 \times 16 \times 16 supercells comprising almost 50,000 atoms. The predicted band gap of MAPbBr3 decreases by approximately 0.15 eV as the supercell size increases, in line with a recent report on FAPbI345. Agreement between HAMSTER and DFT is observed up to 4 \times 4 \times 4 supercells (768 atoms), beyond which direct DFT calculations become impractical.

Figure 4 presents the variation of the band gap with temperature as predicted by HAMSTER across different supercell sizes, together with experimental data for comparison. The objective here is not to reproduce the experimental trend with exact quantitative accuracy, but to employ it as a benchmark for assessing the reliability of the model. Although minor deviations between results obtained with HAMSTER and experimental data remain, the predicted temperature dependence captures the experimental behavior well, with the agreement becoming more favorable at larger supercell sizes. This consistency between experiment and computation supports the robustness and overall validity of the approach.

To assess the computational performance of HAMSTER, we benchmark the total runtime for three calculation modes: (i) full diagonalization of the Hamiltonian, (ii) Lanczos-based partial diagonalization targeting six eigenvalues near the VBM, and (iii) construction and output of the Hamiltonian matrix in Eq. (1) without solving the eigenvalue problem. For each benchmark, the model is evaluated for 10 independent structures, using the same configurations as in Fig. 4, and all calculations are performed on a single compute node using two Intel Xeon Platinum 8168 CPUs. Full diagonalization of the largest system size (16 \times 16 \times 16) is not tractable on the chosen, rather modest hardware and is therefore omitted. The resulting runtimes are shown in Fig. 4c. Since the benchmarks include the first invocation of each function, we fit a function, T(x) = T_0 + ax^n, to the timings and subtract the constant time T_0 = 44 s, i.e., the initial Julia compilation overhead, which disproportionately affects benchmarks with short execution times.

The observed scaling behavior is consistent with theoretical expectations. The construction of the Hamiltonian exhibits nearly ideal linear scaling with system size. Partial diagonalization remains significantly faster than full diagonalization, but deviates from linear behavior, showing an effective scaling exponent of approximately 1.8. This deviation originates from the shift-invert transformation, which can reduce the sparsity of the Hamiltonian. In contrast, full diagonalization yields a scaling exponent of 2.8, close to the cubic scaling expected for conventional DFT implementations. It is worth noting that the absolute computational cost of DFT calculations is substantially higher, primarily due to the overhead associated with repeated self-consistent field iterations and the much larger underlying basis sets. Additional speedups for diagonalization could be achieved through more advanced algorithms, highly parallel implementations, or GPU-accelerated solvers.

Within our approach, the wall time for the largest system considered remains below one hour, demonstrating that HAMSTER can efficiently handle large-scale systems even on modest hardware. When only the Hamiltonian is required, and eigenvalues are not computed, genuinely linear scaling can be achieved. Moreover, many important physical quantities, such as transport properties or the density of states, can be evaluated with linear-scaling methods that exploit the sparsity of the Hamiltonian46,47. Combining our Hamiltonian learning framework with such techniques, therefore, provides a promising route toward efficient simulations of large-scale systems.

Methods
Dataset of electronic-structure calculations
All DFT calculations were performed using the VASP code48 employing the PBE functional49 and accounting for SOC in the case of the halide perovskites. Specifically, for GaAs, we used data from DFT calculations as reported previously26 to obtain the TB model. We randomly selected MD snapshots at a temperature of 400 K to calculate the electronic structure of a 2 \times 2 \times 2 supercell of the primitive structure, comprising 64 atoms, with a 3 \times 3 \times 3 k-point grid, for training and validating the ML model with 10 and 48 calculations, respectively. The DFT electronic structure calculations for comparison with HAMSTER were performed for a 4 \times 4 \times 4 supercell of the diamond structure, comprising 128 atoms, using a 3 \times 3 \times 3 k-point grid, and randomly selecting 100 MD snapshots at a temperature of 400 K.

For CsPbBr3, we set the plane-wave kinetic energy cutoff to 300 eV, and sampled the electronic structure of the unit cell on a 12 \times 12 \times 12 k-point grid to generate data for the TB model. We selected 24 random snapshots from an MD trajectory at a temperature of 425 K to calculate the electronic structure of a 2 \times 2 \times 2 supercell using a 3 \times 3 \times 3 k-point grid for training and validating the ML model with 12 calculations each. For the DFT band gap calculations, we selected 100 random snapshots from MD runs at each temperature and performed self-consistent calculations for all of them using a 3 \times 3 \times 3 k-point grid.

For MAPbBr3, we used a cut-off energy of 400 eV and sampled the electronic structure of the unit cell on an 8 \times 8 \times 8 k-point grid to obtain the data for the TB model. We selected 24 random snapshots from an MD trajectory of a 2 \times 2 \times 2 supercell at a temperature of 300 K to calculate the electronic structure, using a 2 \times 2 \times 2 k-point grid, for training and validating the ML model with 12 calculations each. For the DFT band gap calculations, we selected 100 random snapshots from MD runs at a temperature of 300 K and performed self-consistent calculations for all of them using a 3 \times 3 \times 3 k-point grid. The DFT electronic structure results for the 4 \times 4 \times 4 supercell were obtained with a sampling at the \Gamma-point only, for a total of 100 MD snapshots, at a temperature of 300 K.

For the calculations in HAMSTER, we used the same snapshots as in the DFT calculations when comparing the two. For the 8 \times 8 \times 8 and 16 \times 16 \times 16 cells, we chose a sample size of 20 snapshots each, as we found that the electronic-structure calculations showed a faster statistical convergence with respect to sample size compared to the smaller supercells.

Molecular dynamics calculations
All DFT-based MD calculations were performed in VASP48 using an NVT ensemble. Specifically, for GaAs, we employed a 4 \times 4 \times 4 supercell, used the \Gamma-point, and a timestep of 8 fs to perform a DFT-MD production run containing 4000 steps after equilibration at a temperature of 400 K.

For CsPbBr3, we employed a 2 \times 2 \times 2 supercell, 3 \times 3 \times 3 k-point grid, Tkatchenko-Scheffler van der Waals corrections50, and a timestep of 8 fs to perform a DFT-MD production run containing 2000 steps at each temperature.

For large-scale MD calculations of CsPbBr3, we first trained three MLFFs in VASP5,51,52, at 425 K, 525 K, and 625 K, respectively, using the Bayesian on-the-fly scheme and a timestep of 6 fs for a 4 \times 4 \times 4 supercell of the cubic structure. In order to ensure robustness and convergence of the resulting MLFFs, the training runs lasted 100 ps. During training, the default set of hyperparameters was first employed, specifying the following tags in VASP: ML_RCUT1=8 \AA, ML_RCUT2=5 \AA, ML_SION1=0.5, ML_SION2=0.5. Using the DFT data collected during on-the-fly training, we trained three MACE models9, one per temperature, using 2 message passing layers, messages with 128 channels and max_L of 2, correlation order of 3, spherical harmonics up to l=3, and a radial cutoff of 10 \AA. The weighted mean squared error loss function was implemented with initial energy and force weights being 1 and 100, respectively, increasing to 1000 and 100 in the last 20% of the epochs, with stochastic weight averaging enabled. The MACE models were trained with a batch size of 5 for a total of 600 epochs. The validation sets were chosen to be 10% of the configurations in each dataset. The reference energies for isolated atoms were obtained through spin-polarized DFT calculations. The resulting MACE force fields were used with the LAMMPS code53, and applied to a 16 \times 16 \times 16 supercell of the cubic structure for 5000 time steps of 0.6 fs.

For MAPbBr3, we applied a classical force field from the literature54 using the LAMMPS code53. We equilibrated the system, including its lattice, in successive NVT and NpT runs lasting approximately 22 ps in total at each temperature, before running 5 ps long production runs from which we extracted one snapshot every 50 fs. We repeated this procedure for the different supercell sizes, namely the 2 \times 2 \times 2, 4 \times 4 \times 4, 8 \times 8 \times 8, and 16 \times 16 \times 16 supercell of cubic MAPbBr3.

Contribution of spin-orbit coupling in tight binding model
One can compute the matrix elements of the SOC Hamiltonian as

H_{soc,ij} = \gamma \langle \chi_i | bL \cdot bS | \chi_j \rangle, (7)

where \gamma is a learnable parameter that depends on the atomic species, and bL and bS denote the orbital and spin-angular momentum operators, respectively. In the case of a basis of AOs given in spherical harmonics, we can analytically evaluate Eq. (7) (see, e.g., ref. 55). For systems described using hybrid orbitals (HOs), as employed in our previous work on GaAs26, the SOC Hamiltonian in the HO basis is computed via the transformation

H_{soc}^{ho} = F^y H_{soc} F, (8)

where the matrix, F, contains the coefficients that expand the HOs in terms of AOs (see Tab. I in ref. 26). We note that this model only includes on-site SOC interactions; off-site contributions and spin-flip terms between different atomic centers are neglected.

Principal component analysis
To gain deeper insight into the descriptor space, we performed a PCA on the samples selected by the ML model. Prior to PCA, each descriptor was transformed to the range [0, 1], and centered to ensure comparability. For calculations with an increasing number of training structures, new structures were added incrementally to the set while retaining all previously included ones.

Sampling of kernel support points
From the entire set of descriptor vectors for all training structures, a subset of N_p kernel support points is sampled using k-means clustering with N_{cluster} clusters. To account for the varying importance of each cluster, we weighted the number of points sampled from the i-th cluster by

w_i = \alpha s_i + (1 - \alpha) v_i, (9)

where s_i, v_i are the normalized size and variance of the i-th cluster, respectively, and \alpha = 0.5 for all our calculations. At least one point was sampled from each cluster to ensure that each type of interaction was sampled at least once.

Parameter optimization
We perform separate training runs to optimize each set of model parameters, thereby avoiding numerical instability. Specifically, we optimize V_v for the TB model (Eq. (3)), h_n for the ML model (Eq. (4)), and the \gamma parameters for the SOC model (Eq. (7)) using gradient descent with the ADAM optimizer56. For the TB parameters, learning rates between 0.1 and 0.3 were found to work well, and all TB parameters were initialized to 1. The ML corrections are expected to be small, so a smaller learning rate of 0.01 or 0.005 is used, with all ML parameters initialized to 0. SOC parameters are also initialized to 1. Since the SOC model has only a few parameters, it is largely insensitive to the learning rate and converges within just a few iterations.

We optimize the TB model using the non-distorted primitive unit cell and using 12 randomly selected snapshots from MD as validation data to prevent overfitting; the ML and SOC models are optimized using 12 randomly selected snapshots from an MD run at a single temperature, as described above.

The TB model is optimized for up to 1000 training steps, with early stopping to prevent overfitting. For the ML and SOC models, we first run a combined training of 200 optimization steps, during which both ML and SOC parameters are fitted simultaneously. Next, the ML parameters are reinitialized to 0 and refitted for an additional 200 steps while keeping the SOC parameters fixed. This two-stage procedure is practical because the ML and SOC models strongly influence the band gap and are interdependent.

At each iteration, the gradient of the loss function with respect to each parameter is required. Using the chain rule, this can be expressed as

\frac{\partial L}{\partial p} = \sum_{n,k,i,j,R} \frac{\partial L}{\partial E_{nk}} \frac{\partial E_{nk}}{\partial H_{ij}^R} \frac{\partial H_{ij}^R}{\partial p}, (10)

where p denotes a model parameter and R are the lattice translation vectors. When computing the loss, increased weight was assigned to the \Gamma-point and to the highest valence and lowest conduction bands, while reducing the weight of the two uppermost bands. The term on the right of Eq. (10), \frac{\partial H_{ij}^R}{\partial p}, is straightforward to compute since the model parameters enter the Hamiltonian matrix elements linearly. The term, \frac{\partial E_{nk}}{\partial H_{ij}^R}, can be evaluated using the Hellmann-Feynman theorem and the eigenvectors obtained from the diagonalization of the Hamiltonian. Using Eq. (1), we can identify the derivative, \frac{\partial H_{ij}^R}{\partial p}, simply as e^{ik \cdot R}. Finally, the first term, \frac{\partial L}{\partial E_{nk}}, depends on the chosen loss function. In our work, we adopted a weighted MAE to quantify the discrepancy between predicted and target energy eigenvalues:

L = \frac{1}{\sum_n w_n \sum_k w_k} \sum_{n,k} w_n w_k |E_{nk}^{target} - E_{nk}^{pred}|, (11)

where w_n, w_k are the weights of the n-th band and the k-point k, respectively, and, \sum_n w_n, \sum_k w_k are the sum of the respective weights. To further accelerate the gradient computation, which is the main bottleneck during optimization, we explicitly made use of the sparsity of the Hamiltonian to only compute elements of the sum over i, j in Eq. (10) for which the respective matrix element H_{ij}^R is non-zero.

Nevertheless, our optimization strategy has some limitations. First, the TB, ML, and SOC parameters are trained in separate stages to avoid numerical instability, which prevents fully capturing their mutual interdependence. Consequently, subtle correlations between these components may not be entirely reflected in the final model. Second, gradient-based optimization may be less effective if the Hamiltonian exhibits strong nonlinear dependence on its parameters, and it cannot guarantee convergence to the global minimum. Despite these limitations, our approach provides a stable and computationally efficient framework for fitting effective Hamiltonians.

Hyperparameter optimization
For parameters that cannot be efficiently optimized using gradient descent, we perform hyperparameter optimization using a Tree-structured Parzen Estimator (TPE)57. In each hyperparameter iteration, the TPE proposes a new candidate using a Bayesian estimator. For each candidate, a full parameter optimization is performed, with the model parameters reinitialized at the start of the run. Once the optimal hyperparameter setting is identified, a final optimization is carried out to obtain the corresponding optimal model parameters.

For the TB model, the effective core charge parameter (see Sec. C in ref. 26) is treated as a hyperparameter, since computing its gradient is impractical. For the ML model, the hyperparameters include the number of kernel support points N_p, the number of clusters N_{cluster}, and the kernel width \sigma. As described in the main text, the structural entries in the descriptor can either be fixed or dynamically updated based on lattice distortions. For CsPbBr3 and GaAs, the structural entries are kept fixed. For MAPbBr3, they are updated dynamically to account for thermal lattice changes. In this case, the dynamic entries are normalized by the cutoff radius or the maximum angle 2\pi, and the environment defined in Eq. (5) is scaled by a fixed factor, env_scale, which is treated as an additional hyperparameter. The final ML models comprise 250 parameters for GaAs, 910 for CsPbBr3, and 420 for MAPbBr3.

Band window selection and weighting
As described in detail above, the training and validation datasets consisted of energy eigenvalues obtained from DFT. Here, the choice of the energy window containing the bands used in the optimization must be made carefully, as it depends on the pseudopotential used in the DFT calculations, and the number of Kohn-Sham orbitals is generally large. For CsPbBr3 and MAPbBr3, we included the valence bands in the TB basis, together with additional conduction bands to account for contributions of the remaining orbitals (five for CsPbBr3, four for MAPbBr3). When computing the loss, increased weight was assigned to the \Gamma-point and to the highest valence and lowest conduction bands, while reducing the weight of the two uppermost bands. These numbers were doubled when SOC is included in the calculation and further scaled with the size of the supercell. The reason for this is that these bands typically exhibit larger errors, which can disproportionately influence the fit as the model attempts to minimize their contributions. Moreover, they may include components from higher-lying conduction bands due to band entanglement, making them difficult for the model to predict with sufficient accuracy. Since these states lie far from the band gap region, they are also less critical for most target applications, further justifying their reduced influence during training.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript aims to contribute to the prediction of optoelectronic properties in large-scale materials. It develops a machine learning model based on tight-binding density functional theory, with the parameters of an effective Hamiltonian fitted using a kernel regression approach. The manuscript claims improved efficiency in computing eigenvalues.\nHowever, the text is difficult to follow for non-specialist readers. Several equations or derivations are omitted and need further discussion or proper citation. For example, Equation 4 (the descriptor x_i) is unclear, and its relationship to Equation 5 is not explained. The optimization and fitting of parameters and hyperparameters are confusing. Authors should discuss the limitations of these method.\nSeveral figure captions are incomplete; for instance, Figure 2a does not clearly state which property is studied. Clusters in Fig 2b are not clear. Figure 3 omits values for the Hamster(16) model. Fig 4 captions has typos.\nCode repository misses step on how to form the required h5 files, and how to form the configuration and main input files.\nThe scope of the study may be better suited to more specialized journals, such as npj Computational Materials. I do not recommend the current version for publication in Nature Communications.", "reviewer_2": "The authors present an interesting study on data drive ∆ML learning of Hamiltonians of complex systems. The article is well motivated. I have the following queries :\n1. If domain of application is optoelectronics, then would you not say modeling excited states are crucial ?.\n2. Referring to “However, ENNs typically require large datasets for training, which can be a limiting factor from a computational perspective.”\nOne major reason why equivariant NN method gained prominence was that they reduced the training sample size requirement and made neural networks more data efficient. So if this sentence is written with in a specialistic context, then it should be amended.\n3. It looks like you are using ∆ML to learn the different between tight binding and ground truth of a Hamiltonian matrix element. Can you pre-train you model Hamiltonian to some specific chemical systems and fine tune on completely different chemical system ? . Kindly discuss/ provide clarity on this fundamental transferability challenges related to addressing this question.\n4. The TB Hamiltonian is defined in Fourier Space to account for periodic boundary condition (equation 1). Is the machine learned part of the delta learning (Equation 4) agnostic toward real space or periodic formulation of equation 1 ? .\n5. Would you not want equation (5) also to have Fourier space formulation to account for periodic effects given equation (1) (also a local descriptor will not be able to capture non local effect by the virtue of fixed cutoff) ?. And this does not look like an equivariant formulation.\n6. Additional comments:\n- Please give clear explanations with respect to variables in Figure 1b.\n- I would also like to see Figure 1c plotted in RMSE including a curve with out the delta ML correction (i.e. pure TB) and see how much RMSE improves relatively when ML correction is added. Would you also kindly provide a similar result for MAPbBr3 ?\n-How are the residuals estimated (general comment)?.\n7. Minor comment :\nThe quality of the plotting could be better (e.g. parts of Fig 2, 4 are not legible ).", "reviewer_3": "The manuscript “Physics-informed Hamiltonian learning for large-scale optoelectronic property prediction” by Schwade et al. presents a model that follows a “learn-the-rules, not the answers” approach. Instead of directly predicting target properties, the authors learn the parameters of a physically motivated Hamiltonian model and then compute observables from it. While such approaches are usually more challenging to design than end-to-end deep neural networks, they offer important advantages. The resulting models encode physics beyond standard rotational, translational, and permutational symmetries, are more interpretable, and typically generalize better from limited data.\nThe authors demonstrate the effectiveness of their HAMSTER model for GaAs, CsPbBr3, and MAPbBr3, achieving mean absolute errors in band-gap prediction on the order of 50 meV with only a small number of reference structures. They also show good scaling to larger system sizes, up to approximately 20k atoms, which is particularly relevant for realistic device-level simulations where electronic properties must be evaluated on large configurations.\nThe manuscript is clearly written, well structured, and the results are convincingly presented. I have only one minor suggestion: it would be very helpful if the authors could provide estimates of computational cost as a function of system size (e.g., wall-clock time for representative system sizes). Since the scalability of the approach is one of its key strengths, these practical performance details support readers in assessing applicability to their own problems.\nWith this minor remark, I strongly support publication in Nature Communications. The approach is a significant and timely contribution to physics-informed machine learning for materials modeling, and it will be of broad interest to both academic and applied research communities."}, "author_response": "Comment 1: Themanuscriptaimstocontributetothe predictionofoptoelectronicproperties\nin large-scale materials. It develops a machine learning model based on tight-binding density\nfunctionaltheory,withtheparametersofane!ectiveHamiltonianfittedusingakernelregression\napproach. The manuscript claims improved e\"ciency in computing eigenvalues.\nResponse: We thank the reviewer for the review of our work. We addressed the reviewers’\nconcerns below.\nComment 2: However, the text is di\"cult to follow for non-specialist readers. Several equa-\ntions or derivations are omitted and need further discussion or proper citation. For example,\nEquation 4 (the descriptor x ) is unclear, and its relationship to Equation 5 is not explained.\ni\nTheoptimizationandfittingofparametersandhyperparametersareconfusing. Authorsshould\ndiscuss the limitations of these method. Several figure captions are incomplete; for instance,\nFigure 2a does not clearly state which property is studied. Clusters in Fig 2b are not clear.\nFigure 3 omits values for the Hamster(16) model. Fig 4 captions has typos. Code repository\nmisses step on how to form the required h5 files, and how to form the configuration and main\ninput files.\nThe scope of the study may be better suited to more specialized journals, such as npj Compu-\ntational Materials. I do not recommend the current version for publication in Nature Commu-\nnications.\nResponse: We thank the reviewer for this valuable feedback. We have carefully revised the\nmanuscript to improve clarity and accessibility. In particular, we have added additional expla-\nnations where equations or derivations were previously omitted, clarified the definition of the\ndescriptor in Eq. 4 and its relationship to Eq. 5, and improved the description of the parameter\noptimization and hyperparameter fitting procedure. We now also include an explicit discussion\nof the limitations of our method.\nWe have further revised the parameter optimization for CsPbBr after identifying a numeri-\n3\ncal instability in large supercells at high temperatures following the transition to the MACE\npotential from VASP’s MLFF. For consistency, we have also re-evaluated all calculations for\nMAPbBr , although no analogous instability was observed for that system.\nthat the substantial improvements to clarity, documentation, and methodological transparency\nsignificantly strengthen the manuscript, and we hope that the revised version now meets the\nstandards for publication in Nature Communications.\nChanges:\n• We have augmented the explaining paragraph for Eq. 4 (p. 3, left column), which now\nreads:\nTo this end, each non-zero matrix element in the Hamiltonian of Eq. 3, determined by a\ncuto! radius r , is associated with a descriptor vector, xR, whose specific construction is\ncut ij\ndetailed below. From the full set of descriptor vectors, we select a subset of size N using\np\nthe k-means clustering algorithm. These serve as kernel support points and are used to\npredict the corrections for remaining or unseen matrix elements (see Methods for details\non the sampling procedure).\n• We have added additional information on how the descriptor vectors are constructed on\np. 4, left column, which reads:\nxR =(Z ,Z ,!rR,ωR,ωR,εR,h ,h )T,\nij i j ij i j ij env,i env,j\nwhere !rR is the distance between atoms i and j with the translation vector R applied to\nij\natom j, ωR (x=i,j) is the angle between the orientation of the orbital on atom x and the\nx\nbonding axis; and εR is the angle between the orbital orientations on atoms i and j.\nij\n• Werewroteandrestructuredthe(hyper-)parameteroptimization,pleaserefertothehigh-\nlighted section in the revised manuscript on page 9-11.\n• We added discussion of limitations of our approach at appropriate places:\nOn p. 8, left column as well as right column (see comments 2 and 4 by reviewer 2)\nOn p. 10, right column:\nNevertheless, ouroptimizationstrategyhassomelimitations. First, theTB,ML,andSOC\nparameters are trained in separate stages to avoid numerical instability, which prevents\nfully capturing their mutual interdependence. Consequently, subtle correlations between\nthese components may not be entirely reflected in the final model. Second, gradient-based\noptimization may be less e!ective if the Hamiltonian exhibits strong nonlinear dependence\non its parameters, and it cannot guarantee convergence to the global minimum. Despite\nthese limitations, our approach provides a stable and computationally e\"cient framework\nfor fitting e!ective Hamiltonians.\n• We improved visual clarity of Fig. 2b and revised and expanded the text describing the\nprocedure and results on page 5.\n• We revised the captions of Figs. 3 and 4 and corrected typos.\n• We have added data for large-scale calculations for CsPbBr after switching to MACE\n3\nmolecular dynamics (see updated figure below). We have reevaluated all calculations to\nensurenumericalstabilityandconsistency. Wehaveaddedadiscussionofthetemperature\ndependence of the band gap in comparison to experimental data on p. 6, left and right\ncolumn:\nThisimprovedstatisticalconvergenceallowsustoassessthetemperaturedependenceofthe\nband gap. While the absolute value of the band gap is a!ected by the known inaccuracies of\nPBE-level DFT, the resulting temperature dependence can still be meaningfully compared\nto experiment. For the 16 16 16 supercell, we obtain a slope of 3.95 10 → 4eV/K,\n→ → →\nin good agreement with the experimental value of 3.41 10 → 4eV/K reported in Ref. 40.\n→\nWe note, however, that thermal expansion is not included in the simulation and that the\ntemperature range considered here extends beyond that of the experimental study, both of\nwhich are expected to a!ect the band gap values.\nWe have also added a description how the MACE models are trained to the methods\nsection on p. 9, left column, which reads:\nUsing the DFT data collected during on-the-fly training, we trained three MACE models\n[9], one per temperature, using 2 message passing layers, messages with 128 channels and\nmax_L of 2, correlation order of 3, spherical harmonics up to l=3, and a radial cuto! of\nforce weights being 1 and 100, respectively, increasing to 1000 and 100 in the last 20% of\nthe epochs, with stochastic weight averaging enabled. The MACE models were trained with\nabatchsizeof5foratotal of600epochs. Thevalidationsetswerechosentobe10%ofthe\nconfigurations in each data set. The reference energies for isolated atoms were obtained\nthrough spin-polarized DFT calculations. The resulting MACE force fields were used with\nthe Lammps code [53], and applied to a 16 16 16 supercell of the cubic structure for\n→ →\n5000 time steps of 0.6fs.\n3\na b\nDFT HAMSTER (2) VBM VBM\nHAMSTER (16)\nExp. slope\nTemperature (K) Band index (-)\nFig. 3 Transferability of the model across temperatures and large-scale\ncalculations.\nComment 3: Code repository misses step on how to form the required h5 files, and how to\nform the configuration and main input files.\nResponse: We have expanded the documentation to include the requested information and\nalso included links within the README itself. For your convenience, we provide links here:\nmain input file, structure data, eigenvalue data\nReviewer 2:\nComment 1: The authors present an interesting study on data drive !ML learning of Hamil-\ntonians of complex systems. The article is well motivated. I have the following queries :\nResponse: We thank the reviewer for their positive assessment of our work and constructive\nfeedback. We gladly address the remaining comments below.\nComment 2: If domain of application is optoelectronics, then would you not say modeling\nexcited states are crucial ?.\nResponse: Whileweagreewiththereviewerthatexcited-statepropertiescanbeanimportant\naspect in the modeling of optoelectronic materials, oftentimes researchers are studying merely\nthe ground-state properties of such materials with methods like DFT. Here, our framework\nprovides a very interesting alternative for large-scale computations. That said, in principle, our\nframework can also be extended to construct e!ective Hamiltonians for excited states, provided\nappropriate reference data (e.g., excited-state energy eigenvalues) are available. Likewise, the\nmethod is flexible with respect to the level of theory used in the underlying electronic-structure\ncalculations that generate the training data.\nHowever, extending the present study to explicitly model excited states lies beyond the scope\nof the current work.\nChanges: We now explicitly mention excited states in the discussion on p. 8, right column\ncolumn which reads:\nThese could include hybrid functionals or many-body perturbation theory, potentially extending\nthe framework to including excited-state e!ects that can play an important role in the modeling\nof optoelectronic properties.\nComment 3: Referring to “However, ENNs typically require large datasets for training, which\n4\ncan be a limiting factor from a computational perspective.” One major reason why equivariant\nNNmethod gained prominencewasthattheyreducedthe trainingsamplesize requirementand\nmade neural networks more data e\"cient. So if this sentence is written with in a specialistic\ncontext, then it should be amended.\nResponse: Wethankthereviewerforpointingthisout. Weagreethatakeyreasonforthesuc-\ncess of equivariant neural networks is their improved data e\"ciency, as enforcing physical sym-\nmetries significantly reduces the amount of training data required compared to non-equivariant\nmodels.\nOuroriginalwordingwasmisleadinginaspecialistcontextanddidnotsu\"cientlyacknowledge\nthis advantage. Our intent was to contrast ENNs with physics-based and hybrid approaches\nthat, in some cases, can achieve comparable accuracy with even fewer training samples by\nincorporating stronger physical priors. We have revised the text accordingly to clarify this\ndistinction and to better reflect the established strengths of ENNs.\nChanges: The statement about data e\"ciency of ENNs on p. 1, right column now reads:\nAlthough ENNs improve data e\"ciency by encoding symmetry and therefore require fewer sam-\nples than comparable non-equivariant models, they may still demand substantial amounts of\ntraining data to capture complex dynamical behavior, such as that exhibited in multi-component\nmaterials at finite temperature.\nComment 4: It looks like you are using !ML to learn the di!erent between tight binding and\nground truth of a Hamiltonian matrix element. Can you pre-train you model Hamiltonian to\nsomespecificchemicalsystemsandfinetuneoncompletelydi!erentchemicalsystem? . Kindly\ndiscuss/provideclarityonthisfundamentaltransferabilitychallengesrelatedtoaddressingthis\nquestion.\nResponse: The reviewer raises an important point about chemical transferability in our !ML\nframework. In principle, transfer learning is possible when the source and target systems share\nsimilarchemicalenvironmentsandorbitalcharacter. Forexample,onecouldpre-trainthemodel\non a simple perovskite, such as CsPbBr , and then fine-tune it on a mixed-halide system (e.g.,\nx 3 x\n→\nlearned parameters should serve as a meaningful starting point for the more complex system.\nIt is also conceivable to initialize parameters for one halogen (e.g., Br) as a starting point for\nanother (e.g., I), provided they have comparable valence electron configurations. However, we\ndo not expect such transferability to hold in general when the chemical environments, orbital\nsymmetries, or electronic properties di!er substantially. In such cases, the parameters learned\nduring pre-training may not remain valid, and the model would require extensive re-training.\nThus, chemical transferability is feasible but we expect it to be inherently limited by the degree\nof similarity between the source and target material systems.\nChanges: We have added a paragraph that discusses chemical transferability on p. 8, left col-\numn, which reads:\nAn important consideration that follows from this capability is chemical transferability. In prin-\nciple, transfer learning is feasible when the source and target systems share similar local en-\nvironments, bonding motifs, and orbital character, allowing pre-trained parameters to serve as\nmeaningful initializations. This may hold for materials that di!er only by modest compositional\nchanges or chemically similar substitutions. However, when electronic structure or orbital sym-\nmetries di!er substantially, previously learned parameters are unlikely to remain valid, and\nretraining becomes necessary. At the same time, these considerations highlight the longer-term\n5\npotential for developing a more general model that spans broader classes of materials and can be\ne\"ciently adapted to new systems with minimal additional training.\nComment5: TheTBHamiltonianisdefinedinFourierSpacetoaccountforperiodicboundary\ncondition (equation 1). Is the machine learned part of the delta learning (Equation 4) agnostic\ntoward real space or periodic formulation of equation 1 ? .\nResponse: While the TB Hamiltonian is expressed in Fourier space to incorporate periodic\nboundary conditions, the underlying matrix elements are computed in real space. In our !ML\nscheme,theMLmodeloperatesdirectlyonthesereal-spaceatomicenvironments. Translational\nsymmetryisautomaticallyrespectedbecauseperiodicimagesofatomsshareidenticallocalenvi-\nronments,leadingtoidenticalTBandMLterms. Consequently,themodelise!ectivelyagnostic\nto whether the Hamiltonian is later represented in real or reciprocal space; the periodicity is\nfully encoded through the real-space atomic positions that serve as input.\nChanges: We have added a short discussion of the translational symmetry of the environment\ncomponent of the descriptor on p. 4, left column which reads:\nThe local environment is identical for all lattice translation vectors, R, as a consequence of\ntranslational symmetry. Note that the integrals in Eq. 5 reduce to those in Eq. 3 when V =1,\nv\nand therefore need to be evaluated only once for both the TB and ML contributions.\nComment 6: Would you not want equation (5) also to have Fourier space formulation to\naccountforperiodice!ectsgivenequation(1)(alsoalocaldescriptorwillnotbeabletocapture\nnon local e!ect by the virtue of fixed cuto!) ?. And this does not look like an equivariant\nformulation.\nResponse: Weagreewiththereviewerthatourcurrentapproachisstrictlylocalandtherefore\ncannot capture long-range or non-local e!ects. A reciprocal-space formulation could, in prin-\nciple, incorporate such interactions; however, it would also destroy the sparse structure of the\nHamiltonian and significantly increase computational cost. Incorporating message-passing or\ngraph-based neural network architectures is indeed a promising future direction for addressing\nnon-localcontributions,althoughsuchmodelsmayrequiresignificantlylargertrainingdatasets.\nRegarding equivariance, our descriptor and TB-based construction enforce the required sym-\nmetry constraints by design through the analytic structure of the Hamiltonian. This approach\ndi!ers from typical equivariant formulations, which embed equivariance directly into the net-\nwork architecture. Instead of transforming equivariantly under input changes, the Hamiltonian\ngeneratedbyourmodelisfullyinvariantundertransformationsthatshouldleaveitunchanged.\nWe have revised the manuscript to clarify this distinction.\nChanges: We revised the paragraph about equivariance which now reads\nFinally, weemphasizethattheorderingofvalueswithinthedescriptorvectorinEq. 6iscritical,\nas inconsistent arrangements can violate the symmetry properties of the Hamiltonian. To pre-\nserve these symmetries, the descriptor is constructed such that its structure remains invariant\nunderoperationsthatshouldproduceequalmatrixelements,forexampletherelationH i R j =H j→i R\n(see, e.g., Ref. 16 for a detailed discussion of Hamiltonian symmetry properties). This should\nbe distinguished from equivariance, which instead ensures that the model output transforms con-\nsistently under symmetry operations, as employed in equivariant models [13-15, 18, 31].\nComment 7: Additional comments: - Please give clear explanations with respect to variables\n6\nin Figure 1b.\n- I would also like to see Figure 1c plotted in RMSE including a curve with out the delta ML\ncorrection (i.e. pure TB) and see how much RMSE improves relatively when ML correction is\nadded. Would you also kindly provide a similar result for MAPbBr3 ?\n-How are the residuals estimated (general comment)?.\nResponse: We gladly provide Fig. 1c plotted in RMSE (excluding the uppermost conduction\nbands which carry reduced weight during optimization, see below) as well as an analysis of the\ntrainingandvalidationerrorforMAPbBr versusthenumberoftrainingstructures(seebelow).\nk-points and structures.\nChanges: We have added a horizontal line in Fig. 1c to indicate the error of the TB model.\nFor Fig. 2a, we have added a data point for zero training structures that indicates the error\nwithout ML corrections.\nWe revised the caption of Fig. 1b which now reads:\nSchematic visualization of the environment descriptor for the matrix element between atoms i\nand j, which treats local environments of the two atoms separately. Di!erent atomic species are\nindicated by purple and gray colors. Atoms within a distance r , which is indicated by dashed\ncut\ncircles, are labeled as k , with x = 1,2,3,... The s and p orbitals of selected atoms are shown\nx\nschematically as red circles and red-blue handles, respectively.\nRoot-mean-squared error versus number of training structures for GaAs\nTraining and validation mean-absolute error versus number of training structures\nfor MAPbBr .\n).\nResponse: Wethankthereviewerforpointingoutthatthevisualclarityofsomeofourfigures\n7\ncan be improved. We revised the figures accordingly.\nChanges: We improved the presentation of Fig. 2a, Fig. 2b, Fig. 3a, Fig. 4a, and Fig. 4b.\nComment 9: There is a README section in the GitHub detailing installation instruction\nResponse: We thank the reviewer for assessing the availability of our code repository.\nReviewer 3:\nComment 1: Themanuscript“Physics-informedHamiltonianlearningforlarge-scaleoptoelec-\ntronic property prediction” by Schwade et al. presents a model that follows a “learn-the-rules,\nnot the answers” approach. Instead of directly predicting target properties, the authors learn\nthe parameters of a physically motivated Hamiltonian model and then compute observables\nfrom it. While such approaches are usually more challenging to design than end-to-end deep\nneural networks, they o!er important advantages. The resulting models encode physics beyond\nstandard rotational, translational, and permutational symmetries, are more interpretable, and\ntypically generalize better from limited data.\nThe authors demonstrate the e!ectiveness of their HAMSTER model for GaAs, CsPbBr3, and\nMAPbBr3, achieving mean absolute errors in band-gap prediction on the order of 50 meV with\nonly a small number of reference structures. They also show good scaling to larger system\nsizes, up to approximately 20k atoms, which is particularly relevant for realistic device-level\nsimulations where electronic properties must be evaluated on large configurations.\nThe manuscript is clearly written, well structured, and the results are convincingly presented.\nI have only one minor suggestion: it would be very helpful if the authors could provide esti-\nmatesofcomputationalcostasafunctionofsystemsize(e.g., wall-clocktimeforrepresentative\nsystem sizes). Since the scalability of the approach is one of its key strengths, these practical\nperformance details support readers in assessing applicability to their own problems.\nWith this minor remark, I strongly support publication in Nature Communications. The ap-\nproachisasignificantandtimelycontributiontophysics-informedmachinelearningformaterials\nmodeling, and it will be of broad interest to both academic and applied research communities.\nResponse: Wethankthereviewerfortheirkindandconstructivefeedbackandgladlyincluded\nnew data and further discussion about the scaling of our approach in the revised manuscript.\nChanges: We have added Fig. 4c (see updated figure below) that shows the runtime of Ham-\nster for full, partial, and no diagonalization of the Hamiltonian. We added 2 paragraphs on\np. 7 that read:\nTo assess the computational performance of Hamster, we benchmark the total runtime for\nthree calculation modes: (i) full diagonalization of the Hamiltonian, (ii) Lanczos-based partial\ndiagonalization targeting six eigenvalues near the VBM, and (iii) construction and output of the\nHamiltonian matrix in Eq. 1 without solving the eigenvalue problem. For each benchmark, the\nmodel is evaluated for 10 independent structures, using the same configurations as in Fig. 4, and\nall calculations are performed on a single compute node using two Intel Xeon Platinum 8168\nCPUs. Full diagonalization of the largest system size (16 16 16) is not tractable on the\n→ →\nchosen, rather modest hardware and is therefore omitted. The resulting runtimes are shown in\nFig. 4c. Since the benchmarks include the first invocation of each function, we fit a function,\nT(x)=T +axn, to the timings and subtract the constant time T =44s, i.e., the initial Julia\n0 0\ncompilation overhead, which disproportionately a!ects benchmarks with short execution times.\nthe Hamiltonian exhibits nearly ideal linear scaling with system size. Partial diagonalization re-\nmainssignificantlyfasterthanfulldiagonalization,butdeviatesfromlinearbehavior,showingan\ne!ective scaling exponent of approximately 1.8. This deviation originates from the shift–invert\ntransformation, which can reduce the sparsity of the Hamiltonian. In contrast, full diagonaliza-\ntion yields a scaling exponent of 2.8, close to the cubic scaling expected for conventional DFT\nimplementations. It is worth noting that the absolute computational cost of DFT calculations\nis substantially higher, primarily due to the overhead associated with repeated self-consistent\nfield iterations and the much larger underlying basis sets. Additional speedups for diagonaliza-\ntion could be achieved through more advanced algorithms, highly parallel implementations, or\nGPU-accelerated solvers.\nWithin our approach, the wall time for the largest system considered remains below one hour,\ndemonstrating that Hamster can e\"ciently handle large-scale systems even on modest hard-\nware. WhenonlytheHamiltonianisrequiredandeigenvaluesarenotcomputed, genuinelylinear\nscaling can be achieved. Moreover, many important physical quantities, such as transport prop-\nertiesorthedensityofstates, canbeevaluatedwithuptolinearscalingbyexploitingthesparsity\nof the Hamiltonian [46, 47]. Combining our Hamiltonian learning framework with such linear-\nscaling methods therefore provides a promising route toward e\"cient simulations of large-scale\nsystems.\na\nb\nc\nCubic scaling\nLinear scaling\nFig. 4 Band gap of MAPbBr as a function of supercell size and temperature,\n3\nwith associated computational scaling.\n• Added competing interests statement.\n• Changed data availability statement. Source data is now available in the Zenodo repo\nhttps://doi.org/10.5281/zenodo.18485403.\n• Added one sentence to acknowledgments.\n10"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for their very careful consideration and thorough responses to my comments. All of my suggestions have been appropriately addressed and incorporated into the revised manuscript and/or the Supporting Information. I agree with all of the authors’ responses and have no further comments. I recommend the manuscript for publication.", "reviewer_2": "My authors have addressed my queries", "reviewer_3": "As in my previous assessment, I regard the manuscript titled “Physics-informed Hamiltonian learning for large-scale optoelectronic property prediction” as a valuable and timely contribution to the field. The work addresses a significant challenge in computational materials science: the accurate and scalable prediction of optoelectronic properties beyond the system sizes accessible to conventional ab initio methods.\nI concur with the other referees on the importance of clarifying certain methodological aspects and of more explicitly discussing the limitations of the present framework in capturing excited-state phenomena. Addressing these points was important for strengthening the manuscript's clarity and impact. That said, I consider the methodological progress demonstrated in this work to be both substantial and of broad interest. The proposed physics-informed Hamiltonian learning framework represents a significant step toward overcoming the computational bottlenecks in large-scale simulations of optoelectronic properties. Given the central role of such properties across a wide range of practical applications, from photovoltaics to semiconductor device design, extending predictive modeling to realistic system sizes is utterly important.\nOverall, I find the manuscript well written, clearly structured, and technically sound. In my view, it merits publication in Nature Communications."}, "author_response": ""}]

Decision letter