diff --git a/README.md b/README.md
index 0f890f5..5654120 100644
--- a/README.md
+++ b/README.md
@@ -99,5 +99,6 @@ which is also a good introduction to the software.
### Variational basis state encoder
An efficient algorithm to encode phonon states in electron-phonon systems for quantum computation.
-See [examples](https://github.com/tencent-quantum-lab/TenCirChem/tree/master/example).
+See [examples](https://github.com/tencent-quantum-lab/TenCirChem/tree/master/example)
+and the [tutorial](https://tencent-quantum-lab.github.io/TenCirChem/tutorial_jupyter/vbe_tutorial_td.html).
Reference paper: https://arxiv.org/pdf/2301.01442.pdf (published in PRR).
diff --git a/README_CN.md b/README_CN.md
index 52e5c84..831e9c3 100644
--- a/README_CN.md
+++ b/README_CN.md
@@ -96,6 +96,7 @@ TenCirChem 使用 Academic Public License 发布。
### 变分基矢编码器
变分基矢编码器是编码电声子系统中的声子用以量子计算的高效算法。
-请参考[例子](https://github.com/tencent-quantum-lab/TenCirChem/tree/master/example).
+请参考[例子](https://github.com/tencent-quantum-lab/TenCirChem/tree/master/example)
+及[教程](https://tencent-quantum-lab.github.io/TenCirChem/tutorial_jupyter/vbe_tutorial_td.html)。
参考论文: https://arxiv.org/pdf/2301.01442.pdf (发表于PRR).
diff --git a/docs/source/tutorial_jupyter/figs/vbe_gs_Fig1.svg b/docs/source/tutorial_jupyter/figs/vbe_gs_Fig1.svg
new file mode 100644
index 0000000..072fbb3
--- /dev/null
+++ b/docs/source/tutorial_jupyter/figs/vbe_gs_Fig1.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/tutorial_jupyter/figs/vbe_gs_Fig2.svg b/docs/source/tutorial_jupyter/figs/vbe_gs_Fig2.svg
new file mode 100644
index 0000000..aa8218e
--- /dev/null
+++ b/docs/source/tutorial_jupyter/figs/vbe_gs_Fig2.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/tutorial_jupyter/figs/vbe_td_Fig1.svg b/docs/source/tutorial_jupyter/figs/vbe_td_Fig1.svg
new file mode 100644
index 0000000..73bc659
--- /dev/null
+++ b/docs/source/tutorial_jupyter/figs/vbe_td_Fig1.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/tutorial_jupyter/figs/vbe_td_Fig2.svg b/docs/source/tutorial_jupyter/figs/vbe_td_Fig2.svg
new file mode 100644
index 0000000..8ceb0c1
--- /dev/null
+++ b/docs/source/tutorial_jupyter/figs/vbe_td_Fig2.svg
@@ -0,0 +1 @@
+
\ No newline at end of file
diff --git a/docs/source/tutorial_jupyter/vbe_tutorial_groundstate.ipynb b/docs/source/tutorial_jupyter/vbe_tutorial_groundstate.ipynb
new file mode 100644
index 0000000..ce8cceb
--- /dev/null
+++ b/docs/source/tutorial_jupyter/vbe_tutorial_groundstate.ipynb
@@ -0,0 +1,656 @@
+{
+ "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",
+ "\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",
+ "\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
+}
\ No newline at end of file
diff --git a/docs/source/tutorial_jupyter/vbe_tutorial_td.ipynb b/docs/source/tutorial_jupyter/vbe_tutorial_td.ipynb
new file mode 100644
index 0000000..475547b
--- /dev/null
+++ b/docs/source/tutorial_jupyter/vbe_tutorial_td.ipynb
@@ -0,0 +1,903 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "959d157b-cd6c-46f5-8762-e30dba930d20",
+ "metadata": {},
+ "source": [
+ "# Variational Basis State Encoder (Time Dependent)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cf2da1ef",
+ "metadata": {},
+ "source": [
+ "## 1 Background\n",
+ "This tutorial is for the time evolution of variational basis state encoder (VBE). Here, spin-boson model is taken as an example. As shall be presented below, the key is to calculate $\\dot{\\theta}_j$ and $\\dot{C}[l]^*$. The algorithm realization are presented in section 2. For more theoretical derivations, see https://doi.org/10.1103/PhysRevResearch.5.023046."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "78152079",
+ "metadata": {},
+ "source": [
+ "## 2 Algorithm Realization"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4b13a28e",
+ "metadata": {},
+ "source": [
+ "### 2.1 Imports"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "3fdceaa4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from scipy.integrate import solve_ivp\n",
+ "from opt_einsum import contract\n",
+ "import tensorcircuit as tc\n",
+ "\n",
+ "from tencirchem import set_backend, Op, Mpo, Model, OpSum\n",
+ "from tencirchem.dynamic import get_ansatz, get_deriv, get_jacobian_func, qubit_encode_basis, sbm\n",
+ "from tencirchem.applications.vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f90cb726",
+ "metadata": {},
+ "source": [
+ "### 2.2 Initialize\n",
+ "In this tutorial, we study the time evolution of the following spin-boson model:\n",
+ "\n",
+ "$$\\hat{H}=\\frac{\\epsilon}{2} \\hat{\\sigma}_z + \\Delta\\hat{\\sigma}_x + \\sum_j g_j \\omega_j \\hat{\\sigma}_z (\\hat{b}^\\dagger_j + \\hat{b}_j)\n",
+ "+ \\sum_j \\omega_j \\hat{b}^\\dagger_j \\hat{b}_j$$\n",
+ "\n",
+ "Here, `epsilon`, `delta`, `omega_list` and `g_list` correspond to $\\epsilon$, $\\Delta$, $\\omega_j$ and $g_j$ in the Hamiltonian, respectively.\n",
+ "In this section, the parameters are defined and the circuit is initialized. The schematic diagram of the circuit is plotted in Fig. 1, where blue square corresponds to qubit representing spin and green circles correspond to qubits representing vibrations. \n",
+ "\n",
+ "Fig. 1 Schematic diagram of the circuit."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "766e26d2",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "2023-08-08 16:12:43.927000: 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', 'v-2-0-bottom', 'v-2-1-bottom'] \n",
+ "psi_index_bottom: ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'v-2-0-bottom', 'v-2-1-bottom'] \n",
+ "b_dof_vidx: [array([1, 2]), array([3, 4])] \n",
+ "psi_shape2: [2, 2, 2, 2, 2]\n"
+ ]
+ }
+ ],
+ "source": [
+ "set_backend(\"jax\")\n",
+ "\n",
+ "epsilon = 0\n",
+ "delta = 1\n",
+ "omega_list = [0.5, 1]\n",
+ "g_list = [0.25, 1]\n",
+ "\n",
+ "nmode = len(omega_list) # number of phonon modes\n",
+ "# make sure correct input\n",
+ "assert nmode == len(g_list)\n",
+ "\n",
+ "# two qubit for each mode\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] * n_qubit_per_mode + [1] * n_qubit_per_mode)\n",
+ "b_dof_pidx = np.array([1, 2]) # index of basis that need VBE\n",
+ "\n",
+ "n_dof = len(dof_nature)\n",
+ "psi_shape2 = [2] * n_dof\n",
+ "\n",
+ "psi_idx_top, psi_idx_bottom, b_dof_vidx = get_psi_indices(dof_nature, b_dof_pidx, n_qubit_per_mode)\n",
+ "print(\n",
+ " \"psi_index_top: \",\n",
+ " psi_idx_bottom,\n",
+ " \"\\npsi_index_bottom: \",\n",
+ " psi_idx_bottom,\n",
+ " \"\\nb_dof_vidx: \",\n",
+ " b_dof_vidx,\n",
+ " \"\\npsi_shape2: \",\n",
+ " psi_shape2,\n",
+ ")\n",
+ "\n",
+ "\n",
+ "# spin boson model has been define in tencirchem.dynamic.model\n",
+ "def get_model(epsilon, delta, nmode, omega_list, g_list, nlevels):\n",
+ " ham_terms = sbm.get_ham_terms(epsilon, delta, nmode, omega_list, g_list)\n",
+ " basis = sbm.get_basis(omega_list, nlevels)\n",
+ " return Model(basis, ham_terms)\n",
+ "\n",
+ "\n",
+ "nbas = 16 # number of phonon levels (basis)\n",
+ "\n",
+ "b_shape = tuple([2] * n_qubit_per_mode + [nbas])\n",
+ "\n",
+ "assert len(omega_list) == nmode\n",
+ "assert len(g_list) == nmode\n",
+ "model = get_model(epsilon, delta, nmode, omega_list, g_list, [nbas] * nmode)\n",
+ "\n",
+ "h_mpo = Mpo(model)\n",
+ "\n",
+ "# generate the quantum circuit with defined qubits\n",
+ "circuit = tc.Circuit(1 + nmode * n_qubit_per_mode)\n",
+ "psi0 = circuit.state()\n",
+ "n_layers = 3 # layers of ansatz"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "41235844",
+ "metadata": {},
+ "source": [
+ "### 2.3 Get variational Hamiltonian ansatz terms\n",
+ "In this section, we will generate variational hamiltonian ansatz terms. The following ansatz is adopted:\n",
+ "\n",
+ "$$\\ket{\\phi} = \\prod_n e^{i \\theta_n \\hat {h}_{a_n}} \\ket{\\phi_0}$$\n",
+ "\n",
+ "if\n",
+ "\n",
+ "$$\\hat{H}=\\sum_a J_a \\hat{h}_a$$\n",
+ "\n",
+ "Note that the phonon operators have been transformed to qubit operators based on Gray Encoding."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "d8bb10fa",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[BasisHalfSpin(dof: spin, nbas: 2), BasisHalfSpin(dof: ('v0', 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ('v0', 'TCCQUBIT-1'), nbas: 2), BasisHalfSpin(dof: ('v1', 'TCCQUBIT-0'), nbas: 2), BasisHalfSpin(dof: ('v1', 'TCCQUBIT-1'), nbas: 2)]\n",
+ "Op('X', ['spin'], 1.0)\n",
+ "Op('X', [('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y', [('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z', [('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X', [('v0', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('X X', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X Y', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X Z', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y', [('v0', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Y X', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y Y', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y Z', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z', [('v0', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z X', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z', [('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z', ['spin'], 2.0)\n",
+ "Op('Z X', ['spin', ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y', ['spin', ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z', ['spin', ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X', ['spin', ('v0', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z X X', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X Y', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X Z', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y', ['spin', ('v0', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z Y X', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y Y', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y Z', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z', ['spin', ('v0', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z Z X', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z Y', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z Z', ['spin', ('v0', 'TCCQUBIT-0'), ('v0', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X', [('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y', [('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z', [('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X', [('v1', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('X X', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X Y', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('X Z', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y', [('v1', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Y X', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y Y', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Y Z', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z', [('v1', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z X', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z', [('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X', ['spin', ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y', ['spin', ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z', ['spin', ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X', ['spin', ('v1', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z X X', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X Y', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z X Z', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y', ['spin', ('v1', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z Y X', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y Y', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Y Z', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z', ['spin', ('v1', 'TCCQUBIT-0')], 1.0)\n",
+ "Op('Z Z X', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z Y', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n",
+ "Op('Z Z Z', ['spin', ('v1', 'TCCQUBIT-0'), ('v1', 'TCCQUBIT-1')], 1.0)\n"
+ ]
+ }
+ ],
+ "source": [
+ "def get_vha_terms():\n",
+ " basis = sbm.get_basis(omega_list, [nbas_v] * nmode)\n",
+ " spin_basis = qubit_encode_basis(basis, \"gray\")\n",
+ "\n",
+ " spin_ham_terms = OpSum([Op(\"X\", [\"spin\"], 1.0)])\n",
+ " for i in range(nmode):\n",
+ " complete_list = []\n",
+ " for j in range(n_qubit_per_mode):\n",
+ " complete = OpSum()\n",
+ " dof = (f\"v{i}\", f\"TCCQUBIT-{j}\")\n",
+ " for symbol in \"IXYZ\":\n",
+ " complete += Op(symbol, dof)\n",
+ " complete_list.append(complete)\n",
+ " complete_real = complete_list[0]\n",
+ " for c in complete_list[1:]:\n",
+ " complete_real = complete_real * c\n",
+ " spin_ham_terms.extend(complete_real)\n",
+ " spin_ham_terms.extend(Op(\"Z\", \"spin\") * complete_real)\n",
+ " spin_ham_terms = OpSum([op.squeeze_identity() for op in spin_ham_terms.simplify() if not op.is_identity]).simplify()\n",
+ " return spin_ham_terms, spin_basis\n",
+ "\n",
+ "\n",
+ "spin_ham_terms, spin_basis = get_vha_terms()\n",
+ "\n",
+ "print(spin_basis)\n",
+ "for i in range(len(spin_ham_terms)):\n",
+ " print(spin_ham_terms[i])\n",
+ "\n",
+ "\n",
+ "theta0 = np.zeros(n_layers * len(spin_ham_terms), dtype=np.float64)\n",
+ "ansatz = get_ansatz(spin_ham_terms, spin_basis, n_layers, psi0)\n",
+ "jacobian_func = get_jacobian_func(ansatz)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f6a646d6",
+ "metadata": {},
+ "source": [
+ "### 2.4 Time Evolution\n",
+ "This sections defines the function involved in time evolution. The derivations $\\dot{\\theta}_j$ are calculated by following equation:\n",
+ "\n",
+ "$$\\sum_j \\textrm{Re} \\frac{\\partial \\bra{\\phi}}{\\partial \\theta_k} \\frac{\\partial \\ket{\\phi}}{\\partial \\theta_j}\\dot{\\theta}_j = \n",
+ "\\textrm{Im} \\frac{\\partial \\bra{\\phi}}{\\partial \\theta_k} \\hat {H}_1$$\n",
+ "\n",
+ "where\n",
+ "\n",
+ "$$\\hat{H}_1=\\prod_{l}\\hat{B}[l] \\hat{H} \\prod_{l}\\hat{B}[l]^\\dagger$$\n",
+ "\n",
+ ", and $\\frac{\\partial \\ket{\\phi}}{\\partial \\theta_j}$ is calculated via `theta_deriv`.\n",
+ "The time derivations of $B[l]$ are calculated by following equations:\n",
+ "\n",
+ "$$i\\rho[l]\\dot{B}[l]^* = (1-\\hat{P}[l])\\bra{\\phi}\\tilde{H}'[l]\\ket{\\phi}$$\n",
+ "$$\\dot{B}[l]^* = -i \\rho[l]^{-1}(1-\\hat{P}[l])\\bra{\\phi}\\tilde{H}'[l]\\ket{\\phi}$$\n",
+ "\n",
+ "where\n",
+ "\n",
+ "$$\\tilde{H}'=\\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 projection operator\n",
+ "\n",
+ "$$\\hat{P}[l]=\\hat{B}[l]^\\dagger\\hat{B}[l]$$\n",
+ "\n",
+ "and the reduced density matrix\n",
+ "\n",
+ "$$\\rho[l]_{nn'}=\\textrm{Tr} \\left[ \\bra{\\phi}\\ket{n}_l \\bra{n'}_l \\ket{\\phi} \\right ]$$\n",
+ "\n",
+ "An example that illustrates the function in detail is also presented below.\n",
+ "\n",
+ "Fig. 2 Graphic illustration of (a) $\\rho[l]$, (b) $P[l]$, (c) $\\bra\\phi\\tilde{H}'[l]\\ket\\phi$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "7597768e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def deriv_fun(t, theta_and_b):\n",
+ " # split \\thera and b to independent arrays\n",
+ " theta = theta_and_b[: len(theta0)]\n",
+ " psi = ansatz(theta)\n",
+ " b_array = theta_and_b[len(theta0) :].reshape(nmode, nbas_v, nbas)\n",
+ "\n",
+ " # get contracted H\n",
+ " h_contracted = get_contracted_mpo(h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom)\n",
+ " # get the derivation of \\theta\n",
+ " theta_deriv = get_deriv(ansatz, jacobian_func, theta, h_contracted)\n",
+ "\n",
+ " psi = psi.reshape(psi_shape2)\n",
+ " b_deriv_list = []\n",
+ " for i in range(nmode):\n",
+ " b = b_array[i]\n",
+ " # calculate rho\n",
+ " indices_base = [(\"contract\", ii) for ii in range(n_dof)]\n",
+ " psi_top_indices = indices_base.copy()\n",
+ " psi_bottom_indices = indices_base.copy()\n",
+ " for j in b_dof_vidx[i]:\n",
+ " psi_top_indices[j] = (\"top\", j)\n",
+ " psi_bottom_indices[j] = (\"bottom\", j)\n",
+ " out_indices = [(\"top\", j) for j in b_dof_vidx[i]] + [(\"bottom\", j) for j in b_dof_vidx[i]]\n",
+ " args = [psi.conj(), psi_top_indices, psi, psi_bottom_indices, out_indices]\n",
+ " rho = contract(*args).reshape(1 << n_qubit_per_mode, 1 << n_qubit_per_mode)\n",
+ "\n",
+ " from scipy.linalg import pinv\n",
+ "\n",
+ " rho += np.eye(len(rho)) * 1e-5\n",
+ " rho_inv = pinv(rho)\n",
+ "\n",
+ " b = b.reshape(nbas_v, nbas)\n",
+ " # calculate projector\n",
+ " proj = b.conj().T @ b\n",
+ "\n",
+ " # derivative\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",
+ " args.append(b_array[i].reshape(b_shape))\n",
+ " args.append([f\"v-{k}-{l}-bottom\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-bottom\"])\n",
+ " # output indices\n",
+ " args.append([f\"v-{k}-{l}-top\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-top\", \"mpo-0\", f\"mpo-{len(h_mpo)}\"])\n",
+ "\n",
+ " # take transpose to be compatible with previous code\n",
+ " b_deriv = contract(*args).squeeze().reshape(nbas_v, nbas).T\n",
+ " b_deriv = np.einsum(\"bf, bg -> fg\", b_deriv, np.eye(nbas) - proj)\n",
+ " b_deriv = -1j * np.einsum(\"fg, fh -> hg\", b_deriv, rho_inv.T)\n",
+ " b_deriv_list.append(b_deriv)\n",
+ " return np.concatenate([theta_deriv, np.array(b_deriv_list).ravel()])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c552000f",
+ "metadata": {},
+ "source": [
+ "Let's illustrate the code in detail.\n",
+ "Code that calculates $\\rho[l]$ (see Fig. 2(a))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "1477e876",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "mode 0\n",
+ "out_indices: [('top', 1), ('top', 2), ('bottom', 1), ('bottom', 2)]\n",
+ "args: \n",
+ " [DeviceArray([[[[[1.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]]],\n",
+ "\n",
+ "\n",
+ "\n",
+ " [[[[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]]]], dtype=complex128), [('contract', 0), ('top', 1), ('top', 2), ('contract', 3), ('contract', 4)], DeviceArray([[[[[1.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]]],\n",
+ "\n",
+ "\n",
+ "\n",
+ " [[[[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]]]], dtype=complex128), [('contract', 0), ('bottom', 1), ('bottom', 2), ('contract', 3), ('contract', 4)], [('top', 1), ('top', 2), ('bottom', 1), ('bottom', 2)]]\n",
+ "rho: \n",
+ " [[1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# split \\thera and b to independent arrays\n",
+ "\n",
+ "b_list = []\n",
+ "for _ in range(nmode):\n",
+ " b = np.eye(nbas)[:nbas_v] # nbas_v * nbas\n",
+ " b_list.append(b)\n",
+ "theta_and_b = np.concatenate([theta0, np.array(b_list).ravel()]).astype(complex)\n",
+ "\n",
+ "theta = theta_and_b[: len(theta0)]\n",
+ "psi = ansatz(theta)\n",
+ "b_array = theta_and_b[len(theta0) :].reshape(nmode, nbas_v, nbas)\n",
+ "\n",
+ "# get contracted H\n",
+ "h_contracted = get_contracted_mpo(h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom)\n",
+ "# get the derivation of \\theta\n",
+ "theta_deriv = get_deriv(ansatz, jacobian_func, theta, h_contracted)\n",
+ "\n",
+ "psi = psi.reshape(psi_shape2)\n",
+ "b_deriv_list = []\n",
+ "i = 0\n",
+ "print(\"mode \", i)\n",
+ "b = b_array[i]\n",
+ "\n",
+ "# calculate rho[l]\n",
+ "indices_base = [(\"contract\", ii) for ii in range(n_dof)]\n",
+ "psi_top_indices = indices_base.copy()\n",
+ "psi_bottom_indices = indices_base.copy()\n",
+ "for j in b_dof_vidx[i]:\n",
+ " psi_top_indices[j] = (\"top\", j)\n",
+ " psi_bottom_indices[j] = (\"bottom\", j)\n",
+ "out_indices = [(\"top\", j) for j in b_dof_vidx[i]] + [(\"bottom\", j) for j in b_dof_vidx[i]]\n",
+ "print(\"out_indices: \", out_indices)\n",
+ "args = [psi.conj(), psi_top_indices, psi, psi_bottom_indices, out_indices]\n",
+ "print(\"args: \\n\", args)\n",
+ "rho = contract(*args).reshape(1 << n_qubit_per_mode, 1 << n_qubit_per_mode)\n",
+ "print(\"rho: \\n\", rho)\n",
+ "\n",
+ "from scipy.linalg import pinv\n",
+ "\n",
+ "rho += np.eye(len(rho)) * 1e-5\n",
+ "rho_inv = pinv(rho)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1e4ebde6",
+ "metadata": {},
+ "source": [
+ "Code that calculates $\\hat{P}[l]$, (see Fig. 2(b))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "6780faa5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "proj: [[1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+ " [0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j\n",
+ " 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "b = b.reshape(nbas_v, nbas)\n",
+ "# calculate projector P[l]\n",
+ "proj = b.conj().T @ b\n",
+ "print(\"proj: \", proj)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "08928a9f",
+ "metadata": {},
+ "source": [
+ "Code that calculates $\\bra{\\phi}\\tilde{H}'[l]\\ket{\\phi}$ (see Fig. 2(c))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "5115b03a",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "args in derivative: \n",
+ " [DeviceArray([[[[[1.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]]],\n",
+ "\n",
+ "\n",
+ "\n",
+ " [[[[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j]]]]], dtype=complex128), ['p-0-top', 'v-1-0-top', 'v-1-1-top', 'v-2-0-top', 'v-2-1-top'], DeviceArray([[[[[1.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]]],\n",
+ "\n",
+ "\n",
+ "\n",
+ " [[[[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]],\n",
+ "\n",
+ "\n",
+ " [[[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j]]]]], dtype=complex128), ['p-0-bottom', 'v-1-0-bottom', 'v-1-1-bottom', 'v-2-0-bottom', 'v-2-1-bottom'], , ['mpo-0', 'p-0-top', 'p-0-bottom', 'mpo-1'], , ['mpo-1', 'p-1-top', 'p-1-bottom', 'mpo-2'], , ['mpo-2', 'p-2-top', 'p-2-bottom', 'mpo-3'], array([[[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]]]), ['v-2-0-bottom', 'v-2-1-bottom', 'p-2-bottom'], array([[[1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
+ " 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
+ " 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j]],\n",
+ "\n",
+ " [[0.-0.j, 0.-0.j, 1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
+ " 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j],\n",
+ " [0.-0.j, 0.-0.j, 0.-0.j, 1.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j,\n",
+ " 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j, 0.-0.j]]]), ['v-2-0-top', 'v-2-1-top', 'p-2-top'], array([[[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]],\n",
+ "\n",
+ " [[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
+ " [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
+ " 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]]]), ['v-1-0-bottom', 'v-1-1-bottom', 'p-1-bottom'], ['v-1-0-top', 'v-1-1-top', 'p-1-top', 'mpo-0', 'mpo-3']]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# derivative\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",
+ "args.append(b_array[i].reshape(b_shape))\n",
+ "args.append([f\"v-{k}-{l}-bottom\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-bottom\"])\n",
+ "# output indices\n",
+ "args.append([f\"v-{k}-{l}-top\" for l in range(n_qubit_per_mode)] + [f\"p-{k}-top\", \"mpo-0\", f\"mpo-{len(h_mpo)}\"])\n",
+ "print(\"args in derivative: \\n\", args)\n",
+ "\n",
+ "# take transpose to be compatible with previous code\n",
+ "b_deriv = contract(*args).squeeze().reshape(nbas_v, nbas).T"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e5197573",
+ "metadata": {},
+ "source": [
+ "Contract tensors and calculate $\\dot{B}[l]^*$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "1921f511",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "b_deriv = np.einsum(\"bf, bg -> fg\", b_deriv, np.eye(nbas) - proj)\n",
+ "b_deriv = -1j * np.einsum(\"fg, fh -> hg\", b_deriv, rho_inv.T)\n",
+ "b_deriv_list.append(b_deriv)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "075f0e23",
+ "metadata": {},
+ "source": [
+ "## 2.5 Main Structure of the Function\n",
+ "This is the main stucture of the code. We set parameters for time evolution and initialize the system. Then, we update `theta_and_b` in time evolution and calcutate $\\langle \\sigma_z \\rangle$ `z` and $\\langle \\sigma_x \\rangle$ `x` step by step as the time evolves."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "04e5b336",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "0 1.0\n",
+ "1 0.9801368820199942\n",
+ "2 0.9221622297057478\n",
+ "3 0.8307164281767181\n",
+ "4 0.712902625694087\n",
+ "5 0.5774908043700016\n",
+ "6 0.4340489975856556\n",
+ "7 0.29215043357543824\n",
+ "8 0.1607516949764221\n",
+ "9 0.04771089791814444\n",
+ "10 -0.040502642430565255\n",
+ "11 -0.09904787191380995\n",
+ "12 -0.12493534752914795\n",
+ "13 -0.11727456640490394\n",
+ "14 -0.07753204265425653\n",
+ "15 -0.009732009482324864\n",
+ "16 0.0795128404043315\n",
+ "17 0.18126828150880098\n",
+ "18 0.2849420221890886\n",
+ "19 0.37936481084384754\n",
+ "20 0.454192768911805\n",
+ "21 0.5014051186981502\n",
+ "22 0.5165776398039041\n",
+ "23 0.49960262655096854\n",
+ "24 0.454624729908236\n",
+ "25 0.38915769753158064\n",
+ "26 0.31259780970942813\n",
+ "27 0.2345813095925135\n",
+ "28 0.1636633691524335\n",
+ "29 0.10653226239863313\n",
+ "30 0.06766729976756988\n",
+ "31 0.04926876472593326\n",
+ "32 0.05139921686178229\n",
+ "33 0.07232882932406118\n",
+ "34 0.10899831508240451\n",
+ "35 0.15745403958628454\n",
+ "36 0.21315107699434557\n",
+ "37 0.27112011043252915\n",
+ "38 0.3261026901074539\n",
+ "39 0.3728086768061548\n",
+ "40 0.40637666288671287\n",
+ "41 0.4229683862535386\n",
+ "42 0.42032486089352117\n",
+ "43 0.3981218276511652\n",
+ "44 0.358045462603464\n",
+ "45 0.303614837159146\n",
+ "46 0.23991582373561202\n",
+ "47 0.17311347674692437\n",
+ "48 0.10948894653071074\n",
+ "49 0.05440718215103625\n",
+ "50 0.011717807210271824\n",
+ "51 -0.016427087542680092\n",
+ "52 -0.02981545362112685\n",
+ "53 -0.030162355384666367\n",
+ "54 -0.020379375716308275\n",
+ "55 -0.0033909784547891857\n",
+ "56 0.018576266600602798\n",
+ "57 0.04368717824936481\n",
+ "58 0.06971819393413235\n",
+ "59 0.09364096597083588\n",
+ "60 0.11186039266314592\n",
+ "61 0.12076156704155695\n",
+ "62 0.11755115511155295\n",
+ "63 0.10090752714979544\n",
+ "64 0.07138982669111976\n",
+ "65 0.03253624871953725\n",
+ "66 -0.011552300380283664\n",
+ "67 -0.05607764510179323\n",
+ "68 -0.09628855507563477\n",
+ "69 -0.12836308326138923\n",
+ "70 -0.14994039236833612\n",
+ "71 -0.1596458713814789\n",
+ "72 -0.15576794369843444\n",
+ "73 -0.1373071561893098\n",
+ "74 -0.10671249014129201\n",
+ "75 -0.06743879382842288\n",
+ "76 -0.022285878126101024\n",
+ "77 0.024280809221682687\n",
+ "78 0.06810377510953083\n",
+ "79 0.10678673871439198\n",
+ "80 0.13914885388461803\n",
+ "81 0.1611886576331503\n",
+ "82 0.16581373326788112\n",
+ "83 0.1542497273174869\n",
+ "84 0.13048909287272364\n",
+ "85 0.09708896314020657\n",
+ "86 0.057310903820287765\n",
+ "87 0.015533678436258006\n",
+ "88 -0.024192504501293052\n",
+ "89 -0.05790366853755518\n",
+ "90 -0.0797592370408079\n",
+ "91 -0.08582739907891124\n",
+ "92 -0.07396770719727741\n",
+ "93 -0.048080239560664026\n",
+ "94 -0.012411290850183795\n",
+ "95 0.031070370091046573\n",
+ "96 0.08042891268644264\n",
+ "97 0.13037553828085885\n",
+ "98 0.17451072357728148\n",
+ "99 0.20513840470656897\n"
+ ]
+ }
+ ],
+ "source": [
+ "theta0 = np.zeros(n_layers * len(spin_ham_terms), dtype=np.float64)\n",
+ "ansatz = get_ansatz(spin_ham_terms, spin_basis, n_layers, psi0)\n",
+ "jacobian_func = get_jacobian_func(ansatz)\n",
+ "\n",
+ "b_list = []\n",
+ "for _ in range(nmode):\n",
+ " b = np.eye(nbas)[:nbas_v] # nbas_v * nbas\n",
+ " b_list.append(b)\n",
+ "theta_and_b = np.concatenate([theta0, np.array(b_list).ravel()]).astype(complex)\n",
+ "z_list = [1]\n",
+ "x_list = [0]\n",
+ "\n",
+ "tau = 0.1\n",
+ "steps = 100\n",
+ "\n",
+ "dummy_model = get_model(epsilon, delta, nmode, omega_list, g_list, [nbas_v] * nmode)\n",
+ "z_op = Mpo(dummy_model, Op(\"Z\", \"spin\", factor=1)).todense()\n",
+ "x_op = Mpo(dummy_model, Op(\"X\", \"spin\", factor=1)).todense()\n",
+ "\n",
+ "for n in range(steps):\n",
+ " print(n, float(z_list[-1]))\n",
+ " # update `theta_and_b`\n",
+ " sol = solve_ivp(deriv_fun, [n * tau, (n + 1) * tau], theta_and_b)\n",
+ " theta_and_b = sol.y[:, -1]\n",
+ " theta = theta_and_b[: len(theta0)]\n",
+ " psi = ansatz(theta)\n",
+ " z = psi.conj().T @ (z_op @ psi)\n",
+ " x = psi.conj().T @ (x_op @ psi)\n",
+ " z_list.append(z.real)\n",
+ " x_list.append(x.real)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "82642da3",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEeCAYAAABrB7XiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBbElEQVR4nO3de1yb9d0//lcSjqWFEHqmJ0LPra0E2qq1altQp7ZOBbRzurmtibrtvncS7m73d7e7t3sV1p3c4deg29Rtdm3SOut0WlLbarW2QHqwZ0joAWgpEC6gnCHX74/kSoEGyOE6JXk/H48+7hvCdeUTB+988v6835+PgmVZFoQQQsKaUuoBEEIIER4Fe0IIiQAU7AkhJAJQsCeEkAgQJfUApLZo0SKkp6f7fV1tbS1SU1NFu06K56TrIvM6KZ6TruP3OqvVipqamsHfZCOcTqdj161bx7755pt+Xbdu3bqAni/Q66R4TrouMq+T4jnpOn6ue/PNN9l169axSUlJNz0W8TP71NRU7N69W7Tn27Bhg+jXBvOcYj4fvT5+rwtGuL/GcH19GzZswIYNG5CZmXnzgwG93YQRKWZNoSLcXyO9vtAX7q+Rz/hEC7QBkmK2JbZwf430+kJfuL9GPl+fgmWl66C12+2wWCwwmUwoLS316Zri4mKo1WoAAMMwKCgo8OvxodavXy9qGocQQoTmLa5JNrO3Wq2wWCxgGAYOh8Ona7hArtfrodfrodPpYDAYfH6cEEIilaQzewAwm83YvHkzKioqRv3Z5ORkVFdXe2buAKBQKMC9hNEe94Zm9oSQcCOrmb2/7HY7GIYZFMg5Fotl1Mf55nTS/nGEkNARUsHeG7VaDYZhRn2cT6Un6rDqx++jsa2L1/sSQohQQr7OXqPRwOFweJ3RD3x8OLW1tVi/fr3na65OdTi9fU4U/t2Kyiut+NofP8VbL9wDlTJk3jMJIWFo27Zt2LZtm+fr2tram34m5IP9aIu7oz3ub1NVdJQSf/v2nVj94gfYd+oqfrrzBF7Mu9Xn6wkhhG9DJ6kDJ7CckJmSarVar99nGAZarXbUx/m0cJoav//6CgDAL985jXcqLvN6f0II4VtIBXu1Wu01N5+dnT3q43zLu30Wnr93HgDAYDyEhlbK3xNC5EvyYD9cmsVut6O4uHjQ9zZt2jSossZsNkOv1/v8ON9+9kQG5qcmoa2rD6Un6gR7HkIICZZkwd5ut6OwsBBGoxFWqxWFhYUoKSnxPG6xWGA0GgddU1BQAIZhUFJSgpKSEpSVlQ36mdEe51t0lBIP6lzbj+47eVWw5yGEkGBJ3lQltWCbqg6cvoqHXvoQk5LiUPnyI1AoFDyOjhBC/BfSTVVydducCYiPUaG+pQtnalukHg4hhHhFwT5IsdEqrJw3EQDwIaVyCCEyFfHBnmuqGtiQ4K/ViycDAPadvMLXsAghxG/btm3D+vXrw7OpKlh8nFS1epEr2B88ew3dvf2IjVbxMTRCCPEL11wV0k1VcrZomhoTEuPQ0dOPI1WNUg+HEEJuQsGeB0qlAqsXTQIA7DtFeXtCiPxQsOfJPYu4vD0Fe0KI/FCw5wmXt7dWO9DR3SfxaAghZDAK9jyZlpKAlHGxcLIsztW1Sj0cQggZhII9jxZOSwIAnK5hpB0IIYQMQcGeRwtT1QBAnbSEENmhYM+jBe6Z/Rma2RNCZCbigz0fHbScBanuYE8ze0KIBKiDdgR8dNByFkxTAwAuN3WgtbMXifHRvNyXEEJ8QR20IklOiMGU5HgAwFma3RNCZISCPc8WplJFDiFEfijY82w+l7evoZk9IUQ+KNjzbKE7b0+LtIQQOaFgz7MF1FhFCJEhCvY8mz/VFezrW7rQ1NYt8WgIIcSFgj3PxsVHY+b4BACUyiGEyEfEB3s+m6o43CItlV8SQsQ0UlNVxAd7rqlqw4YNvN2TW6SlvD0hREwbNmzA7t27kZqaetNjER/shXBjkZZm9oQQeaBgL4B5UxIBALb6NolHQgghLhTsBTBr4lgAwFWmk06tIoTIAgV7ASQnxCBpjGsTtEuN7RKPhhBCKNgLQqFQYNYE1+z+QsN1iUdDCCEU7AUzkwv21yjYE0KkR8FeIDSzJ4TIScQHeyGaqgAgzb1IW00ze0KISOikqhHweVLVQLMmuLZMoJk9IUQsdFKVBLjyywvXroNlWYlHQwiJdBTsBTI9JQEKBdDR04+G1i6ph0MIiXCSp3GKi4uhVqsBAAzDoKCgYMSfz8vLQ05ODrKysjzXcbRaLSwWC4xGI3JycqDValFaWoply5YhNzdXoFfgXWy0CqnJY1Dj6MCFhnZMTIoX9fkJIWQgSYM9F+j1ej0AwGKxwGAwwGg0DnuN1WqF2Wy+6fu5ubkwmUxgGMbzM1qtFoWFhaIHes6siWNdwf7adSyfPV6SMRBCCCBxGmfz5s3Iz8/3fJ2dnY2SkpIRr8nNzQXLsoP+GY1GmEwmz89UVFSAZVnYbDbPG4kUqPySECIXkgV7u90OhmFuSsUArhn+cAwGw00/m5WVxffweEHll4QQuZAsjWO3271+X61Wg2GYYa/TarWe/59hGNjtdmRnZw/6mR07dkCj0QAAysrKUFRUFPyAAzCTyi8JITIh+QLtUBqNBg6Hw6efLSwsvCm/r9VqodPpPG8KDocDeXl5g9I8A3FNVRyuTpUPs2jLBEKICLZt2zaoMTQkmqp8DfRWq9Xr93U63aCv8/PzYTAYhk0ZCdVUBdxI49Q2d6Cnrx8xUSpBnocQEtmGTlJl1VQ1MB0zEMMwwz42kNFoRHp6+k3fH1qpwwX44dJGQpqQGIcxMSqwLG11TAiRlqTBXq1Wew3CQ3Pw3lgslptm6gzDIC8vb9A9ufy/L28gfFMoFLT7JSFEFiQtvdy0adOgyhuz2TyoVNJut6O4uNjrtXa7/aYArlarYTQaB32/pKQEubm5XlM4YvBsm9BAM3tCiHQkzdkXFBSguLjYU1tvs9kGLbhy3bDeumq1Wq2n4mag/Pz8QW8QTU1Nwy7OiiFtApVfEkKkJ/kC7UjbI+j1+mGbomw2m9fvq9XqUbdcEBNXfnmxkYI9IUQ6tBGawKaluIJ9bVOHxCMhhEQyCvYCm+4O9pebKGdPCJFOxAd7oU6q4kxLGQMAuNbahZ6+fkGegxBCADqpakRCNlUBwPhxsYiNVqK714m65k5PVy0hhPCNTqqSkEKhwDSNa3ZfQ6kcQohEKNiLIFXjytvX0CItIUQiFOxFwOXtL1OwJ4RIhIK9CLg0Tq2D0jiEEGlQsBdBagqlcQgh0qJgL4Lp7jROjYOCPSFEGhTsReBJ41A1DiFEIhEf7IVuqgJupHGYjl60dfYK9jyEkMg2UlNVxAd7rqmKr6MIvUmMj0bSmGgAQC2lcgghAtmwYQN2796N1NTUmx6L+GAvlmmeRVpK5RBCxEfBXiSpGqq1J4RIh4K9SG7U2lOwJ4SIj4K9SDxpHAr2hBAJULAXCbdlAuXsCSFSoGAvEk9jFeXsCSESoGAvEm7ny1pHB1iWlXg0hJBIE/HBXoymKgCYmhwPAOjq7UdjW7egz0UIiUzUVDUCMZqqACA2WoVJSXEAqCKHECIMaqqSiRv72tMiLSFEXBTsRcSVX9bSIi0hRGQU7EXkaaxqpmBPCBEXBXsRTUl2Bfs6ytkTQkRGwV5EqRpXRU6to1PikRBCIg0FexFN5Wb2lMYhhIiMgr2IuJ0v65qpsYoQIq6ID/ZiNVUBwBR3Y1V3rxNN16mxihDCL2qqGoFYTVUAEBOlwoREV2NVHeXtCSE8o6YqGbmxSEt5e0KIeCjYi4wWaQkhUoiSegDFxcVQq9UAAIZhUFBQMOLPWywWGI1G5OTkQKvVorS0FMuWLUNubm7A9xSTZ5E2RGf23b39iI1WST0MQoifJJ3Zc0FZr9dDr9dDp9PBYDCMeA3DMLBarTAYDDAYDEhPT/ca6P25p5imerpoQytn393bj6/+4SCmP2fGjk8vSD0cQoifFKyENYDJycmorq72zMIBQKFQjFiWaDabkZ2dPeiaYO65fv167N69O5DhB+TNg3YYSj7D6kWTsbtwjWjPG4yunn58+Xcf44PjdQCAKJUC279zN+5dOlXikRFCvPEW1ySb2dvtdjAM4zVoWywW2dyTb6khdvB4Z08fnvjNAXxwvA7xMSrcs3AS+vpZfPl3H+NwZYPUwyOE+EjSYO+NWq0GwzAjXrtjxw6YzWaYzWYUFhbyck+xcAu0V0JkgfbX/zqNvSevIiE2Cju/fw92/uAe5CyZgs6efuT96gAaWrukHiIhxAeSL9AOpdFo4HA4hn1cq9VCp9NBq9UCABwOB/Ly8mAymQK6J9dUxdmwYYOgNfdczr6tqw+tnb1IjI8W7LmCxbIsdhy6CADY8nQWVi2YBAD467dXYfWLH+BMbQt2fnYRz947T8phEhLxtm3bNqgx1FtTleyC/UiBHgB0Ot2gr/Pz82EwGEacuY90T66pSiwJsVFITohBc3sP6hwdSExNEu25/XXqMgNbfRtio5V4OGu65/sJsVH46j3pKPy7FebDFOwJkdrQSerACSwn4DTOrl27sGvXrkAv98zMh2IYZtjHANcC7UBcft5utwd8T7Fx2ybIPW//z7LLAIDsW6Zi3JBPII8snwGFAjhc2YhLjXTyFiFy51ewP3bsGJ599lksW7YMdrsdNpsNWVlZeO6553Ds2DG/nlir1UKtVnvNs2dnZ3u9hmEY5OXlDbqGm9FrtdqA7imFUFmk/WfZJQDAF5dNv+mxKcljsGq+K62z8/BFUcdFCPHfqMG+tbUVv/jFL5CVlQWj0QiDwYCysjL84Ac/wAsvvIDy8nLo9Xps3boVy5Ytw5YtW9Da2urTk2/atGlQlYzZbIZer/d8bbfbUVxc7PlarVbDaDQOmqWXlJQgNzfXM8Mf7Z5yEAqLtGdqGJyra0VMlBJfyLh5nw0AeOy2mQAA8yEK9oTI3bA5+507d8JoNEKhUCAvLw/l5eXD3iQjIwNbt24FALzyyitYs2YNUlJSYDAY8Oijjw57XUFBAYqLi1FSUgIAsNlsMBqNnse5btmBHbD5+fmD3gCampoGLc6Odk85SA2BxiouhbNm8WQkjYnx+jMPZ03H998ow4lLzThX14J5U+W7/kBIpPPaVJWfn4/09HTo9XqkpaUFdOPq6moYjUZUV1dj+/btQQ9UKGI3VQHA6wds+NafDuPepVOx8/v3iPrcvlrxw3dxuqYFRv1t+NKdw693PPbL/dhzvA6bvrgYP3x0iYgjJIQMx1tc8zqz37FjR9BPlpaWhpdeeino+4Sjqe4FWrnuj3OurgWna1oQrVLigYxpI/5s3m0zsed4HUyfXaRgT4iM8dJUdeHCBT5uEzHkvkC7/9RVAMCqBROhTvCewuE8qJuGKJUCVVfbqCqHEBnjJdgnJyd7qnGOHj2KTZs28XFbUYh5UhWHW6Btbu9BZ0+faM/rqwq7qy9hxezxo/7suPho3DpLAwD49Nw1QcdFCBmZ4CdVVVdXo6CgAB9++CEyMjKQn5+P++67j49bC07Mk6o4SWOikRDryqDVyXCR1lrdBADQaVN8+vk75k4EAHx6nvbKIURKgp9UtXbtWhQXFyMzMxMXLlxAeXk5bDYbH7cOSwqF4sZWx03ySuW0dvbi/BVX6Wymj8H+9rkTANDMnhA54yXYb968GeXl5UhKSoLNZoPVakVFRQUftw5bqdwircxq7Y9VO8CywIzxCZ7zckfDBftzda1obAu9jdFaOnrAtPdIPQxCBMXL3jh6vR4tLS0AXLN8q9WK5uZmJCVR3fVwpsp0kbaCS+GkaXy+JmVcLOanJuFsbQsOnW/AusybO27l6rPKBjz6i31o6+rD1OR4LJymxvfXLcKd8ydKPTRCeMXbFscDA/sLL7yA5uZmvm4dljzHE8psZm+1+5ev59zhSeWETt7+80vNyP3lfrR1uRbJ65o7Yfn8Ch7/9QFcbLgu8egI4ZdfwT4lxfcAkJGR4fdgIsmN8kt5LdBaq12VOJlpfgb7ea5gf+h8aOTtK6+04uHifWjp6MVtcyag8uVHYPl/OVg+ezxaO3uhLzmEfqdT6mESwhu/gn1zc7PP+96QkXHll3JqrGpo7cKlxnYoFMCtfqRxgBsVOccuNKO9W37lpAP19jmR96v9aGjtwtKZyTB//25MVsdjxZwJePXZOzA2LgqfnmvAb987K/VQCeGN32mc6upqbNmyBZs2bcKuXbso+Afoxv448gn2Fe4UztwpiX4fqjJ9fAKmacag38mirKpRiOHxxvTZBdjqr2NiUhzeemH1oL1/0iaORdGTmQCAn+08gRMXKR1JwoPfwV6n02HPnj2oqKhAQUEBkpOT8cMf/lCIsYlCiqYqAJiqcVXjXGvpQk9fv6jPPRxPvt7PFA6HS+XIuQTT6WTxm3fPAACev3ee14qjp+7S4qHMaejtd6Lo7ZNiD5GQgPHaVLV161bs2bMHe/bsQVVVFSorK+FwOEKmiWooKZqqACBlbCxio13/+a/IpLGqgsvXa/1L4XDumCf/5qo9J+pwprYF4+Ki8PU1c7z+jEKhwP97zLXPz3tHa+icXRIyeG2qevzxxwd9rdVqsXXrVuj1erz66quBjzLCKBQKpCbLp/ySZdmAK3E4t81xzeyt9iY4nTdtpioLv373NADgmdVzRtz3Z+E0NbK0KejrZ/HmwWqxhkeIYPwK9tnZ2di7d6/Xxx577DEqt/TTVBmVX15u6kBjWzeiVArcMj05oHvMm5qI+BgV2rr6UFXfxvMIg3e4sgGfnmtAtEqJb943+rm5T9+dDgB444ANXnYCJySk+BXst27dihdeeAHHjx/3+nhycmBBIlJxi7Q1Mtgy4XQNAwCYNzUJcTGqgO4RpVJiyUzX78BRd3OWnPzufVd1zRMrZ3neaEfy2G0zMSZGhfNXWvFZpbwXnQkZjV/BnkvZZGRk4Pnnnx8U9C9cuICqqireBxjOPOWXMpjZn611dUAvSA2u65nrvD3qzv/LRUtHD/591LVo9fy9o8/qASAxPhqPrHAdvfjGAdrriYQ2v3P22dnZcDgcaGhoQEZGBlQqFVQqFfLy8kK6KkcKqe6KHDk0Vp2tc5XQzp+aGNR9MtyVPFaZBfv3jtaip8+JeVMTsXiG759Av+JO5ew6fBGtnb1CDY8QwQW0XYJarYbJZILT6URVVRWqqqpQVlaGxMTgAkWk8eTsZbBAe849s58f5Mw+w723/fELDll1oL515BIA4JHlM/y67rY54zF3SiI6evqxu/yyEEMjRBRB742TlpYW8Dm1kc5TjSNxGodlWZytcwX7YA8NnzNlHMbGRaGjpx/n6+TRcNfS0YO9n18BADzqZ7BXKBSeN4g9x+t4HxshYvEa7Lds2eI5eSoYx44dw5YtW4K+j5CkaqoCbizQXmU60dsn3Sy4pqkD17v6EKVSIH3SuKDupVIqPSdXySWVMzCFs2Ca2u/rc5ZMAQDsO3kFff3y+bRCyFB+N1WtXbsWP//5zzFnzhxs2bLFry0RWltb8Ytf/AJZWVnYvHkzsrOzAx+5CKRqqgKACYlxiFIpwLJAfYt0eXtuVj97ciKio4LfCDVDZou0gaZwOFnpKUhOiAHT0YsjMt8KgkQ2v5uqMjIysGPHDlRWViIpKQlr1qzBfffdh127dg37JDt37sS9996LtWvXIjk5GeXl5di+fTtuvfVW3l5IuFEqFZ6KHCkbq/iqxOFwFTlWGZRfBpPC4aiUSmTf4prdl564wtvYCBHTqNO4jRs3ory8HFu3bsWRI0cwe/ZsPPfcczh27BiOHTuGZ599FnPmzEF5eTmMRiPKysrwjW98Q4yxh4Ub5ZdSzuz5qcThcBU5n19iJE1PAcGncDg5S6cCoLw9CV0+n1SVlpaGl156CS+99BL27t2Ln//851AoFNDr9di6dauQYwxr01LGAJXSzuzPuBuqgq3E4WgnjkXSmGi0dPTibF0LbvGj1JFv77graAJN4XC4mf2JS824ynRisjo+6LERIqaAjiVcu3Yt1q5dy/dYIpLUaRyWZXGOm9nzFOwVCgUyZmmw/3Q9rNUOyYJ9b58T+09dBQDcf+vNOUx/TEiMgy5NA2u1A6Un6vDUXel8DJEQ0fB2LCEJDNdYJVWt/ZXmTrR29kKlDL4SZyBPc5Vdurz94apGtHX1IWVcrKf+Pxj3ulM5lLcnoYiCvcSkPsSEq8TRThqH2OjA9sTxhqvIOXZBuoqc0hOu/Hr2LVOgVCqCvl/OElewpxJMEooo2EtM6i5avitxOFyt/cnLjGSHs1jcM3Au3x6sTK3GU4JZIeEnFkICEfHBXsqmKuBGF+0VplOS7QXOcNsk8FSJw5k1IQHJCTHo6XPiTE0Lr/f2xVWmEycuNUOhANbyFOxVSiVWLZgEADh4Vr6ncZHIxetJVeFGyqYqAJikjoNKqUBfP4t6RvwTkc7W8rs4y1EoFFjKbXcsQSrH4q6tz5il8Xr0YKBWhsDRiyRy8XpSFeGXSqnE1GTXIm2NyKkcVyUOPxugeXOrhHl7iztfz+XZ+XLnfNfM/tD5Bsrbk5BCwV4GpqUkAABqmtpFfd7Gtm40t/dAoQBmT+avEofDVcCIHez7+p348KSr5DJ7CT8pHM6i6UlIGhONtq4+fH6J4fXehAgpoDp7PhUXF0OtVgMAGIZBQUGBT9cAgM3mOlDCaDR6HrNYLDAajcjJyYFWq0VpaSmWLVuG3Nxc/gfPk+kpY3AIwKVGcYN95ZVW9/MnID6G/1+FoYu0MVH8VfuMpMLehOb2HiQnxCArwPN0h6NSKnH73Al4/1gdDp6t91QdESJ3ks7suUCv1+uh1+uh0+lgMBhGvKawsBAFBQUoKCjwBPmcnBzP4wzDwGq1wmAwwGAwID09XdaBHgCmj+dm9uKmcaquus6JFWJWDwBpE8dCPSYa3b3iLtJye+Hcs2gyolT8/4qvnDcRAPDJuQbe702IUCQN9ps3b0Z+fr7n6+zsbJSUlAz781wgZxjG8z2DwQCLxQK73e75XkVFBViWhc1mg16vF2TsfJruTuNcFjmNU3nVNbOfM1mYQ2cUCgWWumf3Yi7Sfujuml2zeLIg979zvivYf3ruGpxOOoichAbJgr3dbgfDMJ4UzkAWi2XY68rLywcFdq1WCwCD3gBCzbQUV/ml2MFe6Jk9cCOVI1bevqWjB+U2Vw38msX85us5S2dqkBAbheb2Hk/pKiFyJ1nOfmDAHkitVg8buNVqNZqbmwd9j3tj4II+AOzYsQMajSvIlJWVoaioiIcRC2dGijRpHC5nP2eKcMdJir1I+/GZa+h3skifNA4z3OkxvkVHKXHbnPHYe/IqDp6tx6LpakGehxA+Sb5AO5RGo4HD4Xtg2Lx5M4xGo+cTglarhU6n8wR/h8OBvLw8mEwmr9dzTVWcDRs2iF5zn+oO9s3tPWjr7MW4+GjBn7Pf6YS9/joAgWf2aTcWaXv7nLwcjjKSfadc+frVi4RJ4XBWzp+EvSev4pNzDTDkzBP0uQgZzbZt2wY1hnprqpJdsPcn0BcWFsJgMAzKy+t0ukE/k5+fD4PBMGzKiGuqklJifDTUY6LBdPSipqk9qH3XfXWxoR29/U7ERis9aSQhDNzu+ExtC5bMFHYHzH2n6gEAqwXK13O45qqDZ6+BZVkoFMHvvUNIoIZOUgdOYDmS5ewHpl0GYhhm2McGMpvNSE9Pv2kB1mw2D/qaC/DDpY3kgqvIuSxSKofL16dPGgeVUrhfA1cnrTiLtDVN7ai80gqlQoG73NsaCCVTm4LYaCUaWrtgq28T9LkI4YOkwV6tVnsNwqOdW8vl6blAzzCMZ8E3Ly9v0D25/L8vbyBSEruxqspdiTNboEqcgXRaV7AXevMwrpFKp9VAnRAj6HPFRquQ6a7hP3SeSjCJ/Elaerlp06ZBlTdms3nQTN1ut3saqDhWqxVWqxU6nQ52ux1WqxWbN2+GRqOBWq2G0WgcFNhLSkqQm5vrNYUjJ9PdqZRLIgX7yiuu2eicKcLl6znL0scDgOCHde/jSi4FztdzbpvjSuVQsCehQNKcfUFBAYqLiz219TabzWs3LNdVyzAM1q5dC4ZhUFhYOOheXMVNfn7+oDeIpqamYRdn5WSayBU5Ys7sl892BfvTNYxgC9BOJ+s5lWq1QCWXQ90+1xXsP6sU9k2MRCarvQnvVNTgf/KW8nI/yRdoR9oegeus5XgrvRxKrVb7tOWC3MwQubGq0p2znyNgJQ5nsjoeM8Yn4FJjO6zVTbh7If8z75OXGTS2dSMhNgrLZ/O7RcJwuDexyiutaGjt4nV3TRK5mtt78L+m4/jTvkqwLLBizvigj9UEaCM02eAqYsSY2bd393nOvBVjZg8Ay9JdAfhIlTB5e25L4zvnTxRtDx7N2FjPoS+fVVIqhwRv38mr0BW8g1c/dAX6x++Y5WlMDBYFe5ngtkyodXQIvnWuzT2rT06IQcq4WEGfi8Pl7ctswqQ8LJ9zWxqLk8LheFI55ymVQ4LT1dOPZ185hMa2bsxPTcJ7m9bi1WfvwGR1PC/3j/hgL/VJVZzJ6nhEq5Tod7K4ynQK+lxcvl7Iztmhls2+sUjLsvzuJ3O9q9cTbPk6lcpXt83lFmnpMBMSnNf2V6GuuROpmjH4+Cf3e05F8wedVDUCqU+q4iiVCqRqXO/gQtfac9skCNk5O9TSmcmIiVKiqa0b1deu83rvj87Uo7ffibSJY5E+SbzXBNyY2R+70IyO7j5Rn5uEj86ePmx55xQA4IX1ixAXE1gqkk6qChFcRc5lgfe15xqqxJzZx0arPMcU8p3K4Q4WX7t4iuidrDPHJ2CyOh69/U5Yq+kQchKYP31YhfqWLswYn4Cn7hKmJ4iCvYyItdXxja2NxZ0Fc9UrZTwv0nL71/N9KpUvFAqFZ3Z/iPL2JADt3X341b9OAwAKHl4sWIEBBXsZmS5CRQ7LsgMaqsSb2QPCLNLa6ttgv3YdUSrht0gYzu2UtydBeG1fFRpau5A2cSy+tDJNsOeRvM6e3MClcYTsor3W0oXWzl4oFIB2ojQz+xOXmtHZ08fLUYjcrP72uRNE2S3UGy7YH65sRL/TKeheQ2R4juvdKD1Rh8tNHXjijlmevye52/ZJNQDg2/fPF3RXWAr2MjJjvPD743ApnJnjEwJeBArUtJQxmKyOx1WmExV2h+fEp2CUnnCVXK4VqWvWm1tmqJEYH43Wzl58fonhrS6a+OZCw3V860+H8fGZa3C6K71efu8MtupvwwMZ0yQe3cjO1rbg+MVmRKkUeHTFTEGfi6YgMsLtfHmpsZ338kTOeYlSOIArv32Hexb80emrQd+vu7cfH59xpU5ylkwN+n6BUimVnhLMg2cplSOm3j4nvvL7gzhwuh5OlsWi6Wosmq5Gc3sPHv/1R/ivv1eg3yls30owTIcuAACyb5kieM8LBXsZmTk+AQoFcL2rD41t3YI8h+d0KpE6Z4da466D33sy+GC///RVtHf3YUpyPBZLfFrUjUPIKdiLafM/P4e12oHkhBiUbX4Qn/3fAzjw4n345n2uA2X+8ME5/P79cxKP0juWZWH67CIAIP/2WYI/X8QHe7k0VQGu8sTUZNcirV2gPdLPi3AU4Ui4HSnLbU1g2nuCutfbZZcBAOsyp0GplPbwEC4l9clZOoRcLJ+eu4ZfvuOqYnn5meWY7966IjZahZeezMSvv7IMAPCznSc8kxw5Kbc3ofradSTERuEBHT/pJmqqGoFcmqo4WndTkJ3nxiNOlfuXfq5EwX76+ATMnZIIJ8viwOn6gO/T1+/Eu1bXL/TDWTP4Gl7AMmZpMCZGRYeQi6StsxcbjYfgZFk8uUqLLy6/+Xfg62tmY83iyejq7ce3/nxYdm/COz69AAB4KHMaEmL5WT6lpqoQkjZxLACgWoCZfXdvPy40uBZ/xdjHfjhr3EcGcufFBuLTcw1wXO+GZmws7nAfESil6CglVrj3t6dUjvBKLOdxqbEdM8cnoPjLmV5/RqFQ4OVnliMhNgqfnmvAqx9WijzK4fX1O7Hz8CUAQN7twi7McijYy4wn2Asws6++dh1OlsXYuCjeNlcKxBp35cyHQeTt3y53/aE8qEtFlEoev8ZcKocWaYXV0d2H379/FgDwo8eWIHGEktuZE8bifx+/FQDw4+3H0NDaJcYQR/XRmXo0tHYhZVws1iwSp5JMHn8lxEPINM75AYuzUh6QvWrBRESrlKi+dj2gtQmnk8U7FTUAgIeXTed7eAFbOSDYC1VNFYhzdS34xyfVOH7Bgd4++Vam+Or1AzY0tnVj1oQE5N02+qz4G2vmIGOWBu3dfdi6Rx6Lte+Uu35/12dNF7S2fiAK9jKjdc/shVig9VTiSJjCAYCxcdGeBqtAZvfl9iZcae7EuLgo3CPAQSiByky7cQg5dziMlBpau/Dd18qwfNN72Gg8hDt//D6mGkx48uWP0dbZK/XwAtLd24/fvncGAPDdBxf69KlOqVTg++sWAgBe2VuJ613Svnank8W7R13Bfl2meH0AFOxlJs09s29s6+b9D7JSgg3QhrP2FleQ/vCk/3l7rgrn/ltTERstbmPYSOJiVJ4tIaRO5bxrrUGG+xAMJ8ti6cxkJI2JRldvP3aXX8YTv/kIXT39ko4xEG8erEatowNTkuPx5CrfNwx7KHMaZk8eh+b2Hry23ybgCEdnrb4xWRFziw8K9jKTGB/taa7gO29/vs5VJSJVJc5AXN7+wOl6v1IL/U4n/lnmytevz5JPCofD5e0/PhN4pVGwymyN+OofPkFLRy+WzEjGe5vW4uBPv4BLf8zFe5vWYmxcFD46U4+v/vGTkErr9Dud+M27rlLL//jCAr/e6FVKJf7zgQUAgN/9+wx6+qR7o/uX1TWrz1kyVdTJCgV7GUoTIJXDsqwkWxsP59ZZyZiQGIfWzl58cLzO5+s+OF6HS43tSE6Iwb1LpeuaHc49i7hPLFcl6dy82HAdj//6I3T19uO+pVNx4Cf3eQ7BUCoVWLVgEnZ8927ERivxrrUG33m9TPQxBqr0xBXYr11HckIMnlk92+/rN6xMw2R1POqaO7Hj0EUBRugbrmT4IRFTOAAFe1k1VXHSBVikbWzrRrO7iUnsAz68USmV+NKdrh3+/rK/yufrSkrPAwCeuisdY3iqTebT8vTxSIyPhuN6N45WO0R97tbOXuT96gAaWrtwyww1/vL8Sq857VULJuGv31oFpUKBNw7YBDsqkm9/3uf6PXlylTaguvTYaJWns/a3752RZBG98korzta2IFqlFGSyQk1VI5BbUxUwsPySv5k9tzg7PWWMbILkV+9JB+DazOySDwe2nL/Sir0nr0KhADZmzxF6eAGJjlJ6ZvfcIehiKfxbBc7UtmCyOh6m790z4i6gX8hIxQb3m+1PTMfFGmLAapra8cEx1yfAQGb1nK+tmYOE2CicrW3BZ5Xiv8m9607hrFowEUljYni/PzVVhRghau3ltDjLmT05EXcvnASWBd44MPqiGTerv//WVMyaMFbo4QWMO/S89IR4wf7Dk1fwt4/tUCiA17+5EqmaMaNe88NHbkFMlBIHTtdjHw97FQnpjQM2OFkWqxZMDGrNKTE+Go+ucHXbvu7D7xzfuHz9Qzxtj+APCvYylObeZ57PnL3UG6AN52vuWdobH9nQ1z98jru1sxd/P2gHADybM1eUsQUq273ZW7mtCY7rwmxoN9D1rl58+89HAACG7Lm4Y55vW0fPGJ+Ar69x/ff/ifmYrHoDBurrd3oC89eCmNVzvnK36xPlW4cvoqUjuP2Z/HGtpRNHqlyfJvjaC8cfFOxlKH2Sa9Za4+hAdy8/VQPcfi3zpsor2D+UOQ3jx8XiSnMn3j82/ELttoN2XO/qw5wpibKqrfdmWkoCFqQmwcmyosyY/9d8HJca2zFjfAL+J2+pX9e+sH4xEmKjUGF3eBrV5Ob9Y3Woa+5EyrhYrMsMvgJr+ezxmJ+ahI6efuz8TLyF2nettWBZQJem8emTF98o2MvQhMQ4JMRGgWWBizwdPn6mhgEALJR4K+ChYqJUnnrpEst5r7PLS43t+NnOEwAAQ/YcyXe49AW3v36pwHn7MlsjtrrTWy8/sxxj4/w7rWtCYpxn0XLL7lO8j48P3AL+l1dpeSlVVCgUntm9mKmcf1W4+kMe4uENKxAU7GVIoVDwukjb0tGDy+5zbRdOUwd9P749s3o2VEoF9p266jl4mdPb58Qzf/wETEcvsrQpQS3OiYnL21tO1AmWHul3OvG918vAsq6ywrW3BLbHynP3zUNMlBJHLzhEryAaTT3TCYt77YML0Hx4YuUsRKuUsFY7cOJiM2/3HU5rZy/2u3d5FbNrdiAK9jLl2SOnPvhF2tM1rhTO1OR4JCfwXwEQrPRJ4/AL986FL5qO42130xQA/HTnCRypakTSmGj85ZsrERMln47Zkdw+dwLGxKhQ39KFzy8xgjzHa/ttOHahGUljovGzJ24N+D7jx8XhYXeDmj9lsGLYefginCyLrPQUXosLxo+L8wTd1w8I/5pLj9ehp8+JOVMSPfvui42CvUzxObPnUjiLZJbCGWhj9lzPwutG4yF860+Hse6lvfi1u2PyD19fIesKnKFio1W4a6GrmenfR/nPhTe2dXlKJn/06BJMTApuF9OvuRdqTYcuyGrfnB3uY/seF+Akp6fdnxRMhy7ytjY2nHe4FI4EC7OciA/2cmyqAm5siGbjYWZ/yh3sF6Sqg76XkDZ/SYd7l0xBZ08/Xj9g83zsff7eeXh4mfQHlPjri+4x/+PTC7yncn5iOo7m9h4snq7GxrXB9xysnDcRc6Yk4npXn+eoPKlVXW1Fhd0BlVKBx3zY3dJf9yyahKnJ8Whu78G/j97chMSX7t5+7HF3ia/LEjbYU1PVCOTYVAUAc91VM2d5OPXo1GXXPRZNl+bjo6+iVEq89s078f11C1H48GJs3Xgb9r94H156Uif10AKyPms64mNUqLrahnJ7E2/3PVLV6FlY3PJ0Fi/7+SsUCjzjbnJ7bZ88Ujkm95YGaxZPxoTEON7vr1IqPY1lXFmvEA6crkdbl+us5My0FMGeB6CmqpDELaReamxHaxAfq1mWxWkujSPDxdmhxsVH48W8W/Hfjy3Bk6u0yNSmSLr3fjDGxUd7Nmv7xyfVvNyzr9+J//zLEbAs8KU70zwHnfNhw51pslmoZVkW293H9gl5GPeX7nRVgpWeuIJrLZ2CPAeXwnkwQ9qzkinYy5RmbCymJLvysFywDsRVphPN7T1QKhSYN1XeM/tw9MQdswAA5s8u8bLT4h/3nMPJywySE2Lwsycygr7fQOPHxeGL7sNgfOloFpK12gFbfRviY1SCbhg2d0oilqWnoN95482FT/3OG2clC53CGY3km6QUFxdDrVYDABiGQUFBQdDXBHJPOVo8XY0rzZ04fZnBbXMCO2f11GUGAJA+eRziYkKjkiWc3LNoMiYlxaG+pQt7jl8JKnBdbmzH/7n7DX76RIYgqY0nV2mx49BF7DpyCcVfzhTtFKWhuMO4H9RN87t3wF9PrtKizNaEvx+sxrfun8/rJ8kDp13HD2rGxmLVfPH2rvdG0pk9F5T1ej30ej10Oh0MBkNQ1wRyT7niUjmngpjZn3bn/BdNo1m9FKJUSuS7Z/f/+DTwVA7LsvjeG2Xo6OnH7XMn4Ck/Du7wx10LJmFiUhwc17uxN4CDZfjgdLJ4q0y8w7gfXTETsdFKnLrM8F5z/49PLgAAHlsxQ7I3To6kz75582bk5+d7vs7OzkZJSUlQ1wRyT7niSiW5BdZAcDN7OZddhrsNK12LgP8+WuvZZtpfr+234f1jdYhWKfHbry4TLPcbpVIid4UrwJrcZY9i+6yyAVeaO5E0JhprFwt/GHdyQgwezHB94vrbx/wt1LZ393ny9dwbvpQkC/Z2ux0Mw3jSLQNZLJaArgnknnLGLaiermECLt07HSJll+HslhnJuGWGGj19Tvzh/bN+X195pRX/9fcKAMCPc5dggcAL7dxs+l8VNZKc1/rWEdes/oGMaaKd5PTUXa5PSts/vcDbcY3vWWtwvasPsyYkYIX7zGUpSRrsvVGr1WAYJqBrArmnnM2bmgiVUoHm9h7UNftfKdDvdOJsrWu3S7mXXYa7wocXAwBe/vcZ1Dk6fL6ut8+Jb2z9FB09/bh74ST8xxcWCDVEj0xtCrQTx6Kjpx/vWYWrP/fG6WTxT/cZw9xWxGJYvXgypqeMQXN7D3aXX+blntyC7+N3pMmiokx21TgajQYOh39lX6NdM9LjXFMV909OzVWx0SrMnuzaNoFLx/jDXn8dXb39iI9ReTpyiTTWZ03Hijnj0dnTj5/tOuHzdT/ecQzWageSE2Jg1N8uSumeQqHwpB22i5zKOVTZgKuMK4WzZrF4u5uqlEo8dRd/m6M1tHZ5Dq8RI4XDNVNx/0KiqcrfQO/LNSM9zjVVcf/k1ly1KIhFWi6FM39qElRK2f1PHVEUCgV+vsHVHPa3j+04eWn0hcCX/30Gv3enfV5+Zrmo2+LmuWvb935+BQ2tXaI971uHXSmcB3XTRN8H6am7tFAogI/O1MMW5FkSuw5fRL+ThS5NE9RhK77imqm4f7JqqtJqvVcTMAwz7GOjXRPIPeVu8Qw1AOB0ADP7oxccg+5BpLV89ng8snwGWBb44bajcDqHX4d586AdP9p2FADwYt5SfHG5uNtFzJ2SiIxZGvQ7WfzzyKXRL+BBv9OJt90plEdEfr2A6xyCHPfOocH0GbAsi79/7Kq8ErIhzF+SBnu1Wu01z56dnR3QNYHcU+6CKb887D5jc7kMFoeIy4t5SxGtUmLfqat4+g8H0dnTN+jxfqcTxtJzeP7VwwCAb90/H997aKEUQ/Us1Iq1V86h842SpHAGevpu14Zwf//Yjt6+4U9OG8nhqkYcveBAbLRSFlU4HEk/22/atGlQlYzZbIZer/d8bbfbUVxc7Nc1oz0eariSyXN1rX798vX2OVHh3o8l0IYswj/tpHEw6m9DTJQSb5ddxgOb9+LUZQaXGttRbmtE9k9L8YO/VqDfyWLDyjT83xMZki3uPbZiJhQK4ND5Blzm6RCdkbx1xPWmIkUKh/OFjKmYkOhqgnsvwN1K//jBOQCuWb0QjW+BkjTYFxQUgGEYlJSUoKSkBGVlZTAajZ7HLRbLoK99uWa0x0PNjJQEjI2LQk+f06884vGLDnT29CM5IUaUnCHxXd7ts/B2wRokJ8Sg3NaE2370HhZ9722s/skelNuaMC4uClueysT/t3GFpHupTNWM8ey9Yz4s7Oy+r9+Jt464Uji5Auxw6auYKJXnkJTfvnfG75Lny43tnmqe5++dx/v4giH5dgkjbWXAdcH6c40vj4cSpVKBBalJKLM14dRlxueDDz5zp3BWzBkfEsf4RZo750/E3h/fC0PJIZyra0VPnxNOlsUDGakoejITUyU4o9SbvNtn4eDZazAfuojvPihcOumjMze2FZD6jOFnc+bid++fQZmtCYfON/h8gDsAlOw9j34ni7sXTsLiGckCjtJ/kgd7MrpF09UoszXhxKVmn/f1PlzZAABYQSkc2ZozJREf/s99Ug9jRA9nTcf33yjDiUvNOFvbItgpS2b3usAXl02XfFuBSep4fOlOLf6yrwq/ee+Mz8G+vbvPsz30czKb1QMyLL0kN1uW7lpgPXj2mk8/z7KsZ2ZP+XoSjJRxsZ4tC8wCLdR29/bjHXfqQ4y9cHzx7fvnQ6FwbXHh65kSf//YDqajF9qJY3H/rVMFHqH/KNiHgNWLXB9rK+xNaOkYfW+Vi43tuMp0IlqlRKZWI/TwSJjjygdNh/g/cQsA9p68AqajF1OS43H7XHlMTuZMScS6TNd2z79978yoP8+092DzW58DAJ6/b54s+1rkNyKRyfVYwoGmj09A+qRx6HeyPs3uuRTO0lnJiI+hTB0JzgO6VMTHqGC/dh1lNv5O3OLsdH9ieHT5DFkFye886NqaYvunF1B1tXXEn/3ZzhNobOvGvKmJeGb1bDGG5xUdSzgCuR5LOBQ3u99/6uqoP/vZeffiLNXXEx6MjYvGw+4Tt/76Eb+HmnR093kO9xDinNlgLEsfj7WLJ6O337U/0XClz8cvOPDK3koAwC+fzpKsbBSgYwnDwj3uYL/Pl2DvntnL5SMxCX1Pu8sRzZ9d5HUnzHetNWjvdu0MmaUV9nzWQPz+6yuQnBCDCrsDP3enaQZyOll8741yOFkWj62YgbslriQaCQX7EHHXwklQKFzNVSPtmtjS0ePptqXFWcKXO+dPhHbiWFzv6vNsQcyH1/a7Pik8sVIeO0MONS0lAb99ZjkA4Jf/OjUojdrd24/vvl6GI1WNSIiN8ux9JFcU7ENEckIMdGmuxdb9p4ef3e8/dRUsC6RNHItJ6nixhkfCnEKhwFPu2f0bB/g54KPqais+OlMPhQJ42r3jpBw9snwGvrxKC5YF8n61H//5lyPYc7wOOT8txZ/dpZY//5JONr0Rw6FgH0Lu8SFvz23gJMVGUiS8PXlnGpQKBT6rbMC5usBPT+NwWwnn3DIF08cnBH0/If3iqUxkajW43tWHP++rwmO/3I+jFxzQjI3Frh/cg69JuCjrKwr2IWS1J29f77UErqapHaXuPbS5k3cI4cuU5DG4d6mr5v6vHwU3u+/p6/fsDPnVEAiUY+Oi8eGP78O7/7UW+be7zqxdMWc8Pvnp/chZIr+aem+oLi+ErJg9AXHRKlxlOnGurvWmbsa/fmQHywKrFkzE7Mm0Hw7h39N3p+P9Y3V482A1/vvRJYiLCazy5D1rLRpauzApKQ73L725ckSOlEoF7lo4CXctnIStfU5EqRSyXGcYDs3sQ0hcjAp3zne1bv9xz7lBj/U7nZ6yuK/eLf+ZEglN9y9NRapmDBpau4I6nPu1/a5c95fv0kq+PUIgoqOUIRXoAQr2IdFUNdAP1i8C4Ppj4bYwBoB9J6/iclMHkhNisN5dE00I36KjlPiuu9noV/86FdCe75VXWrH3pGvd6Ss0MeEVNVWNIFSaqjgr503EhpWzwLLA914vQ7/T9cf2mnux6/E7ZgX80ZoQXzx9dzomJsXhclMH/vFptd/XF719EgDwhYxUOhuZZ9RUFWZ++ngGEuOjYa124Mfbj+PRLfvwdplrIyluL25ChBIfE4X/+IJrdv/Ld055Jhy+OFvbgh3uQ8x/9MgtQgyPDIOCfQiapI7Hfz+2BIDrUOrSE1egUirw3QcXym4PbRKevr5mNjRjY2Grv45dh31vsnrpn5+DZYF1mdOwdBZt0icmCvYhauPaOVg+ezyiVAo8fXc6rEUP4X8fv1XqYZEIMTYuGt+637Vn+093nkBb5+hbKJy6zGCXu/v2hzSrFx2VXoaoKJUS//7hWvT0OTE2Llrq4ZAIZMiZh7/sq0L1tev4wV/LYdTfPuLP/9+uE2BZV8MffQIVH83sQ1hMlIoCPZFMYnw0Xn32DigVCrx5sBomdy7eG2PpObxTUQOlQoFNX1ws3iCJBwV7QkjA7pg3EYUPu8qBv/NaGez1bTf9TOmJOhT8zQoAeDF/KRZMU4s5ROJGwZ4QEpSChxfjtjkT0NrZi7v+5328urcSTieLfqcTB8/W4yu/Pwgny+Kpu7T4zgMLpB5uxFKwQpwzFkIyMzORmpqKDRs2hEytPSFyU+vowJMvf4QKuwMAMHvyODS0dqGlw7Vwe+f8iXi7YLWkB3tEgm3btmHbtm2ora1FRUXFoMciPtivX78eu3fvlnoYhIS8fqcTr+6txE9Mx9HW1QcAGBcXhZwlU/GrryxDyrhYiUcYObzFNarGIYTwQqVUwpAzD+uzpmP/6XrMm5KIJTOTEaWibLEcULAnhPBqSvIYbFiZJvUwyBD0lksIIRGAgj0hhEQACvaEEBIBKNgTQkgEoGAfoFA57CQY4f4a6fWFvnB/jXy+vogP9oGeVBXuv2RA+L9Gen2hL9xfYyBxiU6qGobYJ1UF88sZ6LVi/0GIPU56ffwL99cYrq+PO6nKm4gP9mKjYC+f6wIV7q8vmOcMldcY7q/P28w+4rdLWLRoEdLT/T/Kr7a21us5j0JdJ8Vz0nWReZ0Uz0nX8Xud1WpFTU3NoO9FfLAnhJBIQGkcQgiJABTsCSEkAlCwJ4SQCEC7XvqpuLgYarUaAMAwDAoKCqQdkACKi4sBADabDQBgNBqlHI6gcnJyUFpaKvUweFdYWOgpPNBoNMjNzZV4RPwqKSkBwzBQq9Ww2WzYtGmT5+8y1NjtdlgsFphMJq+/i7zFHJb4rKioiDUajZ6vS0tLWb1eL+GI+FdQUDDoa71ez2ZnZ0s0GmGZTCY23P4EmpubWZ1OxzY3N7Msy7IVFRVh9xqLioo8r49lXa85NzdXugEFoaKigjUajWxRURGr0+luepzPmBNevwUCU6vVg37JWJYNqz+k5uZmNjs7e9Br5IKFzWaTbmACaG5uZo1GY1j978eyrjfroqKiQd8rLS2VaDTC8BYUQ31CYjKZvL4uPmMO5ex9ZLfbPR8bh7JYLOIPSCDl5eWw2+2er7VaLQDXx8dwsmPHDuTn50s9DN4VFxcjOzvbkxoAgOzsbIlHxS+tVoucnBzP76Tdbvf8noYTvmMOBXsfDQyAA6nV6rAJhGq1Gs3NzdDpdJ7vcb9U4fTHZLFYwi4AAjd+R7kgodVqYTAYwmoyAgCvvPIKHA4HkpOTUVhYCIvFEpbrSnzHHAr2QdJoNHA4HFIPQzCbN2+G0WgM2cUvb7hAGG644KBWq6HT6aDValFUVIS8vDyJR8YvtVoNg8GA3NxcFBcXw2Qyhc2EyxeBxhwK9kEK50BfWFgIg8EAvV4v9VB4U1JSEnaVKRyNRgMAyMrK8nyPmwWG0+y+sLAQWq0WJpMJNpsNDocDmZmZUg9LNIHGHAr2PhpuJhius0Sz2Yz09PSwCvRWq3VQIAw33O/h0GCgVquHTQmEGi5FxaXhtFotKioqoNVqYTabJR4dv/iOOVRn7yOtVuv5oxn6Hzrc8r/cLJAL9AzDwOFwhPybmsPhgNVq9bw+ro+guLgYWq025Gf8arUaWq32pt9RhmHC5k3Obrd7TSmGW6oKECDmBFTDE6GG1ryaTKawq7OvqKhgi4qKWJvNxtpsNraiooItKCi4qfwrHIRjDbrJZBrUK2EymUK+LHGooeXBLMuG/N+h0Wj0qc4+mJhDu176aWA3m81mQ1FRkbQD4hHDMEhLS/O62BVuvyZmsxnbt2+H2WyGXq9HXl5e2HxC47pLAaCpqSmsfkcB1+/p5s2bkZKS4lmT0Ov1IVlEYLfbYTQaYbFYYLVaUVBQcFP6lK+YQ8GeEEIiAC3QEkJIBKBgTwghEYCCPSGERAAK9oQQEgEo2BNCSASgYE8IIRGAgj0hhEQACvaEEBIBKNgTEgSz2SybTcYGds4SMhQFe0ICZLFYZHVKkl6vx8aNG6UeBpEpCvYkogU6K2cYBkajEQUFBTyPKDibNm2CwWCQehhEhijYk4hltVphtVoDupY72EVudDod7Ha7bFJLRD4o2JOIFczpTeXl5bLdJdNgMITdTpckeBTsSUSyWq0oLCwM6Fqz2Szrw0Byc3OxY8cOqYdBZIaCPYk4ZrMZRqMRAGA0GmEwGGAwGHyuZNm+fTtycnKGfbykpMTzz2Aw+HRcntVqRU5ODpKTk1FSUuL5fl5eHpKTk5Genu7T2DgajSbgFBUJUwEdeUJIGADAmkwmv6/T6XSszWbz+pjJZGJzc3MHfU+tVvv8PAAGnUzEsixbUFDAarVav8aYm5t7031IZKOZPSF+stvt0Gg0wz4+dEadnZ2N0tJSn+7t7bQlf2f1gGtmTzX3ZCAK9oT4iWGYYY/Ay83N9RxkzjAMrFar58B2X4z0JuIPtVqNpqYmXu5FwgMFe0J4ZjabkZmZiY0bN8LhcITk2agk/FCwJ8Rt4MLoSLhDroe7x8aNG2EymWAymZCdnR30bD2QdAzDMEhJSQnqeUl4oWBPItbQGbevQVWj0QzbtMTVuA/cQmFgCseXN5Sh4+DSQv5wOByy2caByAMFexKxsrKyUFZWBsC1qKrT6Xy6jutS9QUXuL3938zMzJvKMrOzswcFd4ZhvHbEDnc9R0579hB5oGBPIpbRaITVakVxcTEsFovPHbGPP/74sNU1FRUVKC0tRXFxMUpKSmCxWGAymaBWq1FYWIjc3FzPz9rt9psWbouKiuBwOAZdn5eXBwDIzMwcVOnj7XoOwzA+v3mRyKBgWZaVehCEhJrMzExUVFRIPQyvLBYLSktLacsEMgjN7AkJQFZWVlB76wipqKhIlpu0EWlRsCckAEVFRZ4tF+SEW3ugfD0ZitI4hATIbDbD4XBAr9dLPRSPvLw8mEwmqYdBZIhm9oQEKDc3d8QyTLGVlJTglVdekXoYRKZoZk8IIRGAZvaEEBIBKNgTQkgEoGBPCCERgII9IYREgP8fNeTUNlh7B3cAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# plot the outcome\n",
+ "from matplotlib import pyplot as plt\n",
+ "import mpl_config\n",
+ "\n",
+ "plt.plot(tau * np.arange(101), np.array(z_list))\n",
+ "plt.xlabel(\"t (a.u.)\")\n",
+ "plt.ylabel(r\"$\\langle \\sigma_z \\rangle$\")\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "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
+}
\ No newline at end of file
diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst
index d3b3570..6e2c551 100644
--- a/docs/source/tutorials.rst
+++ b/docs/source/tutorials.rst
@@ -12,3 +12,5 @@ Tutorials
tutorial_jupyter/noisy_simulation.ipynb
tutorial_jupyter/sbm_dynamics.ipynb
tutorial_jupyter/marcus.ipynb
+ tutorial_jupyter/vbe_tutorial_groundstate.ipynb
+ tutorial_jupyter/vbe_tutorial_td.ipynb
diff --git a/example/vbe_td.py b/example/vbe_td.py
index 1880dc9..a957b5d 100644
--- a/example/vbe_td.py
+++ b/example/vbe_td.py
@@ -16,8 +16,7 @@
from tencirchem import set_backend, Op, Mpo, Model, OpSum
from tencirchem.dynamic import get_ansatz, get_deriv, get_jacobian_func, qubit_encode_basis, sbm
-
-from vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args
+from tencirchem.applications.vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args
set_backend("jax")
diff --git a/example/vbe_ti.py b/example/vbe_ti.py
index 6a6a636..5f3be8e 100644
--- a/example/vbe_ti.py
+++ b/example/vbe_ti.py
@@ -17,8 +17,7 @@
from tencirchem import set_backend, Op, BasisSHO, BasisSimpleElectron, Mpo, Model
from tencirchem.dynamic import get_ansatz, qubit_encode_op, qubit_encode_basis
from tencirchem.utils import scipy_opt_wrap
-
-from vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args
+from tencirchem.applications.vbe_lib import get_psi_indices, get_contracted_mpo, get_contract_args
backend = set_backend("jax")
@@ -130,12 +129,12 @@ def f(x):
es = []
for k, new_b in enumerate(sols):
if not np.allclose(new_b @ new_b.T, np.eye(4)):
- print(f"Enforcing orthogonality for the {k}th root of b_array[{i}]")
+ # print(f"Enforcing orthogonality for the {k}th root of b_array[{i}]")
new_b = np.linalg.qr(new_b.T)[0].T
b_array[i] = new_b
e = psi @ get_contracted_mpo(h_mpo, b_array, n_qubit_per_mode, b_dof_pidx, psi_idx_top + psi_idx_bottom) @ psi
es.append(e)
- print(np.array(es))
+ # print(np.array(es))
lowest_id = np.argmin(es)
return sols[lowest_id]
diff --git a/tencirchem/applications/__init__.py b/tencirchem/applications/__init__.py
new file mode 100644
index 0000000..0cf883e
--- /dev/null
+++ b/tencirchem/applications/__init__.py
@@ -0,0 +1,4 @@
+# Copyright (c) 2023. The TenCirChem Developers. All Rights Reserved.
+#
+# This file is distributed under ACADEMIC PUBLIC LICENSE
+# and WITHOUT ANY WARRANTY. See the LICENSE file for details.
diff --git a/example/vbe_lib.py b/tencirchem/applications/vbe_lib.py
similarity index 100%
rename from example/vbe_lib.py
rename to tencirchem/applications/vbe_lib.py