{ "cells": [ { "cell_type": "markdown", "id": "ce0e891d-f149-4005-9dc3-524bd9c8c047", "metadata": {}, "source": [ "# Variational Basis State Encoder (Ground State)" ] }, { "cell_type": "markdown", "id": "9ceab38a", "metadata": {}, "source": [ "## 1 Backgroud\n", "This tutorial is for variational basis state encoder (VBE) for the ground state of the Holstein model.\n", "\n", "$$ \\hat{H} = -\\sum_{i,j} V \\hat{a}_i^\\dagger \\hat{a}_j + \\sum_i \\omega \\hat{b}_i^\\dagger \\hat{b}_i + \\sum_i g\\omega \\hat{a} _i^\\dagger \\hat{a}_i (\\hat{b}_i^\\dagger + \\hat{b}_i) $$\n", "\n", "To calculate the ground state of Holstein model accurately, many levels of phonons are needed, which will cost too many qubits in quantum circuit. There exists another idea that we can view linear conbination of phonons as an effective phonon mode, which is possible to save the qubits in phonon encoding:\n", "\n", "$$\\hat{B}[l]=\\sum_{n=1}^{2^{N_l}} \\ket{n}_l\\sum_m C[l]_{mn} \\bra{m}_l$$\n", "\n", "Here, we transform the original phonon basis $\\ket{m}$ to the encoded basis $\\ket{n}$, via transformation operator $\\hat{B}[l]$. The form of $\\hat{B}[l]$ is the central of VBE and the algorithm are presented in section 2. For more details, see https://doi.org/10.1103/PhysRevResearch.5.023046" ] }, { "cell_type": "markdown", "id": "420ea9d3", "metadata": {}, "source": [ "## 2 Algorithm Realization" ] }, { "cell_type": "markdown", "id": "af57906d", "metadata": {}, "source": [ "### 2.1 Imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "bd24d246", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy\n", "from opt_einsum import contract\n", "import tensorcircuit as tc\n", "\n", "from tencirchem import set_backend, Op, BasisSHO, BasisSimpleElectron, Mpo, Model\n", "from tencirchem.dynamic import get_ansatz, qubit_encode_op, qubit_encode_basis\n", "from tencirchem.utils import scipy_opt_wrap\n", "from tencirchem.applications.vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args" ] }, { "cell_type": "markdown", "id": "1b29e1ed", "metadata": {}, "source": [ "## 2.2 Initial setups\n", "In this section, we set intital parameters for coming sections. Here, `JAX` is used as backend. `nsite`, `omega`, `v` correspond to the site number, phonon frequency $\\omega$, transfer intergral $V$ in Holstein model, respectively:\n", "\n", "$$ \\hat{H} = -\\sum_{i,j} V \\hat{a}_i^\\dagger \\hat{a}_j + \\sum_i \\omega \\hat{b}_i^\\dagger \\hat{b}_i + \\sum_i g\\omega \\hat{a} _i^\\dagger \\hat{a}_i (\\hat{b}_i^\\dagger + \\hat{b}_i) $$\n", "\n", "Each site possesses one phonon mode, which is represented by 2 qubit per phonon (see `n_qubit_per_mode`). Considering gray encoding is adopted, the number of phonon basis (`nbas_v`) is $2^2$. `psi_index_top` and `psi_index_bottom` correspond to the physical index of ket and bra, `b_dof_vidx` correspond to the qubits that need VBE, `psi_shape2` is the physical bond dimension of each qubit state.\n", "The stucture of wavefunction and operator are presented in Fig.1. Note that the related arguments and functions are also marked.\n", "![fig1](../statics/vbe_gs_Fig1.svg)\n", "Fig. 1 The structure of wavefunction and operator. Blue squares correspond to qubit representing spin, green circles correspond to qubits representing vibrations, purple circles correspond to $B[l]$, and orange squares correspond to Matrix Product Operators (MPO)." ] }, { "cell_type": "code", "execution_count": 2, "id": "37646a1e", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2023-08-08 15:44:36.794994: E external/org_tensorflow/tensorflow/compiler/xla/stream_executor/cuda/cuda_driver.cc:267] failed call to cuInit: CUDA_ERROR_NO_DEVICE: no CUDA-capable device is detected\n", "WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "psi_index_top: ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'p-2-bottom', 'v-3-0-bottom', 'v-3-1-bottom', 'p-4-bottom', 'v-5-0-bottom', 'v-5-1-bottom'] \n", " psi_index_bottom: ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'p-2-bottom', 'v-3-0-bottom', 'v-3-1-bottom', 'p-4-bottom', 'v-5-0-bottom', 'v-5-1-bottom'] \n", " b_dof_vidx: [array([1, 2]), array([4, 5]), array([7, 8])] \n", " psi_shape2: [2, 2, 2, 2, 2, 2, 2, 2, 2]\n" ] } ], "source": [ "backend = set_backend(\"jax\")\n", "\n", "nsite = 3\n", "omega = 1\n", "v = 1\n", "# two qubit for each mode\n", "# modify param_ids before modifying this\n", "n_qubit_per_mode = 2\n", "nbas_v = 1 << n_qubit_per_mode\n", "\n", "# -1 for electron dof, natural numbers for phonon dof\n", "dof_nature = np.array([-1, 0, 0, -1, 1, 1, -1, 2, 2])\n", "# physical index for phonon mode\n", "b_dof_pidx = np.array([1, 3, 5])\n", "\n", "psi_idx_top, psi_idx_bottom, b_dof_vidx = get_psi_indices(dof_nature, b_dof_pidx, n_qubit_per_mode)\n", "\n", "n_dof = len(dof_nature)\n", "psi_shape2 = [2] * n_dof\n", "\n", "print(\n", " \"psi_index_top: \",\n", " psi_idx_bottom,\n", " \"\\n psi_index_bottom: \",\n", " psi_idx_bottom,\n", " \"\\n b_dof_vidx: \",\n", " b_dof_vidx,\n", " \"\\n psi_shape2: \",\n", " psi_shape2,\n", ")\n", "\n", "c = tc.Circuit(nsite * 3) # generate quantum circuit\n", "c.X(0) # prepare one-electron initial state\n", "n_layers = 3 # layers of ansatz" ] }, { "cell_type": "markdown", "id": "9df7f0b5", "metadata": {}, "source": [ "## 2.3 Get Variational Hamiltonian Ansatz (VHA) Terms\n", "In this section, we will generate variational hamiltonian ansatz terms. The following ansatz is adopted:\n", "\n", "$$\\ket{\\phi}=\\prod_l^L {\\prod_{} e^{\\theta_{ijk} {(\\hat{a}_j^\\dagger \\hat{a}_k - \\hat{a}_k^\\dagger \\hat{a}_j)}} \n", "\\prod_j {e^{\\theta_{lj}\\hat{a}_j^\\dagger \\hat{a}_j (\\hat{b}_j^\\dagger - \\hat{b}_j)}}} \\ket{\\phi_0} $$\n", "\n", "The anasatz is transformed from electron-phonon basis to qubit basis through `qubit_encode_op()` and `qubit_encode_basis()`" ] }, { "cell_type": "code", "execution_count": 3, "id": "aa59a6bf", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ansatz_terms: \n", " [Op('X Y', [0, 1], 0.5j), Op('Y X', [0, 1], -0.5j), Op('Y', [((0, 0), 'TCCQUBIT-1')], 0.1830127018922193j), Op('Z Y', [((0, 0), 'TCCQUBIT-0'), ((0, 0), 'TCCQUBIT-1')], -0.6830127018922193j), Op('Y', [((0, 0), 'TCCQUBIT-0')], -0.3535533905932738j), Op('Y Z', [((0, 0), 'TCCQUBIT-0'), ((0, 0), 'TCCQUBIT-1')], 0.3535533905932738j), Op('Z Y', [0, ((0, 0), 'TCCQUBIT-1')], -0.1830127018922193j), Op('Z Z Y', [0, ((0, 0), 'TCCQUBIT-0'), ((0, 0), 'TCCQUBIT-1')], 0.6830127018922193j), Op('Z Y', [0, ((0, 0), 'TCCQUBIT-0')], 0.3535533905932738j), Op('Z Y Z', [0, ((0, 0), 'TCCQUBIT-0'), ((0, 0), 'TCCQUBIT-1')], -0.3535533905932738j), Op('X Y', [1, 2], 0.5j), Op('Y X', [1, 2], -0.5j), Op('Y', [((1, 0), 'TCCQUBIT-1')], 0.1830127018922193j), Op('Z Y', [((1, 0), 'TCCQUBIT-0'), ((1, 0), 'TCCQUBIT-1')], -0.6830127018922193j), Op('Y', [((1, 0), 'TCCQUBIT-0')], -0.3535533905932738j), Op('Y Z', [((1, 0), 'TCCQUBIT-0'), ((1, 0), 'TCCQUBIT-1')], 0.3535533905932738j), Op('Z Y', [1, ((1, 0), 'TCCQUBIT-1')], -0.1830127018922193j), Op('Z Z Y', [1, ((1, 0), 'TCCQUBIT-0'), ((1, 0), 'TCCQUBIT-1')], 0.6830127018922193j), Op('Z Y', [1, ((1, 0), 'TCCQUBIT-0')], 0.3535533905932738j), Op('Z Y Z', [1, ((1, 0), 'TCCQUBIT-0'), ((1, 0), 'TCCQUBIT-1')], -0.3535533905932738j), Op('X Y', [0, 2], -0.5j), Op('Y X', [0, 2], 0.5j), Op('Y', [((2, 0), 'TCCQUBIT-1')], 0.1830127018922193j), Op('Z Y', [((2, 0), 'TCCQUBIT-0'), ((2, 0), 'TCCQUBIT-1')], -0.6830127018922193j), Op('Y', [((2, 0), 'TCCQUBIT-0')], -0.3535533905932738j), Op('Y Z', [((2, 0), 'TCCQUBIT-0'), ((2, 0), 'TCCQUBIT-1')], 0.3535533905932738j), Op('Z Y', [2, ((2, 0), 'TCCQUBIT-1')], -0.1830127018922193j), Op('Z Z Y', [2, ((2, 0), 'TCCQUBIT-0'), ((2, 0), 'TCCQUBIT-1')], 0.6830127018922193j), Op('Z Y', [2, ((2, 0), 'TCCQUBIT-0')], 0.3535533905932738j), Op('Z Y Z', [2, ((2, 0), 'TCCQUBIT-0'), ((2, 0), 'TCCQUBIT-1')], -0.3535533905932738j)] \n", "spin_basis: \n", " [BasisHalfSpin(dof: 0, nbas: 2), BasisHalfSpin(dof: ((0, 0), 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ((0, 0), 'TCCQUBIT-1'), nbas: 2), BasisHalfSpin(dof: 1, nbas: 2), BasisHalfSpin(dof: ((1, 0), 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ((1, 0), 'TCCQUBIT-1'), nbas: 2), BasisHalfSpin(dof: 2, nbas: 2), BasisHalfSpin(dof: ((2, 0), 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ((2, 0), 'TCCQUBIT-1'), nbas: 2)] \n", "param_ids: \n", " [1, -1, 0, 2, 3, 4, 5, 6, 7, 8, 9, -9, 10, 11, 12, 13, 14, 15, 16, 17, 18, -18, 19, 20, 21, 22, 23, 24, 25, 26] \n", "ansatz: \n", " .ansatz at 0x7f2fd4ed29e0>\n" ] } ], "source": [ "def get_vha_terms():\n", " # variational Hamiltonian ansatz (vha) terms\n", "\n", " g = 1 # dummy value, doesn't matter\n", " ansatz_terms = []\n", " for i in range(nsite):\n", " j = (i + 1) % nsite\n", " ansatz_terms.append(Op(r\"a^\\dagger a\", [i, j], v))\n", " ansatz_terms.append(Op(r\"a^\\dagger a\", [j, i], -v))\n", " ansatz_terms.append(Op(r\"a^\\dagger a b^\\dagger-b\", [i, i, (i, 0)], g * omega))\n", "\n", " basis = []\n", " for i in range(nsite):\n", " basis.append(BasisSimpleElectron(i))\n", " basis.append(BasisSHO((i, 0), omega, nbas_v))\n", "\n", " ansatz_terms, _ = qubit_encode_op(ansatz_terms, basis, boson_encoding=\"gray\")\n", " spin_basis = qubit_encode_basis(basis, boson_encoding=\"gray\")\n", " # this is currently hard-coded for `n_qubit_per_mode==2`\n", " # if the values of param_ids are opposite to each other, the values of the parameters are forced to be opposite in the optimization.\n", " param_ids = [1, -1, 0, 2, 3, 4, 5, 6, 7, 8] + [9, -9] + list(range(10, 18)) + [18, -18] + list(range(19, 27))\n", " return ansatz_terms, spin_basis, param_ids\n", "\n", "\n", "ansatz_terms, spin_basis, param_ids = get_vha_terms()\n", "ansatz = get_ansatz(ansatz_terms, spin_basis, n_layers, c, param_ids)\n", "\n", "print(\n", " \"ansatz_terms: \\n\",\n", " ansatz_terms,\n", " \"\\nspin_basis: \\n\",\n", " spin_basis,\n", " \"\\nparam_ids: \\n\",\n", " param_ids,\n", " \"\\nansatz: \\n\",\n", " ansatz,\n", ")" ] }, { "cell_type": "markdown", "id": "7ecd5d1a", "metadata": {}, "source": [ "## 2.4 Cost Functions for VQE Part\n", "The VQE parameters $\\theta_k$ are optimized via following equation:\n", "\n", "$$\\partial {\\bra{\\phi}\\hat{H}_1\\ket{\\phi}}/\\partial {\\theta_k} = 0$$\n", "\n", "where\n", "\n", "$$\\hat{H}_1=\\prod_{l}\\hat{B}[l] \\hat{H} \\prod_{l}\\hat{B}[l]^\\dagger$$" ] }, { "cell_type": "code", "execution_count": 4, "id": "16f66ff0", "metadata": {}, "outputs": [], "source": [ "def cost_fn(params, h):\n", " state = ansatz(params)\n", " return (state.conj() @ (h @ state)).squeeze().real\n", "\n", "\n", "vg = backend.jit(backend.value_and_grad(cost_fn))\n", "opt_fn = scipy_opt_wrap(vg)" ] }, { "cell_type": "markdown", "id": "80df3090", "metadata": {}, "source": [ "## 2.5 Get Hamiltonian Terms and Basis\n", "In this section, we generate the operator of the Holstein Hamiltonian presented in Section 2.2. The format of the operator are shown in Fig. 1. Note that the number of phonon levels are controlled by `nbas`." ] }, { "cell_type": "code", "execution_count": 5, "id": "93d727ad", "metadata": {}, "outputs": [], "source": [ "def get_ham_terms_and_basis(g, nbas):\n", " terms = []\n", " for i in range(nsite):\n", " terms.append(Op(r\"b^\\dagger b\", (i, 0), omega))\n", " terms.append(Op(r\"a^\\dagger a b^\\dagger+b\", [i, i, (i, 0)], g * omega))\n", " j = (i + 1) % nsite\n", " terms.append(Op(r\"a^\\dagger a\", [i, j], -v))\n", " terms.append(Op(r\"a^\\dagger a\", [j, i], -v))\n", "\n", " basis = []\n", " for i in range(nsite):\n", " basis.append(BasisSimpleElectron(i))\n", " basis.append(BasisSHO((i, 0), omega, nbas))\n", "\n", " return terms, basis" ] }, { "cell_type": "markdown", "id": "2d065fcd", "metadata": {}, "source": [ "## 2.6 Update $B[l]$ in Iteration\n", "In this section, the function that calculates $B[l]$ are defined:\n", "\n", "$$(1-\\hat{P}[l])\\bra{\\phi}\\tilde{H'}\\ket{\\phi}=0$$\n", "\n", "where\n", "\n", "$$\\hat{P}[l]=\\hat{B}[l]^\\dagger\\hat{B}[l]$$\n", "\n", "and\n", "\n", "$$\\tilde{H'}=\\hat{H}_{contracted}\\hat{B}[l]^\\dagger=\\left( \\prod_{k\\neq l}\\hat{B}[l] \\hat{H} \\prod_{k\\neq l}\\hat{B}[k]^\\dagger \\right) \\hat{B}[l]^\\dagger$$\n", "\n", "The graphic representarion of `h_contracted` is presented in Fig. 2. Obiviously, if $\\hat{H}$ are provided, we can obtain $\\hat{B}[l]$ by solving the equation mentioned above. Considering this is a non-linear equation, several initial guesses are needed to avoid local minimum, which is controlled by `nroot`.\n", "![fig2](../statics/vbe_gs_Fig2.svg)\n", "Fig. 2 Graphic representation of `h_contracted`" ] }, { "cell_type": "code", "execution_count": 6, "id": "8ee8576c", "metadata": {}, "outputs": [], "source": [ "def solve_b_array(psi, h_mpo, b_array, i):\n", " nbas = b_array.shape[-1]\n", " # get the input of tensor contraction function `contract`\n", " args = get_contract_args(psi, h_mpo, b_array, i, n_qubit_per_mode, psi_idx_top, psi_idx_bottom, b_dof_pidx)\n", " k = b_dof_pidx[i]\n", " # output indices\n", " args.append(\n", " [\n", " f\"v-{k}-0-bottom\",\n", " f\"v-{k}-1-bottom\",\n", " f\"p-{k}-bottom\",\n", " f\"v-{k}-0-top\",\n", " f\"v-{k}-1-top\",\n", " f\"p-{k}-top\",\n", " \"mpo-0\",\n", " f\"mpo-{len(h_mpo)}\",\n", " ]\n", " )\n", " # get contracted_h and reshape the dofs named v-{k}-0-bottom(top) and v-{k}-1-bottom(top) to one dof with dimension 4\n", " contracted_h = contract(*args).reshape(4, nbas, 4, nbas)\n", " nroot = 3\n", "\n", " def f(x):\n", " x = x.reshape(nroot, 4, nbas)\n", " # calculate P[l]\n", " p = contract(\"abc, abd -> acd\", x.conj(), x)\n", " return contract(\"abcd, kab, kde -> kce\", contracted_h, x, (np.array([np.eye(nbas)] * nroot) - p)).ravel()\n", "\n", " # solve the equation mentioned above to obtain B[l]\n", " sols = scipy.optimize.root(f, [b_array[i].flatten()] * nroot, method=\"df-sane\").x.reshape(3, 4, nbas)\n", "\n", " sols = list(sols) + [b_array[i].copy()]\n", " b_array = b_array.copy()\n", " es = []\n", " for k, new_b in enumerate(sols):\n", " # ensure the orthomormal constraint of B[l]\n", " if not np.allclose(new_b @ new_b.T, np.eye(4)):\n", " # print(f\"Enforcing orthogonality for the {k}th root of b_array[{i}]\")\n", " new_b = np.linalg.qr(new_b.T)[0].T\n", " b_array[i] = new_b\n", " e = psi @ get_contracted_mpo(h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom) @ psi\n", " es.append(e)\n", " # print(np.array(es))\n", " lowest_id = np.argmin(es)\n", " return sols[lowest_id]" ] }, { "cell_type": "markdown", "id": "859c6642", "metadata": {}, "source": [ "## 2.7 Main Structure of the Function\n", "This section is the main part of the funtion. The codes contain following parts:\n", "(i) Initialize the parameters and functions, some initializations are also preformed in section 2.2\n", "(ii) Search for ground state, where $\\theta_k$ are updated via VQE and $B[l]$ are calculated via functions in Section 2.6." ] }, { "cell_type": "code", "execution_count": 7, "id": "9da5b86b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "g: 1.5, nbas: 4\n", "Iter 0 VQE energy: -3.1462862688238387\n", "Iter 1 VQE energy: -3.1467101393155392\n", "Iter 2 VQE energy: -3.1516072207104755\n", "Iter 3 VQE energy: -3.146809472707097\n", "Iter 4 VQE energy: -3.1436199703762058\n", "Iter 5 VQE energy: -3.1436211792549265\n", "Iter 6 VQE energy: -3.1436161320531255\n", "Iter 7 VQE energy: -3.1436163145946456\n", "Iter 8 VQE energy: -3.143610688816215\n", "Iter 9 VQE energy: -3.143615904093609\n", "g: 1.5, nbas: 8\n", "Iter 0 VQE energy: -3.143618182488605\n", "Iter 1 VQE energy: -3.200840457440427\n", "Iter 2 VQE energy: -3.2122982581916806\n", "Iter 3 VQE energy: -3.2148415941913107\n", "Iter 4 VQE energy: -3.214454540682066\n", "Iter 5 VQE energy: -3.214612010154128\n", "Iter 6 VQE energy: -3.2147110746262615\n", "Iter 7 VQE energy: -3.2147665497225395\n", "Iter 8 VQE energy: -3.2148131940534586\n", "Iter 9 VQE energy: -3.2154252365889784\n", "g: 1.5, nbas: 12\n", "Iter 0 VQE energy: -3.1514261568965116\n", "Iter 1 VQE energy: -3.207712921399021\n", "Iter 2 VQE energy: -3.2143507766993564\n", "Iter 3 VQE energy: -3.2151165867343625\n", "Iter 4 VQE energy: -3.215551620056812\n", "Iter 5 VQE energy: -3.21568849222177\n", "Iter 6 VQE energy: -3.215871719331443\n", "Iter 7 VQE energy: -3.215941510624362\n", "Iter 8 VQE energy: -3.2160628060732046\n", "Iter 9 VQE energy: -3.216149074896171\n", "g: 1.5, nbas: 16\n", "Iter 0 VQE energy: -3.1514406055444124\n", "Iter 1 VQE energy: -3.2076420281600306\n", "Iter 2 VQE energy: -3.2142584040499904\n", "Iter 3 VQE energy: -3.215027472815238\n", "Iter 4 VQE energy: -3.215298468441994\n", "Iter 5 VQE energy: -3.215473370191188\n", "Iter 6 VQE energy: -3.215603547226202\n", "Iter 7 VQE energy: -3.2157277687243084\n", "Iter 8 VQE energy: -3.21585866940053\n", "Iter 9 VQE energy: -3.215945718302776\n", "g: 1.5, nbas: 20\n", "Iter 0 VQE energy: -3.1514390786099793\n", "Iter 1 VQE energy: -3.207490813088342\n", "Iter 2 VQE energy: -3.214102161196929\n", "Iter 3 VQE energy: -3.2148182787029125\n", "Iter 4 VQE energy: -3.215041112156586\n", "Iter 5 VQE energy: -3.2151813495719064\n", "Iter 6 VQE energy: -3.2153433350190888\n", "Iter 7 VQE energy: -3.215488882659097\n", "Iter 8 VQE energy: -3.2156246556276216\n", "Iter 9 VQE energy: -3.2157174737499026\n", "g: 1.5, nbas: 24\n", "Iter 0 VQE energy: -3.1514387774077157\n", "Iter 1 VQE energy: -3.207367390115564\n", "Iter 2 VQE energy: -3.2139409419868894\n", "Iter 3 VQE energy: -3.214681092894311\n", "Iter 4 VQE energy: -3.2148647347163752\n", "Iter 5 VQE energy: -3.2149976397508198\n", "Iter 6 VQE energy: -3.2150898148896068\n", "Iter 7 VQE energy: -3.215209498240965\n", "Iter 8 VQE energy: -3.215350330348768\n", "Iter 9 VQE energy: -3.215489480343017\n", "g: 1.5, nbas: 28\n", "Iter 0 VQE energy: -3.151439525393469\n", "Iter 1 VQE energy: -3.2073726435058725\n", "Iter 2 VQE energy: -3.21394347501351\n", "Iter 3 VQE energy: -3.2146357208139835\n", "Iter 4 VQE energy: -3.2147953021340396\n", "Iter 5 VQE energy: -3.2149183613661005\n", "Iter 6 VQE energy: -3.2150120210391324\n", "Iter 7 VQE energy: -3.2151326321686358\n", "Iter 8 VQE energy: -3.2152609684807225\n", "Iter 9 VQE energy: -3.2153772130446274\n", "g: 1.5, nbas: 32\n", "Iter 0 VQE energy: -3.1514407536942515\n", "Iter 1 VQE energy: -3.2073747584311443\n", "Iter 2 VQE energy: -3.213913960989289\n", "Iter 3 VQE energy: -3.214575999427622\n", "Iter 4 VQE energy: -3.2147908760265476\n", "Iter 5 VQE energy: -3.214908337484859\n", "Iter 6 VQE energy: -3.214995903540079\n", "Iter 7 VQE energy: -3.2151048002575178\n", "Iter 8 VQE energy: -3.2151780903906246\n", "Iter 9 VQE energy: -3.215267190685132\n", "[array(-3.1436159), array(-3.21542524), array(-3.21614907), array(-3.21594572), array(-3.21571747), array(-3.21548948), array(-3.21537721), array(-3.21526719)]\n", "g: 3, nbas: 4\n", "Iter 0 VQE energy: -5.970388609022889\n", "Iter 1 VQE energy: -5.970388750660544\n", "Iter 2 VQE energy: -5.970380747881972\n", "Iter 3 VQE energy: -5.9703823363984405\n", "Iter 4 VQE energy: -5.970361506210726\n", "Iter 5 VQE energy: -5.970387322685639\n", "Iter 6 VQE energy: -5.970395203902614\n", "Iter 7 VQE energy: -5.970391647335198\n", "Iter 8 VQE energy: -5.970399166505385\n", "Iter 9 VQE energy: -5.970381692433662\n", "g: 3, nbas: 8\n", "Iter 0 VQE energy: -5.97038378078057\n", "Iter 1 VQE energy: -7.5564103165274705\n", "Iter 2 VQE energy: -7.867598000881585\n", "Iter 3 VQE energy: -7.886883542215796\n", "Iter 4 VQE energy: -7.8885422575017525\n", "Iter 5 VQE energy: -7.88860270448971\n", "Iter 6 VQE energy: -7.888976697629529\n", "Iter 7 VQE energy: -7.889032282397925\n", "Iter 8 VQE energy: -7.87508588628647\n", "Iter 9 VQE energy: -7.88070265416792\n", "g: 3, nbas: 12\n", "Iter 0 VQE energy: -5.97067566818566\n", "Iter 1 VQE energy: -8.243806185428152\n", "Iter 2 VQE energy: -8.734803721364637\n", "Iter 3 VQE energy: -8.761741304437681\n", "Iter 4 VQE energy: -8.763332322992516\n", "Iter 5 VQE energy: -8.763771458365108\n", "Iter 6 VQE energy: -8.764185047602272\n", "Iter 7 VQE energy: -8.764384932213805\n", "Iter 8 VQE energy: -8.76466850977698\n", "Iter 9 VQE energy: -8.764911481717341\n", "g: 3, nbas: 16\n", "Iter 0 VQE energy: -5.96862038743091\n", "Iter 1 VQE energy: -8.429794930919833\n", "Iter 2 VQE energy: -9.026873692202507\n", "Iter 3 VQE energy: -9.056442107809515\n", "Iter 4 VQE energy: -9.05694061960338\n", "Iter 5 VQE energy: -9.057045669005818\n", "Iter 6 VQE energy: -9.057238086937822\n", "Iter 7 VQE energy: -9.057282106387905\n", "Iter 8 VQE energy: -9.057566284361702\n", "Iter 9 VQE energy: -9.057828042375913\n", "g: 3, nbas: 20\n", "Iter 0 VQE energy: -5.970618486961766\n", "Iter 1 VQE energy: -8.46027786865908\n", "Iter 2 VQE energy: -9.083147656265663\n", "Iter 3 VQE energy: -9.114099808352503\n", "Iter 4 VQE energy: -9.11560775880261\n", "Iter 5 VQE energy: -9.115877265080734\n", "Iter 6 VQE energy: -9.116097672280478\n", "Iter 7 VQE energy: -9.116229576435904\n", "Iter 8 VQE energy: -9.116393057481334\n", "Iter 9 VQE energy: -9.116547264265742\n", "g: 3, nbas: 24\n", "Iter 0 VQE energy: -5.970625270862546\n", "Iter 1 VQE energy: -8.45711233270321\n", "Iter 2 VQE energy: -9.088394950174024\n", "Iter 3 VQE energy: -9.119585611496353\n", "Iter 4 VQE energy: -9.120622157156225\n", "Iter 5 VQE energy: -9.12076874586676\n", "Iter 6 VQE energy: -9.120868946965105\n", "Iter 7 VQE energy: -9.120961208842665\n", "Iter 8 VQE energy: -9.121266704999321\n", "Iter 9 VQE energy: -9.121400092721187\n", "g: 3, nbas: 28\n", "Iter 0 VQE energy: -5.970428599007128\n", "Iter 1 VQE energy: -8.457873653767907\n", "Iter 2 VQE energy: -9.080983847442083\n", "Iter 3 VQE energy: -9.11887015866691\n", "Iter 4 VQE energy: -9.120581625216305\n", "Iter 5 VQE energy: -9.120780165844552\n", "Iter 6 VQE energy: -9.120861448711434\n", "Iter 7 VQE energy: -9.120944883939789\n", "Iter 8 VQE energy: -9.121064294222425\n", "Iter 9 VQE energy: -9.121158073032065\n", "g: 3, nbas: 32\n", "Iter 0 VQE energy: -5.969871735326939\n", "Iter 1 VQE energy: -8.45468700588645\n", "Iter 2 VQE energy: -9.081070798150298\n", "Iter 3 VQE energy: -9.11871782591465\n", "Iter 4 VQE energy: -9.120482174048545\n", "Iter 5 VQE energy: -9.1206758213174\n", "Iter 6 VQE energy: -9.12076495864178\n", "Iter 7 VQE energy: -9.12087305152589\n", "Iter 8 VQE energy: -9.120974560170358\n", "Iter 9 VQE energy: -9.121128971234404\n", "[array(-3.1436159), array(-3.21542524), array(-3.21614907), array(-3.21594572), array(-3.21571747), array(-3.21548948), array(-3.21537721), array(-3.21526719), array(-5.97038169), array(-7.88070265), array(-8.76491148), array(-9.05782804), array(-9.11654726), array(-9.12140009), array(-9.12115807), array(-9.12112897)]\n" ] } ], "source": [ "vqe_e = []\n", "thetas = np.zeros((max(param_ids) + 1) * n_layers)\n", "\n", "for g in [1.5, 3]:\n", " for nbas in [4, 8, 12, 16, 20, 24, 28, 32]:\n", " print(f\"g: {g}, nbas: {nbas}\")\n", "\n", " # take gray encoding as an initial guess for `b_array`\n", " b_list = []\n", " for i in range(max(dof_nature) + 1):\n", " b = np.eye(nbas)[:nbas_v] # nbas_dummy * nbas\n", " b_list.append(b)\n", " b_array = np.array(b_list)\n", "\n", " # initialize, get hamitonians and basis, see section 2.5\n", " terms, basis = get_ham_terms_and_basis(g, nbas)\n", " model = Model(basis, terms)\n", " h_mpo = Mpo(model)\n", "\n", " # searching for the ground state.\n", " for i_iter in range(10):\n", " h_contracted = get_contracted_mpo(\n", " h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom\n", " )\n", " # get \\theta_k via VQE\n", " opt_res = scipy.optimize.minimize(opt_fn, args=(h_contracted,), x0=thetas / 2, jac=True, method=\"L-BFGS-B\")\n", " print(f\"Iter {i_iter} VQE energy: {opt_res.fun}\")\n", " thetas = opt_res.x\n", " psi = ansatz(thetas).real\n", " # Update b[l] via functions in section 2.6\n", " for i in range(len(b_array)):\n", " b_array[i] = solve_b_array(psi, h_mpo, b_array, i)\n", " vqe_e.append(opt_res.fun)\n", "\n", " print(vqe_e)" ] }, { "cell_type": "code", "execution_count": 8, "id": "a1f91b3b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Energy')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGwCAYAAABRgJRuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAts0lEQVR4nO3deXgUVb7G8bcTSYclHZIQspgQCCCI7IsI4gYooGzXBVQQIlxGFFwRTcbRADOSAC53dBjFq8Myg/rghshcA44gyiIoig5bFC4CmoBMkDTL2Cyp+0ekr00INJ2Q6j79/TxPP49Vfbr6d05K+7XqVJXDsixLAAAAIS7C7gIAAACqA6EGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIF9hdQE0qKytTUVGRYmJi5HA47C4HAAD4wbIsHTx4UKmpqYqIqPx4TFiFmqKiIqWnp9tdBgAACMDu3buVlpZW6fthFWpiYmIklQ+Ky+WyuRoAAOAPt9ut9PR07+94ZcIq1Jw85eRyuQg1AACEmLNNHWGiMAAAMAKhBgAAGCFkQs3AgQPVqFEjRUdHKyUlRXfccYeKiorsLgsAAASJkAk111xzjRYsWKDCwkK99dZb2r59u26++Wa7ywIAAEHCYVmWZXcRgVi0aJEGDx4sj8ejWrVq+fUZt9ut2NhYlZaWMlEYAIAQ4e/vd0he/bR//37Nnz9f3bt3P2Og8Xg88ng83mW3210T5QEAABuEzOknSXr00UdVt25dJSQkaNeuXXr33XfP2D4vL0+xsbHeFzfeAwDAXLaGmuzsbDkcjjO+tm7d6m0/ceJEffnll1q6dKkiIyM1YsQInensWU5OjkpLS72v3bt310S3AACADWydU7Nv3z6VlJScsU1mZqaioqIqrP/++++Vnp6u1atXq1u3bn59H3NqAAAIPSExpyYxMVGJiYkBfbasrEySfObM2OFEmaV1O/brx4M/q2FMtC5tEq/ICB6WCQBATQuJicJr167VZ599ph49eiguLk7bt2/X448/rqZNm/p9lOZ8KNhYrMnvbVZx6c/edSmx0cod0Ep9W6fYVhcAAOEoJCYK16lTR2+//bZ69eqlFi1aaPTo0Wrbtq1WrFghp9NpS00FG4t199++8Ak0krSn9Gfd/bcvVLCx2Ja6AAAIVyF7n5pAVNecmhNllnpMW1Yh0JzkkJQcG62Vj/bkVBQAAFUUEnNqQtW6HfsrDTSSZEkqLv1Z63bsV7emCTVXmE3CfV5RuPdfYgzof3j3X2IMgqX/hJoA/Hiw8kATSLtQFu7zisK9/xJjQP/Du/8SYxBM/Q+JOTXBpmFMdLW2C1XhPq8o3PsvMQb0P7z7LzEGwdZ/Qk0ALm0Sr5TYaFV2YM2h8pR6aZP4miyrRp0oszT5vc063YSsk+smv7dZJ8rMnLIV7v2XGAP6H979lxiDYOw/oSYAkREO5Q5oJUkVgs3J5dwBrYw+n3ou84pMFO79lxgD+h/e/ZcYg2DsP6EmQH1bp+iF4R2VHOt7iik5NlovDO9o/HnUcJ9XFO79lxgD+h/e/ZcYg2DsPxOFq6Bv6xRd2yo5KGZ817Rwn1cU7v2XGAP6H979lxiDYOw/R2qqKDLCoW5NEzSo/YXq1jQhLAKNxLyicO+/xBjQ//Duv8QYBGP/CTUISLjPKwr3/kuMAf0P7/5LjEEw9p9Qg4CF+7yicO+/xBjQ//Duv8QYBFv/eUwCqixY7iRpl3Dvv8QY0P/w7r/EGJzv/vv7+02oAQAAQc3f329OPwEAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAI4RcqPF4PGrfvr0cDoc2bNhgdzkAACBIhFyoeeSRR5Sammp3GQAAIMiEVKh5//33tXTpUj311FN2lwIAAILMBXYX4K+9e/dqzJgxWrhwoerUqePXZzwejzwej3fZ7Xafr/IAAIDNQuJIjWVZysrK0tixY9W5c2e/P5eXl6fY2FjvKz09/TxWCQAA7GRrqMnOzpbD4Tjja+vWrXr++ed18OBB5eTknNP2c3JyVFpa6n3t3r37PPUEAADYzWFZlmXXl+/bt08lJSVnbJOZmakhQ4bovffek8Ph8K4/ceKEIiMjNWzYMM2dO9ev73O73YqNjVVpaalcLleVagcAADXD399vW0ONv3bt2uUzH6aoqEh9+vTRm2++qa5duyotLc2v7RBqAAAIPf7+fofEROFGjRr5LNerV0+S1LRpU78DDQAAMFtITBQGAAA4m5A4UnOqxo0bKwTOmgEAgBrEkRoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjhEyoady4sRwOh88rPz/f7rIAAECQuMDuAs7FlClTNGbMGO9yTEyMjdUAAIBgElKhJiYmRsnJyX6393g88ng83mW3230+ygIAAEEgZE4/SVJ+fr4SEhLUoUMHzZgxQ8ePHz9j+7y8PMXGxnpf6enpNVQpAACoaQ7Lsiy7i/DHM888o44dOyo+Pl6rV69WTk6O7rzzTj3zzDOVfuZ0R2rS09NVWloql8tVE2UDAIAqcrvdio2NPevvt62hJjs7W9OmTTtjmy1btqhly5YV1v/lL3/RXXfdpUOHDsnpdPr1ff4OCgAACB4hEWr27dunkpKSM7bJzMxUVFRUhfWbNm1S69attXXrVrVo0cKv7yPUAAAQevz9/bZ1onBiYqISExMD+uyGDRsUERGhhg0bVnNVAAAgFIXE1U9r1qzR2rVrdc011ygmJkZr1qzRgw8+qOHDhysuLs7u8gAAQBAIiVDjdDr1+uuva9KkSfJ4PGrSpIkefPBBPfTQQ3aXBgAAgkRIhJqOHTvq008/tbsMAAAQxELqPjUAAACVIdQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBECCjWHDx+u7joAAACqJKBQk5SUpFGjRmnlypXVXQ8AAEBAAgo1f/vb37R//3717NlTF110kfLz81VUVFTdtQEAAPgtoFAzePBgLVy4UD/88IPGjh2rV199VRkZGerfv7/efvttHT9+vLrrBAAAOCOHZVlWdWzo+eef18SJE3X06FE1aNBAY8eOVXZ2turUqVMdm68WbrdbsbGxKi0tlcvlsrscAADgB39/vy+oypfs3btXc+fO1Zw5c7Rz507dfPPNGj16tL7//ntNmzZNn376qZYuXVqVrwAAAPBLQKHm7bff1uzZs7VkyRK1atVK99xzj4YPH6769et723Tv3l0XX3xxddUJAABwRgGFmjvvvFO33nqrVq1apS5dupy2TWpqqh577LEqFRcSyk5IO1dLh/ZK9ZKkjO5SRKTdVQEAEHYCmlNz5MiRoJor469qn1OzeZFU8Kjk/tWVX65Uqe80qdXAqm8fAACc3zk1x48fl9vtrrDe4XDI6XQqKioqkM2Gls2LpAUjJJ2SCd3F5euHzCPYAABQgwK6pLt+/fqKi4ur8Kpfv75q166tjIwM5ebmqqysrFqL/fvf/66uXbuqdu3aiouL0+DBg6t1+34rO1F+hObUQCP9/7qC7PJ2AACgRgR0pGbOnDl67LHHlJWVpUsvvVSStG7dOs2dO1e/+93vtG/fPj311FNyOp367W9/Wy2FvvXWWxozZoymTp2qnj176vjx49q4cWO1bPuc7Vzte8qpAkty/1DerskVNVYWAADhLKBQM3fuXD399NMaMmSId92AAQPUpk0bzZo1Sx9++KEaNWqkJ598slpCzfHjx3X//fdrxowZGj16tHd9q1atzvg5j8cjj8fjXT7dKbOAHNpbve0AAECVBXT6afXq1erQoUOF9R06dNCaNWskST169NCuXbuqVt0vvvjiC/3www+KiIhQhw4dlJKSon79+p31SE1eXp5iY2O9r/T09GqpR/WSqrcdAACosoBCTXp6ul555ZUK61955RVvcCgpKVFcXFzVqvvF//7v/0qSJk2apN/97ndavHix4uLidPXVV2v//v2Vfi4nJ0elpaXe1+7du6ulHmV0L7/KSY5KGjgk14Xl7QAAQI0I6PTTU089pVtuuUXvv/++9z41n3/+ubZu3ao333xTkvTZZ59p6NChZ9xOdna2pk2bdsY2W7Zs8U44fuyxx3TTTTdJkmbPnq20tDS98cYbuuuuu077WafTKafTeU5980tEZPll2wtGqDzY/HrC8C9Bp28+96sBAKAGBRRqBg4cqMLCQs2aNUuFhYWSpH79+mnhwoVq3LixJOnuu+8+63YmTJigrKysM7bJzMxUcXGxJN85NE6nU5mZmdV2iuuctRpYftn2ae9Tk8/l3AAA1LBzDjXHjh1T37599eKLLyovL69KX56YmKjExMSztuvUqZOcTqcKCwvVo0cPbx3fffedMjIyqlRDlbQaKLW8gTsKAwAQBM451NSqVUtff/31+ailUi6XS2PHjlVubq7S09OVkZGhGTNmSJJuueWWGq2lgohILtsGACAIBDRRePjw4aedKHw+zZgxQ7feeqvuuOMOdenSRTt37tSyZcuqbTIyAAAIbQE9++nee+/VvHnz1Lx5c3Xq1El169b1ef+ZZ56ptgKrU7U/+wkAAJx35/XZTxs3blTHjh0lSd98843Pew5HZZc5AwAAnD8BhZrly5dXdx0AAABVEtCcmpO2bdumJUuW6N///rckKYAzWQAAANUioFBTUlKiXr166aKLLtL111/vvY/M6NGjNWHChGotEAAAwB8BhZoHH3xQtWrV0q5du1SnTh3v+qFDh6qgoKDaigMAAPBXQHNqli5dqiVLligtLc1nffPmzbVz585qKQwAAOBcBHSk5vDhwz5HaE7av3//+XnWEgAAwFkEFGquuOIKzZs3z7vscDhUVlam6dOn65prrqm24gAAAPwV0Omn6dOnq1evXvr888919OhRPfLII9q0aZP279+vVatWVXeNAAAAZxXQkZrWrVvrm2++UY8ePTRo0CAdPnxYN954o7788ks1bdq0umsEAAA4q4AekxCqeEwCAACh57w+JkGSDhw4oHXr1unHH39UWVmZz3sjRowIdLMAAAABCSjUvPfeexo2bJgOHTokl8vl87wnh8NBqAEAADUuoDk1EyZM0KhRo3To0CEdOHBAP/30k/e1f//+6q4RAADgrAIKNT/88IPuu+++096rBgAAwA4BhZo+ffro888/r+5aAAAAAhbQnJobbrhBEydO1ObNm9WmTRvVqlXL5/2BAwdWS3EAAAD+CuiS7oiIyg/wOBwOnThxokpFnS9c0g0AQOg5r5d0n3oJNwAAgN3OaU7N9ddfr9LSUu9yfn6+Dhw44F0uKSlRq1atqq04AAAAf51TqFmyZIk8Ho93eerUqT6XcB8/flyFhYXVVx0AAICfzinUnDr9JoyesAAAAIJcQJd0AwAABJtzCjUOh8PnkQgn1wEAANjtnK5+sixLWVlZcjqdkqSff/5ZY8eOVd26dSXJZ74NAABATTqnUDNy5Eif5eHDh1dow8MsAQCAHc4p1MyePft81QEAAFAlTBQGAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAI4REqPnoo4/kcDhO+/rss8/sLg8AAASBC+wuwB/du3dXcXGxz7rHH39cH374oTp37mxTVQAAIJiERKiJiopScnKyd/nYsWN69913de+998rhcNhYGQAACBYhEWpOtWjRIpWUlOjOO+88YzuPxyOPx+Nddrvd57s0AABgk5CYU3OqV155RX369FFaWtoZ2+Xl5Sk2Ntb7Sk9Pr6EKAQBATbM11GRnZ1c6Afjka+vWrT6f+f7777VkyRKNHj36rNvPyclRaWmp97V79+7z1RUAAGAzW08/TZgwQVlZWWdsk5mZ6bM8e/ZsJSQkaODAgWfdvtPplNPprEqJAAAgRNgaahITE5WYmOh3e8uyNHv2bI0YMUK1atU6j5UBAIBQE1JzapYtW6YdO3boP//zP+0uBQAABJmQCjWvvPKKunfvrpYtW9pdCgAACDIhdUn3q6++ancJAAAgSIXUkRoAAIDKEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMMIFdhcAA5SdkHaulg7tleolSRndpYhIu6sCAIQZQg2qZvMiqeBRyV30/+tcqVLfaVKrgfbVBQAIO5x+QuA2L5IWjPANNJLkLi5fv3mRPXUBAMISoQaBKTtRfoRG1mne/GVdQXZ5OwAAagChBoHZubriERofluT+obwdAAA1gFCDwBzaW73tAACoIkINAlMvqXrbAQBQRYQaBCaje/lVTnJU0sAhuS4sbwcAQA0ImVDzzTffaNCgQWrQoIFcLpd69Oih5cuX211W+IqILL9sW1LFYPPLct987lcDAKgxIRNq+vfvr+PHj2vZsmVav3692rVrp/79+2vPnj12lxa+Wg2UhsyTXCm+612p5eu5Tw0AoAY5LMs63TW5QeVf//qXEhMT9fHHH+uKK66QJB08eFAul0sffPCBevfu7dd23G63YmNjVVpaKpfLdT5LDi/cURgAcB75+/sdEncUTkhIUIsWLTRv3jx17NhRTqdTs2bNUsOGDdWpU6dKP+fxeOTxeLzLbre7JsoNPxGRUpMr7K4CABDmQiLUOBwO/eMf/9DgwYMVExOjiIgINWzYUAUFBYqLi6v0c3l5eZo8eXINVgoAAOxi65ya7OxsORyOM762bt0qy7I0btw4NWzYUJ988onWrVunwYMHa8CAASouLq50+zk5OSotLfW+du/eXYO9AwAANcnWOTX79u1TSUnJGdtkZmbqk08+0XXXXaeffvrJ51xa8+bNNXr0aGVnZ/v1fcypAQAg9ITEnJrExEQlJiaetd2RI0ckSRERvgeWIiIiVFZWdl5qAwAAoSUkLunu1q2b4uLiNHLkSH311Vf65ptvNHHiRO3YsUM33HCD3eUBAIAgEBKhpkGDBiooKNChQ4fUs2dPde7cWStXrtS7776rdu3a2V0eAAAIAiFxn5rqwpwaAABCj7+/3yFxpAYAAOBsCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEa4wO4CgJBXdkLauVo6tFeqlyRldJciIu2uCgDCDqEGqIrNi6SCRyV30f+vc6VKfadJrQbaVxcAhCFOPwGB2rxIWjDCN9BIkru4fP3mRfbUBQBhilADBKLsRPkRGlmnefOXdQXZ5e0AADWCUAMEYufqikdofFiS+4fydgCAGkGoAQJxaG/1tgMAVBmhBghEvaTqbQcAqDJCDRCIjO7lVznJUUkDh+S6sLwdAKBGEGqAQEREll+2LalisPlluW8+96sBgBpEqAEC1WqgNGSe5ErxXe9KLV/PfWoAoEZx8z2gKloNlFrewB2FASAIEGqAqoqIlJpcYXcVABD2OP0EAACMQKgBAABGCJlQ88UXX+jaa69V/fr1lZCQoN/85jc6dOiQ3WUBAIAgERKhpqioSL1791azZs20du1aFRQUaNOmTcrKyrK7NAAAECRCYqLw4sWLVatWLc2cOVMREeU57MUXX1Tbtm21bds2NWvWzOYKgTBXdoIrwADYLiRCjcfjUVRUlDfQSFLt2rUlSStXrqw01Hg8Hnk8Hu+y2+0+v4UC4WjzovInlv/6AZ+u1PKbE3KvHgA1KCROP/Xs2VN79uzRjBkzdPToUf3000/Kzs6WJBUXF1f6uby8PMXGxnpf6enpNVUyEB42L5IWjKj4xHJ3cfn6zYvsqQtAWLI11GRnZ8vhcJzxtXXrVl1yySWaO3eunn76adWpU0fJyclq0qSJkpKSfI7enConJ0elpaXe1+7du2uwd4Dhyk6UH6GRdZo3f1lXkF3eDgBqgMOyrNP9F6lG7Nu3TyUlJWdsk5mZqaioKO/y3r17VbduXTkcDrlcLr3++uu65ZZb/Po+t9ut2NhYlZaWyuVyVal2IOzt+ESa2//s7UYuNv/mhOE+pyjc+y8xBue5//7+fts6pyYxMVGJiYnn9JmkpCRJ0l/+8hdFR0fr2muvPR+lATibQ3urt12oCvc5ReHef4kxCKL+h8ScGkn605/+pC+++ELffPONZs6cqfHjxysvL0/169e3uzQgPNVLqt52oSjc5xSFe/8lxiDI+h8yoWbdunW69tpr1aZNG7300kuaNWuW7rvvPrvLAsJXRvfy/xuTo5IGDsl1YXk7E4X7nKJw77/EGARh/0Mm1MybN08lJSXyeDz66quvdMcdd9hdEhDeIiLLDy9Lqhhsflnum2/uvIKdqyv+36kPS3L/UN7OROHef4kxCML+h0yoARCEWg2UhsyTXCm+612p5etNnk8Q7nOKwr3/EmMQhP0PiZvvAQhirQZKLW8Ivys/wn1OUbj3X2IMgrD/hBoAVRcRaf5l26c6OafIXazTzylwlL9v6pyicO+/xBgEYf85/QQAgQj3OUXh3n+JMQjC/hNqACBQ4TynSKL/EmMQZP239Y7CNY07CgM4L7ibbHj3X2IMguSOwoQaAAAQ1Pz9/eb0EwAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwQlg9pfvkzZPdbrfNlQAAAH+d/N0+20MQwirUHDx4UJKUnp5ucyUAAOBcHTx4ULGxsZW+H1bPfiorK1NRUZFiYmLkcJz6mHT7uN1upaena/fu3TyTKkCMYdUwflXHGFYdY1g1Jo+fZVk6ePCgUlNTFRFR+cyZsDpSExERobS0NLvLqJTL5TJuR6xpjGHVMH5VxxhWHWNYNaaO35mO0JzERGEAAGAEQg0AADACoSYIOJ1O5ebmyul02l1KyGIMq4bxqzrGsOoYw6ph/MJsojAAADAXR2oAAIARCDUAAMAIhBoAAGAEQg0AADACocZGkyZNksPh8Hm1bNnS7rKC2scff6wBAwYoNTVVDodDCxcu9Hnfsiw98cQTSklJUe3atdW7d299++239hQbhM42fllZWRX2yb59+9pTbBDKy8tTly5dFBMTo4YNG2rw4MEqLCz0afPzzz9r3LhxSkhIUL169XTTTTdp7969NlUcfPwZw6uvvrrCfjh27FibKg4uL7zwgtq2beu9wV63bt30/vvve98P9/2PUGOzSy65RMXFxd7XypUr7S4pqB0+fFjt2rXTzJkzT/v+9OnT9dxzz+nFF1/U2rVrVbduXfXp00c///xzDVcanM42fpLUt29fn33ytddeq8EKg9uKFSs0btw4ffrpp/rggw907NgxXXfddTp8+LC3zYMPPqj33ntPb7zxhlasWKGioiLdeOONNlYdXPwZQ0kaM2aMz344ffp0myoOLmlpacrPz9f69ev1+eefq2fPnho0aJA2bdokif1PFmyTm5trtWvXzu4yQpYk65133vEul5WVWcnJydaMGTO86w4cOGA5nU7rtddes6HC4Hbq+FmWZY0cOdIaNGiQLfWEoh9//NGSZK1YscKyrPL9rVatWtYbb7zhbbNlyxZLkrVmzRq7ygxqp46hZVnWVVddZd1///32FRVi4uLirJdffpn9z7IsjtTY7Ntvv1VqaqoyMzM1bNgw7dq1y+6SQtaOHTu0Z88e9e7d27suNjZWXbt21Zo1a2ysLLR89NFHatiwoVq0aKG7775bJSUldpcUtEpLSyVJ8fHxkqT169fr2LFjPvtgy5Yt1ahRI/bBSpw6hifNnz9fDRo0UOvWrZWTk6MjR47YUV5QO3HihF5//XUdPnxY3bp1Y/9TmD3QMth07dpVc+bMUYsWLVRcXKzJkyfriiuu0MaNGxUTE2N3eSFnz549kqSkpCSf9UlJSd73cGZ9+/bVjTfeqCZNmmj79u367W9/q379+mnNmjWKjIy0u7ygUlZWpgceeECXX365WrduLal8H4yKilL9+vV92rIPnt7pxlCSbr/9dmVkZCg1NVVff/21Hn30URUWFurtt9+2sdrg8c9//lPdunXTzz//rHr16umdd95Rq1attGHDhrDf/wg1NurXr5/3n9u2bauuXbsqIyNDCxYs0OjRo22sDOHq1ltv9f5zmzZt1LZtWzVt2lQfffSRevXqZWNlwWfcuHHauHEj8+CqoLIx/M1vfuP95zZt2iglJUW9evXS9u3b1bRp05ouM+i0aNFCGzZsUGlpqd58802NHDlSK1assLusoMDppyBSv359XXTRRdq2bZvdpYSk5ORkSaow03/v3r3e93BuMjMz1aBBA/bJU4wfP16LFy/W8uXLlZaW5l2fnJyso0eP6sCBAz7t2QcrqmwMT6dr166SxH74i6ioKDVr1kydOnVSXl6e2rVrpz/+8Y/sfyLUBJVDhw5p+/btSklJsbuUkNSkSRMlJyfrww8/9K5zu91au3atunXrZmNloev7779XSUkJ++QvLMvS+PHj9c4772jZsmVq0qSJz/udOnVSrVq1fPbBwsJC7dq1i33wF2cbw9PZsGGDJLEfVqKsrEwej4f9T5x+stXDDz+sAQMGKCMjQ0VFRcrNzVVkZKRuu+02u0sLWocOHfL5v7UdO3Zow4YNio+PV6NGjfTAAw/oD3/4g5o3b64mTZro8ccfV2pqqgYPHmxf0UHkTOMXHx+vyZMn66abblJycrK2b9+uRx55RM2aNVOfPn1srDp4jBs3Tq+++qreffddxcTEeOcpxMbGqnbt2oqNjdXo0aP10EMPKT4+Xi6XS/fee6+6deumyy67zObqg8PZxnD79u169dVXdf311yshIUFff/21HnzwQV155ZVq27atzdXbLycnR/369VOjRo108OBBvfrqq/roo4+0ZMkS9j+JS7rtNHToUCslJcWKioqyLrzwQmvo0KHWtm3b7C4rqC1fvtySVOE1cuRIy7LKL+t+/PHHraSkJMvpdFq9evWyCgsL7S06iJxp/I4cOWJdd911VmJiolWrVi0rIyPDGjNmjLVnzx67yw4apxs7Sdbs2bO9bf79739b99xzjxUXF2fVqVPH+o//+A+ruLjYvqKDzNnGcNeuXdaVV15pxcfHW06n02rWrJk1ceJEq7S01N7Cg8SoUaOsjIwMKyoqykpMTLR69eplLV261Pt+uO9/DsuyrJoMUQAAAOcDc2oAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQaoAw8t1338nhcHifpRMMtm7dqssuu0zR0dFq3759QNuYM2eO6tevX611hRuHw6GFCxfaXQZQJYQaoAZlZWXJ4XAoPz/fZ/3ChQvlcDhsqspeubm5qlu3rgoLC30exPdrJ8fN4XB4n1A8ZcoUHT9+vIarPT+ysrJ4PhlQDQg1QA2Ljo7WtGnT9NNPP9ldSrU5evRowJ/dvn27evTooYyMDCUkJFTarm/fviouLta3336rCRMmaNKkSZoxY0bA3wvAPIQaoIb17t1bycnJysvLq7TNpEmTKpyK+a//+i81btzYu3zy/+6nTp2qpKQk1a9f33v0YuLEiYqPj1daWppmz55dYftbt25V9+7dFR0drdatW2vFihU+72/cuFH9+vVTvXr1lJSUpDvuuEP/+te/vO9fffXVGj9+vB544AE1aNCg0qd4l5WVacqUKUpLS5PT6VT79u1VUFDgfd/hcGj9+vWaMmWKHA6HJk2aVOmYOJ1OJScnKyMjQ3fffbd69+6tRYsW+bRZsmSJLr74YtWrV88bgvyt5eSpubffflvXXHON6tSpo3bt2mnNmjU+3/HWW2/pkksukdPpVOPGjfX000/7vN+4cWNNnTpVo0aNUkxMjBo1aqSXXnqp0n7540x/j5deekmpqakqKyvz+cygQYM0atQo7/K7776rjh07Kjo6WpmZmZo8eXKlR7qOHj2q8ePHKyUlRdHR0crIyDjj/goEC0INUMMiIyM1depUPf/88/r++++rtK1ly5apqKhIH3/8sZ555hnl5uaqf//+iouL09q1azV27FjdddddFb5n4sSJmjBhgr788kt169ZNAwYMUElJiSTpwIED6tmzpzp06KDPP/9cBQUF2rt3r4YMGeKzjblz5yoqKkqrVq3Siy++eNr6/vjHP+rpp5/WU089pa+//lp9+vTRwIED9e2330qSiouLdckll2jChAkqLi7Www8/7Hffa9eu7XOE6MiRI3rqqaf017/+VR9//LF27drls72z1XLSY489pocfflgbNmzQRRddpNtuu837479+/XoNGTJEt956q/75z39q0qRJevzxxzVnzhyfbTz99NPq3LmzvvzyS91zzz26++67VVhY6Hfffu1sf49bbrlFJSUlWr58ufcz+/fvV0FBgYYNGyZJ+uSTTzRixAjdf//92rx5s2bNmqU5c+boySefPO13Pvfcc1q0aJEWLFigwsJCzZ8/3ydQA0HL7seEA+Fk5MiR1qBBgyzLsqzLLrvMGjVqlGVZlvXOO+9Yv/7XMTc312rXrp3PZ5999lkrIyPDZ1sZGRnWiRMnvOtatGhhXXHFFd7l48ePW3Xr1rVee+01y7Isa8eOHZYkKz8/39vm2LFjVlpamjVt2jTLsizr97//vXXdddf5fPfu3bstSVZhYaFlWZZ11VVXWR06dDhrf1NTU60nn3zSZ12XLl2se+65x7vcrl07Kzc394zb+fW4lZWVWR988IHldDqthx9+2LIsy5o9e7Ylydq2bZv3MzNnzrSSkpL8ruXk2Lz88sve9zdt2mRJsrZs2WJZlmXdfvvt1rXXXuuzjYkTJ1qtWrXyLmdkZFjDhw/3LpeVlVkNGza0XnjhBb/6dyp//h6DBg3y7kuWZVmzZs2yUlNTvftGr169rKlTp/ps469//auVkpLiXZZkvfPOO5ZlWda9995r9ezZ0yorK6u0ZiAYcaQGsMm0adM0d+5cbdmyJeBtXHLJJYqI+P9/jZOSktSmTRvvcmRkpBISEvTjjz/6fK5bt27ef77gggvUuXNnbx1fffWVli9frnr16nlfLVu2lFQ+/+WkTp06nbE2t9utoqIiXX755T7rL7/88oD6vHjxYtWrV0/R0dHq16+fhg4d6nO6qk6dOmratKl3OSUlxdvvc6mlbdu2PtuQ5N3Oli1bTruNb7/9VidOnDjtNhwOh5KTkyv8Dfzlz99j2LBheuutt+TxeCRJ8+fP16233urdN7766itNmTLFZxtjxoxRcXGxjhw5UuE7s7KytGHDBrVo0UL33Xefli5dGlDtQE27wO4CgHB15ZVXqk+fPsrJyVFWVpbPexEREbIsy2fdsWPHKmyjVq1aPssOh+O0606db3Emhw4d0oABAzRt2rQK7538kZekunXr+r3N6nDNNdfohRdeUFRUlFJTU3XBBb7/+Tpdv08dQ3/8ejsnr0g7l/GrrJZz3cZJ/vw9BgwYIMuy9Pe//11dunTRJ598omeffdZnG5MnT9aNN95YYRvR0dEV1nXs2FE7duzQ+++/r3/84x8aMmSIevfurTfffDOgPgA1hVAD2Cg/P1/t27dXixYtfNYnJiZqz549sizL+8NanfeW+fTTT3XllVdKko4fP67169dr/Pjxksp/0N566y01bty4QnA4Fy6XS6mpqVq1apWuuuoq7/pVq1bp0ksvPeft1a1bV82aNbO1losvvlirVq3yWbdq1SpddNFFioyMDKi2s/Hn7xEdHa0bb7xR8+fP17Zt29SiRQt17NjRZxuFhYXnNH4ul0tDhw7V0KFDdfPNN6tv377av3+/4uPjq9wn4Hwh1AA2atOmjYYNG6bnnnvOZ/3VV1+tffv2afr06br55ptVUFCg999/Xy6Xq1q+d+bMmWrevLkuvvhiPfvss/rpp5+8V8qMGzdO//3f/63bbrtNjzzyiOLj47Vt2za9/vrrevnll8/px3vixInKzc1V06ZN1b59e82ePVsbNmzQ/Pnzq6Uf56I6apkwYYK6dOmi3//+9xo6dKjWrFmjP/3pT/rzn/9c5fpKS0srBNeEhAS//x7Dhg1T//79tWnTJg0fPtxnO0888YT69++vRo0a6eabb1ZERIS++uorbdy4UX/4wx8q1PLMM88oJSVFHTp0UEREhN544w0lJydzg0MEPebUADabMmVKhVMTF198sf785z9r5syZateundatW3dOVwadTX5+vvLz89WuXTutXLlSixYtUoMGDSTJe0TjxIkTuu6669SmTRs98MADql+/vs/8HX/cd999euihhzRhwgS1adNGBQUFWrRokZo3b15tfanJWjp27KgFCxbo9ddfV+vWrfXEE09oypQpFU4fBuKjjz5Shw4dfF6TJ0/2++/Rs2dPxcfHq7CwULfffrvPtvv06aPFixdr6dKl6tKliy677DI9++yzysjIOG0tMTExmj59ujp37qwuXbrou+++0//8z/+c898fqGkOK5CTzgAAAEGG2A0AAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAI/wfgLmHzivvtJQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot the results\n", "from matplotlib import pyplot as plt\n", "\n", "nbas = [4, 8, 12, 16, 20, 24, 28, 32]\n", "plt.scatter(nbas, vqe_e[0:8], label=\"g=1.5\")\n", "plt.scatter(nbas, vqe_e[8:], label=\"g=3.0\")\n", "plt.xlabel(\"Number of Phonon Levels\")\n", "plt.ylabel(\"Energy\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.7.13" } }, "nbformat": 4, "nbformat_minor": 5 }