Solving the Schrödinger equation for systems of many interacting bosons or fermions is classified as an NP-hard problem due to the complexity of the required many-dimensional wave function, resulting in an exponential growth of degrees of freedom. Reducing the dimensionalities of quantum mechanical many-body systems is an important aspect of modern physics, ranging from the development of efficient algorithms for studying many-body systems to exploiting the increase in computing power. To write software that can fully utilize the available resources has long been known to be an important aspect of these endeavors. Despite tremendous progress having been made in this direction, traditional many-particle methods, either quantum mechanical or classical ones, face huge dimensionality problems when applied to studies of systems with many interacting particles.

Over the last 2 decades, quantum computing and machine learning have emerged as some of the most promising approaches for studying complex physical systems where several length and energy scales are involved. Machine learning techniques and in particular neural-network quantum states Goodfellow et al. [1] have recently been applied to studies of many-body systems, see, for example, Refs. Carleo and Troyer [2]; Carrasquilla and Torlai [3]; Pfau et al. [4]; Calcavecchia et al. [5]; Carleo et al. [6]; Boehnlein et al. [7]; Adams et al. [8]; Lovato et al. [9], in various fields of physics and quantum chemistry, with very promising results. In many of these studies, results have aligned well with exact analytical solutions or are in close agreement with state-of-the-art quantum Monte Carlo calculations.

The variational and diffusion Monte Carlo algorithms are among the most popular and successful methods available for ground-state studies of quantum mechanical systems. They both rely on a suitable ansatz for the ground-state of the system, often dubbed the trial wave function, which is defined in terms of a set of variational parameters whose optimal values are found by minimizing the total energy of the system. Devising flexible and accurate functional forms for the trial wave functions requires prior knowledge and physical intuition about the system under investigation. However, for many systems we do not have this intuition, and as a result it is often difficult to define a good ansatz for the state function.

According to the universal approximation theorem, a deep feedforward neural network can represent any continuous function within a certain error Hornik et al. [10] — see also Refs. Murphy [11]; Hastie et al. [12]; Bishop [13]; Goodfellow et al. [1] for further discussions of deep learning methods. Since the variational state in principle can take any functional form, it is natural to replace the trial wave function with a neural network and treat it as a machine learning problem. This approach has been successfully implemented in recent works, see, for example, Refs. Pfau et al. [4]; Carleo and Troyer [2]; Cassella et al. [14]; Adams et al. [8]; Lovato et al. [9], and forms the motivation for the present study. Here, the neural network of choice was derived from so-called Gaussian-binary restricted Boltzmann machines, much inspired by the recent contributions by Carleo et al., see, for example, Refs. Carleo and Troyer [2]; Carleo et al. [6]. Unlike binary-binary restricted Boltzmann machines Sieber and Gehringer [15], the approximation properties of Gaussian-binary restricted Boltzmann aren’t yet well-understood. However, they are considered highly-expressive networks for certain problems, such as dimensionality reduction and continuous probability density estimation. Note that neural-network representations of variational states are more general, as they don’t in principle require prior knowledge of the ground-state wave function, thereby opening the door to systems that have yet to be solved. Particular attention however has to be devoted to the symmetries of the problem, whose inclusion is critical to achieve accurate results.

In this work, we will focus on systems of electrons confined to move in two-dimensional harmonic oscillator systems, so-called quantum dots. These are strongly confined electrons and offer a wide variety of complex and subtle phenomena which pose severe challenges to existing many-body methods. Due to their small size, quantum dots are characterized by discrete quantum levels. For instance, the ground states of circular dots show similar shell structures and magic numbers as seen for atoms and nuclei. These structures are particularly evident in measurements of the change in electrochemical potential due to the addition of one extra electron. Here, these systems will serve as our test of the applicability of restricted Boltzmann machines as artificial neural network variational states.

The theoretical foundation and the methodology are explained in Section 2. The subsequent sections present our results with an analysis of computational methods and resources. In the last section we present our conclusions and perspectives for future work.

For any Hamiltonian H ̂ and trial wave function ψ _T , the variational principle guarantees that the expectation value of the energy E _T is greater than or equal to the true ground state energy E ₀, E 0 ≤ E T = 〈 ψ T | H ̂ | ψ T 〉 〈 ψ T | ψ T 〉 . (1) Thus, approximate solutions to the time-independent Schrödinger equation can be obtained by choosing a careful parameterization of the wave function and minimizing the energy E _T with respect to the parameters. Since the integrals representing E _T are normally high dimensional, it is most efficient to evaluate them by means of Monte Carlo methods E T ≈ ⟨ E L ⟩ = 1 n ∑ i = 1 n E L R i , R i ∼ | ψ T R | 2 . (2) This involves collecting n samples of configurations and averaging over the so-called local energies E L R = 1 ψ T R H ̂ ψ T R . (3)

We apply the variational Monte Carlo (VMC) method to various circular quantum dots systems. These are systems of interacting electrons confined to move in a two-dimensional harmonic oscillator well. The (scaled) ¹ Hamiltonian is given by H ̂ = 1 2 ∑ i − ∇ i 2 + ω 2 r i 2 + ∑ j ≠ i 1 r i j , (4) where ω is the oscillator frequency, r _i is the distance between electron i and the origin, and r _ij is the distance between electrons i and j. We will henceforth assume the total number of electrons N to be even and the total spin of the system to be zero.

A simple ansatz can be built starting from the analytical solutions to the non-interacting case. The harmonic oscillator eigenfunctions are given by ϕ m , n x , y ∝ e − ω x 2 + y 2 H m ω x H n ω y , (5) where H _n are the Hermite polynomials of degree n. To constrain the antisymmetry of the many-body wave function, products of the lowest N/2 spatial states and the two spin states ξ _±(σ) are used as a basis for a Slater determinant ψ SD R = det ϕ m , n x i , y i ξ k σ i , where m, n, k label the single-particle state, i labels the particle, and R contains all coordinates of the N particles. As an aside, we do not include the spin projections σ _i as explicit inputs to the wave function as we will describe how to treat them separately in Section 2.2. We then define a reference state by pulling the common exponential term out of the determinant and inserting a single variational parameter α ψ Ref R ; α = e − α ω ∑ i x i 2 + y i 2 det H m ω x i H n ω y i ξ k σ i . (6) Correlations among electrons can be handled by a Padé-Jastrow factor Drummond et al. [16], g R ; β = exp ∑ i = 1 N ∑ j > i N a i j r i j 1 + β r i j , (7) where β is a variational parameter and a i j = 1 / 3 if  σ i = σ j 1 if  σ i ≠ σ j , in order for the Kato cusp condition to be satisfied Huang et al. [17]. The product of the Slater determinant and the Padé-Jastrow factor is commonly named the Slater-Jastrow ansatz, ψ Slater - Jastrow R ; α , β = ψ Ref R ; α × g R ; β . (8)

In this section we compare the computational complexity, ground state energy, energy convergence, contribution from the different Hamiltonian terms and electron densities for various trial wave function ansatzes. The RBM ansatz consists of a Slater determinant with Hermite polynomials as the basis, multiplied with the RBM marginal distribution, as presented in Eq. 10 We have not made any attempt to include optimized single-particle state functions through mean-field optimizations like Hartree-Fock theory. The second ansatz is the RBM ansatz with a correlation factor described in Eq. 7, abbreviated RBM + PJ. We have also include results obtained using a traditional Slater-Jastrow ansatz (Eq. 8), and as reference we use a plain Slater determinant (without a correlation factor, Eq. 6). The diffusion Monte Carlo (DMC) results obtained by Høgberget [24] are our main reference for the ground state energy. For quantum dots with two electrons we compare our results with corresponding analytical ones from Ref. Taut [25].

Figure 3 displays the average CPU time per iteration as a function of system size for the different wave function ansatzes. To obtain these estimates, we averaged the CPU time per iteration over 10³ iterations. The RBM + PJ and Slater-Jastrow ansatzes are the most computationally demanding due to the Padé-Jastrow factor, which depends on the relative distance between electrons and requires updating of a matrix containing electron-electron distances. Each process has its own memory, and the matrix is not shared across processes.

The impact of parallel processing on computational overhead is minimal compared to other aspects of the code, except for the burn-in period. In fact, the variation in CPU time as a result of noise is much more important than the variation due to parallel processing. Additional information on CPU time per iteration for different wave function ansatzes is presented in Table 1. The reported calculations were performed for 10⁴ Monte Carlo cycles and an oscillator frequency of ω = 1.0, with similar results for other values of ω. Each benchmark simulation was run on a single core without a burn-in period. Production runs require more cycles and consequently longer CPU times per iteration. For instance, for a system size of N = 90 with 2²⁰ Monte Carlo cycles run on 30 nodes (960 threads), the RBM, RBM + PJ, and Slater-Jastrow ansatzes required 0.55960, 2.7345, and 1.3001 s per iteration, respectively, with approximately 50,000 iterations needed for convergence.

The ground state energies of two-dimensional quantum dots of various sizes and frequencies are presented in Table 2. The RBM + PJ and Slater-Jastrow ansatzes provide the lowest ground state energies as expected as the correlation factor does explicit satisfy the Kato’s cusp conditions. The relative errors with respect to diffusion Monte Carlo (DMC) calculations tend to be less than 0.2% for both methods. The various results obtained with the RBM + PJ ansatz show that this ansatz dominates for small quantum dots, but is outperformed by the Slater-Jastrow ansatz for larger systems. We suspect this is due to the fact that the former ansatz is more complex and contains significantly more parameters than the latter, and has therefore a hard time finding the global minimum. The optimization could be improved with a more sophisticated optimization algorithm.

For low frequency dots (ω < 0.28), the RBM produces ground state energies lower than the reference energy, but fails for large high frequency dots (N > 6, ω > 0.28). For low frequency dots the interaction energy dominates, and the RBM manages to learn the correlations. For high frequency dots, it is the one-body part of the Hamiltonian which dominates and the interaction part is less important for the final expectation value of the energy. For these systems, the Slater determinant part of the wave function ansatz is the most important one, meaning in turn that the RBM ansatz we have used may not capture well the structure of the single-particle states that define the Slater determinants.

We also note that the reference values are similar to corresponding Hartree-Fock energies up to 30 particles Mariadason [26], which is expected as both approaches try to find the optimal single Slater determinant without a correlation factor.

Some simulations were also performed for the RBM and RBM + PJ ansatzes with sorted network inputs to enforce anti-symmetry under exchange of two particles, as suggested by Ref. Saito [27], among others. Sorted inputs showed promising results, where the ground state energy typically dropped in the fourth or fifth digit. For example, the RBM + PJ ansatz found the ground state energy of a system with N = 30 electrons and ω = 0.5 to be 187.311(1) Hartree. This is lower than the Slater-Jastrow energy. We will investigate this further in a follow-up paper [28].

For the traditional ansatzes that do not contain artificial neural networks, the learning rate determines to a large extent how fast the trial wave function converges towards the true ground state wave function. However, for the RBM and the RBM + PJ approaches, the results are sensitive to the chosen learning rate values. A too large learning rate can easily cause exploding energies, and with a too small learning rate the obtained energies with a given trial function might not converge at all. Our strategy is to find the highest learning rate that does not lead to exploding energies. In general this requires efficient grid searching methods. Here, when the energies have converged, we decrease the learning rate by a factor of 10 ³ and let it run until it converges again. Also, the RBMs tend to converge step-wise, making it hard to know whether or not they have converged.

In Figure 4, the convergence of the local energy is plotted for our four ansatzes for N = 2 electrons. Here the results are compared with analytical calculations for N = 2. A similar behavior is shown for N = 56 electrons and ω = 0.1 in Figure 5. For the N = 2 electrons case, RBM + PJ turns out to be the most accurate ansatz with small absolute errors. During training, the Padé-Jastrow parameter is rather constant, but the adjustment of weights seems important for the model to reach the energy minimum. For larger systems, the Slater-Jastrow ansatz provides a slightly lower energy than the RBM + PJ ansatz. Since the number of inputs for the RBM increases with number of particles, the network contains additional training parameters as we increase the number of electrons. Therefore, the networks get more complex for large quantum dots, and they are naturally harder to train.

The spatial one-body density plots for ω = 1.0 and the RBM and RBM + PJ ansatzes are presented in Figures 6, 7. The electron densities have a wave shape for both ansatzes, with two nodes for N = 2, three nodes for N = 6 and so on similar to those observed in Refs. Ghosal et al. [29]; Høgberget [24]. The RBM seems to exaggerate the states with more distinct peaks (higher peaks and lower wave valleys) compared to the RBM + PJ results. This can be explained by more localized electrons, which becomes even more apparent in low frequency dots (see Figure 8), as the interactions are modelled differently. The RBM would hardly be able to model the correct electron-electron distances, as the network itself is purely linear. The electron density was found to be shape-invariant for high frequency dots (ω > 0.28), with decreasing spatial range as the frequency increases.

When reducing the frequency further, the interaction energy dominates over the kinetic energy and harmonic oscillator energy (see Figures 10, 11), and the electrons naturally become more localized. The Padé-Jastrow factor solves this by radially localizing the electrons and conserving the circular symmetry, such that the electrons are confined to specific orbitals. This can be seen from Figure 8, where the spatial one-body density is plotted for low frequency dots (ω = 0.1) and system sizes N = 2, 6 and 12. Radial localization is also what we would expect from the Hamiltonian, which is strictly circular symmetric. On the other hand, the RBM seems to localize the electrons both in radial and angular direction, with the number of electrons corresponding to the number of peaks. This is a nonphysical solution to the problem, and shows that the RBM ansatz breaks down for low frequencies. The RBM + PJ ansatz, on the other hand, confines the electrons in orbitals. Distinct peaks in radial direction is what is expected from Wigner crystallization, which we might see indications of with density parameters r _s ≈ 6.7, r _s ≈ 1.2, and r _s = 0.3 respectively for the three system sizes with the RBM + PJ ansatz. Notice for instance the small “pit” on top of the N = 2 plot for RBM + PJ, which isn’t seen for higher oscillator frequencies. In Figure 9, we have decreased the one-body density even further to ω = 0.01 for N = 2. As expected, the electrons become even more localized and are clearly showing Wigner crystallization effects with a density parameter of r _s ≈ 29. For the RBM ansatz, the electrons are strongly localized (in all directions) and the electron densities barely overlap. For the RBM + PJ ansatz, we observe strong orbital confinement where the Wigner crystallization wasn’t the target of this study, but the framework seems capable of a profound study of this phenomenon.

To understand the behavior of the one-body density for various frequencies, we investigate the expectation values of the kinetic energy ( ⟨ T ̂ ⟩ ) , the harmonic oscillator potential energy ( ⟨ V ̂ ext ⟩ ) and the two-body interaction energy ( ⟨ V ̂ ext ⟩ ) . The energy distribution is plotted for the RBM + PJ ansatzes for N = 2 and N = 20 (Figures 10, 11) for frequencies ω ∈ {0.01, 0.1, 0.28, 0.5, 1.0, 2.0, 3.0, 5.0, and 10.0}. For large values of the oscillator frequency it is the one-body part of the Hamiltonian which dominates in absolute value (kinetic energy and harmonic oscillator potential energy) compared with the expectation value of the two-body interaction. For such frequencies we notice also that the results can almost be interpreted in terms of the virial theorem. This theorem provides a useful relation between the kinetic and potential energy Fock [30]. For circular quantum dots without the two-electron interaction, the theorem reads 2 ⟨ T ̂ ⟩ = 2 ⟨ V ̂ ext ⟩ . From our results we notice that, due to the interaction energy, the ratio between kinetic and harmonic oscillator energies are not exactly equal to two for large frequencies. However, for such frequencies the interaction energy plays a less prominent role, resulting thus in a ratio which is close to the ideal value. When we decrease the frequency however, and the system becomes more dilute with an increase of the mean distance between particles, one notices a somewhat counter intuitive behavior. The expectation value of the kinetic energy and the harmonic oscillator potential energy decrease (and the electrons become more localized as seen in the one-body densities above). However, many-body correlations increase in importance with decreasing frequency. This is reflected in the increased role of the expectation value of the two-body Coulomb interaction, an effect which is simply due to the infinite range of the Coulomb interaction. If we were to multiply the Coulomb interaction with a finite range factor, this effect would disappear. Figures 10, 11 show this behavior rather clearly.

The quantum dot systems were studied using a general framework for variational Monte Carlo simulations (The code is available on github.com/evenmn/VMaChine). Importance sampling was used to accelerate the simulations Metropolis et al. [31]. To minimize the local energy, we applied the Adam optimizer with β ₁ = 0.9 and β ₂ = 0.999 as suggested by Ref. Kingma and Ba [32]. We use an adaptive number of cycles, starting from 2²⁰ and then increased to 2²⁴ for the ten last iterations after the energy had converged, and further to 2³⁰ for the very last iteration to reduce the statistical uncertainty and noise in electron density. The particle step length was chosen to get an acceptance ratio close to 99.5%, spawning from 10¹ for the smallest and weakest systems to 10^–3 for the large and narrow oscillators. The optimal learning rate was found by grid searches, and varied from 10¹ to 10^–5. Both the step length and the learning rate depend strongly on the system and the trial wave function.

For the RBM and RBM + PJ ansatzes, only the raw particle positions were input to the RBM. This choice was made for performance reasons, as inputting processed positions like the electron-electron distances would lead to significantly larger computational efforts, with an increased network complexity.

The Gaussian parameter was initialized to α = 1.0, which is the analytical optimum without interaction. We use Xavier initialization for the RBM weights, putting them all close to zero Glorot and Bengio [33]. The special case with all weights set to zero corresponds to our reference ansatz with α = 1.0. For all the RBMs but the most narrow ones (ω = 1.0), the number of hidden nodes was set to 6, giving 40 to 1272 parameters for the various system sizes. For ω = 1.0, we used H = 12 hidden nodes to achieve a lower energy.

All the simulations were run on Intel Xeon E5-2670 CPUs. In total the computational cost of this project was of the order of 10⁶ CPU hours with largest amount of cycles spent, for obvious reasons, on the largest systems.

We found that the RBM ansatz gives a significantly lower ground state energy than the reference ansatz for low frequency quantum dots. This may indicate that the RBM manages to capture some of the electron-electron correlations. Based on the one-body density plots, the RBM found the electrons to be localized both angularly and radially compared to the ansatzes containing a Padé-Jastrow factor, which confine electrons radially only. For high frequency dots, the RBM fails in the sense that the obtained ground state energy is larger than the reference energy. This can be explained by the fact that when the interactions get less important, the reference ansatz is a good guess. The RBM + PJ ansatz gives energies close to the DMC energies for small quantum dots. The ansatz performance needs to be seen in light of the computational cost, as the ansatzes containing a Padé-Jastrow factor are far more computationally intensive.

From the one-body density and energy distribution plots, it is apparent that the RBM ansatz isn’t able to capture the correlations at the same level as the Padé-Jastrow factor. Because of the linearity of the network, it is impossible for it to compute the distance between the particles, which is crucial to model the interactions correctly. One solution to this could be to input the electron-electron distance into the network, as discussed by Pfau et al. [4], but this would increase the computational intensity. Also, despite including a Slater determinant, the trial wave function is not necessarily anti-symmetric when including an RBM. We performed some simulations where anti-symmetry was forced by sorting the network inputs, which showed promising results in terms of the ground state energy. However, although promising results can be obtained, RBMs are less flexible than general neural networks that make fewer assumptions about the specific mathematical forms of the trial functions. We have encoded explicitly the anti-symmetry via a Slater determinant. Furthermore, two-body correlations are constructed using a Jastrow factor. An RBM with Gaussian distributions is capable of capturing the one-body part of the problem, but is less flexible in finding two-body or more complicated many-body correlations. Although the RBM results reported here are promising compared with existing VMC calculations, recent results with neural networks like those presented in, for example, Refs. Pfau et al. [4]; Cassella et al. [14]; Lovato et al. [9]; Adams et al. [8]; Carrasquilla and Torlai [3] offer much more flexible and promising research venues for deep learning methods applied to many-body problems. Results obtained with deep neural networks for these systems will be presented in a future work [28].