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.

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[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-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-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 on N fermionic modes, and one- and two-body integrals t and u. 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 formwhere 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 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 formwhere . For all the spin encodings considered, the final outcome is a spin Hamiltonian with the general formdefined as a linear combination with real coefficients h of , N-fold tensor products of single-qubit Pauli operators I, σ, σ, σ. Additionally, under such mappings, there is a one to one correspondence between spin configuration and the original particle occupations . 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-of-equilibrium[26-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 formHere, are complex-valued network parameters , 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-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 valueThis expectation value can be evaluated using Monte Carlo sampling using the fact that the energy (and, analogously, any other observable) can be written aswhere we have defined the local energy

Given samples drawn from the distribution , the average over the samples 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. In generic QC problems, this cost scales in the worst case with , as compared to the linear scaling in typical condensed matter systems with local interaction.

Sampling from is performed using Markov chain Monte Carlo (MCMC), with a Markov chain constructed using the Metropolis–Hastings algorithm[38]. Specifically, at each iteration, a configuration is proposed and accepted with probabilityThe sample 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 rather than the corresponding spin basis . The proposed configuration 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[39,40], closely related to the natural-gradient method used in machine learning applications[10,41].

Computational complexity

The main computational cost of the procedure arises from the evaluation of the local energy (Eq. (5)) of the samples generated. This gives an overall computational complexity of where Nvar = MN + M + N is the number of parameters in the network, Nop is the number of Pauli strings in the spin Hamiltonian defined by Eq. (4) and Nsamp is sample size.

However, as can be seen in Fig. 2, there are only small number of relevant configurations in the wavefunction, thus each sample only contains a few unique configurations. By caching amplitudes the computational cost can be significantly reduced to where Nunique ≪ Nsamp is the average number of unique configurations in each sample. Typically, for a sample size of 10,000 there are only about few hundred unique samples.

Fig. 2

Electronic correlations.

Probabilities (in logarithmic scale) of the 500 most probable configurations in the FCI (blue), RBM (orange), and CISD (green) wavefunctions for the equilibrium nitrogen N2 molecule in the STO-3G basis.

Potential energy surfaces

We first consider small molecules in a minimal basis set (STO-3G). We show in Fig. 1 the dissociation curves for C2 and N2, 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 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.

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 C2 and b N2, 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.

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 = ∣∣ in Eq. (4), i.e., the number of single-qubit Pauli operators in that are different from the identity I, depends on the specific mapping used. Jordan–Wigner and parity mappings have linear scalings w = O(N), while the Bravyi–Kitaev encoding has a more favorable scaling , 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 , 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 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. 3.

Fig. 3

Comparison of different spin mappings.

Accuracy of the RBM (green star) representations for three different mapping types (Jordan–Wigner, Parity, and Bravyi–Kitaev) and three different molecules (LiH, C2, and N2) in their equilibrium configuration in the STO-3G basis. The geometries used are reported in the Methods section.

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

Fig. 4

Sampling size dependence of the converged energies.

Converged energy of H2O in the 6-31g basis (26 spin-orbitals) as the number of samples used for each VMC iteration is varied. The converged energy for the samples obtained using the Metropolis algorithm (blue circles) matches that obtained using exact sampling (green crosses), beating the accuracy of CISD and approaching chemical accuracy (red line) for the largest sample size. In the inset, we also show the variational energy as the number of hidden units is increased from 2 to 26.

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).

Table 1

Equilibrium energies (in Hartree) as obtained by different methods.

Molecule	RBM	Jastrow	CISD	CCSD	CCSD(T)	FCI
H2	−1.1373	−1.1373	−1.1373	−1.1373	−1.1373	−1.1373
LiH	−7.8826	−7.8814	−7.8827	−7.8828	−7.8828	−7.8828
NH3	−55.5277	−55.4770	−55.5258	−55.5280	−55.5281	−55.5282
H2O	−75.0232	−74.9784	−75.0221	−75.0231	−75.0232	−75.0233
C2	−74.6892	−74.5001	−74.6371	−74.6745	−74.6876	−74.6908
N2	−107.6767	−107.5924	−107.6591	−107.6717	−107.6738	−107.6774

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 C2 and N2 where we use α = 2.

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.

Table 2

Equilibrium configurations used for the ground-state calculations presented in the main text. The coordinates (x, y, z) are given in angstroms (Å).

Molecule	Basis	Geometry
H2	STO-3G	H(0, 0, 0)
		H(0, 0, 0.734)
LiH	STO-3G	Li(0, 0, 0)
		H(0, 0, 1.548)
		N(0, 0, 0.149)
NH3	STO-3G	H(0, 0.947, −0.348)
		H(0.821, −0.474, −0.348)
		H(−0.821, −0.474, −0.348)
C2	STO-3G	C(0, 0, 0)
		C(0, 0, 1.26)
N2	STO-3G	N(0, 0, 0)
	.	N(0, 0, 1.19)
		H(0, 0.769, −0.546)
H2O	STO-3G	H(0, −0.769, −0.546)
		O(0, 0, 0.137)
		H(0, 0.795, −0.454)
H2O	6-31G	H(0, −0.795, −0.454)
		O(0, 0, 0.113)

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 , appearing in the local energy, Eq. (9). Since H is a sum of products of Pauli operators, the goal is to efficiently compute matrix elements of the formwhere 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 thatand the matrix element is readily computed aswhere n is the total number of σ operators in the string of Pauli matrices.

Simulation details

The optimization follows the stochastic reconfiguration scheme as detailed in the supplementary material of ref. [10]. Given a variational ansatz depending on parameters {α}, the parameter update δα is given by solution of the linear equationwhere 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.

Peer Review File