Fermionic neural-network states for ab-initio electronic structure
Kenny Choo Antonio Mezzacapo Giuseppe Carleo




Department of Physics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
kenny.choo@uzh.ch

IBM T.J. Watson Research Center, Yorktown Heights, NY, USA
amezzac@us.ibm.com

gcarleo@flatironinstitute.org
Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA


03.65.Aa,03.67.-a,03.67.Ac


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 the test molecules, we recover almost the entirety of the correlation energy. We systematically improve upon coupled cluster methods and Jastrow wave functions, reaching levels of chemical accuracy or better. Finally, we discuss routes for future developments and improvements of the methods presented.

Fermionic neural-network states for ab-initio electronic structure
    Giuseppe Carleo
    December 30, 2023
==================================================================






  
*Introduction.- 

Predicting 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 <cit.>, 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 <cit.>, correlated geminal wave-functions <cit.>, and matrix product states <cit.>. Stochastic projection methods systematically improving upon variational starting points are for example the fixed-node Green's function Monte Carlo <cit.> and constrained-path auxiliary field Monte Carlo <cit.>. 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 <cit.>. These approaches to fermionic problems are however  comparatively less explored than for lattice spin systems. Two 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 <cit.>. 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 parametrize permutation symmetric many-body fermionic orbitals <cit.>, in the spirit of "backflow" many-body variational wave functions <cit.>, and only very recently applied to electronic structure <cit.>.  

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 (<5 kcal/mol) or better. The current challenges in extending the method to larger basis sets and molecules are also discussed.    



  
* Electronic structure on spin systems.-We consider many-body molecular fermionic Hamiltonians in second quantization formalism, 

    H=∑_i, j t_ij  c^†_i c_j +  ∑_i, j, k, m u_ijkm  c^†_i c^†_k c_m c_j,
 
where we have defined fermionic annihilation and creation operators with the anticommutation relation {c^†_i, c_j} = δ_i,j on N fermionic modes, and one- and two-body integrals t_ij and u_ijkm. The Hamiltonian (<ref>) can be mapped to interacting spin models via the Jordan-Wigner <cit.> mapping, or the more recent parity or Bravyi-Kitaev <cit.> encodings, which have been developed in the context of quantum simulations. These three encodings can all be expressed in the compact form

    c_j   →1/2∏_i∈ U(j)σ^x_i×(σ^x_j∏_i∈ P(j)σ^z_i-iσ^y_j∏_i∈ R(j)σ^z_i)
    
    c^†_j   →1/2∏_i∈ U(j)σ^x_i×(σ^x_j∏_i∈ P(j)σ^z_i+iσ^y_j∏_i∈ R(j)σ^z_i),

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 <cit.>, 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 

    c_j   →( ∏_i=0^j-1σ_i^z )  σ^-_j
    
    c_j^†   →( ∏_i=0^j-1σ_i^z )  σ^+_j,

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

    H_q =∑_j=1^r h_j  σ_j,

defined as a linear combination with real coefficients h_j of σ_j,
N-fold 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 (<ref>) as a starting point for our variational treatment. 



  
* Neural-network quantum states.-
Once the mapping is performed, we use neural-network quantum states (NQS) introduced in <cit.> to parametrize the ground state of the Hamiltonian (<ref>).   
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 <cit.> and out-of-equilibrium <cit.> 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) <cit.>. 
For a system of N spins, the many-body amplitudes take the compact form 


    Ψ_M(σ⃗; )    = e^∑_i a_i σ^z_i∏_j=1^M 2 coshθ_j(σ⃗), where
    θ_j(σ⃗)    = b_j + ∑_i^N W_ijσ^z_i.

Here,  are complex-valued network parameters = {a,b,W}, and the expressivity of the network is determined by the hidden unit density defined by α = M/N where M is number of hidden units. The simple RBM ansatz can efficiently support volume-law entanglement <cit.>, and it has been recently used in several applications <cit.>. 

One can then train the ansatz Eq.(<ref>) with a variational learning approach known as Variational Monte Carlo (VMC), by minimizing the energy expectation value

    E() = ⟨Ψ_M | H_q | Ψ_M ⟩/⟨Ψ_M|Ψ_M ⟩.

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

    E()= ∑_σ⃗ E_loc(σ⃗) |Ψ_M(σ⃗)|^2 /∑_σ⃗ |Ψ_M(σ⃗)|^2 ,

where we have defined the local energy 

    E_loc(σ⃗) = ∑_σ⃗'Ψ_M(σ⃗')/Ψ_M^*(σ⃗)⟨σ⃗' | H_q | σ⃗⟩.


Given samples  drawn from the distribution |Ψ_M(σ⃗)|^2/∑_σ⃗ |Ψ_M(σ⃗)|^2, the average over the samples Ê() = ⟨ E_loc(σ⃗)⟩_
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 a generic QC problems, this cost scales in the worst case with 𝒪(N^4), as compared to the linear scaling in typical condensed matter systems with local interaction.

Sampling from |Ψ_M(σ⃗)|^2 is performed using Markov chain Monte Carlo (MCMC), with a Markov chain σ⃗_0 →σ⃗_1 →σ⃗_2 →… constructed using the Metropolis-Hastings algorithm <cit.>. Specifically, at each iteration, a configuration σ⃗_prop is proposed and accepted with probability 

    P(σ⃗_k+1 = σ⃗_prop) = min(1, | Ψ_M(σ⃗_prop)/Ψ_M(σ⃗_k)|^2 ).

The sample  then corresponds to the configurations of the Markov chain downsampled at an interval K, i.e. {σ⃗_0, σ⃗_K, σ⃗_2K, …}. For the simulations done in this work, we typically use K=10N with sample size of approximately 100000.

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 , the ansatz can be optimized using the stochastic reconfiguration method <cit.>, closely related to the natural-gradient method used in machine learning applications <cit.>.







  
*Potential Energy surfaces.-
We first consider small molecules in a minimal basis set (STO-3G). We show in Fig. <ref> 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. <ref> the probability distribution of the most relevant configurations in the wavefunction. We contrast between the RBM and configuration interaction limited to single 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 (<ref>), 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 (<ref>). Note that one could in principle use generalized superfast mappings <cit.> ,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 (<ref>). 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 <cit.>, 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. <ref>, we analyze the effect of the different encodings on the accuracy of the variational ground-state 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. <ref>. 





  
*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. <ref>, where the largest peaks are several of order of magnitude larger than the distribution tail. As a result of this structure, any uniform sampling scheme drawing states σ⃗ from the VMC distribution |Ψ_M(σ⃗)|^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. <ref>, 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, and that chemical accuracy is reached only for a relatively large number of samples. 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 |Ψ_M(σ⃗)|^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. 



  
*Outlook.- 
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.  

Several future directions can be envisioned. The distinctive peaked structure of the molecular wave function calls for 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 <cit.>. 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. 

 
The Flatiron Institute is supported by the Simons Foundation. A.M. acknowledges support from the IBM Research Frontiers Institute. Neural-network quantum states simulations are based on the open-source software NetKet <cit.>. Coupled cluster and configuration interaction calculations are performed using the PySCF package <cit.>. The mappings from fermions to spins are done using Qiskit Aqua <cit.>. The authors acknowledge discussions with G. Booth, T. Berkelbach, M. Holtzmann, J. E. T. Smith, S. Sorella, J. Stokes, and S. Zhang.  






§ GEOMETRIES FOR DIATOMIC MOLECULES
 

The equilibrium geometries for the molecules presented in this work were obtained from the CCCBDB database  <cit.>. For convenience, we present them in Table <ref>.







§ 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_q | σ⃗⟩, appearing in the local energy, Eq. <ref>. Since H_q is a sum of products of Pauli operators, the goal is to efficiently compute matrix elements of the form 

    ℳ(σ⃗,σ⃗') = ⟨σ⃗' | σ^ν_1_1σ^ν_2_2…σ^ν_N_N | σ⃗⟩,

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 σ⃗^' such that

    σ_i^'=σ_i   ν_i∈(I,Z)
    σ_i^'=-σ_i   ν_i∈(X,Y)

and the matrix element is readily computed as 

    ℳ(σ⃗,σ⃗') = (i^n_y)Π_k:v_k∈(y,z)σ_k,

where n_y is the total number of σ^y operators in the string of Pauli matrices.

