Master’s thesis Electronic structure calculations on quantum computers Hugo A˚stro¨m 2022 Supervisor: Dage Sundholm Examiners: Prof. Dage Sundholm & Dr. Susi Lehtola University of Helsinki Faculty of Science Tiedekunta – Fakultet – Faculty Koulutusohjelma – Utbildningsprogram – Degree programme Opintosuunta – Studieinrikting – Study track Tekijä – Författare – Author Työn nimi – Arbetets titel – Title Työn laji – Arbetets art – Level Aika – Datum – Month and year Sivumäärä – Sidoantal – Number of pages Tiivistelmä – Referat – Abstract Avainsanat – Nyckelord – Keywords Säilytyspaikka – Förvaringställe – Where deposited Muita tietoja – Övriga uppgifter – Additional information Faculty of science Theoretical and computational methods Chemistry Electronic structure calculations on quantum computers Master's thesis I discuss recent work regarding electronic structure calculations on quantum computers. I introduce quantum computing and electronic structure theory, and then discuss different mappings from electrons and excitation operators, to qubits and unitary operators, mainly Jordan–Wigner and Bravyi–Kitaev. I discuss adiabatic quantum computing in connection to state preparation on quantum computers. I introduce the most important algorithms in the field, namely, quantum phase estimation (QPE) and variational quantum eigensolver (VQE). I also mention recent modifications and improvements to these algorithms. Then I take a detour to discuss noise and quantum operations, a model for understanding how quantum computations fail because of noise from the environment. Because of this noise, quantum simulators have risen as a tool for understanding quantum computers and I have used such simulators to do electronic structure calculations on small atoms. The algorithm I have used, QPE, yields the exact result within the employed basis. As a basis I use numerical orbitals, which are very robust due to their flexibility. 11.2022 34 electronic structure, quantum computing, quantum phase estimation, variational quantum eigensolver e-thesis Hugo Åström Contents 1 Introduction 2 1.1 Quantum computers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 The electronic structure problem . . . . . . . . . . . . . . . . . . . . . . 7 2 Mapping of the problem 9 2.1 Electron to qubit mappings . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Hamiltonian Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3 State Preparation 17 4 Quantum Algorithms 18 4.1 Quantum Phase Estimation . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Variational Quantum Eigensolver . . . . . . . . . . . . . . . . . . . . . . 22 5 Noisy Quantum Computation 24 5.1 Quantum Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.2 Simulating Quantum Computers . . . . . . . . . . . . . . . . . . . . . . 27 6 Results 27 6.1 Numerical Orbitals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.2 Quantum Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 7 Discussion 30 2 1 Introduction List of abbreviations FCI Full Configuration Interaction VQE Variational Quantum Eigensolver QPE Quantum Phase Estimation IQPE Iterative Quantum Phase Estimation SE Schro¨dinger equation SCM Symmetry Configuration mapping UCC Unitary Coupled Cluster UCCSD Unitary Coupled Cluster Singles Doubles JW Jordan–Wigner BK Bravyi–Kitaev BKSF Bravyi–Kitaev Superfast STO Slater Type Orbital GTO Gaussian Type Orbital Hartree–Fock HF CC Coupled Cluster 1 Introduction The electronic structure problem is arguably the most fundamental problem in chemistry. Its main task is to find the wave function of a chemical system, such as an atom or a molecule. The wave function determines all the chemical properties of the system. The electronic structure problem is the underlying problem that has to be solved when it comes to all chemically interesting applications such as drug design, finding efficient catalysts or constructing molecules for capturing solar energy. Therefore, finding efficient meth- ods for solving electronic structure problems can be considered of as some kind of ”holy grail” in chemistry. Theoretical and computational chemists have been developing ap- propriate electronic structure methods for decades. The large underlying challenge lies in the fact that with increasing number of particles, the Hilbert space of the problem in- creases exponentially, which means that the Schro¨dinger equation of only the simplest model systems can be solved analytically. Thus, approximations have to be made when one wants to perform useful applications. The ultimate method is the full configuration interaction (FCI) method [1, 2], which provides the exact solution within the employed basis. However, the FCI method scales factorially with system size according to Weyl’s dimension formula [1]. The best ab initio methods, such as second order Møller–Plesset perturbation level (MP2) or coupled cluster singles and doubles with perturbative treat- ment of the triples level (CCSD(T)) [1] are used in large scale applications. They scale as N5 to N7, where N represents the system size. The accurate FCI method can only be 3applied to small systems, and is used mainly for benchmarking. Even though CCSD(T) calculations can be performed on large molecules consisting of several hundreds of atoms [7], the obtained energies differ significantly from total energies that would be obtained at the FCI level due to the approximations made. Advances in quantum physics naturally lead to new approaches for solving electronic structure problems. The most recent and important one is quantum computing. When thinking about computational sciences at a very fundamental level, one might ask what kind of physical systems can be used to perform a computation? For a very long time only classical systems were considered for this task, like the classical bit in an av- erage laptop, which usually makes use of some electromagnetic phenomenon to be either in the state 0 or 1, and together, many bits can perform advanced computations. It wasn’t until the 80s when pioneers like Richard Feynman thought of making the physical sys- tem quantum mechanical [5, 6]. The motivation was that to be able to simulate quantum mechanical systems (such as atoms and molecules) in an efficient manner compared to classical computations, we also need to make the computer quantum mechanical. This seems reasonable since the quantum mechanical world is so different from the classical world. It should be (and is) very difficult to simulate it using classical physics. How- ever, if we make the computer itself quantum mechanical, the story could be different. Quantum chemistry is a field which deals with highly quantum mechanical systems and therefore it has risen as a natural field of application for quantum computers. It is however well known that making quantum computations tractable is an extremely difficult task due to decoherence. Decoherence is the fact that quantum mechanical systems are very sen- sitive and difficult to control. They tend to interact with the environment and instead of staying in the quantum mechanical world, they collapse to some classical state or become entangled with the environment, which causes the computations to fail. For this reason it has taken such a long time for quantum computing to get serious attention from the field of quantum chemistry. More on the concepts of decoherence and entanglement can be found e.g. in refs [3, 4]. Now that quantum computing is becoming more tractable, the need for quantum algo- rithms grows. I have studied a few existing quantum algorithms for solving the electronic structure problem, focusing on the helium atom. First I have computed the one- and two- electron interactions on a classical computer, using a numerical basis set. There are sev- eral advantages to using numerical basis sets instead of analytical functions [45]. They are more robust and flexible, which means that they can be used as exact orbitals. This is an advantage over analytical basis functions such as Gaussian type orbitals (GTO) or Slater 4 1 Introduction type orbitals (STO), which are approximations of exact orbitals. I have then implemented the integrals on a quantum computer simulator and run the quantum phase estimation al- gorithm for different basis set sizes. I conclude the introduction with some basics of the theory behind quantum computers and a short recap of classical calculations. After this I focus on how to map electronic structure problems to a quantum device, i.e., rewriting it in terms of qubits and quantum gates. I consider a few different mappings and transfor- mations. Their task is trying to minimize the number of qubits and gates needed. After that I present two of the most important quantum algorithms for finding the electronic structure, namely quantum phase estimation and variational quantum eigensolver. This is followed by a section on the largest challenge facing quantum computing today, which is noise. Noise is when the quantum computer interacts with the environment, causing the computation to fail. I discuss a theory for modeling noisy quantum channels which produce errors. After this I present the results of my work, and end with a discussion and conclusions. 1.1 Quantum computers There are many differences between how classical computers and quantum computers work in practice. An example is the fact that in quantum computing one measures the result, which makes the quantum computing itself an experimental science. There are also very important differences on the fundamental theoretical level. Two of them are the superposition principle and quantum entanglement. Classical computers use bits as their most basic unit, which can take on values 0 and 1. The corresponding units on quantum computers are called quantum bits, or qubits. A qubit is a two-level quantum mechanical system and they work quite differently from their classical counterpart. A physical qubit can be realized by any quantum mechanical system with two different states, like an atom with two energy levels or a superconducting circuit where the current can go either clockwise or counterclockwise. For the purposes of this thesis however, the qubits will only be abstract mathematical objects, which lets one explore the theory of quantum information. As all other quantum mechanical systems, they are represented as vectors in a Hilbert space, which for a qubit is two-dimensional. Therefore a qubit in the general case is written as a superposition of the two states |0⟩ and |1⟩ as |q⟩ = a|0⟩+ b|1⟩, (1.1) 1.1 Quantum computers 5 where a, b ∈ C are coefficients which satisfy |a|2 + |b|2 = 1. (1.2) Consequently, a qubit can be in the state 0 and 1 simultaneously contrary to the bit, which is always either 0 or 1. A measurement in the basis {|0⟩, |1⟩} will collapse the qubit to the state |0⟩ with probability |a|2 and |1⟩ with probability |b|2. This is the superpo- sition principle. The other important property is quantum entanglement, a phenomenon where two or more qubits are dependent on each other, or entangled. One can not know anything about one specific qubit, without affecting the other qubit(s). It is perhaps best demonstrated by an example. The quantum state |ψ⟩ = |01⟩+ |10⟩√ 2 (1.3) is a highly entangled state. It is impossible to isolate one qubit from the state and write it as a product state. If we measure the value of one qubit, we immediately know the value of the second one. These two properties of quantum mechanical states result in new ways to think of computation both concerning storage of information, and manipulating the information through algorithms. Regarding information storage, one can again think of a string of qubits in a superposition state |ψ⟩ = a0|0⟩+ a1|01⟩+ a3|10⟩+ a4|11⟩. (1.4) This qubit register consisting of two qubits stores four numbers, a0, a1 . . . , and in general, n qubits can store 2n numbers while n classical bits can only store n numbers. This exponential increase in information storage is one of the properties which makes quantum computers desirable. The next question is then: How do we perform computations with these qubits? The qubit being a vector on a Hilbert space, we can of course operate on them with operators on the space. The operators are represented as unitary matrices and are called quantum gates, or just gates. Some of the most common gates are found in table 1.1. These are the quantum equivalents of classical logic gates. From these gates we can put a qubit in a superposition state using the Hadamard gate, put two qubits in an entangled state using the CNOT gate and make up gates which perform 6 1 Introduction Name Symbol Matrix representation Identity I ( 1 0 0 1 ) Pauli-X X ( 0 1 1 0 ) Pauli-Y Y ( 0 −i i 0 ) Pauli-Z Z ( 1 0 0 −1 ) Hadamard H 1√ 2 ( 1 1 1 −1 ) CNOT  1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0  Table 1.1: The most common quantum gates and their symbol and matrix representations. arbitrary rotations of angle θ of the qubit around an axis nˆ Rnˆ(θ) = exp (−iθnˆ · σ⃗ 2 ) = cos ( θ 2 ) I − i sin ( θ 2 ) (nxX + nyY + nzZ), (1.5) where σ⃗ = (X, Y, Z) denotes the vector of Pauli matrices. In general, any unitary trans- formation U of quantum mechanics, which performs the evolution |ψ⟩f = U |ψ⟩i, (1.6) where |ψ⟩f is the final state and |ψ⟩i is the initial state, can be implemented on a quantum computer using a set of universal quantum gates. Finally, at the end of the computation one must perform a measurement in the basis to get the result. The measurement can be represented by a set of measurement operators {Mm}, and the probability of getting the result m is p(m) = ⟨ψ|M †mMm|ψ⟩. (1.7) As an example, in the computational basis one can measure a qubit using the projection 1.2 The electronic structure problem 7 operators |0⟩⟨0| and |1⟩⟨1|, which satisfy the completeness relation∑ m M †mMm = I. (1.8) As mentioned, the implementation of the operators and qubits is an experimental science, and many different methods have been proposed like trapped ion [30]-[34], superconduct- ing quantum computers [35, 36] and photon systems [37, 38]. Now we have everything we need to start building quantum algorithms; the set of qubits, which form a qubit regis- ter, and the quantum gates. 1.2 The electronic structure problem Before completely turning the attention to quantum computations, I will revisit classical quantum chemistry computations. This is necessary for understanding what the electronic structure problem is, why it is so difficult to solve on a classical computer and how it can be implemented on a quantum computer. The electronic structure problem is quite simple to state: solve the Schro¨dinger equation (SE) and you know everything there is to know about the system. However, actually solving the equation is a much harder problem. The time-independent, non-relativistic SE can be stated as follows H|Ψ⟩ = E|Ψ⟩, (1.9) where H is the electronic Hamiltonian, |Ψ⟩ is the quantum state and E the energy corre- sponding to that state. In the formalism of first quantization the electronic Hamiltonian is H = −1 2 N∑ i ∇2i − n∑ j N∑ i Zj rij + N∑ i>j 1 rij , (1.10) where the first term is the kinetic energy operator, the second term gives the electron- nuclear interaction and the third the electron-electron interaction. N is the number of electrons, n the number of nuclei and Zj the charge of nucleus j. The state |Ψ⟩ is an abstract vector in a Hilbert space and can be expanded in some basis {|ψi⟩}ni=1 in the following way |Ψ⟩ = n∑ i=1 cn|ψn⟩, (1.11) where ci ∈ C are the expansion coefficients. The wavefunction we want to deal with is 8 1 Introduction then obtained as an inner product Ψ(x) = ⟨x|Ψ⟩ = n∑ i=1 cn⟨x|ψ⟩ = n∑ i=1 cnψn(x), (1.12) where the functions ψn(x) are the basis functions. Based on this, all the different computational approaches have individual ways of finding the wavefunction that solves the SE and FCI the one I focus on. Here the state is expanded as a Slater determinant |Φ⟩ = 1√ N ! ∣∣∣∣∣∣∣∣∣∣ ϕ1(1) ϕ2(1) . . . ϕN(1) ϕ1(2) ϕ2(2) . . . ϕN(2) ... ... . . . ... ϕ1(N) ϕ2(N) . . . ϕN(N) ∣∣∣∣∣∣∣∣∣∣ , (1.13) which is a product of one-particle spin functions ϕi(j) of electron j (the basis functions). The Slater determinant preserves the correct anti-symmetric properties of the state. By including all possible determinants made up from the functions in the basis set, one makes sure that all possible excitations are considered in the superposition state. Therefore, the solution will be exact within the basis set. The FCI state is |FCI⟩ = ∑ i Ci|Φi⟩, (1.14) where |Φ⟩ are the Slater determinants. The energy of the system is E = ∑ i ∑ j CiCj⟨Φi|H|Φj⟩, (1.15) where we have the normalization condition∑ i |Ci|2⟨Φi|Φi⟩ = 1, (1.16) with the orthogonality condition ⟨Φi|Φj⟩ = 0, i ̸= j. (1.17) 9This allows us to write down the Lagrangian for FCI as L = ∑ i ∑ j CiCj⟨Φi|H|Φj⟩ − E (∑ i |Ci|2⟨Φi|Φi⟩ − 1 ) (1.18) and differentiating with respect to Ci gives an eigenvalue matrix equation HC = EC, (1.19) where the Hamiltonian matrix elements are Hij = ⟨Φi|H|Φj⟩. (1.20) The number of Slater determinants is Ndet = ( M N/2 + S )( M N/2− S ) , (1.21) where N is the number of electrons and M is the number of orbitals. This large number of determinants is what makes classical computations so inefficient. A detailed discussion on classical quantum chemistry methods can be found e.g. in refs [1, 2] or other similar textbooks. It has now been demonstrated that FCI indeed scales factorially, which makes it useful mainly for benchmarking. The more tractable methods mentioned (CCSD(T) and MP2) yield energies which differ significantly from the FCI energy. Therefore, we now turn our attention to quantum computations instead. 2 Mapping of the problem It is necessary to map the problem to the quantum device before the actual calculation. Thus, one has to rewrite the problem in such a way that it can be implemented on a quantum computer. The idea is to map the states of the quantum system onto the qubit register. The Hamiltonian is then mapped to a combination of unitary operators. In this section a few approaches to map the space of states to qubits is presented. The aim is to minimize the number of qubits needed, which is especially important for today’s early quantum computers. Then a few different ways to implement the Hamiltonian will be dis- cussed. The operator is mapped to a sum of Pauli operators which are easily implemented on a quantum device. The methods use as few operators as possible to generate shallow 10 2 Mapping of the problem algorithms, while preserving the antisymmetric properties of the quantum states. 2.1 Electron to qubit mappings The electronic structure problem starts with defining the Hamiltonian. The non-relativistic electronic Hamiltonian for N particles in second quantization formalism can be written as H = ∑ pq hpqa † paq + 1 2 ∑ pqrs vpqrsa † pa † qasar, (2.1) where the weights, hpq and vpqrs, are one- and two electron integrals respectively. The one electron weights are obtained as hpq = − ∫ χ∗p(r) ( 1 2 ∇2 + ∑ i Zi |Ri − r| ) χq(r)dr, (2.2) which is the expectation value of the kinetic energy term and the Coulomb attraction term between the electron and the nuclei, as in equation 1.2. The two-electron interaction terms are given by vpqrs = ∫ χ∗p(r1)χ ∗ q(r2)χr(r2)χs(r1) |r1 − r2| dr1dr2. (2.3) The integrals are calculated on a classical computer and implemented on the quantum device. Since the integrals have to be calculated only once, they do not lead to any over- head to the simulation. In equation 2.1 ai and a † i are fermionic annihilation and creation operators, respectively, i.e., they annihilate/create an electron in the spin orbital i. They obey the anticommutation relations aia † j − a†jai = {ai, a†j} = δij (2.4) {ai, aj} = {a†i , a†j} = 0, (2.5) which is important because it means that the Fermi statistics of fermions is automatically encoded on the operator level. The second quantization formalism is more convenient than the first quantization of equation 1.2 when doing quantum simulations because it reduces the number of qubits needed [15]. The Hamiltonian acts on the space of states, the Fock space. It is defined as F = K⊕ n=0 Hn,K , (2.6) 2.1 Electron to qubit mappings 11 which denotes a direct sum of the Hilbert spacesHn,K forK spin orbitals. n is the number of particles. An arbitrary basis state in the Fock space can be written as |ψ⟩ = |n1, n2, . . . , nN⟩, ni ∈ {0, 1}, (2.7) where each ni is the occupation number of the spin orbital i. It is 0 or 1 if the spin orbital is unoccupied or occupied respectively. Using the creation operators a generic state can be written as |n1, n2, . . . ⟩ = (a†1)n1(a†2)n2 . . . |0⟩. (2.8) Considering the anticommutation relations in equations 2.4 and 2.5, the action of ai and a†i on that state is a†i | . . . , ni, . . . ⟩ = δni,0(−1) ∑ j 0. This operator has a form which can be implemented on a quantum computer. I will not however go further into the implementation of the operator and instead focus on the more used Trotterization method. These methods for implementing the Hamiltonian introduce errors. One source for errors is that the operator in the truncated Taylor series method is not exactly unitary, since the series expansion is cut off. However, it is close to unitary and the error can be controlled. Implementation of Hamiltonians remains an active area of research and other approaches such as qubitization [19] have been suggested. Thus, we now have the tools to prepare Hamiltonians on a quantum computer. 17 3 State Preparation The last step to consider before discussing the quantum algorithms is preparation of the initial state, which is one of the strengths of quantum computers. Preparing a good initial state with sufficient overlap with the true ground state can be a crucial starting point to many algorithms. For regular closed shell singlet molecules the Hartree–Fock (HF) state is usually a good guess. Around the equilibrium geometry the overlap tends to be about |⟨ΨHF|Ψexact⟩|2 ≈ 0.9, but this overlap can decrease for molecules beyond these criteria [10]. My focus will lie on adiabatic state preparation, since it is probably the best known state preparation technique in quantum computing. The idea stems from the adiabatic theorem, which makes the following statement: A physical system remains in the eigenstate of a Hamiltonian when acted upon with a perturbation if the evolution is sufficiently slow and there is a gap between the eigenvalue of the system and the rest of the spectrum of the Hamiltonian. One can now make the Hamiltonian time dependent H(ξ), ξ ∈ [0, 1], (3.1) where the ground state of H(0) is some known state that can be prepared on a quantum computer, and the ground state of H(1) is the desired state. H(1) could e.g. be obtained by slowly turning on some perturbation say H(1) = H0 +H ′, (3.2) where H0 = H(0) and H ′ is a perturbation. Now one can make use of the time-dependent SE i d dt |Ψ(t)⟩ = H(t/tf )|Ψ(t)⟩, t ∈ [0, tf ] (3.3) to evolve the system along a determined path. Given that the spectrum of the Hamiltonian along this path is gapped and the evolution is sufficiently slow, the final state will be the ground state of the H(1). The total time, tf , depends on the gap between the two lowest eigenvalues of the Hamiltonian along the path, and a parameter ε as [18] tf > ε g2min , (3.4) 18 4 Quantum Algorithms where gmin = min(E1(t)− E0(t)) (3.5) is the gap and ε = max ∣∣∣∣⟨ψ1(t)|dH(t/tf )dt |ψ0(t)⟩ ∣∣∣∣, (3.6) where |ψ0(t)⟩ is the ground state and |ψ1(t)⟩ is the first excited state. Thus, if tf is constant or grows polynomially, there will be no excess overhead to the simulation. On the other hand, large tf might require deep circuits, which is still a problem as mentioned. Another challenge could lie in the fact that the spectrum along the whole path might not be known and degeneracies could lead to the wrong final state. An algorithm designed for state preparation for molecules with open shell character has also been suggested [20]. It is based on Serber construction [47, 48] and has the advantage that the circuit depth is only 2. An alternative technique is variational state preparation. This will be discussed in more detail in section 4.2 in connection to VQE. This now leads us to the subject of quantum algorithms. 4 Quantum Algorithms The algorithm of most interest is Quantum Phase Estimation (QPE) and state preparation and implementation of the Hamiltonian are the starting points for that. Now that we have these tools, we can start to explore the algorithm. As mentioned, the exact solution to the SE within the employed basis is classically given by FCI. The corresponding algorithm for quantum computers is QPE [9]. It is an al- gorithm for determining eigenvalues of unitary operators. Like FCI, it gives the exact solution to the SE within the basis. The advantage is that it scales polynomially with system size, compared to factorially for FCI. As mentioned, the starting points are the initial state and the implementation of the Hamiltonian. The task of the algorithm is to project the initial state onto an eigenstate of the Hamiltonian. A measurement then yields the corresponding eigenvalue of the eigenstate exactly. The success rate is proportional to the overlap |⟨ψinitial|ψexact⟩|2 between the initial state and the exact state. Thus, state preparation can be of great value. Given exponentially decreasing overlap between the HF state and multiconfigurational states, an exponential numbers of runs of the algorithm is also needed. For singlet closed-shell molecules around their optimal geometry the HF state is usually sufficient. 4.1 Quantum Phase Estimation 19 The problem with QPE is that it is sensitive to noise. Therefore, VQE is another important algorithm, especially for near time quantum computers because it is more robust. It takes advantage of quantum computers’ ability to prepare and measure a quantum state, while utilizing classical computers to keep the quantum circuits from becoming too deep. It is based on Ritz variational principle [21] for quantum mechanics, which is a common method in many classical algorithms as well. The two algorithms will be discussed in detail in this section. 4.1 Quantum Phase Estimation QPE is an algorithm for finding eigenvalues of unitary operators on quantum computers. The eigenvalue of a unitary operator U can be represented by a point on the unit circle in the following way U |u⟩ = u|u⟩ = eiφ|u⟩. (4.1) QPE estimates the value of φ with very high accuracy. The Hamiltonian 2.1 is an Her- mitian operator. Thus, we again construct a unitary operator by exponentiating H as in equation 2.28 U |u⟩ = e−iH |u⟩ = eiφ|u⟩, (4.2) where φ, the ”phase”, is an eigenvalue of H . The algorithm consists of two qubit registers. One with t qubits, the ancilla register, from which we in the end will read out the eigenvalue, and another with m qubits, representing the initial state |u⟩. The computer starts in the product state |ψ⟩ = |0⟩⊗t ⊗ |u⟩ = |0⟩⊗t|u⟩. (4.3) First, we apply to the first register a sequence of Hadamard gates, H⊗t, to put the computer in the equal superposition state |ψ⟩ = 1√ N N−1∑ j=0 |j⟩|u⟩. (4.4) Next, we apply controlled U j-operations with the j:th qubit in the first register as the control qubit and |u⟩ as the target |ψ⟩ = 1√ N N−1∑ j=0 |j⟩U j|u⟩. (4.5) 20 4 Quantum Algorithms We can write |u⟩ as a superposition of states in the computational basis |ψ⟩ = 1√ N N−1∑ j=0 |j⟩U j ∑ k ck|ϕk⟩ (4.6) = 1√ N ∑ k ck N−1∑ j=0 |j⟩U j|ϕk⟩ (4.7) = 1√ N ∑ k ck N−1∑ j=0 |j⟩(λk)j|ϕk⟩, (4.8) where λk is the k:th eigenvalue of U , with eigenvector |ϕk⟩. If we now write the eigen- values as eiφk according to 4.2, we obtain |ψ⟩ = 1√ N ∑ k ck|ϕk⟩ N−1∑ j=0 eiφkj|j⟩ = 1√ N N−1∑ j=0 eiφkj|j⟩|u⟩. (4.9) We can now make an interesting observation: the first qubit register has the form of quantum Fourier transformation (QFT) of the state |φk⟩. We apply the inverse QFT QFT† 1√ N N−1∑ j=0 eiφkj|j⟩ = |φk⟩ (4.10) to obtain a state |ψ⟩ = |φk⟩|u⟩ (4.11) where the first qubit register contains the t-bit representation of λk. We then measure it in the computational basis to obtain the eigenvalue λk with probability |⟨u|ϕk⟩|2 = |ck|2. The circuit diagram showing the first part of the algorithm, up until the inverse QFT, can be found in figure 4.1. A more detailed description of the QFT subroutine can be found in ref [3]. 4.1 Quantum Phase Estimation 21 . . . . . . . . . . . . . . . first register t qubits |0⟩ H |0⟩+ e2πi2t−1φ|1⟩ |0⟩ H |0⟩+ e2πi22φ|1⟩ |0⟩ H |0⟩+ e2πi21φ|1⟩ |0⟩ H |0⟩+ e2πi20φ|1⟩ second register |u⟩ U20 U21 U22 U2t−1 |u⟩ Figure 4.1: Circuit diagram outlining the first part of the QPE algorithm. The accuracy of the algorithm depends on the number of qubits in the first qubits register. One can modify the algorithm such that the accuracy is independent of the number of qubits used. This procedure is called iterative QPE (IQPE) and it reduces the number of qubits needed [10]. As the name suggests, the algorithm proceeds in an iterative manner. In each iteration one digit of the eigenvalue is measured, starting from the last and working towards the first ones. One problem might be that one has to prepare the trial state in each iteration if one does not want to maintain throughout the whole algorithm, but if the state is easily prepared then it is no problem. Otherwise IQPE works in the same way as QPE and the circuit for iteration k can be found in figure 4.2. The rotation Rz(ωk) and the second Hadamard gate form the inverse QFT for the first qubit. |0⟩ H Rz(ωk) H φk |u⟩ U2k−1 Figure 4.2: The circuit for iteration k of IQPE. This version however does not remove the bottleneck of the algorithm, namely the imple- mentation of the operator U . 22 4 Quantum Algorithms 4.2 Variational Quantum Eigensolver QPE is not useful for doing computations on near-term quantum computers. VQE how- ever, is more robust against noise. It is a hybrid quantum-classical algorithm. It utilizes quantum computers’ ability to prepare and measure quantum states but takes advantage of classical computers to avoid deep circuits. The algorithm starts by finding a good trial state, or ansatz, |Ψ(θ⃗)⟩ and a trial quantum circuit U(θ⃗). They are parameterized by a set θ⃗ = (θ1, . . . , θn) of parameters we want to optimize. Contrary to QPE, this method is iterative and in each iteration we prepare a new state and new circuit. At the end of the iteration the energy is measured. Then the parameters are optimized with respect to the energy on a classical computer. That way a new set of parameters are obtained for the next iteration. The outline of the algorithm is as follows: 1. Prepare the circuit U(θ⃗). 2. Run the circuit to get |Ψ(θ⃗)⟩. 3. Measure expectation value of H . 4. Optimize the parameters θ⃗. 5. If it has not converged go back to step 1. There are different ways of preparing the ansatz and the circuit. The original approach was unitary coupled cluster (UCC) [22]. It is based on the classical coupled cluster (CC) method [42, 43], but modified to be unitary. UCC offers several advantages over the non- unitary version [44] but there is no way to implement it efficiently on a classicl computer. Now the new state is obtained as |Ψ(θ⃗)⟩ = U(θ⃗)|Φ⟩, (4.12) where the ansatz |Φ⟩ is e.g. an HF determinant. To construct the circuit we will use the unitary cluster operator, which is given by |ΨUCC⟩ = exp[T − T †]|Φ⟩, (4.13) where T is the cluster operator when we have n electrons T = n∑ k=1 Tk. (4.14) 4.2 Variational Quantum Eigensolver 23 The cluster operator is a sum of excitation operators T1 = ∑ i∈occ k∈virt tai a † aai (4.15) T2 = 1 4 ∑ i,j∈occ a,b∈virt tabij a † aa † baiaj (4.16) ... Tn = 1 (n!)2 ∑ i1,...,in∈occ a1,...,an∈virt ta1···ani1···in a † a1 · · · a†anai1 · · · ain , (4.17) i.e., Tm excites m electrons from occupied orbitals to virtual orbitals. Here the indices i, j, k, l, . . . indicate occupied orbitals and a, b, c, d, . . . indicate virtual orbitals. We can write the circuit as U(θ⃗) = exp [∑ k θk(Tk − T †k ) ] (4.18) and make a Suzuki–Trotter decomposition U(θ⃗) ≈ (∏ k e θk n (Tk−T †k ) )n . (4.19) Usually it suffices to go as high as single and double excitations (UCCSD), i.e., T ≈ T1 + T2 (4.20) to obtain good results. There have been several modifications made to the originally suggested UCCSD approach [23, 24], which are also based on UCC. A modification suited for strongly correlated systems has been suggested, called ADAPT-VQE [25]. It starts by defining a set {tpq , tpqrs, . . . } (4.21) of excitation operators and the HF state. In each iteration the gradient of the energy with respect to each operator in the set is calculated ∂E(n) ∂θi = 〈 ψ(n) ∣∣[H,Oi]∣∣ψ(n)〉, (4.22) 24 5 Noisy Quantum Computation where n denotes the iteration and i the operator index. The operator with the largest gradient is then applied to the trial state, the parameters are optimized and a new energy is obtained. If the gradient is converged the algorithm is done. The final state at convergence is |ψadapt⟩ = (eτN )(eτN−1) . . . (eτ1)|ψHF⟩, (4.23) where each eτn is computed at iteration n. The idea behind this approach is that the algorithm builds the trial state itself step by step, so that the system itself decides the state. They showed that the performance of the method is very good for both uncorrelated and correlated systems, outperforming the original approach. A hardware-efficient method has also been proposed [13]. It uses alternating single qubit rotations, U q,d(θ⃗) and entangling operations, Uent. Here q denotes the qubit index and d the circuit depth. The initial state is the state |00 . . . 0⟩, which yields the trial state |ψ(θ⃗)⟩ = N∏ q=1 ( U q,d(θ⃗) ) Uent N∏ q=1 ( U q,d−1(θ⃗) ) . . . Uent N∏ q=1 ( U q,0(θ⃗) )|00 . . . 0⟩, (4.24) where the single qubit rotations have the form U q,i(θ⃗) = Rz(θ q,i 1 )Rx(θ q,i 2 )Rz(θ q,i 3 ), (4.25) where the Rn are the rotation operators from equation 1.5. This method is useful because it yields shallow circuits. Kandala et al. got good results for their studied molecules, however, there are no results of the method giving good results for highly correlated systems. Moreover, it requiresN(3d+2) parameters, which is higher than that of ADAPT- VQE, where the number of parameters is equal to the number of operators in the set. An improvement to the ADAPT-VQE has also been suggested [27], referred to as a hardware-efficient ADAPT-VQE, or qubit-ADAPT-VQE. It replaces the excitations op- erator set with Pauli strings and introduces a completeness relation for the operators. They show that this results in the operator set scaling linearly with the number of qubits, resulting in a more hardware-efficient algorithm and lower-depth circuits. 5 Noisy Quantum Computation Before discussing the applications of the quantum algorithms, I will take a small detour into the theory of noise. The model of noise is somewhat mathematical and might not be relevant for quantum chemistry computations, but offers insight as to what goes on when 5.1 Quantum Operations 25 the quantum gates don’t work as planned. This will then lead to the subject of quantum computer simulators, which have been the focus of my experiments. One of the biggest challenges with today’s quantum computers is that they are not closed quantum systems, but open quantum systems. This means that they interact with their environment, and this interaction produces noise. This noise can be e.g. heat coming in to the system from the surroundings, or electromagnetic fields interacting with the qubit. It can cause different problems, like collapsing of the superposition of the qubit or a bit flip. These cause the computation to fail. We can study noise both from a theoretical point of view with a mathematical model, and from a practical point of view using quantum computer simulators. The theoretical viewpoint requires construction of a new formalism, which is discussed in this section. This is followed by a discussion on quantum computer simulators. 5.1 Quantum Operations To understand noise we can construct a model called the quantum operations formalism. It is now easier to work in the density matrix formulation. If a quantum system is in the state |ψi⟩ with probability pi, we have an ensemble of pure states, {pi, |ψi⟩}. Then we can define a density matrix, or density operator as ρ = ∑ i pi|ψi⟩⟨ψi|. (5.1) Just as in the state vector representation, all the information about the quantum system is contained within the density matrix and all the postulates of quantum mechanics can be restated in the density operator formalism. We begin by stating that there is a noisy interaction U which is not necessarily a unitary interaction. This is visualized in figure 5.1 where ρenv is the density matrix of the environment and ρ′ is the final state of the system. We can define a map E ρ′ = E(ρ), (5.2) which we call a quantum operation. This model of noise focuses a lot on how to develop a formalism for these operations. Next we can introduce the state space of the interacting environment |ek⟩. The environ- ment starts in some state |e0⟩. Then to obtain the density matrix of the system we trace out the environment from the composite system Uρ⊗ ρenvU †, using partial trace E(ρ) = Trenv ( Uρ⊗ ρenvU † ) (5.3) 26 5 Noisy Quantum Computation = ∑ k ⟨ek|Uρ⊗ |e0⟩⟨e0|U †|ek⟩ (5.4) = ∑ k EkρE † k. (5.5) The Ek = ⟨ek|U |e0⟩ are operators on the space of states of our system. ρ U ρ′ ρenv Figure 5.1: The density matrix of the system interacts through a black box with the envi- ronment to become a new state. The quantum operation satisfies three properties: 1. The probability of the process represented by E occurring is represented by the trace Tr(E). 2. For a set of probabilities {pi}ni=1 and density matrices {ρi}ni=1, E satisfies E ( n∑ i=1 piρi ) = n∑ i=1 piE(ρi). (5.6) 3. E is completely positive. If these three properties are satisfied, then the quantum operation can be represented as in equation 5.5 [3], where Ek satisfy N∑ k=1 EkE † k ≤ I. (5.7) The quantum operations formalism is good for describing noisy quantum channels. As an example, one can consider the bit flip channel. It flips a qubit from |0⟩ to |1⟩ and vice versa with probability p, and leaves it be with probability 1− p. Then we have E0 = √ 1− pI = √ 1− p ( 1 0 0 1 ) E1 = √ pX = √ p ( 0 1 1 0 ) . (5.8) 5.2 Simulating Quantum Computers 27 From this one can see that the quantum operation E flips a qubit with probability p. More on quantum operations and the effects of noise can be found in textbooks like ref [3]. The effect of noisy quantum channels can be studied using quantum computer simulators. 5.2 Simulating Quantum Computers So far I have discussed quantum computers from a theoretical point of view. We want to study quantum computers from a practical one also. This is impractical, or in some cases not possible, to do using real quantum computers. As an example, the goal of my work is to study electronic structure calculations on quantum computers that correspond to classical FCI calculations. This would firstly be very resource consuming to do on actual quantum computers when algorithms are run many times. Secondly, it would not give many useful results since these algorithms are still very prone to errors. The solutions to these problems is to simulate a quantum computation. This can be done on classical computers and there already exist several different programs for simulating both error free and noisy quantum computations [28, 29]. 6 Results The goal of this work is to test quantum computer simulators by calculating the electronic structure and ground state energy of the helium atom using quantum algorithms corre- sponding to classical FCI calculations. First the one- and two electron integrals needed to be calculated and for this purpose a basis set had to be defined. I started only with the HF level without any unoccupied, or virtual orbitals. I optimized the electronic structure and included more orbitals in the next calculation, again optimizing the electronic structure, and so on, building larger basis sets of top of the previous ones. I used numerical orbitals which I will discuss shortly in the next section. After this I implemented the integrals on the quantum simulator to calculate the ground state energy of the atom. I did the cal- culations for a few different basis sets and parameter combinations, and compared the results. The quantum simulations were done on Atos quantum learning machine Kvasi, which has a Python interface. One can implement noisy quantum channels using the quantum operations formalism from section 5.1 to simulate realistic quantum computers of today, however I have opted to simulate error free computations since the goal was to understand how the algorithms work in a perfect environment. In this section I first discuss the idea 28 6 Results behind the numerical basis set I have used, and then go through the implementation of the quantum algorithm in more detail. 6.1 Numerical Orbitals As mentioned in section 2.2, the one- and two-electron weights are calculated on a clas- sical computer before the quantum simulation. To calculate these, we first need a basis set and I have made use of numerical orbitals since they have several advantages over analytical functions. A detailed discussion on the method can be found in refs [39]-[41]. The orbital functions I have used can be divided into two parts, a radial part and an angular part ϕnlm(r, ϑ, φ) = Rnl(r)Ylm(ϑ, φ), (6.1) where the angular part is an analytical function and the radial part is numerical. The radial function is defined on a grid in such a way, that inside each element r ∈ [xi, xi+1] it is defined by a polynomial. At the gridpoints r = xi continuity between the polynomials at adjacent elements is assured. One big advantage when using numerical orbitals is that they are more flexible than com- mon analytical basis functions like Slater type orbitals (STO) or Gaussian type orbitals (GTO), which means that they can be adapted to represent exact orbitals [45]. Therefore, faster convergence to the complete basis set (CBS) limit is achieved. This is important for current quantum devices especially because it will mean that fewer qubits are required. This will lead to further advantages if one uses VQE because fewer qubits means more shallow circuits. As an example, looking at the calculations in ref [13] where they used operators of the form of equation 4.25, it is clear that fewer qubits will result in more shallow circuits. 6.2 Quantum Simulations To map the electronic states to the qubit register I used DM and since I did a single point calculation, the mapping was done considering only symmetry configuration of the ground state. The JW transformation was used to map the Hamiltonian to Pauli operators and then Trotter–Suzuki decomposition was used to build the circuit. The one- and two-electron parts of the Hamiltonian are separately implemented as tensors of rank two and four. By ensuring that these tensors have the right form, one can obtain the desired symmetry configuration. This is perhaps best demonstrated with an example. Take the helium atom and a basis set of the 1- 2- and 3s-orbitals, which gives six spin 6.2 Quantum Simulations 29 orbitals. The ground state has the term symbol 1S0, i.e., zero spin projection and s- symmetry. The full Fock space can first and foremost be reduced to all states containing only two electrons {|Ψ⟩} = {|110000⟩, |101000⟩, . . . , |000011⟩}. (6.2) Further reductions can be made if we make the constraint that sz = 0 so the electrons must have opposite spin. Then the space of states becomes {|Ψ⟩} = {|110000⟩, |100100⟩, |100001⟩, |011000⟩, . . . , |000011⟩} (6.3) and therefore we only want to consider excitations which preserve the total spin of the system in the one- and two-body tensors. Then the one-body tensors have the form hij , such that the element is zero if one of the indices is odd, corresponding to spin down and one even, corresponding to spin up. This ensures that we only excite spin up to spin up and spin down to spin down. The two-electron integrals have the form ⟨x1x2|x2x1⟩, which means that the two-body tensors will have the form hijkl, where the non-zero elements are the ones where i and l are odd and j and k are even or vice versa. After the integrals were implemented it was time to choose a mapping of the operators to quantum gates and as mentioned, I opted for the JW transformation since it is the simplest and most the commonly used one. Another parameter that can be varied is the number of qubits needed for the ancilla qubit register containing the eigenvalue, and it can be obtained by analysis of performance [3]. To obtain φ accurate to n bits we need t qubits according to the following equation t = n+ ⌈ log ( 2 + 1 2ϵ )⌉ . (6.4) Choosing a success rate of 0.9, i.e, ϵ = 0.1 and n = 6, which should be sufficient accuracy, gives t = 8. The energy given when increasing the number of ancilla qubits can be found in figure 6.1. One can notice that it converges at 8 qubits. There is excellent agreement with the HF energy of helium done with QPE and the clas- sical ”exact” result. Unfortunately, I was not able to reproduce this for larger basis sets. There is still some agreement when performing an FCI computation with two s-functions using QPE, −2.887968 a.u. and the classical calculation, −2.8779968 a.u., but the al- gorithm fails for larger basis sets. I have reason to believe that the disagreement is due to limitations in the quantum simulator when simulating QPE, since diagonalizing the 30 7 Discussion Figure 6.1: Energy of helium atom given by QPE as a function of the number of ancilla qubits. Calculations done on Kvasi with HF basis set. Hamiltonian matrix gives eigenvalues which are in good agreement with the exact en- ergy. The limitations could lie in the fact that the gate count rises quite high when going to larger bases, making it hard to simulate classically. 7 Discussion I have explored progress in the field of electronic structure calculations in the age of quantum computing, discussing the important algorithms QPE and VQE and their mod- ifications. I have computed the electronic structure of the He atom on classical quantum computer simulators. These simulators are an important tool in the age of noisy quantum computers. The computations are in good agreement for minimal basis sets, but as the dimensions of the Hamiltonian matrix grows, the algorithm fails. I have used numerical orbitals, which converge faster to the CBS limit than analytical basis sets, due to their abil- ity to represent exact orbitals. This in turn could result in fewer qubits, especially when using VQE. Thus, one might in the future try to implement these numerical calculations to VQE and see if they have an impact. There of course remain many great challenges to overcome before quantum computers can make new contributions to computational chemistry. The biggest challenge remains to 31 minimize the impact of noise in quantum computations. There is progress in this area with many different approaches being taken, e.g. trapped-ion quantum computers [30, 31, 32], but this remains an engineering problem which is a whole field of its own. Moore’s law, which is the statement that the number of transistors in computers doubles about every two years, has held up quite well but stagnated a bit in recent years. Could one perhaps apply a sort of quantum version of this law, holding for qubits? It might not be at all impossible with such a broad spectrum of research going on in this field and new progress being made all the time. 32 References References [1] T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic-structure Theory, (2000), John Wiley & sons [2] A. Szabo, N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, (1989), McGraw-Hill, New York [3] M. Nielsen, I. Chuang, Quantum Computation and Quantum Information, (2010), Cambridge University Press [4] S. Kais, S. A. Rice, A. R. Dinner, Quantum Information and Computation for Chem- istry, (2014), John Wiley & sons [5] R.P. Feynman, ”Simulating physics with computers”, Int. J. Theor. Phys. 21 (1982) 467–488. [6] D. Deutsch, (1985) Proc. R. Soc. Lond. A 400, 97-117. [7] P. R. Nagy, M. Ka´llay, (2019), J. Chem. Theory Comput., 15, 10, 5275–5298 [8] S. Ficher, D. Gunlycke, (2019), arXiv:1907.01493v2 [9] D. Abrams, S. Lloyd, (1999), Phys. Rev. Lett. 83, 5162 [10] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, M. Head-Gordon, (2005), arXiv:quant- ph/0604193 [11] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, R. D. Somma, (2014), Phys. Rev. Lett. 114, 090502 [12] A. Kay, (2020), arXiv:1809.03842. [13] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, (2017), Nature 549, 242. [14] P. Jordan and E. Wigner, (1928), Z. Phys. 47, 631 [15] B. Bauer, S. B. Bravyi, M. Motta, G. Kin-Lic Chan, (2020), Chem. Rev. 120, 22, 12685–12717 [16] S. B. Bravyi, A. Yu. Kitaev, (2002), Ann. Phys. 298: 210–26. References 33 [17] A. Tranter, P. J. Love, F. Mintert, P. V. Coveney, (2018), J. Chem. Theory Comput., 14, 11, 5617–5630 [18] L. Veis, J. Pittner, (2014), J. Chem. Phys. 140, 214111 [19] G. Hao Low, I. L. Chuang, (2019), Quantum 3, 163 [20] K. Sugisaki, S. Yamamoto, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi, T. Takui, (2019), Chem. Phys. Lett. 737 100002 [21] W. Ritz, (1909), DOI:10.1515/crll.1909.135.1 [22] A. Peruzzo, J. McClean, P. Shadbolt, M. Yung, X. Zhou, P. J. Love, A. Aspuru- Guzik, J. L. O’Brien, (2014), Nat. Commun. 5, 4213 [23] P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, A. Aspuru-Guzik, (2019), Quan- tum Sci. Technol. 4 045005 [24] J. Lee, W. J. Huggins, M. Head-Gordon, K. B. J. Whaley, (2019), Chem. Theory Comput. 15, 311–324 [25] H. R. Grimsley , S. E. Economou, E. Barnes, N. J. Mayhall, (2019), Nat. Comm. 10, 3007 [26] K. Setia, J. D. Whitfield, (2018), J. Chem. Phys. 148, 164104 [27] H. L. Tang, V. O. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, S. E., Economou, (2021), PRX Quantum 2, 020310 [28] J. R. McClean, et al., (2019), arXiv:1710.07629v5 [29] G. Aleksandrowicz et al., (2019), Qiskit: An open-source framework for quantum computing [30] R. Srinivas, S. C. Burd, H. M. Knaack, R. T. Sutherland, A. Kwiatkowski, S. Glancy, E. Knill, D. J. Wineland, D. Leibfried, A. C. Wilson, D. T. C. Allcock, D. H. Slichter, (2021), Nature 597, pages 209–213 [31] J. Hilder, D. Pijn, O. Onishchenko, A. Stahl, M. Orth, B. Lekitsch, A. Rodriguez- Blanco, M. Mu¨ller, F. Schmidt-Kaler, U. G. Poschinger, (2022), Phys. Rev. X 12, 011032 34 References [32] T. P. Harty, D. Allcock, C. J. Ballance, L. Guidoni, H. Janacek, N. Linke, D. Stacey, D. M. Lucas, (2014), Phys. Rev. Lett. 113, 220501 [33] C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, D. M. Lucas, (2016), Phys. Rev. Lett. 117, 060504 [34] T. Monz, P. Schindler, J. T. Barreiro, M. Chwalla, D. Nigg, W. A. Coish, M. Harlan- der, W. Ha¨nsel, M. Hennrich, R. Blatt, (2011), Phys. Rev. Lett. 106, 130506 [35] J. M. Chow, J. M. Gambetta, A. D. Corcoles, S. T. Merkel, J. A. Smolin, C. Rigetti, S. Poletto, G. A. Keefe, M. B. Rothwell, J. R. Rozen, M. B. Ketchen, M. Steffen, (2012), Phys. Rev. Lett. 109, 060501 [36] F. Arute, et al., (2019), Nature 574, 505–510 [37] L.-K. Chen, et al., (2017), Class. Quantum Grav. 29 224011 [38] H.-S. Zhong, et al., (2018), Phys. Rev. Lett. 121, 250505. [39] D. Sundholm, J. Olsen, P. A˚. Malmqvist, B. O. Roos, ”Numerical MCSCF in One and Two Dimensions”, in Numerical Determination of the Electronic Structure of Atoms, Diatomic and Polyatomic Molecules, (1989), Proc. NATO Advanced Research Workshop, Versailles 1988, Edited by M. Defranceschi and J. Delhalle. Dordrecht: Reidel, pp. 329-334. [40] D. Sundholm, J. Olsen, (1990), Phys. Rev. A 42, 2614-2621 [41] D. Sundholm, J. Olsen, (1990), Chem. Phys. Lett. 171, 53-57 [42] R. J. Bartlett, M. Musiał, (2007), Rev. Mod. Phys. 79, 291 [43] T. D. Crawford, H. F. Schaefer, (2000), Rev. Comput. Chem. 14, 33 [44] A. G. Taube, R. J. Bartlett, (2006), Int. J. Quantum Chem. 106: 3393-3401 [45] S. Lehtola, (2019), Int J Quantum Chem. 119:e25968 [46] J. Preskill, (2018), Quantum 2, 79 [47] R. Serber, (1934) Phys. Rev. 45, 461–467. [48] R. Serber, (1934), J. Chem. Phys. 2, 697–710