Infusing Self-Consistency into Density Functional Theory Hamiltonian Prediction via Deep Equilibrium Models

Authors: zun  wang
chang  liu
nianlong  zou
he  zhang
xinran  wei
lin  huang
lijun  wu
bin  shao
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 16:43:37

Abstract: In this study, we introduce a unified neural network architecture, the Deep Equilibrium Density Functional Theory Hamiltonian (DEQH) model, which incorporates Deep Equilibrium Models (DEQs) for predicting Density Functional Theory (DFT) Hamiltonians. The DEQH model inherently captures the self-consistency nature of Hamiltonian, a critical aspect often overlooked by traditional machine learning approaches for Hamiltonian prediction. By employing DEQ within our model architecture, we circumvent the need for DFT calculations during the training phase to introduce the Hamiltonian's self-consistency, thus addressing computational bottlenecks associated with large or complex systems. We propose a versatile framework that combines DEQ with off-the-shelf machine learning models for predicting Hamiltonians. When benchmarked on the MD17 and QH9 datasets, DEQHNet, an instantiation of the DEQH framework, has demonstrated a significant improvement in prediction accuracy. Beyond a predictor, the DEQH model is a Hamiltonian solver, in the sense that it uses the fixed-point solving capability of the deep equilibrium model to iteratively solve for the Hamiltonian. Ablation studies of DEQHNet further elucidate the network's effectiveness, offering insights into the potential of DEQ-integrated networks for Hamiltonian learning. We open source our implementation at https://github.com/Zun-Wang/DEQHNet.

Main paper: 
 introduction: 
 Density Functional Theory (DFT) is a framework that has revolutionized the study of physical, chemical, and material systems. By providing insights into the electronic structure of matter, DFT has become an indispensable tool for researchers in various scientific fields. Despite its widespread adoption, the computational intensity of DFT poses a significant bottleneck, particularly when applied to large or complex systems.To navigate the computational challenges of DFT, machine learning (ML) has emerged as a powerful ally, capable of predicting a myriad of molecular and material properties with reduced computational overhead [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] . However, these ML approaches typically require bespoke models for different properties, necessitating extensive retraining or redesign.The Hamiltonian, a fundamental quantity in DFT, encapsulates the entire energy and dynamics of a system. Properties of molecules and materials can be considered downstream tasks of the Hamiltonian. Therefore, once the Hamiltonian is accurately predicted, a broad spectrum of properties can be subsequently derived. This realization has spurred interest in using ML to predict the Hamiltonian directly [22] . Recent efforts in this direction have seen substantial improvements by incorporating equivariant graph neural networks, which respect the underlying symmetries of physical laws [23, 24, 25, 26] .In the context of DFT, the Hamiltonian is not just another property-it inherently involves selfconsistent iterative processes to achieve convergence. Some studies [27, 28, 29, 30] have incorporated this iterative nature into network training by intertwining DFT computations with the loss function during the training of neural networks, and even leveraging unlabeled data to enhance the model's generalizability. Nonetheless, these methods unavoidably rely on DFT iterations to solve the generalized eigenvalue problem, rendering them impractical for larger systems. This is because they integrate DFT iterations within the training process, melding DFT with the training loss. To fundamentally address this issue, a paradigm shift in network architecture is required.In this work, we introduce DEQH model, combines deep equilibrium models (DEQs) [31] with off-the-shelf neural network models for the prediction of DFT Hamiltonians. DEQs, characterized by their ability to model systems at an equilibrium state, offer a structural innovation that inherently captures the self-consistent nature of the Hamiltonian without the need for iterative DFT calculations during training. By embedding the equilibrium concept into the architecture, our model seeks to directly learn the fixed-point representation of the Hamiltonian, circumventing the computational complexity while maintaining fidelity to the principles of DFT. Our key contributions are:• The Hamiltonian's iterative qualities are often neglected by standard machine learning approaches for its direct prediction. Our approach integrates DEQs with off-the-shelf ML frameworks, leveraging node features derived from the Hamiltonian and overlap matrix to harness these iterative aspects. Considering computational efficiency and practical applicability, we adopt QHNet as the backbone of the architecture, leading to the development of DEQHNet. • Traditional machine learning models primarily serve as Hamiltonian predictors, and while recent self-consistency integrating frameworks aim to refine training, they incur high computational costs. DEQH model distinguishes itself by acting fundamentally as a solver, iteratively determining the Hamiltonian with the deep equilibrium model's fixed-point capabilities. This intrinsic self-consistency and computational efficiency render DEQH model a scalable approach for precise eigenstate prediction without a significant increase in complexity. • We have benchmarked DEQHNet against the MD17 [22] and QH9 [32] datasets, demonstrating that the incorporation of Hamiltonian self-consistency can significantly enhance predictive accuracy. • We conduct an ablation study on DEQHNet, and present conjectures regarding the efficacy of networks that incorporate DEQs for the task of learning Hamiltonians.This architectural leap holds the potential to unlock scalable, efficient, and accurate predictions of DFT Hamiltonians, paving the way for rapid advancements in the understanding and design of new materials and molecules.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: We have provided cases, statistical numbers for our motivation, and also extensive experiments for verification to support our claim.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
related works: 
 To date, more and more studies have tackled the intricate challenge of employing machine learning (ML) techniques to model the wavefunction directly. The first study in this domain was conducted by Hegde and Bowen [33] , who utilized kernel ridge regression to ascertain the Hamiltonian matrix. Building upon this, Schütt et al. introduced the SchNOrb neural network architecture [22] , which assembles the Hamiltonian matrix of molecules using a block-wise approach based on atom-pair features. Deep Hamiltonian [24] was developed for the prediction of Hamiltonian matrices of periodic systems. PhiSNet [23] leverage the principles from a series of SE(3)-equivariant models to ensure that predictions maintain an exact physical consistency with the orientation of the input structures. DeepH-E3 [25] formulates the DFT Hamiltonian as a function of the structural configuration of the material. This representation inherently upholds Euclidean symmetry, maintaining this invariance seamlessly even when spin-orbit coupling effects are taken into account. QHNet [26] has been meticulously constructed to eliminate 92% of the tensor product operations to enhance the computational efficiency of the model. DeepH-2 [34] is an equivariant local-coordinate transformer which reduces SO(3) convolutions to SO(2) for efficient Hamiltonian prediction.Additionally, an increasing number of studies have recognized the inherent self-consistent iterative nature of the Hamiltonian's eigenproperties. These works integrate DFT into the training process to enhance network performance by using this self-consistent iterative characteristic. For instance, Ref. [27] introduced a self-consistent loss that utilizes unlabeled data through amortized DFT, assisting the model training and enabling the exploration of a broader chemical space to improve generalization capabilities. Ref. [30] circumvent the creation of datasets by pretraining with DFT iterations and employing an implicit DFT loss, denoted as E(•), as the training loss function. Simultaneously, a similar approach [28] uses a network predicting the Hamiltonian matrix as an input to DFT to obtain energies, which are then used as optimization targets. Ref. [29] has approached the problem by simultaneously predicting both the Hamiltonian matrix and the density matrix, introducing the direct inversion in the iterative subspace (DIIS) error ϵ = HDS -SDH as a training signal to avoid the costly self-consistent field calculations. Furthermore, this method can estimate the accuracy of the predictions, making it amenable to integrate with active learning strategies.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: We have discussed the limitation in Supplementary Material A.12 and provided the potential solutions.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
preliminary: 
 Equivariance Consider a mapping function L : X → Y that transforms inputs from space X to outputs within space Y. The function L is deemed G-equivariant if it consistently upholds the symmetry induced by a group G in its mappings. Specifically, for every element g belonging to group G, the following relation is satisfied:EQUATIONwhere D X and D Y denote the actions of the group G on the spaces X and Y, respectively. This condition guarantees that the function L faithfully translates the symmetry operations performed on the inputs by G into corresponding transformations in the outputs. For additional concepts relevant to this discussion, please refer to the Supplementary Material A.2-A.3.The primary goal of most electronic structure methodologies is to accurately solve the electronic Schrödinger equation:Ĥel Ψ el = E el Ψ el (2)where Ĥel is the electronic Hamiltonian operator, which encapsulates the interactions and dynamics of electrons, Ψ el represents the electronic wavefunction, and E el denotes the ground state energy. Typically, the electronic wavefunction Ψ el is approximated by an antisymmetric combination of molecular orbitals {ψ i (r)} (r denotes the electronic coordinates). Generally, each molecular orbital ψ i is formulated as a linear combination of atom-centered basis functions {ϕ j (r)}:EQUATIONThis formulation leads to the matrix equation:HC = SCϵ (4)in which the Hamiltonian is represented as a matrix H. The overlap matrix S, with elements S ij = ϕ i (r)ϕ j (r)dr, is introduced to account for the non-orthonormality of the basis functions. To obtain the coefficients C ij that define the wavefunction Ψ el and the orbital energies ε, one must solve a generalized eigenvalue problem. More introduction refers to the Supplementary Material A.1.The deep equilibrium model (DEQ), a class of innovative implicit layer models introduced by Ref. [31] , has garnered considerable attention for its remarkable performance across various large-scale vision and natural language processing tasks. These models often rival or even surpass the state of the art established by traditional explicit models, as demonstrated in subsequent studies by Bai et al. [35] . The fundamental principle underpinning this methodology is the concept of an implicit layer that converges to a fixed point through an iterative process.A typical k-layer deep network h : X → Y is defined by a stack of layers,z 1 = x, z i+1 = σ (W i z i + b i ) , i = 1, • • • , k -1, h(x) = W k z k + b k . (5)To allow for the continuous integration of input information throughout the depth of the model, an input injection mechanism is incorporated across the layers, which involves adding a linear transformation of the input, denoted as U x. Consequently, this augmented model can be characterized asEQUATIONTo model an infinitely deep network (k → ∞), it is observed that the layer values typically converge to a fixed point, z * = σ (W z * + U x + b). The inclusion of input injection U x in the model is critical due to the equilibrium point's independence from any initial z 1 value. Otherwise, the network's output would become paradoxically invariant to the input. Drawing on the insight that deep sequence model layers often converge to a fixed point, the DEQ approach employs root-finding methods to identify these equilibrium points, thereby optimizing the learning process. For a comprehensive explanation of DEQ, please see the Supplementary Material A.7.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: We open source our implementation at https://github.com/Zun-Wang/ DEQHNet.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
methods: 
 Our DEQ-based framework, depicted in Fig. 1 , is designed as a Hamiltonian solver, contrasting with traditional predictors. It inputs the Hamiltonian, overlap matrix, and structural data to yield the Hamiltonian output, which is then recursively processed through DEQ to find its fixed point. Theoretically, the Hamiltonian and structural information alone could suffice, as structural data can function as a conduit within DEQ to maintain fixed point traits, with atomic orbital representations learned from the structure itself (see Sec. 5.3). Nevertheless, we include the overlap matrix, given its role in DFT's iterative equation (Eq. 4), ease of calculation, and similarity in node feature construction.The matrix element between atom i and j of single-electron operator Ô ∈ { Ĥ, Ŝ} represented in the atomic orbitals {Φ} isEQUATIONwhere µ := (n 1 , l 1 , m 1 ) and ν := (n 2 , l 2 , m 2 ) are the basis functions in the set of atomic orbitals centered at each atom, respectively. The n, l, and m are principal, azimuthal, and magnetic quantum numbers, respectively. Refer to OrbNet-Equi [36] , for each diagonal block of T, i.e., T AA , defined for an on-site atom pair (A, A), there exists a set of T -independent coefficients Q µν nlm such that the following linear transformation ψ, is injective and satisfy equivariance. The existence of Q could be proved by Wigner-Eckart theorem.EQUATIONFor the sake of computational feasibility, a physically-motivated scheme is employed to tabulate Q and produce order-1 equivariant embeddings h A , using on-site 3-index overlap integrals Q:Qµ,ν nlm := Qn1,l1,m1;n2,l2,m2EQUATIONwhere Φ A are the atomic orbital basis, and ΦA are auxiliary Gaussian-type basis functions. Q adheres to equivariance constraints due to its relation to SO(3) Clebsch-Gordan coefficients C lm l1m1;l2m2 ∝ r∈S 2 Y l1m1 (r)Y l2m2 (r)Y * lm (r)dr. Consequently, the role of Q is to couple the indices l 1 and l 2 of the Hamiltonian matrix elements, resulting in features of order l. Details could be found in Supplementary Material A.4-A.6.In this paper, we have chosen to use QHNet [26] because of its versatility and adaptability. QHNet is capable of supporting multiple molecules simultaneously. Its matrix prediction module is designed to predict matrices of varying sizes for different molecules [32] , making it a more flexible and robust choice for our study. As illustrated in Fig. 2 , we present the comprehensive architecture of DEQHNet. To effectively incorporate the self-consistent nature of the Hamiltonian through the DEQ mechanism within the network architecture, it is imperative that the network's filters accept the Hamiltonian as input. Additionally, structural information must be used as an injection to establish dependencies at the fixed point. Consequently, this necessitates the design of a Hamiltonian solver tailored to these requirements. This solver will enable the network to iteratively refine its predictions, ensuring that the output Hamiltonian is self-consistent and in line with the structural information provided, thereby harnessing the full potential of the DEQ framework in electronic structure modeling.Injection As indicated in the boxed area of Fig. 2 (a), the section highlighted represents the injection component of DEQHNet. Within this component, the atomic number {Z} is input as the node feature, while the relative position vectors ⃗ r ij are processed through Radial Basis Functions (RBF) and spherical harmonics to construct the edge features that are invariant and equivariant, respectively. This design ensures that the crucial chemical information, such as the type of elements and their spatial relations, are effectively encoded into the network, allowing DEQHNet to accurately capture the nuances of molecular structures and interactions. Filter Our design draws inspiration from QHNet [26] and Equiformer [15] to handle structural information, as well as node features constructed from the Hamiltonian and the overlap matrix. As depicted in Fig. 2 (d), the Hamiltonian and overlap matrix are initially processed using the previously described diagonal reduction method to construct node features. It is noteworthy that the initial value of the Hamiltonian provided to the Hamiltonian solver is set as the identity matrix I, a deliberate choice to preserve the Hermiticity of the Hamiltonian. These features are then normalized as x H/S(l T ) i before being fed into the module illustrated in Fig. 2 (b ). Within this module, the equivariant node features are combined with the equivariant relative position vectors through a tensor product, resulting in the generation of new node features,EQUATIONAs depicted in Fig. 2 (c), the filtering module within the Tensor Field Network (TFN) layer regulates the influence of messages coming from other nodes. In processing structural features, the inner product I l ij of the l-order irreducible representation on a pair is calculated asEQUATIONSubsequently, the pairwise cosine similarity of the irreducible representations, along with the 0-order irreducible expression, are fed into a MLP to compute attention scores a ij . The calculation of the message m ij from node j to node i is constructed through a tensor product of the transformed relative position vector rij and the feature vector xlin jc ,EQUATIONThe output irreducible representation xl i is then obtained by aggregating the messages m l ij and self-connections to form the updated feature:EQUATIONSubsequently, the node features xH/S,l i constructed from the Hamiltonian and overlap matrix are merged with the node features xl i containing structural information. This amalgamation is employed to construct the iteratively refined Hamiltonian for the solver. As shown in Fig. 2 (e), the block is designed to harness the tensor product between nodes to separately construct on-site features and pairwise features. These are then channeled into the module depicted in Fig. 2(f)-(g) , where the inverse operation of the tensor product [26] , i.e. the tensor expansion operation,EQUATIONis applied to construct the Hamiltonian for the next iteration step of the solver. We have developed a complete Hamiltonian solver that seamlessly integrates with the Differentiable Equation (DEQ) method for fixed-point iterations, enforcing the Hamiltonian's inherent self-consistency. This solver iteratively refines the Hamiltonian, mirroring the self-consistency principle in physics. Utilizing the DEQ framework, it dynamically adjusts parameters to achieve convergence, yielding a Hamiltonian that precisely captures the interactions. This approach not only improves accuracy but also aligns with the iterative nature of electronic structure calculations.As illustrated in Fig. 3 , the DEQH model operates as a Hamiltonian solver, learning the iterative processH * = f (H * , Z, R).In contrast, models like QHNet function as predictors, using the formula H = f (Z, R), where Z and R represent atomic numbers and coordinates, respectively. The function f typically involves a neural network, and H denotes the Hamiltonian. The superscript H * indicates the converged Hamiltonian.In practice, if the dataset includes overlap matrices, these can also be incorporated into the network input, transforming them into equivariant node features akin to the Hamiltonian. Even if overlap matrices are not initially present, they can be easily calculated during data preprocessing using Z and R. The overlap matrix provides a wealth of detailed information, enhancing the model's input. Consequently, the equations for the DEQH and the modified Hamiltonian predictor become H * = f (H * , Z, R, S) and H = f (Z, R, S), respectively, where S is the overlap matrix. 
results: 
 We evaluated the performance of DEQHNet by applying it to the MD17 [22] and QH9 [32] datasets.We have further explored the convergence behavior of DEQHNet. Additionally, we conducted an ablation study to isolate and understand the specific impacts of the DEQ mechanism and the architectural design on our method's effectiveness. For descriptions of the metrics and results pertaining to acceleration ratio, please consult the Supplementary Material A.8-A.9.MD17 dataset [22] encompasses an array of molecular properties, including structures, energies, forces, and Hamiltonians-specifically Fock or Kohn-Sham matrices-as well as overlap matrices pertaining to molecules such as water, ethanol, malondialdehyde, and uracil. All computational analyses were executed at the PBE/def2-SVP level of theory [37, 38] , employing the ORCA electronic structure software [39] . Table 1 shows DEQHNet's performance on the MD17 dataset against QHNet. Except for water, DEQHNet achieves lower mean absolute error (MAE) in Hamiltonian predictions for three molecules, a benefit likely from more extensive training data that prevents overfitting. The predicted Hamiltonian improves orbital coefficients and the MAE for orbital energies has risen, which indicates that lower Hamiltonian MAE doesn't necessarily correlate with better orbital energy predictions.Table 1 : Comparison of the MAEs between QHNet and DEQHNet trained on MD17 dataset. The optimal and the second best values are highlighted in bold and underlining, respectively. It should be noted that due to the design of DeepH [24] , the data used is re-labeled by OpenMX [40] , which may lead to different results compared to other models. Furthermore, considering the high training cost of PhiSNet (ori) [23] , PhiSNet (reproduce) results from the QHNet [26] article are provided for a more balanced and comprehensive comparison. The QH9 dataset [32] , crafted using PySCF [41] , is split into QH9-stable and QH9-dynamic. QH9stable contains 130,831 Hamiltonian matrices from QM9 [42, 43] , while QH9-dynamic has 100geometry trajectories for 999 molecules. DFT calculations employ a grid level of 3, SCF convergence with a 10 -13 tolerance, and a 3.16 × 10 -5 gradient limit. The B3LYP functional and def2SVP GTO basis set are used, with DIIS aiding SCF convergence. QH9-dynamic's molecular dynamics were run at 300K in an NVE ensemble. The QH9-dynamic-100k subset features 0.12 fs step trajectories, taken every 10th step for 1,000 steps. Table 2 displays the results of DEQHNet on the QH9 dataset, where a significant reduction in the MAE of the Hamiltonian across both the QH9-stable and QH9-dynamic subsets can be seen. The decrease in MAE is particularly striking for the QH9-dynamic subset, where it nearly halves. The fidelity of the orbital coefficients, derived from the Hamiltonians predicted by DEQHNet, reached over 99% accuracy. This is attributed to the fact that the overlap matrix and the Hamiltonian offer significant insights into the basis set, thereby facilitating the expression of orbital coefficients. Similar to the trend observed with the MD17 dataset, there is no direct positive correlation between the MAE of the Hamiltonian and the orbital energies. Additional experiments was performed to elucidate this situation, with details provided in the Supplementary Material A.10. Dataset Model H [10 -6 E h ] ↓ ϵ [10 -6 E h ] ↓ ψ [10We further explored DEQHNet through additional experiments. As depicted in Fig. 4 (a), we inferred 50 uracil test set configurations using DEQHNet checkpoints from various training phases. Initially, the model required many iterations to achieve convergence but soon stabilized at two or three iterations, showing quick adaptation to the Hamiltonian's iterative solution. In comparison, PySCF calculations for these configurations, which are structurally similar due to their molecular dynamics origin, consistently converged at 15 SCF iterations. Initially, DEQHNet surpassed this number but soon settled at fewer iterations. This suggests that DEQHNet learns a distilled version of the DFT solver, and unlike PySCF, which considers energy for convergence, it uses Hamiltonian differences for this purpose.Furthermore, we randomly selected a molecule from the QH9-dynamic-geo test set for additional analysis. In Fig. 4 (b), we plotted the MAE of the Hamiltonian at each iteration during the PySCF computation and the DEQHNet inference relative to the true Hamiltonian. DEQHNet starts with an identity matrix as the initial Hamiltonian, leading to a high initial error, while PySCF employs MINAO density [41] -the superposition of orbitals of each atom in isolation-for its initial guess, resulting in a lower starting MAE. DEQHNet typically stabilizes within three iterations, unlike PySCF, which can have initial optimization fluctuations. By initiating DEQHNet with PySCF's starting Hamiltonian [44] , we noted improved stability in convergence. This suggests that DEQHNet's iterative refinement can benefit from a more accurate initial guess, enhancing the stability and reliability of convergence, an important aspect in electronic structure computations. A.1 Electronic structuresThe Hamiltonian for a system composed of electrons and nuclei is given by:EQUATIONwhere electrons are represented by lowercase subscripts, while nuclei, characterized by their charge Z I and mass M I , are denoted by uppercase subscripts.In the general Hamiltonian, the inverse mass of the nuclei, denoted as 1/M I , stands out as the only term that can be considered "small." This allows to establish a perturbation series based on this parameter, which holds broad applicability across the fully interacting electron-nuclei system. Initially, by setting the nuclear mass to infinity, we can disregard the kinetic energy of the nuclei. This leads us to the Born-Oppenheimer or adiabatic approximation. This approximation proves highly effective for various applications, such as calculating nuclear vibration modes in most solids.With the omission of nuclear kinetic energy, the fundamental Hamiltonian for the theory of electronic structure simplifies to:EQUATIONAdopting Hartree atomic units where ℏ = m e = e = 4πϵ 0 = 1, we can express the kinetic energy operator for the electrons, T as:EQUATIONThe external potential acting on the electrons due to the nuclei, V ext, is given by:EQUATIONVint encompasses the electron-electron interactions,EQUATIONwhich quantifies the electron-electron Coulomb repulsion, a fundamental aspect of electronic interactions. The last term, E II , represents the classical interactions between nuclei themselves and encompasses any additional contributions to the system's total energy that are not directly related to the description of electron behavior. In this model, the influence of nuclei on the electrons is integrated into a fixed "external" potential. This formulation remains applicable even when the straightforward nuclear Coulomb interaction is substituted with a pseudopotential, which accounts for core electron effects, albeit these potentials are "non-local".There are two basic independent-particle approaches that may be classified as "non-interacting" and "Hartree-Fock" [45] . They are similar in that each assumes the electrons are uncorrelated except that they must obey the exclusion principle. However, they are different in that Hartree-Fock includes the electron-electron Coulomb interaction in the energy, while neglecting the correlation that is introduced in the true wavefunction due to those interactions. In general, "non-interacting" theories have some effective potential that incorporates some effect of the real interaction, but there is no interaction term explicitly included in the effective Hamiltonian. This approach is often referred to as "Hartree" or "Hartree-like," after D. R. Hartree who included an average Coulomb interaction in a rather heuristic way. More to the point of modern calculations, all calculations following the Kohn-Sham method involve a non-interacting Hamiltonian with an effective potential chosen to incorporate exchange and correlation effects approximately.According to our comprehensive definition, calculations for non-interacting electrons necessitate solving a Schrödinger-like equation, represented as:EQUATIONwhere V σ eff (r) denotes an effective potential influencing each electron of spin σ at position r. To determine the ground state for a collection of non-interacting electrons, one must populate the lowest eigenstates, adhering to the exclusion principle. If the Hamiltonian lacks spin-dependence, the spin states, both up and down, are considered degenerate, effectively doubling the count for these states. The occupation of higher energy eigenstates characterizes the excited states. The primary rationale for employing such independent-particle equations for electrons in materials is founded on density functional theory.At finite temperatures, it is straightforward to utilize the general formulas of statistical mechanics to demonstrate that the equilibrium distribution of electrons conforms to the Fermi-Dirac (or Bose-Einstein) distribution for occupation numbers as a function of energy. The expectation value, aggregates over many-body states Ψ, each defined by a set of occupation numbers n σ i for individual particle states with energies ϵ σ i . With each n σ i being either 0 or 1 and the sum i n σ i = N σ , it can be shown that:EQUATIONwhere ⟨ψ σ i | Ô|ψ σ i ⟩ represents the expectation value of the operator Ô for the one-particle state ψ σ i , and f σ i indicates the probability of finding an electron in state i, spin σ. Specifically, the Fermi-Dirac distribution is given by:EQUATION)where µ represents the Fermi energy or the chemical potential of the electrons. The system's energy, for instance, is the weighted sum of the energies of non-interacting particles:EQUATIONSimilar to the general many-body scenario, a single-body density matrix operator can be defined:EQUATIONwhere the expectation value ⟨ Ô⟩ is equal to Tr(ρ Ô). For a specific representation involving spin and position, ρ is expressed as:EQUATIONwhere the density is the diagonal component:EQUATION)The Hartree-Fock method is another fundamental approach in many-particle theory, which involves constructing a properly antisymmetrized determinant wavefunction for a predefined number, N , of electrons. The goal is to find a single determinant that minimizes the total energy of the fully interacting Hamiltonian. In cases where there is no spin-orbit interaction, this wavefunction, denoted as Φ, can be expressed as a Slater determinant:EQUATIONwhere each ϕ i (r j , σ j ) is a single-particle "spin-orbital" comprising a spatial function, ψ i (r j ), and a spin function, α i (σ j ). It's important to note that ψ i (r j ) is spin-independent in closed-shell systems, aligning with the "spin-restricted Hartree-Fock approximation." These spin-orbitals must not only be linearly independent but also orthonormal to simplify the equations significantly. It can be straightforwardly demonstrated that Φ is normalized to 1.Moreover, if the Hamiltonian does not depend on spin or is diagonal in the spin basis σ = | ↑⟩; | ↓⟩, the expected value of the Hamiltonian using Hartree atomic units with the wavefunction Φ is calculated as follows:EQUATIONThe first term in the Hartree-Fock method consolidates the single-body expectation values by summing over the orbitals, while the third and fourth terms address the direct and exchange interactions between electrons, involving double sums. Conventionally, we include the "self-interaction" (where i = j), a concept deemed spurious but necessary as it cancels out in the aggregation of direct and exchange interactions. When this self-interaction is accounted for, the cumulative effect across all orbitals represents the electron density, and the direct term effectively simplifies to the Hartree energy. The "exchange" term specifically interacts between electrons of the same spin, owing to the orthogonality of spin components in orbitals for opposite spins. These interactions, crucial to the energy computation. The essence of the Hartree-Fock approach is to minimize the total energy subject to all constraints imposed by the wavefunction's degrees of freedom. This minimization process leverages orthonormality to streamline the equations and employs Lagrange multipliers to ensure this condition is met throughout.A real basis of spherical harmonics Y lm : S 2 → R can be defined in terms of their complex analogues Y m l : S 2 → C by settingEQUATIONThe Condon-Shortley phase convention is frequently employed to maintain uniformity and consistency across the analysis. The corresponding inverse equations defining the complex spherical harmonics Y m l : S 2 → C in terms of the real spherical harmonics Y lm : S 2 → R areEQUATIONIt is well-established from the analytical solutions of the hydrogen atom that the eigenfunctions corresponding to the angular component of the wave function manifest as spherical harmonics. Intriguingly, when magnetic terms are absent, the solutions to the non-relativistic Schrödinger equation can be rendered as real functions. This characteristic underpins the prevalent use of realform basis functions in electronic structure computations, as it obviates the need for complex algebra in the software implementations. It is essential to recognize that these real-valued functions occupy an identical functional space to that of their complex counterparts, ensuring no loss of generality or completeness in the solutions.Irreps representation The model's equivariant framework harnesses the special orthogonal group SO(3) to represent the 3D rotational symmetries that are fundamental to molecular structures. It utilizes the irreducible representations (irreps) of SO(3), denoted by an integer l, which resonate with spherical harmonic functions Y lm . These harmonics confer rotational attributes upon the feature vectors, granting the model rotational equivariance and allowing for uniform geometric property assessment.Tensor product In pursuit of heightened expressiveness within the model, it orchestrates interactions among irrep features corresponding to different angular momenta l via the tensor product. This operation combines two irreps with angular momenta l 1 and l 2 to form a new irrep characterized by angular momentum l 3 , employing Clebsch-Gordan coefficients to facilitate an expansion that is weighted by w m1,m2 :h l3,m3 = h l1,m1 ⊗ h l2,m2 = m1,m2w m1,m2 C l3,m3 l1,m1,l2,m2 h l1,m1 h l2,m2 .This mathematical formulation enables the synthesis of complex features from simpler ones, capturing the intricate interplay of angular momentum in the molecular description.The Wigner-Eckart theorem [46] , a cornerstone of representation theory, articulates that the matrix elements of spherical tensor operators, when projected onto a basis of angular momentum eigenstates, can be factorized into two distinct components: a term invariant to angular momentum orientation and a Clebsch-Gordan coefficient. It provides a bridge connecting the symmetry transformation groups governing spatial configurations, as applied within the Schrödinger equation framework, to the fundamental conservation principles of energy, momentum, and angular momentum.Mathematically articulated, the Wigner-Eckart theorem posits that for a given tensor operator T (k) and two states with angular momenta j and j ′ , there exists a reduced matrix element ⟨j|T (k) vertj ′ ⟩ that is invariant with respect to the magnetic quantum numbers m, m ′ , and q. The theorem asserts that for all values of these quantum numbers, the matrix elements of the tensor operator satisfy the relation:EQUATIONwhereT (k) qis the q-th component of the spherical tensor operator T (k) of rank k, and |jm⟩ represents an eigenstate of the total angular momentum operator J 2 and its z-component J z . ⟨j ′ m ′ kq|jm⟩ is the Clebsch-Gordan coefficient, which mediates the coupling of angular momentum j ′ with k to yield j. ⟨j ∥ T (k) ∥ j ′ ⟩ symbolizes the reduced matrix element, a scalar that encapsulates the essence of the tensor operator independent of the magnetic quantum numbers. This theorem elegantly decouples the geometric dependencies from the dynamic properties of the matrix elements, thereby simplifying the computation of transition amplitudes.The order-l equivariant embeddings h l A using on-site three-index overlap integrals Q: Qµ,ν nlm = Qn1,l1,m1;n2,l2,m2EQUATIONwherel 1 l 2 l m 1 m 2 mdenotes the Wigner 3j-symbols, which are proportional to the Clebsch-Gordan coefficients.Adhering to the methodology outlined in Ref. [36] , we define the auxiliary Gaussian-type basis functions, and for the sake of brevity, we consider their form at the origin, x A = 0:EQUATION)where c n,l is a normalization constant such that r || Φn,l,m A (r))|| 2 dr = 1 following standard conventions. For numerical experiments considered in this work the scale parameters γ are chosen as (in atomic units):EQUATIONγ n,l=1 := 32EQUATIONγ n,l=2 := 4.0 • (0.25) n-1 where n ∈ {1, 2, 3, 4}Note that the auxiliary basis ΦA is independent of the atomic numbers thus the resulting h A are of equal length for all chemical elements.Here, we adopt a generalized notation for the DEQ [31] function, denoted as f (z, x). This notation encapsulates the earlier specified function f (x, y) = σ(W z + U x + b), where σ represents the activation function, and W , U , and b denote the weight matrix, input weight matrix, and bias vector, respectively. Our objective here is to determine a stable fixed point z * that satisfies the equilibrium condition z * = f (z * , x).Any deep network-irrespective of its depth or the intricacy of its connections-can be succinctly encapsulated within the framework of a single-layer DEQ model. Remarkably, this encapsulation does not succumb to the typical exponential surge in parameters that is often associated with universal function approximation theorems for single-layer models. In essence, a single-layer DEQ model can represent the functional capacity of any network without necessitating an increase in the number of parameters. Consider a conventional deep network described by the composition of two functions, denoted as y = g 2 (g 1 (x)). This construct can be seamlessly transmuted into a single-layer DEQ model by amalgamating all intermediate variables derived from the computation into an extended vector,EQUATION)At an equilibrium point z * of this function,EQUATION)We can concatenate all intermediate results of a computational graph into the vector z, and designate the function f as the operator that applies the "subsequent" computation in the graph to each of these elements. This theoretical framework unequivocally demonstrates the representational strength of a single DEQ layer.To elucidate the principles of implicit backpropagation tailored to DEQ models-and by extension, to any fixed-point iterative layer-we will primarily concentrate on the DEQ model's specific form herein. Our objective is to compute the vector-Jacobian product ( ∂z * (•) ∂(•) ) T y for a given vector y, where (•) symbolizes any variable with respect to which we seek to differentiate the fixed point, all of which influence the determined fixed point z * . Differentiating both sides of the fixed-point equation yields:EQUATIONwhere z * (•) denotes the case where z * is treated as an implicit function of the differentiated quantity, and z * alone signifies the equilibrium value. The second equality arises by applying the multivariate chain rule. Rearranging terms, we arrive at an explicit expression for the Jacobian:EQUATIONwhere the terms on the right-hand side can be computed using conventional automatic differentiation.To compute the vector-Jacobian product,EQUATION)In practice, the focal term is the solution to the linear system (abbreviated as g):EQUATIONwhich can be expressed as a fixed-point equation:EQUATION)The convergence of this equation via forward iteration hinges on the stability of the Jacobian ∂f (z * ,x) ∂z * , which also underlies the local stability of the forward iterative process.The derivation of the vector-Jacobian product for a DEQ layer can be distilled into a two-step procedure:• Solve the Fixed-Point Equation g = ∂f (z * ,x) ∂z * T g + y. This can be done either through direct analytical inversion or, more commonly, by employing an iterative method that necessitates only the multiplications by ∂f (z * ,x) ∂z * T . Such multiplications can be expediently performed using standard automatic differentiation, as they are vector-Jacobian products in their own right.• Calculate the final vector-Jacobian product ∂z * (•)∂(•) T y = ∂f (z * ,x) ∂(•)T g. This step also leverages the capabilities of conventional automatic differentiation, as it involves the computation of another vector-Jacobian product.By executing these steps, one can integrate DEQ layers into the backpropagation algorithm, thereby harnessing the power of implicit differentiation to effectively train deep equilibrium models.To evaluate the accuracy and computational efficiency of the predicted Hamiltonian matrix, we adopt a set of metrics as demonstrated in Ref. [32] to evaluate the model's performance:The MAE against DFT-computed reference data, considering both diagonal and off-diagonal blocks that represent intra-and inter-atomic interactions, respectively. Given the sparsity induced by distant atom pairs in larger molecules, we separately assess the MAE for diagonal and off-diagonal blocks along with the total MAE of the matrix.The MAE for occupied orbital energies, including HOMO and LUMO levels derived from the predicted and reference Hamiltonian matrices, is used to measure the matrix's predictive accuracy for these critical properties.To determine the similarity between predicted and reference electronic wavefunctions, we compute the cosine similarity of the coefficients for occupied molecular orbitals, which are crucial for inferring chemical properties. Error-level ratio ↑ 0.406 ± 0.066 0.436 ± 0.050 Acceleration ratio Acceleration ratios, such as the achieved ratio and the error-level ratio, could be used to evaluate how well the predicted Hamiltonian matrix speeds up DFT calculations. These ratios compare the number of optimization steps needed when using the predicted matrix versus traditional initial guesses, and the steps required to reach a comparable error level in the DFT SCF cycle.The DEQHNet, after being trained on the QH9 dataset, undergoes evaluation on a meticulously selected subset of molecules. This subset consists of 50 molecules that have been randomly chosen from the overlapping elements of the QH9-stable and QH9-dynamic's respective test subsets. The intersection ensures that the evaluation is conducted on a common ground, facilitating a direct comparison of the performance results between the stable and dynamic regimes. As illustrated in Table 3 , the study underscores the feasibility of leveraging DEQHNet as initial seeds for electronic structure computations, which can potentially reduce the time to convergence significantly.In the experiments with MD17 and QH9, we observed a decrease in the Mean Absolute Error (MAE) of the Hamiltonian, but the MAE of the orbital energy did not show a corresponding positive trend. To investigate this, we conducted additional experiments. We randomly selected a molecule, added a Hermitian Gaussian noise matrix to its Hamiltonian, and solved the corresponding generalized eigenvalue equation. As shown in Fig. 6 , it can be seen that as the noise on the Hamiltonian gradually increases, the range of error in the orbital energy becomes larger. The errors in the orbital energy of QHNet and DEQHNet on MD17 and QH9 are all within the range shown in the figure. Only when the MAE of the Hamiltonian is sufficiently small can it be ensured that the corresponding orbital energy error is also small. In the comparison between QHNet and PhiSNet, we observed similar situations. This somewhat suggests that using existing Hamiltonian error measurements cannot definitively reflect the errors in the eigenvalues. These evidence suggests that it is a nontrivial and long-standing challenge in the domain of Hamiltonian prediction to design a proper metric on the matrix space so that reflects the metric of derived quantities e.g. energy, for which we will investigate in future work. The overarching framework is implemented using PyTorch [47] , PyTorch Geometric [48] , e3nn [49] ,and Torchdeq [50] libraries. The computation of on-site three-index overlap integrals is conducted by Psi4 software [51] , and subsequently converted to adhere to the conventions of PySCF [41] . The loss function for training DEQHNet model is the summation of the Frobenius norm and L1 norm of the absolute error of the matrix of Hamiltonians,EQUATIONwhere N is the number of elements in Hamiltonian matrix. H and Ĥ denote the predicted and ground-truth Hamiltonian respectively. AdamW optimizer [52] was adopted. DEQHNet model was trained on single 16G Nvidia Tesla V100 GPU for MD17 dataset and single 80G Nvidia Tesla A100 GPU for QH9 dataset. It is noteworthy that the iterative process inherent in DEQ may impact training efficiency to a certain extent. In our experiments with the MD17 dataset, we observed that the model eventually converges after only 2-3 DEQ iterations. Therefore, to enhance training efficiency, we set the maximum number of iterations for both the forward and backward passes of DEQ to 3 in our experiments with the QH9 dataset. With this configuration, the training time on the QH9 data was approximately 1.5 times that of QHNet. Furthermore, detailed settings of hyperparameters are summarized in the Table 4 .The Hamiltonian exhibits self-consistent iterative properties. Specifically, in the DFT solution process, the information from the previously obtained Hamiltonian is used to construct the current Hamiltonian for the current DFT solution. This process is repeated until certain criteria meet a set threshold. We hypothesize that introducing DEQ enables the network to learn the iterative mechanism of solving the Hamiltonian, rather than directly learning the solution to the Hamiltonian. This shift could allow the network to better capture inherent physical laws, thereby improving generalizability. Compared to models without self-consistency, DEQ introduces additional orbital information to the network input. This, coupled with the introduction of scientific priors, can expedite DFT convergence even when the model's output is simply used as the initial value for DFT. In contrast with methods that introduce self-consistency (which often include DFT computations in the loss function during model training), our method offers an architectural enhancement. The DEQH model does not directly include DFT computations, leading to a more efficient computational process. Empirically, if the dataset and molecular size are relatively small and the cost of DFT computations during training is acceptable, previous methods introducing self-consistency might outperform ours (as suggested by DEQHNet's results on water, where DEQHNet requires more data). However, when the dataset is large enough, or when it includes large molecules where DFT computation cost is prohibitive, the benefits of the DEQH model become more significant.The DEQH model is a versatile approach that can be seamlessly integrated with current machine learning models for Hamiltonian prediction. In this paper, we opted for QHNet [26] over PhiSNet [23] due to the current implementation of PhiSNet being restricted to single-molecule support. This limitation stems from PhiSNet's matrix prediction module, which is tailored to predict matrices of fixed size for identical molecules [32] . Additionally, based on the ablation studies, we also posit that providing the overlap matrix could enable training on data whose DFT labels are associated with varying basis sets, whereas models lacking it may be restricted to consistent DFT-level datasets. Furthermore, the DEQH model employs a method that is entirely orthogonal to the inclusion of iterative DFT processes during training to achieve the self-consistency of the Hamiltonian-an architectural innovation within the model. Incorporating DFT iterations into the existing Hamiltonian training introduces a more explicit physical prior [? 28, 29, 30] compared to the DEQH model, which requires a larger dataset to enable the model to learn the iterative solving process. However, for cases involving large datasets and larger molecules, the DEQH model may be more appropriate. This not only offers the possibility of direct integration with off-the-shelf training methodologies that incorporate self-consistency but may also assist in enhancing aspects such as model generalizability.NeurIPS Paper Checklist 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We have provided all the information needed to reproduce experimental results in Supplementary Material A.11.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: The mean and variance of the experimental results are reported in Fig. 4 (a) and Table . 3.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: There is no societal impact of the work performed. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: All the data can be found publicly, and we open source our implementation at https://github.com/Zun-Wang/DEQHNet..• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA]Justification: The paper does not release new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: or [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. 
 
 Appendix:  


Decision: The paper introduces a novel neural architecture for quantum Hamiltonian prediction based on deep equilibrium models to model iterative algorithms. Benchmarks against competing methods show in general improvements in the MAE on Hamiltonian predictions. 
Some concerns highlighted by the reviews are around clarity of the presentation and the fact that better Hamiltonian prediction does not translate to better prediction of other molecular properties. The authors addressed these providing further experimental results during the discussion phase that helped understanding the performance of the model.

Reviews: review: paper_summary: paper_summary: This paper introduces the DEQH (Deep Equilibrium Quantum Hamiltonian) model, which integrates Deep Equilibrium Models (DEQs) for predicting quantum Hamiltonians. By incorporating DEQs, the model captures the self-consistency of Hamiltonians without needing iterative Density Functional Theory (DFT) calculations during training, enhancing computational efficiency. DEQHNet, a specific implementation, demonstrates improvements in prediction accuracy on the MD17 and QM9 datasets. The model acts as both a predictor and solver, iteratively refining the Hamiltonian to achieve self-consistency. Ablation studies further validate the effectiveness of this approach.
strength_weakness: strength_weakness: 1. The presentation of the paper could be improved. In this paper, the integration of Deep Equilibrium Models (DEQs) with the Hamiltonian solver is presented in a quite technical way. Additional visual illustrations and diagrams could help clarify the workflow of the DEQH model and its components. 

2. This paper primarily compares DEQHNet with QHNet. Benchmarking with additional methods could help better evaluate the results.
questions: questions: What are the runtime and memory usage for DEQHNet compared to previous methods?

What are the theoretical assumptions underlying the DEQH model, and how do they impact the generalizability of the results?
limitations: limitations: The limitations are discussed.
paper_summary: paper_summary: The authors introduce DEQH model, which combines deep equilibrium models with existing ML models to predict quantum Hamiltonian, and the author adopt QHNet as the backbone and further develop DEQHNet.  The authors evaluate the proposed method on benchmarks MD17 and QM9, the results show some effectiveness.
strength_weakness: strength_weakness: 1. Though I am relatively familiar with ML for Hamiltonian prediction and have gained some knowledge of DEQ, it is still not easy for me to understand the main part of this paper.

2. This work only compare with one baseline named QHNet, there are some other works which can also be baselines [1] [2] [3].

3. Source code is not available.

[1] Unke O, Bogojeski M, Gastegger M, et al. SE (3)-equivariant prediction of molecular wavefunctions and electronic densities[J]. Advances in Neural Information Processing Systems, 2021, 34: 14434-14447.

[2] Gong X, Li H, Zou N, et al. General framework for E (3)-equivariant neural network representation of density functional theory Hamiltonian[J]. Nature Communications, 2023, 14(1): 2848.

[3] Wang Y, Li H, Tang Z, et al. DeepH-2: Enhancing deep-learning electronic structure via an equivariant local-coordinate transformer[J]. arXiv preprint arXiv:2401.17015, 2024.
questions: questions: 1. Why is the model performance so much worse than the baseline in terms of energy? Has the author explored the reasons?

2. Is the method proposed by the author extensible and can it be applied to other existing ML models?
limitations: limitations: The authors have discussed the limitations.
paper_summary: paper_summary: The paper introduces a novel neural network architecture DEQH, extending deep equilibrium models (DEQs) to improve predictions of quantum Hamiltonians. The architecture constrains solutions to ensure self-consistency of the Hamiltonian, thereby improve generalization capability and test accuracy.
strength_weakness: strength_weakness: As mentioned, DEQH model acts fundamentally as a solver, iteratively determining the Hamiltonian with fixed point iterations. It is not directly clear to me why this would necessarily be better, or how this compares, to methods that integrating frameworks that also ensure self-consistency. Why would we expect to be better both computationally and in terms of performance, or is there a clear trade-off?
questions: questions: How do we expect DEQH to compare to methods that ensure self-consistency through integration? Is there are trade-off between performance and computational cost, or do we expect the method to be better overall. Why?
limitations: limitations: The paper proposes a novel way to embed self-consistency into neural network architectures for quantum Hamiltonian prediction and provides empirical evidence of benefits.
paper_summary: paper_summary: Paper presents DEQH model which is a Deep Equilibrium Model to solve quantum Hamiltonians. DEQH uses the deep equilibrium model because it converges to a fixed point, which matches the self-consistent nature of Hamiltonian solving. Paper presents architecture based on the QHNet, results and comparison with QHNet for the MD17 and QH9 datasets, and convergence and ablation study.
strength_weakness: strength_weakness: It seems there are differences between DEQH and QHNet, including the input of the predicted H, but the differences are not clearly highlighted in the paper. this makes the comparison harder to interpret.

The result of H, Psi errors being low while epsilon errors are high for DEQH is strange. What is the meaning of H, epsilon, and Psi, and why is the epsilon error higher in many cases? No error bars on tables and ablation study, makes results harder to interpret.
questions: questions: see above.
limitations: limitations: A limitation I see is that DFT ground truth calculation is still used in the loss. Since DEQH model is a Hamiltonian solver, it would be nice to see how it performs using a self-consistency loss. It seems that avoiding having to run many DFT calculations would be a big improvement, whereas the current results of paper show DEQH gets lower error results than previous model but still needs DFT data and requires 1.5x training time.

Comments: comment: We are very grateful for the increased score and deeply appreciate your valuable feedback.
I just raised my score. Thanks!
We are delighted that our response has met with your satisfaction. We will make sure to incorporate all of the latest discussions and experimental results into the manuscript and enhance its overall presentation promptly once the opportunity for a final polish is granted. Your constructive feedback is highly appreciated.
I thank the author’s for the further explanation and clarifications. I keep my recommendation for acceptance with the score of a 7, conditioned on including discussions and additional experiments from the other rebuttals in the final manuscript.
We are pleased that our response has satisfied you and are very grateful for the increased score.

As outlined in our rebuttal, the DEQH functions as a Hamiltonian solver, learning the iterative solving process $H^* = f(H^*, Z, R)$, while the QHNet acts as a predictor with the formula $H = f(Z, R)$. In practice, if the dataset includes overlap matrices, the network input can also incorporate this information, treating it similarly to the Hamiltonian by transforming it into equivariant node features (even if it's absent from the dataset, the calculation of the overlap matrix is quite straightforward and can be generated during data preprocessing). The overlap matrix offers a wealth of detailed information and can be computed easily with $Z$ and $R$. Accordingly, the equations for the DEQH and the modified Hamiltonian predictor become $H^* = f(H^*, Z, R, S)$ and $H = f(Z, R, S)$, respectively, where $S$ is the overlap matrix. In line with this, we conducted several experiments in our ablation study, testing the QHNet ($H = f(Z, R)$), QHNet w/ S ($H = f(Z, R, S)$), DEQHNet ($H^* = f(H^*, Z, R, S)$), and DEQHNet w/o S ($H^* = f(H^*, Z, R)$) on the Uracil dataset. The results depicted in Figure 4 indicate that:
* Incorporating the overlap matrix as an input leads to a lower Mean Absolute Error (MAE) in predicting the Hamiltonian (DEQHNet's Hamiltonian MAE is lower than DEQHNet w/o S, and QHNet's Hamiltonian MAE is higher than QHNet w/ S). Furthermore, the inclusion of the overlap matrix results in higher similarity in the orbital coefficients, suggesting that the overlap matrix is beneficial for the network's learning of the Hamiltonian and also supports our experimental hypothesis that the overlap matrix provides orbital information, which in turn improves the results for the orbital coefficients.
* Regardless of the presence of the overlap matrix in the network's input, the DEQH's Hamiltonian MAE is consistently lower than that of the QHNet-based Hamiltonian MAE, indicating that the DEQH model benefits from the introduction of self-consistency.

We will ensure that all new experimental results are included in the manuscript and will promptly refine our paper as soon as we are permitted a polish for the final version. We deeply appreciate your valuable feedback.
Thank the authors for response. I still think the ablation study is hard to interpret, but I appreciate the additional comparison with other benchmarks. I have increased score accordingly
We are glad to know that our response is satisfactory to you. We will make sure to include all the new experimental results in the manuscript and improve the presentation as soon as we are permitted a polish for the final version. We are grateful for your valuable feedback.
Thanks a lot for taking the time and effort to answer my questions. I would be considering raising the score in the next two days.
Thank you for sharing your further comment. Nevertheless, we are afraid that your comment is based on a factual misunderstanding. The fact is, DEQHNet outperforms QHNet in all cases, including ethanol and malonaldehyde, with the only exception of water due to limited available labels (500; other molecules have 25,000 labels) as we have explained in the paper. Such results indicate that DEQHNet improves upon existing end-to-end Hamiltonian prediction methods universally.

We are not sure about the specific reason why you have the misunderstanding, but we would like to stress that the primary comparing metric should be Hamiltonian MAE, which directly measures the prediction accuracy of the models. Other metrics, i.e., orbital energy and orbital coefficients, are derived quantities from the Hamiltonian. As we mentioned in the rebuttal, the prevailing MAE metric for Hamiltonian has been found not perfectly aligning with metrics of these derived quantities, which has also been observed in the comparison between QHNet and PhiSNet. In the appendix of our paper and in the rebuttal, we have explained why our method performs well in terms of Hamiltonian error and orbital coefficients. We also pointed out that the reported orbital energy falls within the error range of the current Hamiltonian. Specifically, in Section A.9 of the Supplementary Material and the Question 1 section of the rebuttal, we discussed this discrepancy in detail. We randomly selected a molecule, added a symmetrized Gaussian noise matrix (Hamiltonian matrices need to be symmetric), and solved the corresponding generalized eigenvalue equation. As shown in Fig. 5, as the noise on the Hamiltonian increased, the error range of the orbital energy also increased. The orbital energy errors we reported for both QHNet and DEQHNet on the MD17 and QH9 datasets fall within the range shown in the figure, suggesting that the experimental results are reasonable.
 
Additionally, we can provide further evidence from the QHNet paper, which observed a similar phenomenon. In their study, the PhiSNet achieved a Hamiltonian error of only 15.67x10$^{-6}$  Hartree on water, but the orbital energy error was "ten times worse than that by QHNet (33.76)", approximately  330x10$^{-6}$  Hartree. 

This observation implies that only when the Mean Absolute Error (MAE) of the Hamiltonian is sufficiently small can we ensure that the corresponding orbital energy deviation is also small. In the comparison between QHNet and PhiSNet, we observed similar situations. This somewhat suggests that using existing Hamiltonian error measurements cannot definitively reflect the errors in the eigenvalues. These evidence suggests that it is a nontrivial and long-standing challenge  in the domain of Hamiltonian prediction to design a proper metric on the matrix space so that reflects the metric of derived quantities e.g. energy, for which we will investigate in future work.

Finally, concerning your comment that "the presentation of the paper needs further improvement," we have already stated in our rebuttal that "We will also add further explanations about the integration of these two components and how this combination contributes to the overall model. We hope that these revisions will make our paper more comprehensible to readers with various backgrounds." The PDF we provided in the global rebuttal is merely an initial visualization, as the current NeurIPS review process only permits a one-page PDF that includes images and tables but no text.
Thank you for your efforts and detailed response. However, I still have some concerns about the additional experiments: DEQHNet performed poorly on water, not even surpassing QHNet; its performance on ethanol was also not better than QHNet; and on Malonaldehyde, it showed no advantage compared to the baseline. It only had some advantage on uracil.

Additionally, although the authors provided visualizations to demonstrate the method, I still believe that the presentation of the paper needs further improvement (the Reviewer eJGc also mentioned this point).

In summary, although the idea of DEQHNet seems interesting, its presentation and experimental performance are not convincing enough, so I will maintain my original rating.