Fermionic neural-network states for ab-initio electronic structure

Neural-network quantum states have been successfully used to study a variety of lattice and continuous-space problems. Despite a great deal of general methodological developments, representing fermionic matter is however still early research activity. Here we present an extension of neural-network quantum states to model interacting fermionic problems. Borrowing techniques from quantum simulation, we directly map fermionic degrees of freedom to spin ones, and then use neural-network quantum states to perform electronic structure calculations. For several diatomic molecules in a minimal basis set, we benchmark our approach against widely used coupled cluster methods, as well as many-body variational states. On some test molecules, we systematically improve upon coupled cluster methods and Jastrow wave functions, reaching chemical accuracy or better. Finally, we discuss routes for future developments and improvements of the methods presented.

P redicting the physical and chemical properties of matter from the fundamental principles of quantum mechanics is a central problem in modern electronic structure theory. In the context of ab-initio quantum chemistry (QC), a commonly adopted strategy to solve for the electronic wave-function is to discretize the problem on finite basis functions, expanding the full many-body state in a basis of anti-symmetric Slater determinants. Because of the factorial scaling of the determinant space, exact approaches systematically considering all electronic configurations, such as the full configuration interaction (FCI) method, are typically restricted to small molecules and basis sets. A solution routinely adopted in the field is to consider systematic corrections over mean-field states. For example, in the framework of the coupled cluster (CC) method 1,2 , higher level of accuracy can be obtained considering electronic excitations up to doublets, in CCSD, and triplets in CCSD(T). CC techniques are routinely adopted in QC electronic calculations, and they are often considered the "gold standard" in ab-initio electronic structure. Despite this success, the accuracy of CC is intrinsically limited in the presence of strong quantum correlations, in turn restricting the applicability of the method to regimes of relative weak correlations.
For strongly correlated molecules and materials, alternative, non-perturbative approaches have been introduced. Most notably, both stochastic and non-stochastic methods based on variational representations of many-body wave-functions have been developed and constantly improved in the past decades of research. Notable variational classes for QC are Jastrow-Slater wave-functions 3 , correlated geminal wave-functions 4 , and matrix product states [5][6][7] . Stochastic projection methods systematically improving upon variational starting points are for example the fixed-node Green's function Monte Carlo 8 and constrained-path auxiliary field Monte Carlo 9 . Main limitations of these methods stem, directly or indirectly, from the choice of the variational form. For example, matrix-product states are extremely efficient in quasi-one-dimensional systems, but suffer from exponential scaling when applied to larger dimensions. On the other hand, variational forms considered so-far for higher dimensional systems typically rely on rigid variational classes and do not provide a systematic and computationally efficient way to increase their expressive power.
To help overcome some of the limitations of existing variational representations, ideas leveraging the power of artificial neural networks (ANN) have recently emerged in the more general context of interacting many-body quantum matter. These approaches are typically based on compact, variational parameterizations of the many-body wave-function in terms of ANN 10 . These approaches to fermionic problems are however comparatively less explored than for lattice spin systems. Two main conceptually different implementations have been put forward. In the first, fermionic symmetry is encoded directly at the mean field level, and ANNs are used as a positive-definite correlator function 11 . Main limitation of this ansatz is that the nodal structure of the wave function is fixed, and the exact ground state cannot, in principle, be achieved, even in the limit of infinitely large ANN. The second method is to use ANNs to indirectly parameterize and modify the fermionic nodal structure [12][13][14][15] . In this spirit, "backflow" variational wave functions 16,17 with flexible symmetric orbitals have been introduced 13,14 , and only very recently applied to electronic structure 18,19 .
In this article, we provide an alternative representation of fermionic many-body quantum systems based on a direct encoding of electronic configurations. This task is achieved by mapping the fermionic problem onto an equivalent spin problem, and then solving the latter with spin-based neural-network quantum states. Using techniques from quantum information, we analyze different model agnostic fermion-to-spin mappings. We show results for several diatomic molecules in minimal Gaussian basis sets, where our approach reaches chemical accuracy (<1 kcal/mol) or better. The current challenges in extending the method to larger basis sets and molecules are also discussed.

Results
Electronic structure on spin systems. We consider many-body molecular fermionic Hamiltonians in second quantization formalism, where we have defined fermionic annihilation and creation operators with the anticommutation relation fc y i ; c j g ¼ δ i;j on N fermionic modes, and one-and two-body integrals t ij and u ijkm . The Hamiltonian in Eq. (1) can be mapped to interacting spin models via the Jordan-Wigner 20 mapping, or the more recent parity or Bravyi-Kitaev 21 encodings, which have been developed in the context of quantum simulations. These three encodings can all be expressed in the compact form where we have defined an update U(j), parity P(j), and remainder R(j) sets of spins, which depend on the particular mapping considered 22,23 , and σ ðx;y;zÞ i denote Pauli matrices acting on site i. In the familiar case of the Jordan-Wigner transformation, the update, parity, and remainder sets become U(j) = j, P(j) = {0, 1, ... j − 1}, R(j) = P(j), and the mapping takes the simple form where σ þðÀÞ j ¼ ðσ x j þ ðÀÞiσ y j Þ=2. For all the spin encodings considered, the final outcome is a spin Hamiltonian with the general form defined as a linear combination with real coefficients h j of σ j , Nfold tensor products of single-qubit Pauli operators I, σ x , σ y , σ z .
Additionally, under such mappings, there is a one to one correspondence between spin configuration σ ! and the original particle occupations n ! σ . In the following, we will consider the interacting spin Hamiltonian in Eq. (4) as the starting point for our variational treatment.
Neural-network quantum states. Once the mapping is performed, we use neural-network quantum states (NQS) introduced in ref. 10 to parametrize the ground state of the Hamiltonian in Eq. (4). One conceptual interest of NQS is that, because of the flexibility of the underlying non-linear parameterization, they can be adopted to study both equilibrium 24,25 and out-ofequilibrium [26][27][28][29][30][31] properties of diverse many-body quantum systems. In this work, we adopt a simple neural-network parameterization in terms of a complex-valued, shallow restricted Boltzmann machine (RBM) 10,32 . For a system of N spins, the many-body amplitudes take the compact form Here, W are complex-valued network parameters W ¼ fa; b; Wg, and the expressivity of the network is determined by the hidden unit density defined by α = M/N where M is the number of hidden units. The simple RBM ansatz can efficiently support volume-law entanglement [33][34][35][36] , and it has been recently used in several applications 37 .
One can then train the ansatz given in Eq. (5) with a variational learning approach known as variational Monte Carlo (VMC), by minimizing the energy expectation value This expectation value can be evaluated using Monte Carlo sampling using the fact that the energy (and, analogously, any other observable) can be written as where we have defined the local energy Given samples M drawn from the distribution jΨ M ð σ the average over the samplesÊðWÞ ¼ E loc ð σ ! Þ M gives an unbiased estimator of the energy. Note that the computational cost of evaluating the local energy depends largely on the sparsity of the Hamiltonian H q . In generic QC problems, this cost scales in the worst case with OðN 4 Þ, as compared to the linear scaling in typical condensed matter systems with local interaction.
Sampling from jΨ M ð σ ! Þj 2 is performed using Markov chain constructed using the Metropolis-Hastings algorithm 38 . Specifically, at each iteration, a configuration σ ! prop is proposed and accepted with probability The sample M then corresponds to the configurations of the Markov chain downsampled at an interval K, i.e., For the simulations done in this work, we typically use K = 10N with a sample size of approximately 100,000.
Since the Hamiltonians we are interested in have an underlying particle conservation law, it is helpful to perform this sampling in the particle basis n ! σ rather than the corresponding spin basis σ ! .
The proposed configuration σ ! prop at each iteration, then corresponds to a particle hopping between orbitals. Once a stochastic estimate of the expectation values is available, as well as its derivatives w.r.t. the parameters W, the ansatz can be optimized using the stochastic reconfiguration method 39,40 , closely related to the natural-gradient method used in machine learning applications 10,41 . Potential energy surfaces. We first consider small molecules in a minimal basis set (STO-3G). We show in Fig. 1 the dissociation curves for C 2 and N 2 , compared to the CCSD and CCSD(T). It can be seen that on these small molecules in their minimal basis, the RBM is able to generate accurate representations of the ground states, and remarkably achieve an accuracy better than standard QC methods. To further illustrate the expressiveness of the RBM, we show in Fig. 2 the probability distribution of the most relevant configurations in the wavefunction. We contrast between the RBM and configuration interaction limited to single a b Fig. 1 Dissociation profiles. The accuracy of fermionic neural-network quantum states compared with other quantum chemistry approaches. Shown here are dissociation curves for a C 2 and b N 2 , in the STO-3G basis with 20 spin-orbitals. The RBM used has 40 hidden units, and it is compared to both coupled-cluster approaches (CCSD, CCSD(T)) and FCI energies. and double excitations (CISD). In CISD, the Hilbert space is truncated to include only states which are up to two excitations away from the Hartree-Fock configuration. It is clear from the histogram that the RBM is able to capture correlations beyond double excitations.
Alternative encodings. The above computations were done using the Jordan-Wigner mapping. To investigate the effect of the mapping choice on the performance of the RBM, we also performed select calculations using the parity and Bravyi-Kitaev mappings. All the aforementioned transformations require a number of spins equal to the number of fermionic modes in the model. However, the support of the Pauli operators w j = |σ j | in Eq. (4), i.e., the number of single-qubit Pauli operators in σ j that are different from the identity I, depends on the specific mapping used. Jordan-Wigner and parity mappings have linear scalings w j = O(N), while the Bravyi-Kitaev encoding has a more favorable scaling w j ¼ Oðlog ðNÞÞ, due to the logarithmic spin support of the update, parity, and remainder sets in Eq. (2). Note that one could in principle use generalized superfast mappings 42 , which have a support scaling as good as w j ¼ Oðlog ðdÞÞ, where d is the maximum degree of the fermionic interaction graph defined by Eq. (1). However, such a mapping is not practical for the models considered here because the typical large degree of molecular interactions graphs makes the number of spins required for the simulation too large compared to the other model-agnostic mappings. While these encodings are routinely used as tools to study fermionic problems on quantum hardware 43 , their use in classical computing has not been systematically explored so far. Since they yield different structured many-body wave functions, it is then worth analyzing whether more local mappings can be beneficial for specific NQS representations. In Fig. 3, we analyze the effect of the different encodings on the accuracy of the variational groundstate energy for a few representative diatomic molecules. At fixed computational resources and network expressivity, we typically find that the RBM ansatz can achieve consistent levels of accuracy, independent of the nature of the mapping type. While the Jordan-Wigner allows to achieve the lowest energies in those examples, the RBM is nonetheless able to efficiently learn the ground state also in other representations, and chemical accuracy is achieved in all cases reported in Fig. 3.
Sampling larger basis sets. The spin-based simulations of the QC problems studied here show a distinctive MCMC sampling behavior that is not usually found in lattice model simulations of pure spin models. Specifically, the ground-state wave function of the diatomic molecules considered is typically sharply peaked around the Hartree-Fock state, and neighboring excited states. This behavior is prominently shown also in Fig. 2, where the largest peaks are several order of magnitude larger than the distribution tail. As a result of this structure, any uniform sampling scheme drawing states σ ! from the VMC distribution jΨ M ð σ ! Þj 2 , is bound to repeatedly draw the most dominant states, while only rarely sampling less likely configurations. To exemplify this peculiarity, we study the behavior of the ground state energy as a function of the number of MCMC samples used at each step of the VMC optimization. We concentrate on the water molecule in the larger 6-31g basis. In this case, the Metropolis sampling scheme exhibits acceptance rates as low as 0.1% or less, as a  consequence of the presence of dominating states previously discussed.
In Fig. 4, we vary the sample size and also compare MCMC sampling with exact sampling. We can see that the accuracy of the simulation depends quite significantly on the sample size. The large number of samples needed in this case, together with a very low acceptance probability for the Metropolis-Hasting algorithm, directly points to the inefficiency of uniform sampling from jΨ M ð σ ! Þj 2 . At present, this represents the most significant bottleneck in the application of our approach to larger molecules and basis sets. This issue however is not a fundamental limitation, and alternatives to the standard VMC uniform sampling can be envisioned to efficiently sample less likely-yet important for chemical accuracy-states. Beyond sampling issues, representability is also a factor as can be seen from the inset of Fig. 4. Enough hidden units are required to capture the wavefunction accurately, however, with more hidden units optimization also becomes more challenging, thus finding an appropriate network architecture is also crucial.

Discussion
In this work, we have shown that relatively simple shallow neural networks can be used to compactly encode, with high precision, the electronic wave function of model molecular problems in quantum chemistry. Our approach is based on the mapping between the fermionic quantum chemistry molecular Hamiltonian and corresponding spin Hamiltonians. In turn, the ground state of the spin models can be conveniently modeled with standard variational neural-network quantum states. On model diatomic molecules, we show that a RBM state is able to capture almost the entirety of the electronic excitations, improving on routinely used approaches as CCSD(T) and the Jastrow ansatz (Table 1). Several future directions can be envisioned. The distinctive peaked structure of the molecular wave function calls for the development of alternatives to uniform sampling from the Born probability. These developments will allow to efficiently handle larger basis sets than the ones considered here. Second, our study has explored only a very limited subset of possible neural-network architectures. Most notably, the use of deeper networks might prove beneficial for complex molecular complexes. Another very interesting matter for future research is the comparison of different neural-network-based approaches to quantum chemistry. Contemporary to this work, approaches based on antisymmetric wave-functions in continuous space have been presented 18,19 . These have the advantage that they already feature a full basis set limit. However, the discrete basis approach has the advantage that boundary conditions and fermionic symmetry are much more easily enforced. As a consequence, simple-minded shallow networks can already achieve comparatively higher accuracy than the deeper and substantially more complex networks so-far adopted in the continuum case. On a different note, in a recent article 44 , the use of a unitary-coupled RBM applicable for noisy intermediate-scale quantum devices has been proposed and is also worth exploring.

Methods
Geometries for diatomic molecules. The equilibrium geometries for the molecules presented in this work were obtained from the CCCBDB database 45 . For convenience, we present them in Table 2.
Computing matrix elements. A crucial requirement for the efficient implementation of the stochastic variational Monte Carlo procedure to minimize the ground-state energy, is the ability to efficiently compute the matrix elements of the spin Hamiltonian h σ ! 0 jH q j σ ! i, appearing in the local energy, Eq. (9). Since H q is a sum of products of Pauli operators, the goal is to efficiently compute matrix elements of the form where σ ν i i denotes a Pauli matrix with ν = I, x, y, z acting on site i. Because of the structure of the Pauli operators, these matrix elements are non-zero only for a specific σ ! 0 such that where n y is the total number of σ y operators in the string of Pauli matrices. The basis set considered here is STO-3G, and the corresponding geometries are reported in the Methods section. Energies are reported in Hartrees and statistical uncertainty on RBM and Jastrow states energies are on the last reported digits. The RBM used has a hidden unit density α = 1 for all the molecules apart from C 2 and N 2 where we use α = 2.  Simulation details. The optimization follows the stochastic reconfiguration scheme as detailed in the supplementary material of ref. 10 . Given a variational ansatz Ψðfα k gÞ 2 C 2 n depending on parameters {α k }, the parameter update δα k is given by solution of the linear equation X k 0 hO y k O k 0 i À hO y k ihO k 0 i þ λδ kk 0 h i δα k 0 ¼ Àϵ hO y kĤ i À hO y k ihĤi where O k ¼ ∂ ∂α k log Ψðfα 0 k gÞ Â Ã are the logarithmic derivatives, ϵ is the step size and λ is the regularization parameter. For the simulations done in this paper, we take ϵ = 0.05 and λ = 0.01. The expectation values 〈⋯〉 are estimated with Markov chain Monte Carlo sampling as described in the main text.
The parameters of the RBM are initialized from a random normal distribution with a zero mean and a standard deviation of 0.05.

Data availability
The datasets generated during and/or analyzed during the current study are available from the authors on reasonable request.

Code availability
The code used in the current study is largely based on the open-sourced software NetKet 46 with some custom modifications, which will be made available from the authors upon reasonable request.