{'fraction_non_alphanumeric': 0.029410714923082414, 'fraction_numerical': 0.006317592890031053, 'mean_word_length': 4.542631058358062, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
QUBO-based density matrix electronic structure method

Density matrix electronic structure theory is used in many quantum chemistry methods to"alleviate"the computational cost that arises from directly using wave functions. Although density matrix based methods are computationally more efficient than wave functions based methods, yet significant computational effort is involved. Since the Schr\"odinger equation needs to be solved as an eigenvalue problem, the time-to-solution scales cubically with the system size, and is solved as many times in order to reach charge or field self-consistency. We hereby propose and study a method to compute the density matrix by using a quadratic unconstrained binary optimization (QUBO) solver. This method could be useful to solve the problem with quantum computers, and more specifically, quantum annealers. The method hereby proposed is based on a direct construction of the density matrix using a QUBO eigensolver. We explore the main parameters of the algorithm focusing on precision and efficiency. We show that, while direct construction of the density matrix using a QUBO formulation is possible, the efficiency and precision have room for improvement. Moreover, calculations performing Quantum Annealing with the D-Wave's new Advantage quantum processing units is compared with classical Simulated annealing, further highlighting some problems of the proposed method. We also show some alternative methods that could lead to a better performance of the density matrix construction.


I. INTRODUCTION
The correct choice of the model or method used to represent the electronic structure of a chemical system is crucial to get a good description of any observable (forces, energies, etc.) of such a system. The quality of a model can be estimated as the ratio between its predictive power and its complexity (Quality = Predictive Power /Complexity), the predictive power being related to the amount of information that can be extracted from a calculation, and the complexity mainly referring to how the time-to-solution scales with the system size 1 . Although this definition is vague, it gives us the general understanding that, in order to optimize a method (increase its quality) it is necessary to work on approaches that could decrease the complexity, and, at the same time, avoid losing predictive power; In many cases the predictive power translates into accuracy. Reducing the complexity is the typical approach that has been taken given that computational resources have always been limited due to the ever increasing need of addressing larger systems. However, with the arrival of quantum computers (an alternative computational paradigm) and more specifically, quantum annealers (QAs) 2 , it is possible that more complex formulations (in terms of traditional computation) of a given problem end up having shorter or even instantaneous time-to-solution without accuracy resignation. This has been the case of graph problems including graph partitioning and community detection solved with QAs that was made possible due to a quadratic unconstrained binary optimization (QUBO) reformulation of the problem 3,4 . For the case of quantum chemistry methods, the long term vision to reach is the one in which all the computational burden of the eigenvalue problem and self-consistent process could be shifted towards an "extremely complex," yet easily solvable QUBO problem (See Figure 1). In this paper we take an initial step towards the aforementioned vision by On the left we represent the traditional self-consistent charge (SCC) approach with steps that include the correction of the Hamiltonian elements H; the construction of the density matrix ρ from a Fermi function fF of the Hamiltonian; and the calculation of the effective potential V as a Coulombic function fC of ρ to be used to reconstruct the Hamiltonian in the subsequent SCC step. On the right, the alternative QUBO idea is proposed.
exploring how feasible it is to construct the density matrix (DM) using linear algebra QUBO algorithms. Previous work on QAs used in quantum chemistry has been recently applied to determine the ground and firsts excited states on the Full Configuration Interaction (CI) method by diagonalizing the Full CI (FCI) matrix 5 ; the downfolding of the FCI matrix to reduce the size of the effective matrix 6 ; and Time Dependent Density Functional Theory (DFT) to compute excitations 7 . The density matrix (DM) formalism is ubiquitous to several methods in computational chemistry that compute the electronic structure of a molecular system. Among these methods, the Khon-Sham method within DFT is arguably the most frequently used of the quantum chemistry methods, preferred over post Hartree-Fock methods, for its performance in terms of predictive power relative to its computational cost 8, 9 . The foundations of this theory will be briefly stated below.
The many-body wave function of interacting electrons in a chemical system resulting in the solution of the Schrödinger equation depends on 3×N e spatial coordinates, where N e is the number of electrons in the system. DFT can be used to overcome this complexity, allowing for the calculation of observables just by knowing the electronic density of the system. This theory is based on using the electronic density n(r) (r being the electron coordinates) as a variable of the total energy functional E[n(r)] of the system, and the variational principle to minimize the latter. The theorems of P. Hohenberg and W. Kohn 10 ensure that the energy, wave function, and all other electronic properties are uniquely determined by the electron density n(r). Moreover, for any external potential v ext (r) generated by the nuclear configuration of a molecular system, there is a bijective relation with a unique density function n(r). To solve for n(r) the Kohn-Sham method relies on auxiliary wave functions ψ i (r) that are used to transform the minimization of the energy functional into the following eigenvalue problem: where : Here, v ext (r) is the external potential given by the nuclei; It is the part of the energy functional that characterizes the chemical system in question. E xc is the functional that takes into account the energy of exchangecorrelation. It is written as the sum of two components: The kinetic energy correlation component, which comes from correcting the error by considering the kinetic energy of non-interacting electrons, and the inter-electronic interaction containing the non-kinetic correlation and exchange. The exact form of this functional is not known, but many approximations exists nowadays with varying accuracies and computational costs 11 . In this work we will use the most basic one where the functional is represented by a local density approximation (LDA) of the form: E xc [n(r)] =´ xc (n)n(r)d 3 r, where xc (n) is the correlation and exchange energy per electron for a uniform gas of density n. Equation 2 is solved iterativelly using equation 3 until reaching self-consistency between the density n(r) and the effective potential v ef f (r). When written using a basis set of functions {φ i } i∈1,..,N with S ij =´r φ * i (r)φ j (r)dr and H ij =´r φ * i (r)Ĥφ j (r)dr equation 2 turns into the following generalized eigenvalue problem: where, c k is the k-th component of the expansion of the wavefunction ψ i (r) in the chosen basis set. In matrix form equation 4 becomes: which is a generalized matrix eigenvalue problem, where is a diagonal matrix with entries ij = δ ij i containing all the eigenvalues of matrix H. Finally, C is a matrix containing all the expansion coefficients for each of the N eigenfunctions ψ i (r).
Given the expansion coefficients C, the DM which is our object of interest here, is constructed as: where we have used the Fermi distribution function f n for a system in thermal equilibrium defined as: where k b is the Boltzman constant, F is the Fermi level (or electrochemical potential) and T is the system's temperature. For practical purposes, when temperature is low, f n ( ) = 2Θ( − F ), leading to: The DM can be hence directly constructed by summing the self outer products of each eigenvector (in order of increasing energy i ) up to the number of occupied states (occ = N e /2). Given the DM, many properties can be computed since the expectation value of any operator A can be calculated as A = tr(ρA). For instance, the trace of the DM gives the total number of electrons of the system, and the trace of ρH provides the total electronic energy also known as band energy of the system 12 . Another important property is that we can easily compute n(r) as follows: which can be reused in equation 3 to recomputeĤ. Other methods such as Hartree-Fock, Tight-Binding, and Semiempirical are also based on constructing the singleparticle DM 1, [13][14][15] . From what was explained above, matrices and C arise from diagonalizing the matrix representation ofĤ in the aforementioned basis set. This diagonalization step in all the density matrix based methods is typically the bottleneck of the whole calculation, and scales as O(N 3 ), where N is the number of elements in the basis set (which scales linearly with the system size). Over the years, computational chemists have dedicated lots of effort to develop algorithms that could solve the density matrix with O(N ) scaling [16][17][18] ; leading to the so called "order n methods." Some of them have also adapted algorithms to the ever evolving classical computer architectures 19,20 . The problem with all these order n approaches is that they require a trade-off between accuracy and efficiency. The faster we want the method to be, the lower will be the accuracy of our results. In this paper we analyze the tractibility of quantum annealing as an alternative computational paradigm to get accurate results at a lower computational cost for DM based quantum chemistry methods.

II. METHODS
In this section we will describe the methods used, and in particular we will focus on our recently developed high precision QUBO eigensolver 21 .

A. Quantum annealers
Quantum Annealers (QAs) are specific adiabatic quantum computing devices particularly suitable to solve combinatorial optimization problems 22 . In these devices the exploration of the combinatorial space is given by the delocalization and degree of superposition of the initial quantum state that decays to the state of minimum energy during a process called quantum annealing. The D-Wave QA is an architecture that uses this quantum annealing technique to solve for the ground state of an Ising Model Hamiltonian that is derived from a combinatorial optimization problem. In other words, any problem to be solved with the D-Wave QA will need to be reformulated to be solved as an Ising Model Hamiltonian whose energy is evaluated as follows: for which s i ∈ {−1, 1} are magnetic spin variables; h i are local programmable fields; and J ij are programmable field coupling interactions. An Ising model Hamiltonian can be easily transformed into its boolean equivalent using the following transformation: s = 2q − 1, where 1 = (1, . . . , 1). The boolean equivalent problem that is obtained is referred to as a QUBO problem and it is expressed as follows: with q i ∈ {0, 1}.

B. High-precision QUBO-based eigensolver
The solution of the extremal eigenvectors with high accuracy (eigenvector corresponding to minimal or maximal eigenvalue) was recently possible thanks to a new QUBO formulation of the problem. 21 . This method builds on top of the quantum annealer eigensolver (QAE) work of Teplukhin et. al. 5 . The main difference between this high-precision QUBO eigensolver method and the previous QAE, is that it solves the problem as a fixed point method instead of using binary search 23 . The new method also involves a steepest descent search together with a systematic adjustment of the length of the precision vector to guarantee high accuracy convergence. The precision vector defined as follows is used to convert numbers from binary to real representation: where m is the number of bits used in the representation. In order to transform a number from binary q to real x representation we need to perform x = v T q, where q is the vectorized form of q. Note that with this representation we cover the [−1, 1) segment where the binary representation q of x will satisfy: −1 ≡ (1, 0, 0, 0, ..., 0) ≤ q < (0, 1, 1, 1, ..., 1) ≡ 1. A real field multiplication of x × y has the binary equivalent q T x (v ⊗ v T )q y . A particular example of this multiplication is given in the Appendix.
With the help of this precision vector and the previously introduced multiplication, any quadratic equation x T Bx for which we need to find an optimal x with B ∈ R N ×N , will have binary equivalent binary field vector space. Any linear term r T x with r ∈ R N for which we need to find an optimal x, will have a binary equivalent q T Diag(rI ⊗ v))q, where I is the N × N identity matrix, and Diag() is an operator that, given a in B mN builds an operator Given the Hamiltonian (H) and overlap (S) matrices, every eigenpair (λ, x) will satisfy: (Hx − λSx) = 0, and, in particular, the lowest eigenpair (λ min , x min ) will be an absolute minimum of: Since H − λS is positive semidefinite, we can instead minimize: Since λ min is not known, equation 12 needs to be solved iterativelly as a fixed point. A good initial value for λ happens to be the average eigenvalue of H computed as tr (H) /N 21 . The method converges once the changes in λ min are below a certain preset tolerance. In order to achieve a higher precision in the calculation of both x min and λ min , we have introduced a steepest descent phase also solved as a QUBO problem. Setting the previously converged x min as an initial guess x 0 min , we Taylor expand f (x) = x T (H − λS)x around x 0 min to obtain: where ∇f and H(f ) are the gradient and Hessian of f respectivelly. For the case of f (x), the gradient and Hessian becomes ∇f = 2(H−λS)(x 0 min ) T , and H(f ) = 2(H−λS) respectivelly. Hence, we can compute a high-precision approximation to the lowest eigenpair as: In this case, δ converges to a fixed point, and, once δ does not change anymore, the length of the precision vector is shortened to 10% of its previous value. This procedure guarantees convergence to any arbitrary desired precision. For every "next" eigenpair we want to compute, the previous eigenvalue x min ≡ x i−1 of H needs to be "pushed out" of the eigenspectrum by using the following transformation: where w is an estimation for the "eigenspectrum width." A good estimation of w can be calculated as w = h − l where h and l are the highest and lowest Gershgorin bounds of H respectively. Using the previously explained procedure to compute eigenvectors, the density matrix is constructed as explained before; using a straightforward summation of eigenvector outer products up to the number of occupied states. This is: A python like pseudocode for the full algorithm is detailed in Pseudocode 1. The inputs are the Hamiltonian H, Overlap S, the occupation nocc, and the target precision prec T . The output is the DM ρ. x =solve(x T (H − λS)x) #Send to Annealer Pseudocode 1 implements a straightforward way of computing the density matrix using the high-precision QUBO eigensolver. From this pseudocode we can immediately identify some important properties. We see that the matrix operations in the inner loops are O(N 2 ) unless they can be solved using sparse matrix formats to get linear scaling. The other important point is that the two inner loops (initial and descent phases) have to be solved nocc times, where nocc is the number of occupied states that depends on the number of electrons which scales linearly with the number of orbitals N . From this, we can deduce that even if the solution of the quadratic problem takes no time, the overall scaling of the algorithm gives back the "undesired" O(N 3 ) scaling. Provided the matrix operations in the inner loop could be solved with O(N ) scaling methods, and the quadratic problem takes no time, we will have at best an O(N 2 ) scaling. Even though this seems to be an undesired result, it is a first direct approach to solve the DM using a QUBO eigensolver and deserves some further study.

C. DFT calculations
We used the DFT-based code SIESTA 24,25 to obtain all the molecules' Hamiltonians in a tight-binding like representation. SIESTA is a first-principles electronic structure code for molecules and solids which represents wavefunctions as a combination of atomic-like orbitals 26 . The Kohn-Sham equations are solved self-consistently in a linear combination of atomic orbitals basis set and they are restricted to a fix-size cell condition. This framework allows us to obtain a description of the operator associated with the system energy employing a reduced number of orbitals in the basis set with complete multiple-zeta and polarized bases, depending on the required accuracy. In this study, all first-principles calculations were performed with a minimal single-ζ basis for each atom. Description of the interaction between atoms was conducted within the local density approximation approach 27 for the exchange-correlation functional. The integration over the Brillouin zone (BZ) was performed using a Monkhorst sampling in Γ point. The radial extension of the orbitals had a finite range with a kinetic energy cutoff of 50 meV. A separation of 20Å in a cubic simulation box prevents virtual periodic molecules from interacting. We used different numbers of molecules to get the Hamiltonians (H) and overlap (S) matrices needed to test our algorithm. We also computed H and S for benzene in order to test for resiliency regarding degeneracy of eigenstates.

III. RESULTS AND DISCUSSION
We first studied the propagation of the errors with the number of occupied states. The error in the calculated DM is computed as follows: where ρ ex is the DM obtained using regular diagonalization in double precision which we hereby set as our standard. In order to reach chemical accuracy, the error needs to satisfy error ≤ 1.0[kCal/mol]/ max lm (H lm ) as explained in the Appendix. Results showing scaling of error and iteration with system size, occupation, and degeneracy were performed using classical SA algorithm as implemented in the dwave-neal code for the Ocean API 28 . A comparison with QA running on the D-Wave Advantage 4.1 system is offered at the end of this section. From the Pseudocode 1, we can deduce that the error committed on the calculations of the very first computed eigenpairs is susceptible to be propagated towards the higher eigenpairs, hence, making the calculations of DM with higher occupations less accurate. Figure 2 shows the error as a function of the occupation number for a single water molecule with six total atomic orbitals. Three different precision values were used to show that this problem happens regardless of the precision used. In general, we notice that the higher we set the precision for the computation of the eigenvectors, the less error is committed in the construction of the DM across different values of occupations. The larger the number of occupied orbitals, the larger the error that is committed in the computation of the DM. The problem of error propagation worsens with system size. As it is shown in Figure 3, the error gets larger when the total number of orbitals increases, regardless of the precision and number of bits that are used. We notice that, the error increases with the system size which is a sign that a higher precision will be required for larger systems to reach chemical accuracy. Figure 4 shows the number of iterations (total QUBO solutions) that are needed to construct the DM for two and four bits. The total number of iterations is computed by adding up the number of iterations needed to achieve the desired precision for the calculation of each eigenvector. From Figure 3 it does not seem that much is gained by increasing the number of bits, however, as it is shown in Figure 4, we notice that less iterations are needed to reach convergence when using a larger number of bits. Finally, we computed the number of iterations needed to achieve the desired precision to compute every single state for the benzene molecule. This particular case has been selected since, due to its high symmetry (see molecular representation on the inset of figure 5), the electronic structure of benzene is characterized by having several groups of degenerate molecular orbitals (orbitals that are very close in energy). In general, the closer two molecular orbitals are in energy, the more the degeneracy between them. Figure 5 shows the number of iterations for computing a single eigenvector as a function of the proximity in energy to the next eigenvector. We notice that a group of molecular orbitals that are close in energy to the next (< 0.1 eV) need significantly more iterations to be computed at the desired precision. This point is a significantly limiting factor since the electronic structure of several chemical systems including semiconductors and metals will have a high degree of degeneracy. For semiconductors, degeneracy will most likely happen around the band edges, while for metals it will happen around the Fermi level.

IV. QUANTUM VS SIMULATED ANNEALING
In this section we present a comparison between the results obtained using SA with different sweeps and QA performed with the new D-Wave's Advantage 4.1 machine 29 . The chain strength for QA was calculated as 30% of the largest absolute value in the QUBO matrix 30 . This helps ensure minimal chain breaks contributing to erroneous solutions. The default anneal time of 20 microseconds was used. The system DM calculation was performed over samples Hamiltonians taken from different configurations over a 100 Molecular Dynamics (MD) steps. We have extracted H from 20 configurations evenly spaced over the course of the MD trajectory. The DM was computed using an eigenpair precision of 10 −5 . From Figure 6 we can see that, on average, the results using QA have a larger error. Moreover, we have observed that more iterations are needed to converge. In the case of SA a larger number of sweeps tend to converge to results with lower errors. Means for the error between SA and QA are statistically different even for small values of sweeps. One of the mayor issues is the fact that the QUBO matrices have a high dynamic range which causes a "resolution problem." When the number of bits is large there are many small QUBO entries that could have large relative error at the moment of embedding but that might have a large influence in the results. One could scale the QUBO matrix Q so that the minimun entry has a value of 1.0, but this introduces the need of having very strong chain strengths that could introduce yet another source of error. Moreover, we have observed that some chain length are too long due to the high connectivity or matrix density which also depends on the number of bits. Other types of optimization problems have also revealed no advantage using QA vs SA 31 .

V. ALTERNATIVE APPROACHES
In this section we will briefly describe some alternative approaches that could potentially overcome the issues that were found by using a direct method to construct the DM using a QUBO eigensolver. Three different approaches are hereby suggested together with a discussion of possible issues they might have. 1) Using a QUBO based linear solver on the O(N ) method recently proposed by L. Anderson 32 . This algorithm shifts the computational cost of solving ρ to the solution of a linear system (Ax = b) using Conjugate Gradient method which can be replaced by: and can be easily formulated as a QUBO problem and solved on a QA. The problem that this method could have is the fact that there might be several iterations (QUBO solves) needed to reach the desired precision of x. This will add yet another inner loop into the algorithm.
2) We can use Green Functions (G(z)) method and find ρ by integrating over the complex contour (with variable z ∈ C) containing the eigenspectrum of H. For this proposition we would need to have a method to invert a complex matrix using a QUBO formulation. Provided we can compute the inverse of a matrix using a QUBO formulation, ρ would be determined as: where G(z) = (H − zI) −1 and C is the semicircle above the real axis. Inverting a matrix A can be simply done by solving N linear systems of the form Ax i − e i = 0, with e i being the i-th canonical vector, and x j the jth column of A −1 . This can be easily formulated as a QUBO provided both, the solution and the matrix to be inverted are real. If this is not the case, as for the Green Function method here proposed, the formulation becomes more complicated and, as far as we know, no previous work has addressed this issue.
3) Finally, we could think about solving for ρ "all at once" from a QUBO problem where ρ would be written as a long vector x containing N × N entries, such that: where . is the operator that takes the "floor" value of a real number. In this case, the function to be minimized will be the total energy where in this case 1 k = δ ij , where i = k/N + 1, and j = k − N k/N . Although equation 17 can be easily formulated as a QUBO problem, the energy needs to be minimized under certain constraints given by the properties of ρ; some of which can be extremely difficult to formulate as a QUBO. Moreover, the size of the QUBO problem to be solved will scale as N 2 × N 2 , probably opening up other bottlenecks in the algorithm even if the annealing process require no time.

VI. CONCLUSION
We have demonstrated that the DM can be computed by using a QUBO based eigensolver and summing up the self outer products of the eigenvectors up to the number of occupied states. Although the results obtained with this direct approach are encouraging, many issues still needs to be addressed. We reported on the problem of the O(N 2 ) operation inside loops that compute the egenpairs which leads back to an O(N 3 ) scaling. By experimenting with different occupation numbers, we have seen that the error rapidly increases with the occupation number regardless of the precision at which we compute each eigenvector. We have also seen that the error increases with the system size regardless of the number of bits used to compute the eigenpairs. We observed that by increasing the number of bits the number of iterations to converge decreases significantly. We have identified a problem where the number of iterations to reach precision increases for degenerate states. Finally we showed that QA does not show any advantage when compared to SA for computing DM within this QUBO formulation and more research needs to be done to fully understand the QA vs SA comparison.
We have also proposed other alternative methods to compute the DM and analyze some possible issues. These methods include using linear scaling algorithms based on linear systems that could be solved using QUBO solvers; using Green functions obtained from a QUBO solver that then could be used to compute the DM by a close integral surrounding the eigenspectrum through the complex plane; and solving the full density matrix "all at once" from a QUBO problem using a vectorized form of the DM. Even if all these alternative methods have challenges there is still an opportunity for further improvements. An in depth study of these other methods will be the subject of future work. Multiplying x and y in the binary representation leads to the following operation: x × y = q T x (v T ⊗ v)q y An example with 3-bits follows. Given we have: x = 0.5 ≡ (0, 1, 0), and y = 0.25 ≡ (0, 0, 1) The error in ρ is computed as: ∆ρ = ||ρ − ρ ex ||; where ρ ex is the exact DM constructed using regular diagonalization methods. To ensure chemical accuracy we need to ensure that the total energy error remains below 1.0 kCal/mol. This means: ∆E < 1.0kCal/mol.