{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "16e3036d-b1be-43b8-8838-95dfc14cad60",
   "metadata": {},
   "source": [
    "# Multi-reference Quantum Krylov Algorithm - $H_2$ Molecule\n",
    "\n",
    "The multireference selected quantum Krylov (MRSQK) algorithm is defined in [this paper](https://arxiv.org/pdf/1911.05163) and was developed as a low-cost alternative to quantum phase estimation. This tutorial will demonstrate how this algorithm can be implemented in CUDA-Q and accelerated using multiple GPUs. The [CUDA-Q Hadamard test tutorial](https://nvidia.github.io/cuda-quantum/latest/applications/python/hadamard_test.html) might provide helpful background information for understanding this tutorial.\n",
    "\n",
    "The algorithm works by preparing an initial state, and then defining this state in a smaller subspace constructed with a basis that corresponds to Trotter steps of the initial state. This subspace can be diagonalized to produce an approximate energy for the system without variational optimization of any parameters.\n",
    "\n",
    "In the example below, the initial guess is the ground state of the diagonalized Hamiltonian for demonstration purposes.  In practice one could use a number of heuristics to prepare the ground state such as Hartree Fock or CISD. A very promising ground state preparation method which can leverage quantum computers is the linear combination of unitaries (LCU). LCU would allow for the state preparation to occur completely on the quantum computer and avoid storing an inputting the exponentially large state vector.\n",
    "\n",
    "\n",
    "Regardless of the method used for state preparation, the procedure begins by selecting a $d$-dimensional basis of reference states ${\\Phi_0 \\cdots \\Phi_d},$ where each is a linear combination of Slater determinants: \n",
    "\n",
    "$$ \\ket{\\Phi_I}  =  \\sum_{\\mu} d_{\\mu I}\\ket{\\phi_{\\mu}}. $$\n",
    "\n",
    "\n",
    "From this, a non-orthogonal Krylov Space $\\mathcal{K} = \\{\\psi_{0} \\cdots \\psi_{N}\\}$ is constructed by applying a family of $s$ unitary operators on each of the $d$ reference states resulting in $d*s = N$ elements in the Krylov space where \n",
    "$$ \\ket{\\psi_{\\alpha}} \\equiv \\ket{\\psi_I^{(n)}} = \\hat{U}_n\\ket{\\Phi_I}, $$\n",
    "\n",
    "Therefore, the general quantum state that we originally set out to describe is\n",
    "\n",
    "$$ \\ket{\\Psi} = \\sum_{\\alpha} c_{\\alpha}\\ket{\\psi_{\\alpha}} = \\sum_{I=0}^d \\sum_{n=0}^s c_I^{(n)}\\hat{U}_n\\ket{\\Phi_I}.  $$\n",
    "\n",
    "The energy of this state can be obtained by solving the generalized eigenvalue problem\n",
    "$$ \\boldsymbol{Hc}=\\boldsymbol{Sc}E, $$\n",
    "\n",
    "where the elements of the overlap are\n",
    "\n",
    "$$S_{\\alpha \\beta} = \\braket{\\psi_{\\alpha}|\\psi_{\\beta}} =   \\braket{\\Phi_I|\\hat{U}_m^{\\dagger}\\hat{U}_n|\\Phi_J}$$ \n",
    "\n",
    "and Hamiltonian matrix is\n",
    "\n",
    "$$H_{\\alpha \\beta} = \\braket{\\psi_{\\alpha}|\\hat{H}|\\psi_{\\beta}} =   \\braket{\\Phi_I|\\hat{U}_m^{\\dagger}\\hat{H}\\hat{U}_n|\\Phi_J}.$$\n",
    "\n",
    "The matrix elements for $S$ are computed with the Hadamard test with a circuit shown below for the case of the overlap matrix elements. \n",
    "\n",
    "![Htest](./images/krylovcircuit.png)\n",
    "\n",
    "The $2\\sigma_+$ term refers to measurement of the expectation value of this circuit with the $X+iY$ operator.\n",
    "\n",
    "The Hamiltonian matrix elements are computed with a circuit that includes controlled application of the Hamiltonian. Once the $H$ and $S$ matrices are constructed, the diagonalization is performed classically to produce an estimate for the ground state in question.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b90de056-4dc5-473f-bf50-29c7d8b2ad05",
   "metadata": {},
   "source": [
    "### Setup\n",
    "\n",
    "This cell installs the necessary packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba61665c-dc3b-4e43-b1cf-340855ea68fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install pyscf\n",
    "import cudaq\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "# Single-node, single gpu\n",
    "if cudaq.num_available_gpus() > 0 and cudaq.has_target(\"nvidia\"):\n",
    "    cudaq.set_target(\"nvidia\")\n",
    "else:\n",
    "    print(\"CUDA or GPU support is unavailable. Running with CPU simulator. Performance may be significantly reduced.\")\n",
    "    cudaq.set_target(\"qpp-cpu\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba6e94e-a3e6-48aa-8392-e598429b0be7",
   "metadata": {},
   "source": [
    "The molecule is defined below and its Hamiltonian is extracted as a matrix. The matrix is diagonalized to produce the ground state. The corresponding state vector will be used as the initial state to ensure good results for this demonstration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2fb4192c-054a-4bff-b86f-a8d58cf8bac6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwrite output file: H 0-pyscf.log\n",
      "[pyscf] Total number of orbitals =  2\n",
      "[pyscf] Total number of electrons =  2\n",
      "[pyscf] HF energy =  -1.116325564486115\n",
      "[pyscf] Total R-CCSD energy =  -1.1371758844013327\n",
      "Ground state energy (classical simulation)=  (-1.1371757102406845+0j) , index=  3\n"
     ]
    }
   ],
   "source": [
    "from qchem.classical_pyscf import get_mol_hamiltonian\n",
    "from qchem.hamiltonian import jordan_wigner_fermion\n",
    "\n",
    "# Define H2 molecule\n",
    "\n",
    "# Generate the spin Hamiltonian\n",
    "\n",
    "geometry = 'H 0.0 0.0 0.0; H 0.0 0.0 0.7474'\n",
    "#geometry = 'Li 0.3925 0.0 0.0; H -1.1774 0.0 0.0'\n",
    "#geometry = \"qchem/H2O.xyz\"\n",
    "\n",
    "\n",
    "# Run HF, ccsd and compute the spin molecular hamiltonian using the HF molecular orbitals.\n",
    "molecular_data = get_mol_hamiltonian(xyz=geometry, spin=0, charge=0, basis='sto3g', ccsd=True, verbose=True)\n",
    "\n",
    "# For active space calculations, we can specify the number of active electrons and orbitals.\n",
    "#molecular_data = get_mol_hamiltonian(xyz=geometry, spin=0, charge=0, basis='631g', nele_cas=4, norb_cas=4, \n",
    "#                                     ccsd=True, casci=True, verbose=True)\n",
    "\n",
    "obi = molecular_data[0]\n",
    "tbi = molecular_data[1]\n",
    "e_nn = molecular_data[2]\n",
    "nelectrons = molecular_data[3]\n",
    "norbitals = molecular_data[4]\n",
    "\n",
    "qubits_num = 2 * norbitals\n",
    "\n",
    "# Jordan-Wigner transformation to convert the fermionic Hamiltonian to a spin Hamiltonian\n",
    "hamiltonian = jordan_wigner_fermion(obi, tbi, e_nn, tolerance = 1e-12)\n",
    "\n",
    "# Diagonalize Hamiltonian\n",
    "spin_ham_matrix = hamiltonian.to_matrix()\n",
    "e, c = np.linalg.eig(spin_ham_matrix)\n",
    "\n",
    "# Find the ground state energy and the corresponding eigenvector\n",
    "print('Ground state energy (classical simulation)= ', np.min(e), ', index= ',\n",
    "      np.argmin(e))\n",
    "min_indices = np.argmin(e)\n",
    "\n",
    "# Eigenvector can be used to initialize the qubits\n",
    "vec = c[:, min_indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "562f795e-cf4c-4e43-a0d9-3a3ae346277d",
   "metadata": {},
   "source": [
    "The functions below take the original Hamiltonian defined above, and strip it into a list of Pauli words and a list of coefficients which will be uses later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "014e72eb-ee66-48f5-bb58-73210e621fad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.17028010135220506+0j), (0.17028010135220503+0j), (-0.2200413002242175+0j), (-0.2200413002242175+0j), (0.1683359862516207+0j), (0.12020049071260122+0j), (0.1656068235817425+0j), (0.1656068235817425+0j), (0.12020049071260122+0j), (0.17407289249680213+0j), (-0.04540633286914128+0j), (0.04540633286914128+0j), (0.04540633286914128+0j), (-0.04540633286914128+0j)]\n",
      "['ZIII', 'IZII', 'IIZI', 'IIIZ', 'ZZII', 'ZIZI', 'ZIIZ', 'IZZI', 'IZIZ', 'IIZZ', 'XXYY', 'XYYX', 'YXXY', 'YYXX']\n"
     ]
    }
   ],
   "source": [
    "# Collect coefficients from a spin operator so we can pass them to a kernel.\n",
    "# The identity term is excluded. Its contribution is added back to the \n",
    "# Hamiltonian matrix classically below.\n",
    "def term_coefficients(ham: cudaq.SpinOperator) -> list[complex]:\n",
    "    result = []\n",
    "    for term in ham:\n",
    "        if term.is_identity():\n",
    "            continue\n",
    "        result.append(term.evaluate_coefficient())\n",
    "    return result\n",
    "\n",
    "# Collect Pauli words from a spin operator so we can pass them to a kernel\n",
    "def term_words(ham: cudaq.SpinOperator) -> list[str]:\n",
    "    # Our kernel uses these words to apply exp_pauli to the entire state.\n",
    "    # we hence ensure that each pauli word covers the entire space.\n",
    "    result = []\n",
    "    for term in ham:\n",
    "        if term.is_identity():\n",
    "            continue\n",
    "        result.append(term.get_pauli_word(qubits_num))\n",
    "    return result\n",
    "\n",
    "# Build the lists of coefficients and Pauli Words from the H2 Hamiltonian\n",
    "coefficient = term_coefficients(hamiltonian)\n",
    "pauli_string = term_words(hamiltonian)\n",
    "\n",
    "# Sum of identity-term coefficients\n",
    "# The identity contributes `identity_coef * S` to the Hamiltonian matrix.\n",
    "identity_coef = sum(\n",
    "    term.evaluate_coefficient().real\n",
    "    for term in hamiltonian\n",
    "    if term.is_identity())\n",
    "\n",
    "print(coefficient)\n",
    "print(pauli_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74b4c079-8837-47ae-a484-52ec5f4a160e",
   "metadata": {},
   "source": [
    "In this example, the unitary operators that build the Krylov subspace are first-order Trotter operations at different time steps. The performance here could potentially be improved by increasing the size of the time step, using a higher order Trotter approximation, or using other sorts of approximations. The CUDA-Q kernels below define the unitary operations that construct the $\\psi$ basis.  Each receives the target qubits, the time step, and components of the Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c98db2b7-eaaf-4e3a-83bb-4f7ad0730471",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Applies Unitary operation corresponding to Um\n",
    "@cudaq.kernel\n",
    "def U_m(qubits: cudaq.qview, dt: float, coefficients: list[complex],\n",
    "        words: list[cudaq.pauli_word]):\n",
    "    # Compute U_m = exp(-i m dt H)\n",
    "    for i in range(len(coefficients)):\n",
    "        exp_pauli(dt * coefficients[i].real, qubits, words[i])\n",
    "\n",
    "\n",
    "# Applies Unitary operation corresponding to Un\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_n(qubits: cudaq.qview, dt: float, coefficients: list[complex],\n",
    "        words: list[cudaq.pauli_word]):\n",
    "    # Compute U_n = exp(-i n dt H)\n",
    "    for i in range(len(coefficients)):\n",
    "        exp_pauli(dt * coefficients[i].real, qubits, words[i])\n",
    "\n",
    "\n",
    "# Applies the gate operations for a Hamiltonian Pauli Word\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def apply_pauli(qubits: cudaq.qview, word: list[int]):\n",
    "    # Add H (Hamiltonian operator)\n",
    "    for i in range(len(word)):\n",
    "        if word[i] == 1:\n",
    "            x(qubits[i])\n",
    "        if word[i] == 2:\n",
    "            y(qubits[i])\n",
    "        if word[i] == 3:\n",
    "            z(qubits[i])\n",
    "\n",
    "\n",
    "# Performs Hadamard test circuit which determines matrix elements of S and H of subspace\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def qfd_kernel(dt_alpha: float, dt_beta: float, coefficients: list[complex],\n",
    "               words: list[cudaq.pauli_word], word_list: list[int],\n",
    "               vec: list[complex]):\n",
    "\n",
    "    ancilla = cudaq.qubit()\n",
    "    qreg = cudaq.qvector(vec)\n",
    "\n",
    "    h(ancilla)\n",
    "\n",
    "    x(ancilla)\n",
    "    cudaq.control(U_m, ancilla, qreg, dt_alpha, coefficients, words)\n",
    "    x(ancilla)\n",
    "\n",
    "    cudaq.control(apply_pauli, ancilla, qreg, word_list)\n",
    "    cudaq.control(U_n, ancilla, qreg, dt_beta, coefficients, words)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c8f9fb4-4974-44bf-aa8c-c9e86ffb0dd4",
   "metadata": {},
   "source": [
    "The auxillary function below takes a Pauli word, and converts it to a list of integers which informs applications of this word to a circuit with the `apply_pauli` kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5bed072a-9fb6-4394-9759-b45ba292d0a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_str(pauli_word, qubits_num):\n",
    "\n",
    "    my_list = []\n",
    "    for i in range(qubits_num):\n",
    "        if str(pauli_word[i]) == 'I':\n",
    "            my_list.append(0)\n",
    "        if str(pauli_word[i]) == 'X':\n",
    "            my_list.append(1)\n",
    "        if str(pauli_word[i]) == 'Y':\n",
    "            my_list.append(2)\n",
    "        if str(pauli_word[i]) == 'Z':\n",
    "            my_list.append(3)\n",
    "    return my_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "999675e2-a5b0-43d1-ba80-827f1c81b3dc",
   "metadata": {},
   "source": [
    "The spin operators necessary for the Hadamard test are defined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "297dbc29-e48a-411b-bc42-4fabe180a641",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the spin-op x for real component and y for the imaginary component.\n",
    "x_0 = cudaq.spin.x(0)\n",
    "y_0 = cudaq.spin.y(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1880e50d-6d27-44e4-b27e-329d8c499bfb",
   "metadata": {},
   "source": [
    "Finally, the time step for the unitary operations that build the Krylov space is defined as well as the dimension of the Krylov space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3ae6b210-13c4-40ba-be73-4fb906a1f7ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define parameters for the quantum Krylov space\n",
    "dt = 0.5\n",
    "\n",
    "# Dimension of the Krylov space\n",
    "m_qfd = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d59f13b-e9f8-44e6-b270-be4fc6a80c01",
   "metadata": {},
   "source": [
    "### Computing the matrix elements \n",
    "\n",
    "The cell below computes the overlap matrix. This can be done in serial or in parallel, depending on the `multi_gpu` specification.  First, an operator is built to apply the identity to the overlap matrix circuit when `apply_pauli` is called. Next, the `wf_overlap` array is constructed which will hold the matrix elements. \n",
    "\n",
    "Next, a pair of nested loops iterate over the time steps defined by the dimension of the subspace.  Each m,n combination corresponds to computation of an off-diagonal matrix element of the overlap matrix $S$ using the Hadamard test.  This is accomplished by calling the CUDA-Q `observe` function with the X and Y operators, along with the time steps, the components of the Hamiltonian matrix, and the initial state vector `vec`.\n",
    "\n",
    "The observe function broadcasts over the two provided operators $X$ and $Y$ and returns a list of results. The `expectation` function returns the expectation values which are summed and stored in the matrix.\n",
    "\n",
    "The multi-gpu case completes the same steps, except the `observe_async` command is used. This allows for the $X$ and $Y$ observables to be evaluated at the same time on two different simulated QPUs. In this case, the results are stored in lists corresponding to the real and imaginary parts. These are then accessed later with the `get` command to build $S$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "875124e8-dba8-4ace-828b-d4b03352f9d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the identity operator\n",
    "identity_op = cudaq.SpinOperator.from_word('I' * qubits_num)\n",
    "# Get the Pauli word and convert it to a list of integers\n",
    "# where I=0, X=1, Y=2, Z=3\n",
    "# This is used to apply the identity operator in the kernel.\n",
    "identity_word = identity_op.get_pauli_word()\n",
    "pauli_list = pauli_str(identity_word, qubits_num)\n",
    "\n",
    "# Empty overlap matrix S\n",
    "wf_overlap = np.zeros((m_qfd, m_qfd), dtype=complex)\n",
    "# Loop to solve for S matrix elements\n",
    "for m in range(m_qfd):\n",
    "    dt_m = dt * m\n",
    "    for n in range(m, m_qfd):\n",
    "        dt_n = dt * n\n",
    "        results = cudaq.observe(qfd_kernel, [x_0, y_0], dt_m, dt_n,\n",
    "                                coefficient, pauli_string, pauli_list, vec)\n",
    "        temp = [result.expectation() for result in results]\n",
    "        wf_overlap[m, n] = temp[0] + temp[1] * 1j\n",
    "        if n != m:\n",
    "            wf_overlap[n, m] = np.conj(wf_overlap[m, n])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fa0e0f1-0705-4dd6-a354-30db1526e266",
   "metadata": {},
   "source": [
    "The Hamiltonian matrix elements are computed in the same way, except this time the Hamiltonian is applied as part of the circuit. This is accomplished with the extra for loop, which iterates through the terms in the Hamiltonian, computing an expectation value for each one and then summing the results to produce one matrix element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "562a6063-b0fe-4d69-a4f4-c7d5309436da",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Single GPU:\n",
    "\n",
    "# Empty H matrix\n",
    "# This matrix will be used to store the Hamiltonian matrix elements\n",
    "# for the Krylov subspace.\n",
    "# The matrix is symmetric, so we only need to compute the upper triangle.\n",
    "# The matrix is of size m_qfd x m_qfd, where m_qfd is the dimension of the Krylov subspace.\n",
    "ham_matrx = np.zeros((m_qfd, m_qfd), dtype=complex)\n",
    "\n",
    "# Loops over H matrix terms\n",
    "for m in range(m_qfd):\n",
    "    dt_m = dt * m\n",
    "    for n in range(m, m_qfd):\n",
    "        dt_n = dt * n\n",
    "\n",
    "        # 2 entry array that stores real and imaginary part of matrix element\n",
    "        tot_e = np.zeros(2)\n",
    "\n",
    "        # Loops over the (non-identity) terms in the Hamiltonian, computing expectation values\n",
    "        for coef, word in zip(coefficient, pauli_string):\n",
    "            pauli_list = pauli_str(word, qubits_num)\n",
    "            \n",
    "            results = cudaq.observe(qfd_kernel, [x_0, y_0], dt_m, dt_n,\n",
    "                                        coefficient, pauli_string, pauli_list,\n",
    "                                        vec)\n",
    "\n",
    "            temp = [result.expectation() for result in results]\n",
    "\n",
    "            # Multiplies result by coefficient corresponding to Pauli Word\n",
    "            temp[0] = coef.real * temp[0]\n",
    "            temp[1] = coef.real * temp[1]\n",
    "\n",
    "            # Accumulates results for each Pauli Word\n",
    "            tot_e[0] += temp[0]\n",
    "            tot_e[1] += temp[1]\n",
    "\n",
    "        # Adds back the identity-term contribution.\n",
    "        ham_matrx[m, n] = tot_e[0] + tot_e[1] * 1j + identity_coef * wf_overlap[m, n]\n",
    "        if n != m:\n",
    "            ham_matrx[n, m] = np.conj(ham_matrx[m, n])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1917bd58-4ce3-4e3e-9ccf-096577c0da14",
   "metadata": {},
   "source": [
    "### Determining the ground state energy of the subspace\n",
    "\n",
    "The final step is to solve the generalized eigenvaulue problem with the overlap and Hamiltonian matrices constructed using the quantum computer. The procedure begins by diagonalizing  $S$ with the transform $$S = U\\Sigma U^{\\dagger}$$\n",
    "\n",
    "The eigenvectors $v$ and eigenvalues $s$ are used to construct a new matrix $X':$\n",
    "\n",
    "$$ X' = S ^{\\frac{-1}{2}} = \\sum_k v_{ki} \\frac{1}{\\sqrt{s_k}} v_{kj}.$$\n",
    "\n",
    "The matrix $X'$ diagonalizes $H:$\n",
    "\n",
    "$$ X'^{\\dagger}HX' = ES^{\\frac{1}{2}}C.$$\n",
    "\n",
    "Using the eigenvectors of $H'$, ($^{\\frac{1}{2}}C$), the original eigenvectors to the problem can be found by left multiplying by $S^{\\frac{-1}{2}}C.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8b281533-bf3d-4aed-823b-9d2f926dfe3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def eig(H, s):\n",
    "    #Solver for generalized eigenvalue problem\n",
    "    # HC = SCE\n",
    "\n",
    "    THRESHOLD = 1e-20\n",
    "    s_diag, u = scipy.linalg.eig(s)\n",
    "    s_prime = []\n",
    "    for sii in s_diag:\n",
    "        if np.imag(sii) > 1e-7:\n",
    "            raise ValueError(\n",
    "                \"S may not be hermitian, large imag. eval component.\")\n",
    "        if np.real(sii) > THRESHOLD:\n",
    "            s_prime.append(np.real(sii))\n",
    "\n",
    "    X_prime = np.zeros((len(s_diag), len(s_prime)), dtype=complex)\n",
    "\n",
    "    for i in range(len(s_diag)):\n",
    "        for j in range(len(s_prime)):\n",
    "            X_prime[i][j] = u[i][j] / np.sqrt(s_prime[j])\n",
    "\n",
    "    H_prime = (((X_prime.conjugate()).transpose()).dot(H)).dot(X_prime)\n",
    "    e_prime, C_prime = scipy.linalg.eig(H_prime)\n",
    "    C = X_prime.dot(C_prime)\n",
    "\n",
    "    return e_prime, C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "74c1fa51-4fae-412b-b90d-262e0b3aaf53",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy from QFD:\n",
      "(-1.1359686811350462-4.497484607599205e-09j)\n"
     ]
    }
   ],
   "source": [
    "eigen_value, eigen_vect = eig(ham_matrx[0:m_qfd, 0:m_qfd], wf_overlap[0:m_qfd,\n",
    "                                                                      0:m_qfd])\n",
    "print('Energy from QFD:')\n",
    "print(np.min(eigen_value))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
