Towards a Foundation Model for Neural Network Wavefunctions
Deep neural networks have become a highly accurate and powerful wavefunction
ansatz in combination with variational Monte Carlo methods for solving the
electronic Schr\"odinger equation. However, despite their success and favorable
scaling, these methods are still computationally too costly for wide adoption.
A significant obstacle is the requirement to optimize the wavefunction from
scratch for each new system, thus requiring long optimization. In this work, we
propose a novel neural network ansatz, which effectively maps uncorrelated,
computationally cheap Hartree-Fock orbitals, to correlated, high-accuracy
neural network orbitals. This ansatz is inherently capable of learning a single
wavefunction across multiple compounds and geometries, as we demonstrate by
successfully transferring a wavefunction model pre-trained on smaller fragments
to larger compounds. Furthermore, we provide ample experimental evidence to
support the idea that extensive pre-training of a such a generalized
wavefunction model across different compounds and geometries could lead to a
foundation wavefunction model. Such a model could yield high-accuracy ab-initio
energies using only minimal computational effort for fine-tuning and evaluation
of observables.
Michael Scherbela, Leon Gerard, Philipp Grohs
# Towards a Foundation Model for Neural Network Wavefunctions

###### Abstract

Deep neural networks have become a highly accurate and powerful wavefunction ansatz in combination with variational Monte Carlo methods for solving the electronic Schrodinger equation. However, despite their success and favorable scaling, these methods are still computationally too costly for wide adoption. A significant obstacle is the requirement to optimize the wavefunction from scratch for each new system, thus requiring long optimization. In this work, we propose a novel neural network ansatz, which effectively maps uncorrelated, computationally cheap Hartree-Fock orbitals, to correlated, high-accuracy neural network orbitals. This ansatz is inherently capable of learning a single wavefunction across multiple compounds and geometries, as we demonstrate by successfully transferring a wavefunction model pre-trained on smaller fragments to larger compounds. Furthermore, we provide ample experimental evidence to support the idea that extensive pre-training of a such a generalized wavefunction model across different compounds and geometries could lead to a foundation wavefunction model. Such a model could yield high-accuracy ab-initio energies using only minimal computational effort for fine-tuning and evaluation of observables.
## 1 Introduction

Accurate predictions of quantum mechanical properties for molecules is of utmost importance for the development of new compounds, such as catalysts, or pharmaceuticals. For each molecule the solution to the Schrodinger equation yields the wavefunction and electron density, and thus in principle gives complete access to its chemical properties. However, due to the curse of dimensionality, computing accurate approximations to the Schrodinger equation quickly becomes computationally intractable with increasing number of particles. Recently, deep-learning-based Variational Monte Carlo (DL-VMC) methods have emerged as a high-accuracy approach with favorable scaling \(\mathcal{O}(N^{4})\) in the number of particles \(N\)[1]. These methods use a deep neural network as ansatz for the high-dimensional wavefunction, and minimize the energy of this ansatz to obtain the ground-state wavefunction. Based on two major architectures for the treatment of molecules in first quantization, PauliNet [1] and FermiNet [2], several improvements and applications have emerged. On the one hand, enhancements of architecture, optimization and overall approach have led to substantial improvements in accuracy or computational cost [3, 4, 5, 6, 7]. On the other hand, these methods have been adapted to many different systems and observables: model systems of solids [8, 9], real solids [10], energies and properties of individual molecules [1, 2, 11, 5], forces [12, 13], excited states [14] and potential energy surfaces [13, 15, 16]. Furthermore, similar methods have been developed and successfully applied to Hamiltonians in second quantization [17, 18].

We want to emphasize that DL-VMC is an ab-initio method, that does not require any input beyond the Hamiltonian, which is defined by the molecular geometry. This differentiates it from surrogate models, which are trained on results from ab-initio methods to either predict wavefunctions [19, 20] or observables [21].

Despite the improvements in DL-VMC, it has not yet been widely adopted, in part due to the high computational cost. While DL-VMC offers favorable scaling, the method suffers from a large prefactor, caused by an expensive optimization with potentially slow convergence towards accurate approximations. Furthermore this optimization needs to be repeated for every new system, leading to prohibitively high computational cost for large-scale use. This can be partially overcome by sharing a single ansatz with identical parameters across different geometries of a compound, allowing more efficient computation of Potential Energy Surfaces (PES) [13, 15, 16]. However, these approaches have been limited to different geometries of a single compound and do not allow successful transfer to new compounds. A key reason for this limitation is that current architectures explicitly depend on the number of orbitals (and thus electrons) in a molecule. Besides potential generalization issues, this prevents a transfer of weights between different compounds already by the fact that the shapes of weight matrices are different for compounds of different size.

In this work we propose a novel neural network ansatz, which does not depend explicitly on the number of particles, allowing to optimize wavefunctions across multiple compounds with multiple different geometric conformations. We find, that our model exhibits strong generalization when transferring weights from small molecules to larger, similar molecules. In particular we find that our method achieves high accuracy for the important task of relative energies. Inspired by the success of foundation models in language [22] or vision [23, 24] - which achieve
high accuracy with minimal fine-tuning of an extensively pre-trained base-model - we train a first base-model for neural network wavefunctions.

We evaluate our pre-trained wavefunction model by performing few-shot predictions on chemically similar molecules (in-distribution) and disparate molecules (out-of-distribution). We find that our ansatz outperforms conventional high-accuracy methods such as CCSD(T)-ccpVTZ and that fine-tuning our pre-trained model reaches this accuracy \(\approx\)20x faster, than optimizing a new model. When analyzing the accuracy as a function of pre-training resources, we find that results systematically and substantially improve by scaling up either the model size, data size or number of pre-training steps. These results could pave the way towards a foundation wavefunction model, to obtain high-accuracy ab-initio results of quantum mechanical properties using only minimal computational effort for fine-tuning and evaluation of observables.

Additionally we compare our results to GLOBE, a concurrent work [25], which proposes reparameterization of the wavefunction based on machine-learned, localized molecular orbitals. We find that our method in comparison achieves lower absolute energies, higher accuracy of relative energies and is better able to generalize across chemically different compounds.

## 2 Results

In the following, we briefly outline our approach and how it extends existing work in Sec. 2.1. We show the fundamental properties of our ansatz such as extensivity (Sec. 2.3) and equivariance with respect to the sign of reference orbitals (Sec. 2.4). We demonstrate the transferability of the ansatz when pre-training on small molecules and re-using it on larger, chemically similar molecules. We also compare its performance against GLOBE, a concurrent pre-print [25] in Sec. 2.5. Lastly, we present a first wavefunction base-model pre-trained on a large diverse dataset of 360 geometries and evaluate its downstream performance in Sec. 2.6.

### A multi-compound wavefunction ansatz

Existing high-accuracy ansatze for neural network wavefunctions all exhibit the following structure:

\[\mathbf{h}_{i}=h_{\theta}(\mathbf{r}_{i},\{\mathbf{r}\},\mathbf{R},\mathbf{Z}) \tag{1}\] \[\Phi_{ik}^{d}=\varphi_{dk}(\mathbf{r}_{i})\sum_{\nu=1}^{D_{\rm emb}}F _{k\nu}^{d}h_{i\nu}\] (2) \[\psi=\sum_{d=1}^{N_{\rm det}}\det\left[\Phi_{ik}^{d}\right]_{i,k= 1\ldots n_{\rm el}} \tag{3}\]

Eq. 1 computes a \(D_{\rm emb}\)-dimensional embedding of electron \(i\), by taking in information of all other particles, e.g. by using attention or message passing. Eq. 2 maps these high-dimensional embeddings onto \(n_{\rm el}\times N_{\rm det}\) orbitals (indexed by \(k\)), using trainable backflow matrices \(\mathbf{F}^{d}\) and typically trainable envelope functions \(\varphi_{dk}\). Eq. 3 evaluates the final wavefunction \(\psi\) as a sum of (Slater-)determinants of these orbitals, to ensure antisymmetry with respect to permutation of electrons.

While this approach works well for the wavefunctions of a single compound, it is fundamentally unsuited to represent wavefunctions for multiple different compounds at once. The main problem lies in the matrices \(\mathbf{F}^{d}\), which are of shape \([n_{\rm el}\times D_{\rm emb}]\), and thus explicitly depend on the number of electrons. There are several potential options, how this challenge could be overcome. A
naive approach would be to generate a fixed number of \(N_{\mathrm{orb}}>n_{\mathrm{el}}\) orbitals and truncate the output to the required number of orbitals \(n_{\mathrm{el}}\), which may differ across molecules. While simple to implement, this approach is however fundamentally limited to molecules with \(n_{\mathrm{el}}\leq N_{\mathrm{orb}}\). Another approach is to use separate matrices \(\mathbf{F}_{\mathcal{G}}^{d}\) for each molecule or geometry \(\mathcal{G}\), as was done in [13], but also this approach can fundamentally not represent wavefunctions for molecules that are larger than the ones found in the training set. A third approach would be to not generate all orbitals in a single pass, but generate the orbitals sequentially in an auto-regressive manner, by conditioning each orbital on the previously generated orbitals. While this approach has been successful in other domains such as language processing, it suffers from inherently poor parallelization due to its sequential nature. A final approach - chosen in this work - is to replace the matrix \(\mathbf{F}\) with a trainable function \(f_{\theta}^{a}(\mathbf{c}_{Ik})\), which computes the backflows based on some descriptor \(\mathbf{c}_{Ik}\) of the orbital \(k\) to be generated:

\[h_{i\nu}=h_{\theta}(\mathbf{r}_{i},\{\mathbf{r}\},\{\mathbf{R}\},\{\mathbf{Z} \})_{\nu}\] \[\varphi_{\theta}^{d}(\mathbf{r}_{i},\mathbf{R}_{I},\mathbf{c}_{Ik})=\exp{(-| \mathbf{r}_{i}-\mathbf{R}_{I}|\,g_{\theta}^{s}(\mathbf{c}_{Ik})_{d})}\] \[\Phi_{ik}^{d}=\sum_{I=1}^{N_{\mathrm{muc}}}\varphi_{\theta}^{d}( \mathbf{r}_{i},\mathbf{R}_{I},\mathbf{c}_{Ik})\sum_{\nu=1}^{D_{\mathrm{emb}}}f_{\theta}^{ a}(\mathbf{c}_{Ik})_{d\nu}h_{i\nu}\]

While there are several potential descriptors \(\mathbf{c}_{Ik}\) for orbitals, one particularly natural choice is to use outputs of computationally cheap, conventional quantum chemistry methods such as Density Functional Theory or Hartree-Fock. We compute orbital features based on the expansion coefficients of a Hartree-Fock calculation, by using orbital localization and a graph convolutional network (GCN), as outlined in Sec 4.2. We then map these features to _transferable atomic orbitals (TAOs)_\(\Phi_{ik}^{d}\), using (anti-)symmetric functions \(f_{\theta}^{a}\) and \(g_{\theta}^{s}\) as illustrated in Fig. 1.

### Properties of our ansatz

These TAOs fulfil many properties, which are desirable for a wavefunction ansatz:

* **Constant parameter count**: The number of parameters in the ansatz is independent of system size. In previous approaches [1, 2, 13] the number of parameters grows with the number of

Figure 1: Illustration of the Transferrable Atomic Orbitals, demonstrated on the C=C-bond of Ethene.
particles, making it impossible to use a single ansatz across systems of different sizes. In particular backflows and envelope exponents have typically been chosen as trainable parameters of shape [\(N_{\mathrm{orb}}\times N_{\mathrm{det}}\)]. In our ansatz the backflows \(\mathbf{F}\) are instead computed by a single function \(f_{\theta}\) from multiple inputs \(\mathbf{c}_{Ik}\).
* **Equivariant to sign of HF-orbital**: Orbitals of a HF-calculation are obtained as eigenvectors of a matrix and are thus determined only up to their sign (or their phase in the case of complex eigenvectors). We enforce that the functions \(f_{\theta}^{a}\), \(g_{\theta}^{s}\) are (anti-)symmetric with respect to \(\mathbf{c}_{Ik}\). Therefore our orbitals \(\Phi_{ik}^{d}\) are equivariant to a flip in the sign of the HF-orbitals used as inputs: \(\Phi(-\mathbf{c}_{Ik})=-\Phi(\mathbf{c}_{Ik})\). Therefore during supervised pre-training, the undetermined sign of the reference orbitals becomes irrelevant, leading to faster convergence as demonstrated in Sec. 2.4.
* **Locality**: When using localized HF-orbitals as input, the resulting TAOs are also localized. Localized HF-orbitals are orbitals which have non-zero orbital features \(\widetilde{\mathbf{c}}_{Ik}\) only on some subset of atoms. Since we enforce the backflow \(f_{\theta}^{a}\) to be antisymmetric (and thus \(f^{a}(\mathbf{0})=\mathbf{0}\)), the resulting TAOs have zero contribution from atoms \(I\) with \(\mathbf{c}_{Ik}=\mathbf{0}\).
* **High expressivity**: We empirically find that our ansatz is sufficiently expressive to model ground-state wavefunctions to high accuracy. This is in contrast to previous approaches which were based on incorporating ab-initio orbitals [1], which could not reach chemical accuracy even for small molecules.

### Size consistency of the ansatz

One design goal of the ansatz is to allow transfer of weights from small systems to larger systems. In particular, if a large system consists of many small previously seen fragments, one would hope to obtain an energy which corresponds approximately to the sum

Figure 2: Zero-shot transferability of the ansatz to chemically similar, larger systems. While Moon cannot successfully transfer to larger chains, our ansatz successfully predicts zero-shot energies for up to 2x larger molecules and achieves very high accuracy with little pre-training.
of the fragment energies. One simple test case, are chains of equally spaced Hydrogen atoms of increasing lengths. These systems have been studied extensively using high-accuracy methods [26], because they are small systems which already show strong correlation and are thus challenging to solve. We test our method by pre-training our ansatz on chains of length 6 and 10, and then evaluating the model (with and without subsequent fine-tuning) for chain lengths between 2 and 28. Fig. 2 shows that our ansatz achieves very high zero-shot-accuracy in the interpolation regime (\(n_{\mathrm{atoms}}\) = 10) and for extrapolation to slightly larger or slightly smaller chains (\(n_{\mathrm{atoms}}\) = 4, 12). Even when extrapolating to systems of twice the size (\(n_{\mathrm{atoms}}\) = 20), our method still outperforms a Hartree-Fock calculation and eventually converges to an energy close the the Hartree-Fock solution. Fine-tuning the pre-trained model for only 500 steps, yields near perfect agreement with the specialized MRCI+Q method.

This good performance stands in stark contrast to other approaches such as GLOBE+ FermiNet or GLOBE+Moon, studied in [25]: Both GLOBE-variants yield 5-6x higher errors in the interpolation regime and both converge to much higher energies for larger chains. While our approach yields Hartree-Fock-like energies for very long chains, GLOBE+ FermiNet and GLOBE+Moon yield results that are outperformed even by assuming a chain of non-interacting H-atoms, which would yield an energy per atom of -0.5 Ha. For modest extrapolations (\(n\) = 12 to \(n\) = 20) our zero-shot results yield 3 - 20x lower errors than GLOBE+Moon.

### Equivariance with respect to HF-phase

Due to the (anti-)symmetrization of the TAOs, our orbitals are equivariant with respect to a change of sign of the Hartree Fock orbitals. Therefore, a sign change of the HF-orbitals during HF-pre-training has no effect on the

Figure 3: Accuracy when HF-pre-training against rotated H\({}_{2}\)O molecules, which contain a change of sign in the Hartree-Fock-p-orbitals of the Oxygen atom. Comparing a shared optimization of a backflow-based neural network wavefunction (Standard backflow) against TAOs. Top: HF-pre-training loss averaged over 20 geometries and the last 100 samples for each rotation angle. Bottom: Mean energy error vs. reference, averaged across all geometries.
optimization of the wavefunction. One test case to assure this behaviour is the rotation of a H\({}_{2}\)O molecule, where we consider a set of 20 rotations of the same geometry, leading to a change of sign in the p-orbitals of the Oxygen atom (cf. Fig. 3). We evaluate our proposed architecture and compare it against a naive approach, where we use a standard backflow matrix \(\boldsymbol{F}\), instead of a trainable, anti-symmetrized function \(f_{\theta}^{a}\). In Fig. 3 we can see a clear spike in the HF-pre-training loss at the position of the sign flip for the standard backflow-type architecture, causing slower convergence during the subsequent variational optimization. Although, in this specific instance the orbital sign problem could also be overcome, without our approach by correcting the phase of each orbital to align them across geometries, phase alignment is not possible in all circumstances. For example, there are geometry trajectories, where the Berry phase prevents such solutions [27].

### Transfer to larger, chemically similar compounds

To test the generalization and transferability of our approach, we perform the following experiment: First, we train our ansatz on a dataset of multiple geometries of a single, small compound (e.g. 20 distorted geometries of Methane). For this training, we follow the usual procedure of supervised HF-pre-training and subsequent variational optimization as outlined in Sec. 4.1. After 64k variational optimization steps, we then re-use the weights for different geometries of a larger compound (e.g. distorted geometries of Etiene). We fine-tune the model on this new geometry dataset for a variable number of steps and plot the resulting energy errors in Fig. 4. We do not require supervised HF-pre-training on the new, larger dataset. We perform this experiment for 3 pairs of test systems: Transferring from geometries of Hydrogen-chains with 6 atoms each, to chains with 10 atoms each, transferring from Methane to Etiene, and transferring from Etiene to Cyclobutadiene.

We compare our results to the earlier DeepErwin approach [13], which only partially reused weights, and GLOBE, a concurrent preprint [25] which reuses all weights. To measure accuracy we compare two important metrics: First, the mean energy error (averaged across all geometries \(g\) of the test dataset) \(\frac{1}{N}\sum_{g}(E_{g}-E_{g}^{\mathrm{ref}})\), which reflects the method's accuracy for absolute energies. Second, the maximum relative energy error \(\max_{g}(E_{g}-E_{g}^{\mathrm{ref}})-\min_{g}(E_{g}-E_{g}^{\mathrm{ref}})\), which reflects the method's consistency across a potential energy surface. Since different studies use different batch-sizes and different definitions of an epoch, we plot all results against the number of MCMC-samples used for variational optimization, which is very closely linked to computational cost.

Compared to other approaches, we find that our method yields substantially lower and more consistent energies. On the toy problem of H\({}_{6}\) to H\({}_{10}\) our approach and GLOBE reach the same accuracy, while DeepErwin converges to higher energies. For the actual molecules Etiene (C\({}_{2}\)H\({}_{4}\)) and Cyclobutadiene (C\({}_{4}\)H\({}_{4}\)) our approach reaches substantially lower energies and much more consistent potential energy surfaces. When inspecting the resulting Potential Energy Surface for Etiene, we find that we obtain qualitatively similar results as DeepErwin, but obtain energies that are \(\approx\) 6 mHa lower (and thus more accurate). GLOBE on the other hand does not yield a qualitatively correct PES for this electronically challenging problem. It overestimates the energy barrier at 90\({}^{\circ}\) twist angle by \(\approx\) 50 mHa and yields a spurious local minimum at 10\({}^{\circ}\). We
observe similar results on the Cyclobutadiene geometries, where our approach yields energy differences that are in close agreement to the reference energies, while the GLOBE-results overestimate the energy difference by \(\approx\) 20 mHa.

### Towards a first foundation model for neural network wavefunctions

While the experiments in Sec. 2.5 demonstrate the ability to pre-train our model and fine-tune it on a new system, the resulting pre-trained models are of little practical use, since they are only pre-trained on a single compound each and can thus not be expected to generalize to chemically different systems. To obtain a more diverse pre-training dataset, we compiled

Figure 4: Accuracy when pre-training the model on small compounds and reusing it for larger compounds. Top: Mean energy error vs reference, averaged across all geometries of the test set. Middle: Error of relative energies, measured as \(\max_{g}(E_{g}-E_{g}^{\mathrm{ref}})-\min_{g}(E_{g}-E_{g}^{\mathrm{ref}})\). Bottom: Final PES for the Ethene molecule for each method.
a dataset of 360 distorted geometries, spread across 18 different compounds. The dataset effectively enumerates all chemically plausible molecules with up to 18 electrons containing the elements H, C, N, and O. For details on the data generation see Appendix A. We pre-train a base-model for 500,000 steps on this diverse dataset and subsequently evaluate its performance, when computing Potential Energy Surfaces. We evaluate its error both for compounds that were in the pre-training dataset (with different geometries), as well as for new, larger, out-of-distribution compounds which were not present in the pre-training dataset. We compare the results against a baseline model, which uses the default method of supervised HF-pre-training and subsequent variational optimization.

Fig. 5 shows that fine-tuning this pre-trained model yields substantially lower energies than the usual optimization from a HF-pre-trained model. For example, when optimizing for 8k steps, we obtain 8x lower energy errors for large out-of-distribution compounds, and 12x lower energies for small in-distribution compounds. When evaluating the model for up to 32k steps, we find that for small molecules both approaches converge to the same energy. For large molecules the final energy error obtained by fine-tuning the base-model is 3x lower than optimization of a HF-pre-trained model.

### Scaling behaviour

In many domains, increasing the amount of pre-training, has led to substantially better results, even without qualitative changes to the architecture [28]. To investigate the scalability of our approach, we vary the three key choices, along which one could increase the scale of pre-training: The size of the wavefunction model, the number of compounds and geometries present in the pre-training-dataset, and the number of pre-training steps. Starting from a large model, trained on 18x20 geometries, for 256k pre-training steps, we independently vary each parameter. We test 3 different architectures sizes, with decreasing layer width and depth for the networks \(f_{\theta}\), \(g_{\theta}\), and GCN\({}_{\theta}\) (cf. Appendix C). We test 3 different training sets, with decreasing number of compounds in the training set, with 20 geometries each (cf. Appendix A). Finally, we evaluate model-checkpoints at different amounts of pre-training, ranging from 64k steps to 512k steps. Fig. 6 depicts the accuracy obtained by subsequently fine-tuning the resulting model for just 4000 steps on the evaluation set. In each case, increasing the scale of pre-training clearly

Figure 5: Fine-tuning of pre-trained foundation model (solid lines) vs. fine-tuning Hartree-Fock-pre-trained models (dashed lines) for 70 different geometries. Blue: Small compounds, with geometries similar to geometries in pre-training dataset. Orange: Larger compounds outside pre-training dataset.
improves evaluation results - both for the small in-distribution compounds, as well as the larger out-of-distribution compounds. We find a strong dependence of the accuracy on the model size and number of compounds in the pre-training dataset, and a weaker dependency on the number of pre-training steps. While our computational resources, currently prohibit us from training at larger scale, the results indicate that our approach may already be sufficient to train an accurate multi-compound, multi-geometry foundation model for wavefunctions.

## 3 Discussion

This work presents an ansatz for deep-learning-based VMC, which can in principle be applied to molecules of arbitrary size. We demonstrate the favourable properties of our ansatz, such as extensivity, zero-shot prediction of wavefunctions for similar molecules (Sec. 2.3), invariance to the phase of orbitals (Sec. 2.4) and fast finetuning for larger, new molecules (Sec. 2.5). Most importantly, Sec. 2.6 is, to our knowledge, the first successful demonstration of a general wavefunction, which has successfully been trained on a diverse dataset of compounds and geometries. We demonstrate that the dominating deep-learning paradigm of the last years - pre-training on large data and fine-tuning on specific problems - can also be applied to the difficult problem of wavefunctions. While previous attempts [13, 25] have failed to obtain high-accuracy energies from pre-trained neural network wavefunctions, we find that our approach yields accurate energies and does so at a fraction of the cost needed without pre-training. We furthermore demonstrate in Sec. 2.7 that results can be improved systematically by scaling up any aspect of the pre-training: Model size, data-size, or pre-training-steps.

Figure 6: Error when fine-tuning the pre-trained model for 4000 steps on small in-distribution geometries and larger out-of-distribution geometries.
Despite these promising results, there are many open questions and limitations which should be addressed in future work. First, we find that our ansatz currently does not fully match the accuracy of state-of-the-art single-geometry DL-VMC ansatze. While our approach consistently outperforms conventional variational methods such as MRCI or CCSD at finite basis set, larger, computationally more expensive DL-VMC models can reach even lower energies. Exchanging our message-passing-based electron-embedding, with recent attention based approaches [4] should lead to higher accuracy. Furthermore we have made several deliberate design choices, which each trade-off expressivity (and thus potentially accuracy) for computational cost: We do not exchange information across orbitals and we base our orbitals on computationally cheap HF-calculations. Including attention or message passing across orbitals (e.g. similar to [25]), and substituting HF for a trainable, deep-learning-based model should further increase expressivity. While we currently use HF-orbitals due to their widespread use and low computational cost, our method does not rely on a specific orbital descriptor. We could substitute HF for a separate model such as PhysNet [20] or SchnOrb [29] to compute orbital descriptors \(\mathbf{c}_{Ik}\), leading to a fully end-to-end machine-learned wavefunction. Second, while we include useful physical priors such as locality, we do not yet currently use the invariance of the Hamiltonian with respect to rotations, inversions or spin-flip. E3-equivariant networks have been highly successful for neural network force-fields, but have not yet been applied to wavefunctions due to the hitherto unsolved problem of symmetry breaking [15]. Using HF-orbitals as symmetry breakers, could open a direct avenue towards E3-equivariant neural network wavefunctions. Third, while we use locality of our orbitals as a useful prior, we do not yet use it to reduce computational cost. By enforcing sparsity of the localized HF-coefficients, one could limit the evaluation of orbitals to a few participating atoms, instead of all atoms in the molecule. While the concurrent GLOBE approach enforces its orbitals to be localized at a single position, our approach naturally lends itself to force localization at a given number of atoms, allowing for a deliberate trade-off of accuracy vs. computational cost. Lastly, we observe that our method performs substantially better, when dedicating more computational resources to the pre-training. While we are currently constrained in the amount of computational resources, we hope that future work will be able to scale up our approach. To facilitate this effort we open source our code, dataset as well as model parameters.

## 4 Methods

### Variational Monte Carlo

Considering the Born-Oppenheimer approximation, a molecule with \(n_{\mathrm{el}}\) electrons and \(N_{\mathrm{nuc}}\) nuclei can be described by the time-independent Schrodinger equation

\[\hat{H}\psi=E\psi \tag{4}\]

with the Hamiltonian

\[\hat{H}= -\frac{1}{2}\sum_{i}\nabla_{\mathbf{r}_{i}}^{2}+\sum_{i>j}\frac{1}{| \mathbf{r}_{i}-\mathbf{r}_{j}|}\] \[+\sum_{I>J}\frac{Z_{I}Z_{J}}{|\mathbf{R}_{I}-\mathbf{R}_{J}|}-\sum_{i,I} \frac{Z_{I}}{|\mathbf{r}_{i}-\mathbf{R}_{I}|} \tag{5}\]

By \(\mathbf{r}=(\mathbf{r}_{1},\dots,\mathbf{r}_{n_{\uparrow}},\dots,\mathbf{r}_{n_{\mathrm{el}}} )\in\mathbb{R}^{3\times n_{\mathrm{el}}}\) we denote the set of electron positions divided into \(n_{\uparrow}\) spin-up and \(n_{\downarrow}\) spin-down electrons. For the coordinates and charges of the nuclei we
write \(\mathbf{R}_{I}\), \(Z_{I}\), with \(I\in\{1,\ldots,N_{\mathrm{nuc}}\}\). The solution to the electronic Schrodinger equation \(\psi\) needs to fulfill the anti-symmetry property, i.e. \(\psi(\mathcal{P}\mathbf{r})=-\psi(\mathbf{r})\) for any permutation \(\mathcal{P}\) of two electrons of the same spin. Finding the groundstate wavefunction of a system, corresponds to finding the solution to Eq. 4, with the lowest eigenvalue \(E_{0}\). Using the Rayleigh-Ritz principle, an approximate solution can be found through minimization of the loss

\[\mathcal{L}(\psi_{\theta})=\mathbb{E}_{\mathbf{r}\sim\psi_{\theta}^{2}(\mathbf{r})} \left[\frac{(\hat{H}\psi_{\theta})(\mathbf{r})}{\psi_{\theta}(\mathbf{r})}\right]\geq E _{0}, \tag{6}\]

using a parameterized trial wavefunction \(\psi_{\theta}\). The expectation value in Eq. 6 is computed by drawing samples \(\mathbf{r}\) from the unnormalized probability distribution \(\psi_{\theta}^{2}(\mathbf{r})\) using Markov Chain Monte Carlo (MCMC). The effect of the Hamiltonian on the wavefunction can be computed using automatic differentiation and the loss is minimized using gradient based minimization. A full calculation typically consists of three steps:

1. **Supervised HF-pre-training**: Minimization of the difference between the neural network ansatz and a reference wavefunction (e.g. a Hartree-Fock calculation) \(||\psi_{\theta}-\psi^{\mathrm{HF}}||\). This is the only part of the procedure which requires reference data, and ensures that the initial wavefunction roughly resembles the true groundstate. While this step is in principle not required, it substantially improves the stability of the subsequent variational optimization.
2. **Variational optimization**: Minimization of the energy (Eq. 6) by drawing samples from the wavefunction using MCMC, and optimizing the parameters \(\theta\) of the ansatz using gradient based optimization.
3. **Evaluation**: Evaluation of the energy by evaluating Eq. 6 without updating the parameters \(\theta\), to obtain unbiased estimates of the energy.

To obtain a single wavefunction for a dataset of multiple geometries or compounds, only minimal changes are required. During supervised and variational optimization, for each gradient step we pick one geometry from the dataset. We pick geometries either in a round-robin fashion, or based on the last computed energy variance for that geometry. We run the Metropolis Hastings algorithm [30] for that geometry to draw electron positions \(\mathbf{r}\) and then evaluate energies and gradients. For each geometry we keep a distinct set of electron samples \(\mathbf{r}\).

### Obtaining orbital descriptors from Hartree-Fock

As discussed in Sec. 2.1, our ansatz effectively maps uncorrelated, low-accuracy Hartree-Fock orbitals, to correlated, high-accuracy neural network orbitals. The first step in this approach is to obtain orbital descriptors \(\mathbf{c}_{k}\) for each orbital \(k\), based on a Hartree-Fock calculation.

The Hartree-Fock method uses a single determinant as ansatz, composed of single-particle orbitals \(\phi_{k}\):

\[\psi^{\mathrm{HF}}(\mathbf{r}_{1},\ldots,\mathbf{r}_{n_{\mathrm{el}}})= \det\left[\Phi_{ik}^{\mathrm{HF}}\right]_{i,k=1\ldots n_{\mathrm{ el}}} \tag{7}\] \[\Phi_{ik}^{\mathrm{HF}}:= \phi_{k}^{\mathrm{HF}}(\mathbf{r}_{i}) \tag{8}\]

For molecules, these orbitals are typically expanded in atom-centered basis-functions \(\mu(\mathbf{r})\), with \(N_{\mathrm{basis}}\) functions centered on each atom \(I\):

\[\phi_{k}^{\mathrm{HF}}(\mathbf{r})=\sum_{I=1}^{N_{\mathrm{nuc}}}\sum_{b=1}^{N_{ \mathrm{basis}}}\alpha_{k,Ib}\;\mu_{b}(\mathbf{r}-\mathbf{R}_{I}), \tag{9}\]

The coefficients \(\mathbf{\alpha}_{k}\) and the corresponding orbitals \(\phi_{k}^{\mathrm{HF}}(\mathbf{r})\) are obtained as solutions of an
eigenvalue problem and are typically delocalized, i.e. they have non-zero contributions from many atoms. However, since \(\det[U\Phi]=\det[U]\det[\Phi]\), the wavefunction is invariant under linear combination of orbitals by a matrix \(U\) with \(\det[U]=1\). One can thus choose orbital expansion coefficients

\[\widetilde{\alpha}_{k,Ib}=\sum_{k^{\prime}=1}^{N_{\mathrm{orb}}}\alpha_{k,Ib}U_ {kk^{\prime}} \tag{10}\]

corresponding to orbitals which are maximally localized according to some metric. Several different metrics and corresponding localization schemes, such as Foster-Boys [31] or Pipek-Mezey [32], have been proposed and are easily available as post-processing options in quantum chemistry codes. We use the Foster-Boys method as implemented in pySCF [33].

Due to the fundamentally local nature of atom-wise orbital coefficients \(\widetilde{\boldsymbol{\alpha}}_{Ik}\), which can be insufficient to distinguish orbitals, we use a fully connected graph convolutional neural network (GCN) to add context about the surrounding atoms. We interpret each atom as a node (with node features \(\widetilde{\boldsymbol{\alpha}}_{Ik}\)) and the 3D inter-atomic distance vector \(\boldsymbol{R}_{IJ}\) as edge features:

\[\boldsymbol{c}_{k}=\mathrm{GCN}_{\theta}\left(\{\widetilde{\boldsymbol{\alpha }}_{Ik}\}_{I=1\ldots N_{\mathrm{nuc}}},\{\boldsymbol{R}_{IJ}\}_{I,J=1\ldots N _{\mathrm{nuc}}}\right)\]

We embed the edge features using a cartesian product of Gaussian basis functions of the distance \(R_{IJ}\) and the concatenation of the 3D-distance vector with the constant 1:

\[\widetilde{\boldsymbol{c}}_{IJ} =\exp\left(-\frac{(R_{IJ}-\mu_{n})^{2}}{2\sigma_{n}^{2}}\right) \otimes[1|\boldsymbol{R}_{IJ}]\] \[\boldsymbol{e}_{IJ}^{0} =\mathrm{MLP}(\widetilde{\boldsymbol{c}}_{IJ})\] \[\boldsymbol{c}_{I}^{0} =\boldsymbol{\alpha}_{I}\]

Each layer \(l\) of the GCN consist of the following update rules

\[\boldsymbol{u}_{Ik}^{l} =\sum_{J}\boldsymbol{c}_{Jk}^{l}\odot\left(\boldsymbol{W}_{g}^{ l}\boldsymbol{e}_{IJ}^{l}\right),\] \[\boldsymbol{c}_{Ik}^{l+1} =\sigma(\boldsymbol{W}_{v}^{l}\boldsymbol{c}_{Ik}+\boldsymbol{W}_ {u}^{l}\boldsymbol{u}_{Ik}^{l}),\]

with trainable weight matrices \(\boldsymbol{W}_{g}^{l}\), \(\boldsymbol{W}_{v}^{l}\), \(\boldsymbol{W}_{u}^{l}\) and the SiLU activation function \(\sigma\). After \(L\) iterations we use the final outputs as orbitals features:

\[\boldsymbol{c}_{Ik}:=\boldsymbol{c}_{Ik}^{L}\]

### Mapping orbital descriptors to wavefunctions

To obtain entry \(\Phi_{ik}\) of the Slater determinant, we combine a high-dimensional electron embedding \(\boldsymbol{h}_{i}\) with a function of the orbital descriptor \(\boldsymbol{c}_{Ik}\):

\[h_{i\nu}=h_{\theta}(\boldsymbol{r}_{i},\{\boldsymbol{r}\},\{ \boldsymbol{R}\},\{\boldsymbol{Z}\})_{\nu}\] \[\varphi_{\theta}^{d}(\boldsymbol{r}_{i},\boldsymbol{R}_{I}, \boldsymbol{c}_{Ik})=\exp\left(-|\boldsymbol{r}_{i}-\boldsymbol{R}_{I}|\,g_{ \theta}^{s}(\boldsymbol{c}_{Ik})_{d}\right)\] \[\Phi_{ik}^{d}=\sum_{I=1}^{N_{\mathrm{nuc}}}\varphi_{\theta}^{d}( \boldsymbol{r}_{i},\boldsymbol{R}_{I},\boldsymbol{c}_{Ik})\sum_{\nu=1}^{D_{ \mathrm{amb}}}f_{\theta}^{a}(\boldsymbol{c}_{Ik})_{d\nu}h_{i\nu}\]

The functions \(\mathrm{GCN}_{\theta}^{a}\), \(f_{\theta}^{a}\), and \(g_{\theta}^{s}\) are trainable functions, which are enforced to be (anti-)symmetric with respect to change in sign of their argument \(\boldsymbol{c}\):

\[\text{Symmetric }g_{\theta}^{s}\text{:}\] \[g_{\theta}^{s}(\boldsymbol{c}):=g_{\theta}(\boldsymbol{c})+g_{ \theta}(-\boldsymbol{c})\] \[\text{Antisymm. }f_{\theta}^{a}\text{:}\] \[f_{\theta}^{a}(\boldsymbol{c}):=f_{\theta}(\boldsymbol{c})-f_{ \theta}(-\boldsymbol{c})\] \[\text{Antisymm. }\mathrm{GCN}_{\theta}^{a}\text{:}\] \[\mathrm{GCN}_{\theta}^{a}(\boldsymbol{\alpha},\boldsymbol{R}):= \mathrm{GCN}_{\theta}(\boldsymbol{\alpha},\boldsymbol{R})-\mathrm{GCN}_{\theta }(-\boldsymbol{\alpha},\boldsymbol{R})\]

To obtain electron embeddings \(\boldsymbol{h}_{i}\) we use the message-passing architecture outlined in [5], which is invariant with respect to permutation of electrons of the same spin, or the permutation of ions.

\[\mathcal{G}=\{(\boldsymbol{R}_{I},Z_{I})\}_{I=1\ldots N_{\mathrm{ nuc}}}\] \[\mathcal{E}_{\uparrow}=\{\boldsymbol{r}_{i}\}_{i=1\ldots n_{ \uparrow}}\] \[\mathcal{E}_{\downarrow}=\{\boldsymbol{r}_{i}\}_{i=n_{\uparrow}+1 \ldots n_{\mathrm{el}}}\] \[\boldsymbol{h}_{i}=h_{\theta}^{\mathrm{embed}}(\boldsymbol{r}_{i },\mathcal{G},\mathcal{E}_{\uparrow},\mathcal{E}_{\downarrow})\]
Note that during training, all samples in a batch come from the same geometry, and thus have the same values for \(\mathbf{R}\), \(\mathbf{Z}\), and \(\widetilde{\mathbf{\alpha}}\). While the embedding network \(h^{\text{embed}}_{\theta}\), needs to be re-evaluated for every sample, the networks GCN\({}_{\theta}\), \(f_{\theta}\), and \(g_{\theta}\) only need to be evaluated once per batch, substantially reducing their impact on computational cost.

## Acknowledgements

We gratefully acknowledge financial support from the following grants: Austrian Science Fund FWF Project I 3403 (P.G.), WWTF-ICT19-041 (L.G.). The computational results have been achieved using the Vienna Scientific Cluster (VSC). The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. Additionally, we thank Nicholas Gao for providing his results and data, Ruard van Workum for initial work on the python implementation for multi-compound optimization and Jan Hermann for fruitful discussions.

## Author contributions

MS, LG, and PG conceived the overall idea. MS conceived and implemented the ansatz, built the dataset and designed the experiments. LG gave input on the ansatz and worked on implementation. MS and LG performed the experiments. MS and LG wrote the manuscript with input, supervision and funding from PG.
